# Smoothness of integrated density of states of the Anderson model on Bethe
lattice in high disorder

Dhriti Ranjan Dolai
Indian Institute of Technology Dharwad,
Dharwad - 580011, India.
email<EMAIL_ADDRESS>

M. Krishna
Ashoka University, Plot 2, Rajiv Gandhi Education City,
Rai, Haryana 131029, India
email<EMAIL_ADDRESS>
Abstract: In this work we consider the Anderson model on Bethe lattice and
prove that the integrated density of states (IDS) is as smooth as the single
site distribution (SSD), in high disorder.

MSC (2020): 81Q10, 47B80, 35J10.
Keywords: Anderson Model, random Schrödinger operators, integrated density of
states, Bethe lattice.

##  Introduction

The Bethe lattice $\mathbb{B}$ is an infinite connected graph with no closed
loops and a fixed degree $K+1$ (number of nearest neighbours) at each vertex.
The degree is called the coordination number and the connectivity $K$ is one
less the coordination number. Distance between any two vertex $n_{1}$ and
$n_{2}$ is denoted by $d(n_{1},n_{2})$ and it is the length of shortest path
connecting $n_{1}$ and $n_{2}$. We denote $\ell^{2}(\mathbb{B})$ to be the
Hilbert space
$\bigg{\\{}u:\mathbb{B}\to\mathbb{C};\displaystyle\sum_{n\in\mathbb{B}}|u(n)|^{2}<\infty\bigg{\\}}$.
The Anderson Model on Bethe lattice is a random Schrödinger operator
$H^{\omega}$ on the Hilbert space $\ell^{2}(\mathbb{B})$ and it is given by

$\displaystyle H^{\omega}$ $\displaystyle=\Delta+\lambda
V^{\omega},~{}\lambda>0,~{}\omega\in\Omega,$ (1.1) $\displaystyle(\Delta
u)(n)$
$\displaystyle=\sum_{k:d(k,n)=1}u(k),~{}u=\\{u(n)\\}_{n\in\mathbb{B}}\in\ell^{2}(\mathbb{B}),$
$\displaystyle(V^{\omega}u)(n)$ $\displaystyle=\omega_{n}u(n),$

where $\\{\omega_{n}\\}_{n\in\mathbb{B}}$ are i.i.d real random variables with
common distribution $\mu$ and it has a compact support. Consider the
probability space
$\big{(}\mathbb{R}^{\mathbb{B}},\mathcal{B}_{\mathbb{R}^{\mathbb{B}}},\mathbb{P}\big{)}$,
where $\mathbb{P}=\underset{n\in\mathbb{B}}{\otimes}\mu$ is constructed via
the Kolmogorov theorem. We refer to this probability space as
$\big{(}\Omega,\mathcal{B}_{\Omega},\mathbb{P}\big{)}$ and denote
$\omega=(\omega_{n})_{n\in\mathbb{B}}\in\Omega$. The operator $\Delta$ is
known as the graph Laplacian with diagonal term removed and the potential
$V^{\omega}$ is the multiplication operator on $\ell^{2}(\mathbb{B})$ by the
sequence $\\{\omega_{n}\\}_{n\in\mathbb{B}^{d}}$. We note that the operators
$\\{H^{\omega}\\}_{\omega\in\Omega}$ are bounded self-adjoint. It follows from
the ergodicity (see [12, appendix] and [5, 15]) that the spectrum of the
random operators $H^{\omega}$ is given by
$\sigma(H^{\omega})=[-2\sqrt{K},2\sqrt{K}]+\lambda~{}supp~{}\mu$, a.e
$\omega$.

First, we fix a vertex (we can choose any one) as the origin and denote it by
$0$. Now we define the finite sub-graph $\Lambda_{L}$ and its boundary
$\partial\Lambda_{L}$ as

$\Lambda_{L}=\big{\\{}n\in\mathbb{B}:d(n,0)\leq
L\big{\\}}~{}\&~{}\partial\Lambda_{L}=\big{\\{}(n,m)\in\Lambda_{L}\times\big{(}\mathbb{B}\setminus\Lambda_{L}\big{)}:d(n,m)=1\big{\\}}.$

Let $\chi_{{}_{L}}$ and $P_{n}$ denote the orthogonal projection onto the
subspace $\ell^{2}(\Lambda_{L})$ and $\ell^{2}(\\{\delta_{n}\\})$,
respectively. Here $\\{\delta_{n}\\}_{n\in\mathbb{B}}$ is the standard basis
for the Hilbert space $\ell^{2}(\mathbb{B})$ . Define the finite matrix
$H^{\omega}_{L}$ of size $|\Lambda_{L}|$ as

$H^{\omega}_{L}=\Delta_{L}+V^{\omega}_{L},~{}\Delta_{L}=\chi_{{}_{L}}\Delta\chi_{{}_{L}},~{}V^{\omega}_{L}=\displaystyle\sum_{n\in\Lambda_{L}}\omega_{n}P_{n}.$

We consider the integrated density of states (IDS) $\mathcal{N}(\cdot)$, as it
is given in [12, equation (1.4)], namely

$\displaystyle\mathcal{N}(E):=\mathbb{E}\big{(}\langle\delta_{0},E_{H^{\omega}}(-\infty,E]\delta_{0}\big{)},$
(1.2)

where $E_{A}(\cdot)$ denote the spectral measure of a self-adjoint operator
$A$. The probability measure
$\nu(\cdot):=\mathbb{E}\big{(}\langle\delta_{0},E_{H^{\omega}}(\cdot)\delta_{0}\big{)}$
is known as the density of states measure (DOSm). In the Hypothesis 1.1
(below), we assume the absolute continuity of $\mu$, therefore the Wegner
estimate (see [11]) will give the absolute continuity of the DOSm and in that
case we have $d\nu(x)=\rho(x)dx$ and $\mathcal{N}^{\prime}(x)=\rho(x)$.

In this paper we are interested in the degree of smoothness of $\rho$, for
which we make the following assumptions on $H^{\omega}$ (1.1) :

######  Hypothesis 1.1.

  1. (a)

The single site distribution (SSD) $\mu$ is absolutely continuous with respect
to the Lebesgue measure on $\mathbb{R}$, i.e $d\mu(x)=g(x)dx$ and we also
assume the function $g\in C^{m}_{c}(\mathbb{R})$, set of all continuously
differentiable function up to the order $m$, with compact support.

  2. (b)

Let the compact interval $J\subset\mathbb{R}$ be in the region of exponential
localization. Stated precisely, there exists a constant $C_{s}>0$ such that

$\sup_{\Re(z)\in
J,~{}\Im(z)>0}\mathbb{E}\bigg{(}\bigg{|}\big{\langle}\delta_{n_{1}},\big{(}H^{\omega}_{L}-z\big{)}^{-1}\delta_{n_{2}}\big{\rangle}\bigg{|}^{s}\bigg{)}\leq
C_{s}e^{-\xi_{s}d(n_{1},n_{2})},~{}C_{s}>0,$ (1.3)

for some $\xi_{s}>0,~{}0<s<1$, for any $n_{1},n_{2}\in\mathbb{B}$ and
$\forall~{}L>0$. The decay rate $\xi_{s}$ satisfies the condition
$\xi_{s}>2\ln(K+1)\big{(}3+p\big{)}$ for some $1\leq p\leq m$. We also assume
that the constants $C_{s},~{}\xi_{s}$ do not change if we replace the
distribution $g$ with its derivatives at finitely many vertices.

######  Remark 1.2.

Exponential decay in (1.3) can be obtained with established techniques
developed to prove the localization of Anderson model via fractional moments
methods of the Green’s function. We refer to Aizenman [3], Aizenman-Molchanov
[2] (see also Aizenman et al [1]) for more details. With the same method one
can make $\xi_{s}$ (the decay rate in the exponential (1.3)) arbitrary large
by taking the disorder parameter $\lambda$ large enough, in the definition of
$H^{\omega}$ (1.1). Going through the proofs given in [3, 2] and keep track of
the constants, one can also realize that the assumption about changing the
distribution at finitely many sites is also valid.

Now we are ready to state our main result:

######  Theorem 1.3.

Consider the operator $H^{\omega}$ given in (1.1) satisfying the hypothesis
1.1, then the integrated density of states $\mathcal{N}(\cdot)$ is
continuously differentiable up to the order $p$ in the interval $J$. In other
words,

$\mathcal{N}(E)=\mathbb{E}\big{(}\langle\delta_{0},E_{H^{\omega}}(-\infty,E]\delta_{0}\big{)}\in
C^{p}(J),~{}\text{for some}~{}1\leq p\leq m.$ (1.4)

Several results are known about the smoothness of the integrated density of
states (IDS) of the Anderson model on $\ell^{2}(\mathbb{Z}^{d})$, under
different types of assumptions on the single site distribution (SSD). We refer
to Constantinescu et al. [13], Bovier et al [16], Kaminaga et al [14],
Bellissard-Hislop [17], Carmona-Lacroix [5, Corollary VI.3.2] for more
details. It is also possible to calculate the exact expression of the IDS for
different models, when the SSD is Cauchy, we refer to Lloyd [10], Carmona-
Lacroix [5, problem VI.5.5] and Kirsch-Krishna [18] for more details.
For the one-dimensional random Schrödinger operators, much more rich
literature can be found in Companino-Klein [19], March-Sznitman [25], Simon-
Taylor [8], Speis [20], Klein-Speis [21, 22], Klein et al [23] and Glaffig
[24]. Brodie, in his doctoral thesis [9] investigated the smoothness of the
density of states (DOS) for the random band matrix.
The higher order differentiability of the integrated density of states (IDS)
for continuum model on $\L^{2}(\mathbb{R}^{d})$ was obtained by Dolai et al
[7]. Although, their methods can also be applied to a large class of discrete
model, but as $|\Lambda_{L}|$, the volume of the finite sub-graph (of Bethe
lattice) increases exponentially as $L$ gets larger, the differentiability of
IDS for Bethe lattice model is not immediate from [7].
In [12], Acosta-Klein showed the analyticity of the IDS for the Bethe lattice
model when the single site distribution (SSD) is Cauchy or close to Cauchy (in
some function space). The result in [12] is also valid in the region of
absolutely continuous spectrum.
In this work, we prove that for the Bethe lattice model, the integrated
density of states (IDS) and the single site distribution (SSD) have the same
degree of smoothness, in high disorder.

##  Proof of the results

In this section we will give the proof of our main result Theorem 1.3. Before
going to the proof of this theorem, we need some pre requisite results. First,
we show that the fractional moments of the difference of two resolvents of
finite volume restriction of $H^{\omega}$ to the sub-graphs $\Lambda_{L+1}$
and $\Lambda_{L}$ is decay exponentially, as $L$ increases.
Let $G^{\omega}_{L}(z;n_{1},n_{2})$ denote the Green’s function of the
operator $H^{\omega}_{L}$ and it is defined by
$G^{\omega}_{L}(z;n_{1},n_{2}):=\big{\langle}\delta_{n_{1}},\big{(}H^{\omega}_{L}-z\big{)}^{-1}\delta_{n_{2}}\big{\rangle}$.

######  Lemma 2.1.

Under the condition [(b), Hypothesis 1.1], the expected value of the
difference of the resolvents can be bounded as, for some $B_{s}>0$

$\sup_{\Re(z)\in
J,~{}\Im(z)>0}\mathbb{E}\bigg{(}\big{|}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\big{|}^{\frac{s}{2}}\bigg{)}\leq
B_{s}e^{-(L+1)\big{[}\frac{\xi_{s}}{2}-\frac{2L+3}{L+1}\ln(K+1)\big{]}}.$
(2.1)

###### Proof.

Let’s begin with the resolvent identity

$\displaystyle(H^{\omega}_{L+1}-z)^{-1}-(H^{\omega}_{L}-z)^{-1}$
$\displaystyle\qquad=(H^{\omega}_{L+1}-z)^{-1}\bigg{[}(H^{\omega}_{L}-z)-(H^{\omega}_{L+1}-z)\bigg{]}(H^{\omega}_{L}-z)^{-1}$
$\displaystyle\qquad=(H^{\omega}_{L+1}-z)^{-1}\bigg{[}\chi_{{}_{L}}\Delta\chi_{{}_{L}}-\chi_{{}_{L+1}}\Delta\chi_{{}_{L+1}}-\sum_{n\in\Lambda_{L+1}\setminus\Lambda_{L}}\omega_{n}P_{n}\bigg{]}(H^{\omega}_{L}-z)^{-1}$
$\displaystyle\qquad=(H^{\omega}_{L+1}-z)^{-1}\bigg{(}\sum_{n\in\Lambda_{L+1}\setminus\Lambda_{L}}P_{n}\bigg{)}\Delta\chi_{{}_{L}}(H^{\omega}_{L}-z)^{-1}.$
(2.2)

In the last line, we wrote
$\chi_{{}_{L+1}}=\chi_{{}_{L}}+\sum_{n\in\Lambda_{L+1}\setminus\Lambda_{L}}P_{n}$
and also we have used the fact that the operator
$\bigg{(}\displaystyle\sum_{n\in\Lambda_{L+1}\setminus\Lambda_{L}}P_{n}\bigg{)}(H^{\omega}_{L}-z)^{-1}\equiv
0$.
Now using the definition of the Laplacian $\Delta$ in (2), we write

$\displaystyle G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)$
$\displaystyle\qquad=-\sum_{m\in\Lambda_{L}}\sum_{\begin{subarray}{c}k:d(k,m)=1\\\
k\in\Lambda_{L+1}\setminus\Lambda_{L}\end{subarray}}\big{\langle}\delta_{m},(H^{\omega}_{L}-z)^{-1}\delta_{0}\big{\rangle}\big{\langle}\delta_{0},(H^{\omega}_{L+1}-z)^{-1}\delta_{k}\big{\rangle}.$
(2.3)

A bound for the number of vertex in $\Lambda_{L}$ can be given by
$|\Lambda_{L}|\leq D(K+1)^{L+1}$, similarly
$\big{|}\Lambda_{L+1}\setminus\Lambda_{L}\big{|}\leq D(K+1)^{L+2}$, $D$ is
independent of $L$.
Now take the fractional power of the both side of the equation (2),

$\displaystyle\big{|}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\big{|}^{\frac{s}{2}}$
$\displaystyle\qquad\leq\sum_{m\in\Lambda_{L}}\sum_{\begin{subarray}{c}k:d(k,m)=1\\\
k\in\Lambda_{L+1}\setminus\Lambda_{L}\end{subarray}}\big{|}\big{\langle}\delta_{m},(H^{\omega}_{L}-z)^{-1}\delta_{0}\big{\rangle}\big{|}^{\frac{s}{2}}\big{|}\big{\langle}\delta_{0},(H^{\omega}_{L+1}-z)^{-1}\delta_{k}\big{\rangle}\big{|}^{\frac{s}{2}}.$
(2.4)

After taking the expectation both side of the above inequality (2), we use
Cauchy-Schwartz inequality to write

$\displaystyle\mathbb{E}\bigg{(}\big{|}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\big{|}^{\frac{s}{2}}\bigg{)}$
$\displaystyle\qquad\leq\sum_{\begin{subarray}{c}m\in\Lambda_{L}\\\
k\in\Lambda_{L+1}\setminus\Lambda_{L}\\\
k:d(k,m)=1\end{subarray}}\mathbb{E}\bigg{(}\big{|}\big{\langle}\delta_{m},(H^{\omega}_{L}-z)^{-1}\delta_{0}\big{\rangle}\big{|}^{\frac{s}{2}}\big{|}\big{\langle}\delta_{0},(H^{\omega}_{L+1}-z)^{-1}\delta_{k}\big{\rangle}\big{|}^{\frac{s}{2}}\bigg{)}$
$\displaystyle\qquad\leq\sum_{\begin{subarray}{c}m\in\Lambda_{L}\\\
k\in\Lambda_{L+1}\setminus\Lambda_{L}\\\
k:d(k,m)=1\end{subarray}}\mathbb{E}\bigg{(}\big{|}\big{\langle}\delta_{m},(H^{\omega}_{L}-z)^{-1}\delta_{0}\big{\rangle}\big{|}^{s}\bigg{)}^{\frac{1}{2}}\mathbb{E}\bigg{(}\big{|}\big{\langle}\delta_{0},(H^{\omega}_{L+1}-z)^{-1}\delta_{k}\big{\rangle}\big{|}^{s}\bigg{)}^{\frac{1}{2}}$
$\displaystyle\qquad\leq
D^{2}C_{s}^{2}(K+1)^{2L+3}e^{-\frac{\xi_{s}}{2}(L+1)}.$ (2.5)

In the last line we have used the volume estimation of
$\Lambda_{L},~{}\Lambda_{L+1}\setminus\Lambda_{L}$ and the fractional moments
given in (1.3). Now (2.1) is immediate. ∎

We can always conclude the differentiability of the density of an absolutely
continuous measure from the derivatives of its Borel transform.

######  Lemma 2.2.

Let $d\tilde{\mu}(x)=f(x)dx$ be an absolutely continuous measure on the real
line and $I\subset\mathbb{R}$ an interval. Denote
$F(z):=\int_{\mathbb{R}}\frac{1}{x-z}d\tilde{\mu}(x)$ to be the Borel
transformation of $\tilde{\mu}(\cdot)$. Then for any $m\in\mathbb{N}$

$\underset{x\in
I}{ess~{}sup}~{}\bigg{|}\frac{d^{m}}{dx^{m}}f(x)\bigg{|}<\infty~{}\text{whenever}~{}\sup_{\Re(z)\in
I,~{}\Im(z)>0}\bigg{|}\Im{\bigg{(}\frac{d^{m}}{dz^{m}}F(z)\bigg{)}}\bigg{|}<\infty.$
(2.6)

The proof can be found in [7, Lemma A.1].

The above Lemma will give the higher order differentiability of $\rho$, the
density of states function (DOS), once we show that the imaginary part of the
higher order derivatives of the expected value of the Green’s function of
$H^{\omega}$ is bounded on the upper half plane.

######  Proposition 2.3.

Under the condition [(a) Hypothesis 1.1], the imaginary part of the higher
order derivative of the Green’s function $G^{\omega}_{M}(z;0,0)$ is bounded on
the upper half of the complex plane

$\max_{1\leq\ell\leq
p}\sup_{\Im(z)>0}\bigg{|}\Im{\bigg{(}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{M}(z;0,0)\big{)}\bigg{)}}\bigg{|}<\infty~{}~{}\text{for}~{}~{}1\leq
p\leq m,$ (2.7)

here $M$ is a fixed positive integer.

###### Proof.

Let $z=E+i\epsilon,~{}\epsilon>0$, we can re-write the expression for
$H^{\omega}_{M}-z$ as

$\displaystyle H^{\omega}_{M}-z$
$\displaystyle=\Delta_{M}+\sum_{n\in\Lambda_{M}}(\omega_{n}-E)P_{n}-i\epsilon$
(2.8)
$\displaystyle=\Delta_{M}+\sum_{n\in\Lambda_{M}}\tilde{\omega}_{n}P_{n}-i\epsilon,~{}~{}\tilde{\omega}_{n}=\omega_{n}-E$
$\displaystyle=H^{\tilde{\omega}}_{M}-i\epsilon,~{}~{}\tilde{\omega}=(\tilde{\omega}_{n})_{n\in\Lambda_{M}}.$
(2.9)

We write the expected value of the Green’s function as the integration

$\displaystyle\mathbb{E}\big{(}G^{\omega}_{M}(z;0,0)\big{)}$
$\displaystyle\qquad=\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{M}}|}}}\langle\delta_{0},\big{(}H^{\omega}_{M}-z\big{)}^{-1}\delta_{0}\big{\rangle}\prod_{n\in\Lambda_{M}}g(\omega_{n})d\omega_{n}$
$\displaystyle\qquad=\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{M}}|}}}\langle\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}\prod_{n\in\Lambda_{M}}g(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}.$
(2.10)

Using analyticity of $G^{\omega}_{M}(z;0,0)$ and the representation (2), we
can write its higher order derivatives in terms of the product of the
derivatives of $g$,

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{M}(z;0,0)\big{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle=\frac{d^{\ell}}{dE^{\ell}}\mathbb{E}\big{(}G^{\omega}_{M}(z;0,0)\big{)},~{}~{}z=E+i\epsilon,~{}\epsilon>0$
$\displaystyle=\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{M}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{M}|-1}}\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{M}}|}}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n=0}^{|\Lambda_{M}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$
(2.11)

Use resolvent identity for one dimensional perturbation to write

$\displaystyle\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}$
$\displaystyle=G^{\tilde{\omega}}_{M}(i\epsilon;0,0)$
$\displaystyle=\frac{1}{\lambda\tilde{\omega}_{0}+\big{(}G^{\tilde{\omega}\setminus\tilde{\omega}_{0}}_{M}(i\epsilon;0,0)\big{)}^{-1}},~{}$
$\displaystyle=\frac{1}{\lambda\tilde{\omega}_{0}+A-iB}.$

here
$G^{\tilde{\omega}\setminus\tilde{\omega}_{0}}_{M}(i\epsilon;0,0)=(A-iB)^{-1}$
denote the Green’s function for the matrix
$H^{\tilde{\omega}\setminus\tilde{\omega}_{0}}_{M}:=H^{\tilde{\omega}}_{M}-\lambda\omega_{0}P_{0}$.
It is easy check that $B>0$ for $\Im(z)>0$.
Since $g\in C^{m}_{c}(\mathbb{R})$, therefore we estimate the integral

$\displaystyle\bigg{|}\Im\bigg{(}\int_{\mathbb{R}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}g^{(k_{0})}(\tilde{\omega}_{0}+E)d\omega_{0}\bigg{)}\bigg{|}$
$\displaystyle\qquad=\int_{\mathbb{R}}\frac{B}{(\lambda\omega_{0}+A)^{2}+B^{2}}\big{|}g^{(k_{0})}(\tilde{\omega}_{0}+E)\big{|}d\omega_{0},$
$\displaystyle\qquad\leq\big{\|}g^{(k_{0})}\big{\|}_{\infty}\frac{\pi}{\lambda^{2}},~{}~{}~{}0\leq
k_{0}\leq m.$ (2.12)

The integral in (2) can be written as

$\displaystyle\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{M}}|}}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}\bigg{(}\prod_{n=0}^{|\Lambda_{M}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}$
$\displaystyle\qquad\qquad=\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{M}}|-1}}}\bigg{(}\int_{\mathbb{R}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{M}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}g^{(k_{0})}(\omega_{0}+E)d\omega_{0}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n\neq
0}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$ (2.13)

Now (2.7) can be easily obtained, once we use (2) and (2) in (2). ∎

Using [7, Theorem 2.2] we can estimate the average of the difference of the
resolvents by its fractional moments.

######  Proposition 2.4.

Under the assumption [(a) Hypothesis 1.1], we have

$\displaystyle\bigg{|}\int\langle\delta_{0},\big{(}H^{\omega}_{L+1}-z\big{)}^{-1}-\big{(}H^{\omega}_{L}-z\big{)}^{-1}\delta_{0}\rangle
g(\omega_{0})d\omega_{0}\bigg{|},~{}~{}~{}\Im(z)>0$
$\displaystyle\qquad\leq\Xi\int\bigg{|}\langle\delta_{0},\big{(}H^{\omega}_{L+1}-z\big{)}^{-1}-\big{(}H^{\omega}_{L}-z\big{)}^{-1}\delta_{0}\rangle\bigg{|}^{\frac{s}{2}}\phi_{R}(\omega_{0})d\omega_{0}.$
(2.14)

Here $0<s<2$ and $\phi_{R}$ is a translation of a smooth indicator function
$\chi_{R}$ on a bounded interval which contain the support of the density $g$.
The constant $\Xi>0$ is independent of $\omega$, $z$ and $L$ but depends on
$s$ and $g$.

###### Proof.

Using the notations of [7, Theorem 2.2] we define
$F:=|\delta_{0}\rangle\langle\delta_{0}|$, $\rho_{1}:=g$,
$A:=H_{L+1}^{\omega\setminus\omega_{0}}=H_{L+1}^{\omega}-\omega_{0}~{}|\delta_{0}\rangle\langle\delta_{0}|$
and
$B:=H_{L}^{\omega\setminus\omega_{0}}=H_{L}^{\omega}-\omega_{0}~{}|\delta_{0}\rangle\langle\delta_{0}|$.
Now (2.4) is immediate from [7, (2) of Theorem 2.2] . ∎

######  Remark 2.5.

For more details on the above result and the properties of $\phi_{R}$,
$\chi_{R}$ and $\Xi$, we refer to [7, Theorem 2.2].

Now the above proposition will give the exponential decay of the higher order
derivatives of the expected value of the difference of the resolvents
$G^{\omega}_{L+1}(z;0,0)$ and $G^{\omega}_{L}(z;0,0)$.

######  Proposition 2.6.

Under the Hypothesis 1.1, we have

$\displaystyle\sup_{\Re(z)\in
J,~{}\Im(z)>0}\bigg{|}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}\bigg{|},~{}~{}1\leq\ell\leq
m$
$\displaystyle\qquad\qquad\qquad\leq\tilde{C}_{s}e^{-(L+1)\big{[}\frac{\xi_{s}}{2}-\frac{\ln(K+1)}{L+1}\big{(}(2L+3)+(L+2)(\ell+1)\big{)}\big{]}},~{}\tilde{C}_{s}>0.$
(2.15)

###### Proof.

As in the proof of the inequality (2), we write

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L+1}(z;0,0)\big{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle=\frac{d^{\ell}}{dE^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L+1}(z;0,0)\big{)},~{}~{}z=E+i\epsilon,~{}\epsilon>0$
$\displaystyle=\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|}}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n=0}^{|\Lambda_{L+1}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$
(2.16)

We also write

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle=\frac{d^{\ell}}{dE^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)},~{}~{}z=E+i\epsilon,~{}\epsilon>0$
$\displaystyle=\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L}|-1}}\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L}}|}}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n=0}^{|\Lambda_{L}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$
(2.17)

The above integrand is independent of the random variables
$\\{\omega_{n}\\}_{n\in\Lambda_{L+1}\setminus\Lambda_{L}}$ and under the
assumption [(a), Hypothesis] $\int g^{(j)}(x)dx=\delta_{0,j},~{}1\leq j\leq
m$. We use this two facts to write the (2) as

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle=\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|}}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n=0}^{|\Lambda_{L+1}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$
(2.18)

Now the derivatives of the difference can be written as

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle=\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}$
$\displaystyle\qquad\qquad\qquad\qquad\times\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|}}}\bigg{(}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}-\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n=0}^{|\Lambda_{L+1}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\bigg{)}.$

For $\ell<<L$, in the product
$\displaystyle\prod_{n=0}^{|\Lambda_{L+1}|-1}g^{(k_{n})}(\tilde{\omega}_{n}+E)$
at most $\ell$ number of places $g^{(k_{n})}$ is different from $g$ and all
the other remaining places we have $g^{(k_{n})}=g$. Therefore, w.l.o.g we can
always assume $g^{(k_{0})}=g$. Now the above expression can be rewritten as

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)},~{}~{}0\leq\ell\leq
m$
$\displaystyle\leq\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}$
$\displaystyle\qquad\times\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|-1}}}\bigg{(}\int_{\mathbb{R}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}-\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}g^{(k_{0})}(\tilde{\omega}_{0}+E)d\omega_{0}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n\neq
0}g^{(k_{n})}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\big{)}.$ (2.19)

To estimate the above difference, we are going to use the Proposition 2.4. As
the function $g^{(k_{0})}=g$ is the density of $\mu$, the single site
distribution (SSD) we have,

$\displaystyle\bigg{|}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}\bigg{|},~{}~{}0\leq\ell\leq
m$
$\displaystyle\leq\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}\times$
$\displaystyle\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|-1}}}\bigg{|}\int_{\mathbb{R}}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}-\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}g^{(k_{0})}(\tilde{\omega}_{0}+E)d\omega_{0}\bigg{|}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n\neq
0}\big{|}g^{(k_{n})}(\tilde{\omega}_{n}+E)\big{|}d\tilde{\omega}_{n}\big{)}$
$\displaystyle\leq\Xi\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}\times$
$\displaystyle\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|-1}}}\bigg{(}\int_{\mathbb{R}}\bigg{|}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}-\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}\bigg{|}^{\frac{s}{2}}\phi_{R}(\tilde{\omega}_{0}+E)d\omega_{0}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n\neq
0}\big{|}g^{(k_{n})}(\tilde{\omega}_{n}+E)\big{|}d\tilde{\omega}_{n}\big{)}.$
$\displaystyle\leq\Xi~{}C^{\ell}\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}\times$
$\displaystyle\int_{{}_{\mathbb{R}^{|\Lambda_{{}_{L+1}}|-1}}}\bigg{(}\int_{\mathbb{R}}\bigg{|}\big{\langle}\delta_{0},\big{(}H^{\tilde{\omega}}_{L+1}-i\epsilon\big{)}^{-1}-\big{(}H^{\tilde{\omega}}_{L}-i\epsilon\big{)}^{-1}\delta_{0}\big{\rangle}\bigg{|}^{\frac{s}{2}}\tilde{g}_{0}(\tilde{\omega}_{0}+E)d\omega_{0}\bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times\bigg{(}\prod_{n\neq
0}\tilde{g}_{n}(\tilde{\omega}_{n}+E)d\tilde{\omega}_{n}\big{)}.$ (2.20)

In the above (2), we have used the notation
$\tilde{g}_{n}(x)=\frac{|g^{(k_{n})}(x)|}{\|g^{(k_{n})}\|_{1}},~{}n\neq 0$ and
$\tilde{g}_{0}(x)=\frac{\phi_{R}(x)}{\|\phi\|_{1}}$. As discussed earlier that
in the product $\displaystyle\prod_{n\neq 0}\tilde{g}_{n}$, at most $\ell$
number of places $\tilde{g}_{n}\neq g$ and all the other remaining places
$\tilde{g}_{n}$ will be $g$, with $\|g\|_{1}=1$. Since $\phi_{R}$ is a
compactly supported smooth function and $g$, the density of $\mu$ (SSD) is in
$C^{m}_{c}(\mathbb{R})$, then there exist a constant $C>0$ such that
$\displaystyle\max_{n:~{}0\leq k_{n}\leq
m}\big{\\{}\|g^{(k_{n})}\|_{1},~{}\|\phi_{R}\|_{1}\big{\\}}\leq C$ holds.
Assume $\tilde{\mathbb{E}}(\cdot)$ denote the expectation with respect to the
product probability measure
$\displaystyle\prod_{n}\tilde{g}_{n}(\omega_{n})d\omega_{n}$. Now the (2) can
be rewritten as

$\displaystyle\bigg{|}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}\bigg{|},~{}~{}0\leq\ell\leq
m$
$\displaystyle\leq\Xi~{}C^{\ell}\sum_{\begin{subarray}{c}k_{0}+\cdots+k_{|\Lambda_{L+1}|-1}=\ell\\\
\forall j,~{}k_{j}\geq
0\end{subarray}}\binom{\ell}{k_{0},\cdots,k_{|\Lambda_{L+1}|-1}}$
$\displaystyle\qquad\qquad\qquad\qquad\times\tilde{\mathbb{E}}\bigg{(}\bigg{|}\big{\langle}\delta_{0},\big{(}H^{\omega}_{L+1}-z\big{)}^{-1}-\big{(}H^{\omega}_{L}-z\big{)}^{-1}\delta_{0}\big{\rangle}\bigg{|}^{\frac{s}{2}}\bigg{)},$
(2.21)

here we have used the notation $z=E+i\epsilon$ and
$\omega=\big{(}\tilde{\omega}_{n}+E\big{)}_{n\in\Lambda_{L+1}}$.
In view of the assumption [(b), Hypothesis 1.1] and using the Lemma 2.1, we
estimate the r.h.s of (2)

$\displaystyle\bigg{|}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}\bigg{|}$
$\displaystyle\qquad\qquad\qquad\leq\Xi~{}C^{\ell}\tilde{C}\big{|}\Lambda_{L+1}\big{|}^{\ell+1}B_{s}e^{-(L+1)\big{[}\frac{\xi_{s}}{2}-\frac{2L+3}{L+1}\ln(K+1)\big{]}}$
$\displaystyle\qquad\qquad\qquad\leq\Xi~{}C^{\ell}\tilde{C}B_{s}e^{-(L+1)\big{[}\frac{\xi_{s}}{2}-\frac{\ln(K+1)}{L+1}\big{(}(2L+3)+(L+2)(\ell+1)\big{)}\big{]}},$
(2.22)

here we estimate the number of terms in the sum (2) by
$\tilde{C}|\Lambda_{L+1}|^{\ell+1}$ and $|\Lambda_{L+1}|$, the number of
vertex in $\Lambda_{L+1}$ can be bounded by a constant multiplication of
$(K+1)^{L+2}$. Hence the proposition. ∎

We also require the convergence of $G^{\omega}_{L}(z;0,0)$, the Green’s
function of finite volume approximation $H^{\omega}_{L}$ to
$G^{\omega}(z;0,0)$, the Green’s function of the full operator $H^{\omega}$.
Although, for $\mathbb{Z}^{d}$ model, this convergence of Green’s function is
relatively straightforward, for the Bethe lattice model, it demands proof, and
it was given in [12, Proposition 1.2]. Using the same result [12, Proposition
1.2], we can also conclude the convergence of the derivatives of
$G^{\omega}_{L}(z;0,0)$ to the derivatives of $G^{\omega}(z;0,0)$.

######  Proposition 2.7.

For $\Im(z)>0$ and $\ell\in\mathbb{N}$, we have

$\lim_{L\to\infty}\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)}=\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}$
(2.23)

and

$\lim_{L\to\infty}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)}=\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}.$
(2.24)

###### Proof.

For the proof of (2.23), we refer to [12, Proposition 1.2].
Define $g_{L}(z):=\mathbb{E}\big{(}G^{\omega}_{L}(z;0,0)\big{)}$ and it is
analytic function on the upper half of the complex plane. Now we use the
Cauchy integral formula to write $g^{(\ell)}_{L}$, the $\ell$th order
derivative of $g_{L}$ as

$g_{L}^{(\ell)}(z)=\frac{\ell!}{2\pi
i}\int_{\gamma}\frac{g_{L}(w)}{(w-z)^{\ell+1}}dw,$ (2.25)

here $\gamma$ be a circle with center at $z$, oriented counter clockwise and
it is lie inside the upper half of the complex plane. Now the convergence
(2.24) is immediate once we use (2.23) and dominated convergence theorem in
(2.25). ∎

In view of the Lemma 2.2, to prove the higher order differentiability of
$\rho$, the density of states function (DOSf), we only need to show the
imaginary part of the higher order derivatives of
$\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}$, expected value of the Green’s
function of the full operator $H^{\omega}$ is uniformly bounded in the upper
half of the complex plane.

Proof of the Theorem 1.3: Since $G^{\omega}_{L+1}(z;0,0)$ converges to
$G^{\omega}(z;0,0)$ (see (2.23)), therefore for a fix $M\in\mathbb{N}$ (large
enough), using telescoping series we write

$\displaystyle\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}$
$\displaystyle=\mathbb{E}\big{(}\big{\langle}\delta_{0},(H^{\omega}-z)^{-1}\delta_{0}\big{\rangle}\big{)}$
$\displaystyle=\lim_{L\to\infty}\mathbb{E}\big{(}G^{\omega}_{L+1}(z;0,0)\big{)}$
$\displaystyle=\sum_{L=M}^{\infty}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}+G^{\omega}_{M}(z;0,0).$
(2.26)

Using the convergence of the derivatives of $G_{L}^{\omega}(z;0,0)$, as in
(2.24), we can write its derivative as the infinite series of the derivatives
of its finite volume approximation

$\displaystyle\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}$
$\displaystyle=\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}\big{\langle}\delta_{0},(H^{\omega}-z)^{-1}\delta_{0}\big{\rangle}\big{)},~{}1\leq\ell\leq
p$
$\displaystyle=\lim_{L\to\infty}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}_{L+1}(z;0,0)\big{)}$
$\displaystyle=\sum_{L=M}^{\infty}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\bigg{(}G^{\omega}_{L+1}(z;0,0)-G^{\omega}_{L}(z;0,0)\bigg{)}+\frac{d^{\ell}}{dz^{\ell}}G^{\omega}_{M}(z;0,0).$
(2.27)

Since $\xi_{s}>2\ln(K+1)\big{(}3+p\big{)}$, assumed in [(b), Hypothesis 1.1],
then use of (2.6) and (2.7) in (2) will give

$\displaystyle\sup_{\Re(z)\in
J,~{}\Im(z)>0}\bigg{|}\Im\bigg{(}\frac{d^{\ell}}{dz^{\ell}}\mathbb{E}\big{(}G^{\omega}(z;0,0)\big{)}\bigg{)}\bigg{|}<\infty~{}~{}\text{for}~{}~{}1\leq\ell\leq
p.$ (2.28)

The Green’s function $G^{\omega}(z;0,0)$ is nothing but the Borel
transformation of $\nu$, the density of states measure (DOSm), which is
absolutely continuous in our case and $\mathcal{N}$, the IDS is its
distribution function. Hence an application of the Lemma 2.2 to the estimate
(2.28) will give the result. ∎

Acknowledgment

Dhriti Ranjan Dolai is partially supported by the INSPIRE Faculty Fellowship
grant of DST, Govt of India.

## References

  * [1] Aizenman. M, Schenker. J. H, Friedrich. R. M, Hundertmark. D: Finite volume fractional-moment criteria for Anderson localization, Commun. Math. Phys. 224(1), 219$-$253, 2001.
  * [2] Aizenman. M, Molchanov. S: Localization at large disorder and at extreme energies: an elementary derivation, Commun. Math. Phys, 157(2), 245$-$278, 1993.
  * [3] Aizenman. M: Localization at weak disorder: some elementary bounds, Rev. Math. Phys, 6(5A), 1163$-$1182 1994.
  * [4] Kirsch. W: An Invitation to Random Schrödinger Operators, (With an appendix by Frédéric Klopp) Panor. Synthèses, 25, Random Schrödinger operators, 1, Soc. Math. France, Paris, 1$-$119, 2008.
  * [5] Carmona. C, Lacroix. J: Spectral Theory of Random Schrödinger Operators, Birkhäuser, Boston, 1990.
  * [6] Cycon. H. L, Froese. R. G, Kirsch. W, Simon. B: Schrödinger operators, with application to quantum mechanics and global geometry, Texts and Monographs in Physics, Springer-Verlag, Berlin, 1987.
  * [7] Dolai. D, Krishna. M, Mallick. A: Regularity of Density of States of Random Schrödinger Operators, Commun. Math. Phys, 378(1), 299$-$328, 2020.
  * [8] Simon. B, Taylor. M: Harmonic analysis on $SL(2,R)$ and smoothness of the density of states in the one-dimensional Anderson model. Commun. Math. Phys. 101(1), 1$-$19, 1985.
  * [9] Brodie. Benjamin : Eigenvalue Statistics and Localization for Random Band Matrices with Fixed Width and Wegner Orbital Model, theses and dissertations, University of Kentucky, 2020.
  * [10] P. Lloyd:  Exactly solvable model of electronic states in a three dimensional disordered hamiltonian: non existence of localized states, J. Physics (C), 2, 1717$-$1725, 1969.
  * [11] Combes. J M, Germinet. F, Klein. A: Generalized Eigenvalue-Counting Estimates for the Anderson Model, J. Stat. Phys. 135(2), 201$-$216, 2009.
  * [12] Acosta. V, Klein. A: Analyticity of the density of states in the Anderson model on the Bethe lattice, J. Stat. Phys. 69(1-2), 277$-$305, 1992.
  * [13] Constantinescu. F, Fröhlich. J, Spencer, T: Analyticity of the density of states and replica method for random Schrödinger operators on a lattice, J. Stat. Phys. 34(3$-$4), 571$-$596, 1984.
  * [14] Kaminaga. M, Krishna. M, Nakamura, S: A note on the analyticity of density of states, J. Stat. Phys. 149(3), 496$-$504, 2012.
  * [15] Pastur. L: Spectra of random selfadjoint operators, Russ. Math. Surv. 28, 1$-$73, 1973.
  * [16] Bovier. A, Campanino. M, Klein. A, Perez. J.F: Smoothness of the density of states in the Anderson model at high disorder, Commun. Math. Phys. 114(3), 439$-$461, 1988.
  * [17] Bellissard. J, Hislop. P: Smoothness of Correlations in the Anderson Model at high disorder, Ann. Henri Poincaré 8, 1$-$26, 2007.
  * [18] Kirsch. W, Krishna. M: Analyticity of the density of states for the Cauchy distribution, Preprint, 2020.
  * [19] Campanino. M, Klein. A: A supersymmetric transfer matrix and differentiability of the density of states in the one-dimensional Anderson model, Commun. Math. Phys. 104(2), 227$-$241, 1986 .
  * [20] Speis. A: Weak disorder expansions for the Anderson model ona one-dimensional strip at the center of the band. Commun. Math. Phys. 149(3), 549$-$571, 1992.
  * [21] Klein. A, Speis. A.: Smoothness of the density of states in the Anderson model on a one-dimensional strip, Ann. Phys. 183(2), 352$-$398, 1988.
  * [22] Klein. A, Speis, A: Regularity of the invariant measure and of the density of states in the one-dimensional Anderson model, J. Funct. Anal. 88(1), 211$-$227, 1990.
  * [23] Klein. A, Lacroix. J, Speis. A: Regularity of the density of states in the Anderson model on a strip for potentials with singular continuous distributions. J. Stat. Phys. 57(1$-$2), 65$-$88, 1989.
  * [24] Glaffig. C: Smoothness of the integrated density of states on strips, J. Funct. Anal. 92(2), 509$-$534, 1990.
  * [25] March. P, Sznitman. A: Some connections between excursion theory and the discrete Schrödinger equation with random potentials, Probab. Theory Relat. Fields. 75(1), 11$-$53, 1987. 1$-$31, 2017.
  * [26] Reed. M, Simon. B. Methods of modern mathematical physics: Functional analysis, volume 1, Academic Press Inc. 1980.


# Analytical Characterization of Epileptic Dynamics in a Bistable System

Yuzhen Qin, Ahmed El-Gazzar, Danielle S. Bassett, Fabio Pasqualetti, and
Marcel van Gerven This work was supported in part by the project Dutch Brain
Interface Initiative (DBI2) with project number 024.005.022 of the research
programme Gravitation which is (partly) financed by the Dutch Research Council
(NWO). It is also supported in part by NSF CMMI-2308639. Y. Qin, A. El-Gazzar,
and M. van Gerven are with the Donders Institute for Brain, Cognition and
Behaviour, Radboud University, Nijmegen, the Netherlands. F. Pasqualetti is
with the Department of Mechanical Engineering, University of California at
Riverside. D. S. Bassett is with the Department of Bioengineering, the
Department of Electrical & Systems Engineering, the Department of Physics &
Astronomy, the Department of Psychiatry, and the Department of Neurology,
University of Pennsylvania, and The Santa Fe Institute. (Contact:
yuzhen.qin@donders.ru.nl).

###### Abstract

Epilepsy is one of the most common neurological disorders globally, affecting
millions of individuals. Despite significant advancements, the precise
mechanisms underlying this condition remain largely unknown, making accurately
predicting and preventing epileptic seizures challenging. In this paper, we
employ a bistable model, where a stable equilibrium and a stable limit cycle
coexist, to describe epileptic dynamics. The equilibrium captures normal
steady-state neural activity, while the stable limit cycle signifies seizure-
like oscillations. The noise-driven switch from the equilibrium to the limit
cycle characterizes the onset of seizures. The differences in the regions of
attraction of these two stable states distinguish epileptic brain dynamics
from healthy ones. We analytically construct the regions of attraction for
both states. Further, using the notion of input-to-state stability, we
theoretically show how the regions of attraction influence the stability of
the system subject to external perturbations. Generalizing the bistable system
into coupled networks, we also find the role of network parameters in shaping
the regions of attraction. Our findings shed light on the intricate interplay
between brain networks and epileptic activity, offering mechanistic insights
into potential avenues for more predictable treatments.

##  I INTRODUCTION

Epilepsy affects over 65 million individuals worldwide [1]. Patients with
epilepsy usually experience recurring seizures characterized by the sudden
emergence of excessively synchronous neural activity [2]. Anti-seizure
medications, relying on stabilizing electrical activity in the brain to reduce
the likelihood of seizures, are the most widely used treatment for epilepsy.
However, one-third of the patients remain resistant to such drugs [3].
Neurostimulation and surgery are alternative treatments. Yet, only 50% of
patients after surgical resection have long-term freedom of seizures. When it
comes to neurostimulation, the success rate is even lower. Nearly 95% of
individuals under vagus nerve stimulation therapies continue to have seizures
[4].

One of the main factors for the limited success of current treatments is the
lack of understanding of the precise mechanisms underlying epilepsy as a
network disease [5]. Dynamical system modeling is promising since it offers
the potential to uncover the mechanisms underlying emerging behaviors
associated with epileptic seizures and subsequently identify specific targets
for intervention [6].

Related work. Various models have been proposed to capture epileptic dynamics
[7]. Many of them use two fundamental states to describe epileptic activity in
the brain or brain regions: an interictal one characterized by normal steady-
state neural activity, and an ictal one marked by excessively synchronous
oscillations [8, 9]. In the realm of dynamical systems, these states
correspond to equilibrium and limit cycle, respectively. Some studies focus on
characterizing how parameter alterations destabilize the equilibrium, leading
systems to seizure-like oscillations. In [10, 11, 12], researchers consider
networks of linear-threshold units and construct conditions under which the
limit cycle is the only stable state. Using the same model, the work [13]
analytically investigates how bifurcations induced by external inputs engender
pathological oscillations. Subsequent research endeavors propose network
control strategies to contain the spread of such oscillations [14] or restore
desired oscillations [15, 16, 17]. Additionally, other investigations link
epileptic dynamics with instability in discrete-time linear fractional-order
systems [18].

The studies above primarily focus on exploiting the instability of equilibrium
within systems to describe epileptic activity. Alternative approaches propose
models that allow for bistability. This bistability is ideal for capturing the
transition from a stable equilibrium to a stable limit cycle, which signifies
the onset of epileptic seizures, wherein normal steady-state neural activity
gives way to excessive oscillations [19, 20]. These transitions can be induced
by endogenous or exogenous perturbations. A simple model with two stable
states was introduced to capture similar transitions [21]. Subsequent research
efforts have extended this model to networks of bistable units, exploring the
influence of network structure on transition possibilities through numerical
investigations [22]. Inspired by these studies, we utilize this bistable model
as a mesoscale building block for describing firing-rate neural dynamics.

Figure 1: Nullclines and vector fields of the system (1) in different
parameter regimes. In all cases, $a=1,b=1$, and $\omega=0.4$.

Contribution. Our contributions are threefold. Firstly, employing Lyapunov
methods and LaSalle’s invariance principle, we rigorously identify the regions
of attraction for both the stable equilibrium and limit cycle within the
bistable system. We put forth that disparities in these regions of attraction
serve as distinguishing factors between epileptic and healthy brains. Systems
with smaller regions of attraction towards the equilibrium are more prone to
transitioning into pathological oscillations in the presence of noise.
Secondly, by leveraging the concept of input-to-state stability, we construct
conditions for external perturbations to ensure the system remains in
proximity to the equilibrium. This also formally demonstrates that systems
with smaller regions of attraction are more vulnerable to perturbations.
Finally, we investigate a network comprising bistable units and analytically
estimate the overall network’s regions of attraction. This analysis sheds
light on the interplay between network parameters and the network’s
susceptibility to perturbations.

In addition to our theoretical results, we discover that networks composed of
bistable units can accurately replicate EEG recordings from epileptic brains
through parameter training. Notably, this approach, utilizing bistable
dynamics as a generative model, holds significant promise in understanding the
network pathology of epilepsy. Together with our analytical findings, this
method could offer valuable insights into potential objectives for controlling
epileptic dynamics, thereby informing the design of targeted, predictive, and
closed-loop treatments.

Notation: Given $x\in\mathbb{R}^{n},y\in\mathbb{R}^{m}$, define ${\rm
col}(x,y)\coloneq[x^{\top},y^{\top}]^{\top}$. A function
$\beta\colon[0,a)\times[0,\infty)\to[0,\infty)$ is called a class
$\mathcal{KL}$-function if for any fixed $s$, $\beta(a,s)$ is increasing and,
for any fixed $a$, $\beta(a,s)$ is decreasing and
$\lim_{s\to\infty}\beta(a,s)=0$. Given a set
$\mathcal{C}\subset\mathbb{R}^{n}$, the distance from a point
$x\in\mathbb{R}^{n}$ to $\mathcal{C}$ is denoted as ${\rm
dist}(x,\mathcal{C}):=\inf_{y\in\mathcal{C}}\|x-y\|$.

##  II Problem Formulation

We first introduce the following model to describe the activity of individual
neuronal populations:

$\displaystyle\dot{x}$ $\displaystyle=-\omega
y+x\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)\,,$ (1a)
$\displaystyle\dot{y}$ $\displaystyle=\omega
x+y\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)\,,$ (1b)

where $x,y\in\mathbb{R}$, $\omega,a,b>0$, and $\sigma\in\mathbb{R}$.

This model is an equivalent variation of the one found in [21]. Previous work
has demonstrated that special versions of the system exhibit different
behaviors as the parameters $\omega,\sigma,a,b$ vary [23, 22]. We start by
presenting a comprehensive characterization of the possible behaviors
exhibited by the system (1) across various parameter regimes.

First, let us rewrite the system into polar coordinates. Let
$re^{i\theta}=x+iy$, then we arrive at

$\displaystyle\dot{r}=r(\sigma+2abr^{2}-br^{4})\,,$ (2a)
$\displaystyle\dot{\theta}=\omega\,.$ (2b)

The system (2) has at most three possible equilibria

$r=0\,,\hskip 5.0ptr=\sqrt{a+\gamma_{0}}\,,\hskip
5.0ptr=\sqrt{a-\gamma_{0}}\,,$

with $\gamma_{\mu}:=\sqrt{a^{2}+(1-\mu)\sigma/b}$ for any $\mu\in[0,1]$, which
respectively correspond to the equilibrium ${\rm col}(x,y)=0$, denoted as
$e_{0}$, and the limit cycles,

$\displaystyle\mathcal{C}_{1}\coloneq\\{x,y\in\mathbb{R}\colon
x^{2}+y^{2}={a+\gamma_{0}}\\}\,,$
$\displaystyle\mathcal{C}_{2}\coloneq\\{x,y\in\mathbb{R}\colon
x^{2}+y^{2}={a-\gamma_{0}}\\}$

of the original system (1). Note that not all of these potential equilibria
coexist and are stable simultaneously. We next characterize their stability at
different parameter regimes.

######  Lemma 1.

For the system (1), the following statements hold:

  1. (i)

if $\sigma<-a^{2}b$, the system has a unique stable equilibrium ${\rm
col}(x,y)=0$.

  2. (ii)

if $\sigma=-a^{2}b$, the system has a stable equilibrium ${\rm col}(x,y)=0$
and a semi-stable limit cycle $\mathcal{C}_{1}$.

  3. (iii)

if $-a^{2}b<\sigma<0$, the system has a stable equilibrium ${\rm col}(x,y)=0$,
a stable limit cycle $\mathcal{C}_{1}$, and an unstable limit cycle
$\mathcal{C}_{2}$.

  4. (iii)

if $\sigma\geq 0$, the system has a stable limit cycle $\mathcal{C}_{1}$ and
unstable equilibrium ${\rm col}(x,y)=0$.

The nullclines and vector fields of the system (1) in Fig. 1 provide a more
intuitive illustration of the situations in this lemma. Most interestingly, in
Case (iii), the system exhibits a bistable behavior, where a stable limit
cycle and a stable equilibrium coexist. Around the equilibrium, the system
behaves as a damped oscillator; on the limit cycle, the system is a self-
sustained oscillator with amplitude $\sqrt{a+\gamma_{0}}$ and frequency
${2\pi}/{\omega}$. Note that $\sigma$ serves as the bifurcation parameter,
whereas $a$ primarily governs the oscillation amplitude of the limit cycles
and $\omega$ determines the oscillation frequency. We refer to the system (1)
as a bistable oscillator in the remainder of this paper.

In the bistable regime of $-a^{2}b<\sigma<0$, the system is particularly
suitable for describing epileptic dynamics. The stable equilibrium at the
origin describes ‘background’ neural activity at a normal state without
epileptic seizures. Conversely, the stable limit cycle characterizes the
oscillating neural activity observed during pathological conditions in the
presence of seizures (known as ictal states [22]). The shift from equilibrium
to a limit cycle captures the transition from a normal steady state to seizure
activity. External inputs or internal noise serves as the driving force behind
the switch between these two stable states.

In this paper, we aim to investigate the chance of such transitions by
rigorously studying the regions of attractions of the above two states and
investigating how the system responds to external inputs. Additionally, we
delve into networks comprising these bistable units, exploring the influence
of network parameters on the regions of attraction.

##  III Individual Bistable Oscillators

###  III-A Regions of attraction

######  Definition 1.

Given a system $\dot{z}=f(z,t),z\in\mathbb{R}^{n}$ that has a stable
equilibrium $e_{0}$ and a stable limit cycle $\mathcal{C}$, let $\phi(t,z)$ be
the solution of this system that starts at the initial state $z$. The basin of
attraction of $e_{0}$ is defined by

$\mathcal{A}_{e}\coloneq\\{z\in\mathbb{R}^{n}\colon\lim_{t\to\infty}\phi(t,z)=e_{0}\\}\,.$

The basin of attraction of $\mathcal{C}$ is defined by

$\mathcal{A}_{\ell}\coloneq\\{z\in\mathbb{R}^{n}\colon\lim_{t\to\infty}{\rm
dist}(\phi(t,z),\mathcal{C})\to 0\\}\,.$

Next, we study the regions of attraction of the equilibrium ${\rm col}(x,y)=0$
and the limit cycle $\mathcal{C}_{1}$ for the system (1).

######  Theorem 1 (Regions of Attraction).

For the system (1), consider the bistable regime where $-a^{2}b<\sigma<0$.
Then, the following statements hold:
(i) The region of attraction of ${\rm col}(x,y)=0$ is

$\mathcal{A}_{e}\coloneq\\{x,y\in\mathbb{R}\colon
x^{2}+y^{2}<{a-\gamma_{0}}\\}\,.$

(ii) The region of attraction of the limit cycle $\mathcal{C}_{1}$ is

$\mathcal{A}_{\ell}\coloneq\\{x,y\in\mathbb{R}\colon
x^{2}+y^{2}>{a-\gamma_{0}}\\}\,.$

$x$$y$$(a)$$(b)$$x$$y$noise Figure 2: Illustration of bistability when
$-a^{2}b<\sigma<0$. (a) A stable equilibrium (the origin) and a stable limit
cycle (solid circle) coexist. The region of attraction for the origin is
depicted by the shaded open area. Starting from inside, the system converges
to the origin; starting from outside, it converges to the limit cycle. The
separatrix is the dashed circle. (b) Different regions of attractions provide
a possible interpretation of the discrepancy between healthy and epileptic
brains.

Theorem 1 provides a rigorous characterization of the regions of attraction
for limit cycles and equilibria. An intuitive illustration is shown in Fig.
2a. The unstable limit cycle $\mathcal{C}_{2}$, depicted by the dashed circle
in Fig. 2a, is the separatrix of the two stable modes of behavior. The radius
of the separatrix is given by $a-\gamma_{0}=a-\sqrt{a^{2}+\sigma/b}$. For
fixed $a$ and $b$, as $\sigma$ approaches 0, the radius decreases; as $\sigma$
approaches $a^{2}b$, the radius approaches $a$.

From a dynamical systems perspective, we posit that the different sizes or
shapes of regions of attraction distinguish an epileptic brain from a healthy
one. In a healthy brain, the region of attraction for the equilibrium is
sufficiently large, thus making it unlikely for neural activity to transition
into oscillations in response to noise. In contrast, the epileptic brain
possesses a smaller region of attraction, rendering it more susceptible to
noise. Even under the influence of the same level of noise, the epileptic
brain is more prone to transitioning into oscillations (see Fig. 2b).

###  III-B Forced Bistable Oscillators

Given the noisy nature of the brain, fluctuations often exert continuous
influence on the neural activity. Motivated by this observation, we consider
the forced model given by

$\displaystyle\dot{x}=-\omega
y+x\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)+u_{1}(t)\,,$ (3a)
$\displaystyle\dot{y}=\omega
x+y\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)+u_{2}(t)\,,$ (3b)

where the inputs $u_{1}(t)$ and $u_{2}(t)$ are used to model perturbations to
the system, such as noise and external stimuli. Denote $u\coloneq{\rm
col}(u_{1},u_{2})$. We assume that the system (1) operates in the bistable
regime where $-a^{2}b<\sigma<0$.

Ideally, we want the system to stay around the equilibrium, even when
subjected to external perturbations, as the equilibrium corresponds to the
normal state in the brain. Next, we rigorously explore the criteria that
initial conditions and inputs should meet to achieve this objective.

######  Theorem 2 (Input-to-State Stability).

For any $\mu\in(0,1)$, let

$\mathcal{B}_{\mu}\coloneq\\{x,y\colon x^{2}+y^{2}<a-\gamma_{\mu}\\}\,,$ (4)

where $\gamma_{\mu}:=\sqrt{a^{2}+(1-\mu)\sigma/b}$. For any initial condition
${\rm col}(x(0),y(0))\in\mathcal{B}_{\mu}$ and any input $u$ satisfying

$\displaystyle\sup_{t\geq
0}\|u(t)\|<(1-\varepsilon)\mu|\sigma|\sqrt{a-\gamma_{\mu}}$ (5)

for any $\varepsilon\in(0,1)$, the solution to the system (3) satisfies

$\displaystyle\|{\rm col}(x(t),y(t))\|$
$\displaystyle\leq\beta\left(\left\|{\rm
col}\big{(}x(0),y(0)\big{)}\right\|,t\right)+\frac{1}{(1-\varepsilon)\mu|\sigma|}\sup_{t\geq
0}\|u(t)\|\,,$ (6)

where $\beta(\cdot)$ is a $\mathcal{KL}$-function.

Figure 3: Evolution of the system (3) at different regimes with the same
initial condition and identical continuous perturbations ((a): $\sigma=-0.8$,
(b): $\sigma=-0.2$). In both cases, $a=1,b=1,\omega=4$. Perturbations are
Gaussian with mean $0.2$ and STD $0.5$ truncated at $[-0.3,0.7]$.

###### Proof.

Consider a continuously differentiable function
$V(x,y)=\frac{1}{2}(x^{2}+y^{2})$. Its time derivative satisfies

$\displaystyle\dot{V}(x,y)$
$\displaystyle=(x^{2}+y^{2})\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)$
$\displaystyle\hskip 113.81102pt+xu_{1}(t)+yu_{2}(t)$
$\displaystyle\leq(x^{2}+y^{2})\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right)$
$\displaystyle\hskip 113.81102pt+\|u\|\sqrt{x^{2}+y^{2}}\,,$

where the Cauchy–Schwarz inequality has been used.

Next, we bound the term $\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}$. To do
that, we define a function $f(\eta)\coloneq\sigma+2ab\eta-b\eta^{2}$ with
$\eta=x^{2}+y^{2}\geq 0$. It can be derived that $f(\eta)$ is an increasing
function for $\eta\in[0,a-\gamma_{\mu})$, which implies that

$\displaystyle f(\eta)<f(a-\gamma_{\mu})=\mu\sigma\,.$

Therefore, it follows that

$\displaystyle\dot{V}(x,y)\leq\mu\sigma(x^{2}+y^{2})+\|u\|\sqrt{x^{2}+y^{2}}$
$\displaystyle\leq\varepsilon\mu\sigma(x^{2}+y^{2})+(1-\varepsilon)\mu\sigma(x^{2}+y^{2})+\|u\|\sqrt{x^{2}+y^{2}}\,.$

For any ${\rm col}(x,y)$ that belongs to $\mathcal{B}_{\mu}$ and satisfies

$\displaystyle\sqrt{x^{2}+y^{2}}\geq\frac{1}{(1-\varepsilon)\mu|\sigma|}\,,$

one can derive that

$\displaystyle\dot{V}(x,y)\leq\varepsilon\mu\sigma(x^{2}+y^{2})\leq 0\,,$

where the equality holds only when $x^{2}+y^{2}=0$. By applying the local
version of Theorem 4.18 in [24], one can show that the solution satisfies the
inequality (2). ∎

According to Theorem 2, the parameter $\sigma$ plays a pivotal role in
determining the input-to-state stability of the system (3). A smaller $\sigma$
results in a larger value for $a-\gamma_{\mu}$, consequently allowing for a
greater deviation of the initial condition from the origin and permitting
larger perturbations (see Eqs. (4) and (5)). In Fig. 3, we compare the
behavior of system (3) at different parameter regimes. Different values of
$\sigma$ result in distinct regions of attraction within the system. Even when
starting from the same initial condition and being influenced by identical
continuous perturbations, contrasting outcomes emerge. When the region of
attraction is larger, the system remains in the neighborhood of the origin
(see Fig. 3a); if the region of attraction is small, perturbations drive the
system into undesired oscillations (see Fig. 3b).

##  IV Networks of Bistable Oscillators

The brain is a complex network of interacting neuronal populations. The
underlying connections between the populations have a fundamental influence on
the overall dynamics of the entire brain. To study the network’s influence on
epileptic activity, we consider a network of bistable oscillators. To simplify
the analysis, in this paper we restrict our attention to homogeneous
oscillators coupled by a complete network, with dynamics governed by

$\displaystyle\dot{x}_{k}=$ $\displaystyle-\omega
y_{k}+x_{k}\left(\sigma+2ab(x_{k}^{2}+y_{k}^{2})-b(x_{k}^{2}+y_{k}^{2})^{2}\right)$
$\displaystyle\hskip 113.81102pt+\frac{C}{n}\sum_{\ell=1,\ell\neq
k}^{n}x_{\ell}\,,$ (7a) $\displaystyle\dot{y}_{k}=$ $\displaystyle\omega
x_{k}+y_{k}\left(\sigma+2ab(x_{k}^{2}+y_{k}^{2})-b(x_{k}^{2}+y_{k}^{2})^{2}\right)$
$\displaystyle\hskip 113.81102pt+\frac{C}{n}\sum_{\ell=1,\ell\neq
k}^{n}y_{\ell}\,,$ (7b)

where $k=1,2,\dots,n$, and $C\in\mathbb{R}$ is the coupling strength. For
notational simplicity, we denote
$z\coloneq[x_{1},y_{1},x_{2},y_{2},\dots,x_{n},y_{n}]^{\top}$.

The parameter $\sigma$ is assumed to satisfy $-a^{2}/b<\sigma<0$ to ensure
that each oscillator operates in the bistable regime. Note that our choice of
additive coupling follows from recent work, which demonstrates the advantages
of additive coupling over diffusive coupling [25].

One observes that the origin $z=0$ is still the equilibrium point of the
network system (7). Consistent with previous sections, this origin continues
to represent the normal steady state in the absence of epileptic seizures.
Next, we estimate the region of attraction to the origin and study how it is
influenced by the coupling strength in the network.

######  Theorem 3 (Estimate of Region of Attraction).

For the system (7), the following statements hold:
(i) If $|C|<|\sigma|/b$, the set defined by

$\displaystyle\mathcal{Z}\coloneq\left\\{z\in\mathbb{R}^{2n}\colon\|z\|^{2}\leq
a-\sqrt{a^{2}+\sigma/b+|C|}\right\\}$

belongs to the region of attraction for the equilibrium $z=0$. That is,
starting from any initial condition $z(0)\in\mathcal{Z}$, the solution $z(t)$
to (7) converges to $z=0$ asymptotically.
(ii) If $|C|>|\sigma|/b$, the equilibrium $z=0$ becomes unstable.

###### Proof.

Case (i): We construct the proof using LaSalle’s invariance principle. Let
$\nu\coloneq a-\sqrt{a^{2}+\sigma/b+|C|}$ and
$z_{k}\coloneq[x_{k},y_{k}]^{\top}$.

We first show that the set $\mathcal{Z}$ is invariant. To this end, consider a
continuously differentiable function $V(z)=\frac{1}{2}\|z\|^{2}$. Then, the
set $\mathcal{Z}$ can be equivalently rewritten as

$\mathcal{Z}=\left\\{z\in\mathbb{R}^{2n}\colon V(x)\leq
2a-2\sqrt{a^{2}+\sigma/b+|C|}\right\\}.$

The time derivative of $V(x)$ satisfies

$\displaystyle\dot{V}(z)=\sum_{k=1}^{n}\Big{(}(x_{k}^{2}+y_{k}^{2})\big{(}\sigma+2ab(x_{k}^{2}+y_{k}^{2})-b(x_{k}^{2}+y_{k}^{2})^{2}\big{)}$
$\displaystyle\hskip 113.81102pt+\frac{C}{n}\sum_{\ell=1,\ell\neq
k}^{n}(x_{\ell}x_{k}+y_{\ell}y_{k})\Big{)}$
$\displaystyle=\sum_{k=1}^{n}\|z_{k}\|^{2}(\sigma+2ab\|z_{k}\|^{2}-b\|z_{k}\|^{4})+\frac{C}{n}\sum_{k=1}^{n}\sum_{\ell=1,\ell\neq
k}^{n}z_{k}^{\top}z_{\ell}.$

Observe that

$\displaystyle\sum_{k=1}^{n}\sum_{\ell=1,\ell\neq
k}^{n}z_{k}^{\top}z_{\ell}\leq\left\|\sum_{k=1}^{n}z_{k}\right\|^{2}\leq
n\sum_{k=1}^{n}\|z_{k}\|^{2}\,.$

Then, we arrive at

$\displaystyle\dot{V}(z)\leq\sum_{k=1}^{n}\|z_{k}\|^{2}(\sigma+2ab\|z_{k}\|^{2}-b\|z_{k}\|^{4})+|C|\sum_{k=1}^{n}\|z_{k}\|^{2}\,.$

For any $z\in\mathcal{Z}$, it holds that $\|z_{k}\|^{2}\leq\|z\|^{2}\leq\nu$.
Applying the fact in the proof of Theorem 2 that $f(\eta)$ is increasing for
$\eta\in[0,\nu]$, we have

$\displaystyle\dot{V}(z)$
$\displaystyle\leq\sum_{k=1}^{n}\|z_{k}\|^{2}(\sigma+2ab\nu^{2}-b\nu^{4})+|C|\sum_{k=1}^{n}\|z_{k}\|^{2}$
$\displaystyle=\sum_{k=1}^{n}\|z_{k}\|^{2}(\sigma+2ab\nu^{2}-b\nu^{4}+|C|)\leq
0\,.$

where the last inequality is obtained by substituting
$\nu=a-\sqrt{a^{2}+\sigma/b+|C|}$. Then, the set $\mathcal{Z}$ is invariant.

Finally, observe that $\dot{V}(z)=0$ only when $z=0$, implying that starting
in $\mathcal{Z}$, the solution $z(t)$ to (7) converges to $z=0$
asymptotically.

Case (ii): Linearization the system at $z=0$ leads to a Jacobian matrix, of
which eigenvalues have positive real parts. Then, $z=0$ becomes unstable. ∎

We emphasize that Theorem 3 provides an approximation of the region of
attraction to the origin within the network. While the precise region of
attraction may be larger than the one we have identified, finding it
analytically is always challenging. Nevertheless, our estimation reveals an
intriguing trend: a larger coupling strength $C$ (in amplitude) corresponds to
a smaller region of attraction, indicating that an ictal state that captures
seizures is more likely to emerge. It is worth noting, however, that this
observation may be specific to the complete network scenario with homogeneous
individual oscillators. In more complex and heterogeneous networks, the region
of attraction’s dependence on network structure and connection strengths can
be more intricate.

It is also worth noting that networks of bistable units can exhibit a rich
repository of multi-stable dynamics. This multi-stability provides a promising
framework for future exploration into partial seizures, generalized seizures,
as well as the initiation and propagation of localized seizures.

##  V From Phenomenological Model to Real Data

Figure 4: (a) Illustration of networks of bistable oscillators as generative
models. Recorded EEG signals (left) are used to train the network parameters
and the readout function so that the EEG recordings can be reproduced (right).
(b) A network of 64 oscillators reproduces an EEG recording. The trained model
also exhibits seizure-like activity in the presence of higher levels of noise,
indicating a noise-driven transition. Shaded areas show STD.

When recording neural activity using EEG techniques, each electrode measures a
composite signal generated by millions of neurons. Motivated by this
understanding, we propose a framework to reproduce EEG data (see Fig. 4a).
Within this framework, we consider a network of bistable oscillators, each
describing the dynamics of one neuronal ensemble. Additionally, we incorporate
a nonlinear mapping to simulate the intricate journey of neural electrical
activity traversing brain tissues, the skull, and ultimately reaching the EEG
electrodes. Both the network parameters and the readout mapping are trainable.
In Fig. 4b, we train a 64-node network using seizure-free EEG recording (data
source [26]), which accurately reproduces the real data. Remarkably, the
trained model also demonstrates seizure-like activity when subjected to higher
levels of noise.

We emphasize that this framework empowers us to delve into the hidden dynamics
driving epileptic activity. Analyzing the trained network can potentially
provide insights into the mechanisms of seizure onsets, which can help
identify specific objectives for intervention strategies.

##  VI Concluding Remarks

In this paper, we employ a phenomenological model to capture epileptic
phenomena. At certain parameter regimes, this model exhibits bistability,
where a stable equilibrium and a stable limit cycle coexist. The equilibrium
signifies the normal steady-state neural activity, while the limit cycle
represents pathological neural oscillations. The transitions between these
states effectively capture the onset of seizures, which can arise from
internal or external perturbations.

We have analyzed this model, demonstrating that differences in regions of
attraction for both stable states distinguish epileptic brains from healthy
ones. In addition, by considering a complete network, we provide analytical
insights into how different connection strengths impact the region of
attraction of the entire network. This analysis highlights the profound
influence of network structure and connection weights on modifying the
susceptibility of brain regions to noise, ultimately influencing the onset of
seizures. Complementing our theoretical findings, we also discover that
networks consisting of bistable units can accurately reproduce EEG recordings
from epileptic brains through training parameters and network weights. This
approach provides a framework to analyze the trained model and gain deeper
insights into the network pathology of epilepsy.

In addition, our findings provide valuable insights into potential strategies
for controlling epileptic dynamics. For instance, one approach involves
designing closed-loop or event-based controllers to drive the system’s states
back to the region of attraction, effectively preventing them from converging
into oscillatory patterns. Additionally, implementing controllers that expand
the region of attraction can enhance the robustness of normal states against
perturbations, thereby offering another avenue for intervention.

Proof of Theorem 1: We construct the proof using LaSalle’s theorem.

For Case (i), consider the Lyapunov function candidate
$V_{e}(x,y)=\frac{1}{2}(x^{2}+y^{2}).$ The time derivative of $V_{e}(x,y)$ is

$\displaystyle\dot{V}_{e}(x,y)=(x^{2}+y^{2})\left(\sigma+2ab(x^{2}+y^{2})-b(x^{2}+y^{2})^{2}\right),$

which satisfies $\dot{V}_{e}(x,y)\leq 0$ for any ${\rm
col}(x,y)\in\mathcal{A}_{e}$. Observe that the set $\mathcal{A}_{e}$ can be
equivalently written as $\mathcal{A}_{e}=\\{x,y\in\mathbb{R}\colon
V_{e}(x,y)<2{a-2\gamma_{0}}\\}$. Therefore, $\mathcal{A}_{e}$ is an invariant
set. In $\mathcal{A}_{e}$, the equality $\dot{V}_{e}(x,y)=0$ holds only when
${\rm col}(x,y)=0$, which implies that starting at any initial condition in
$\mathcal{A}_{e}$, the solution to the system (1) converges to the origin
asymptotically.

For Case (ii), the set $\mathcal{A}_{\ell}$ can be divided into two sets
$\mathcal{A}_{\ell}^{\prime}$ and $\mathcal{A}_{\ell}^{\prime\prime}$, where
$\mathcal{A}_{\ell}^{\prime}=\\{x,y\in\mathbb{R}\colon
a-\gamma_{0}<x^{2}+y^{2}\leq a+\gamma_{0}\\},\\\
\mathcal{A}_{\ell}^{\prime\prime}=\\{x,y\in\mathbb{R}\colon x^{2}+y^{2}\geq
a+\gamma_{0}\\}$ with $\gamma_{0}$ defined as before. Consider the Lyapunov
function candidates

$\displaystyle V_{\ell}^{\prime}(x,y)=a+\gamma_{0}-(x^{2}+y^{2}),$
$\displaystyle V_{\ell}^{\prime\prime}(x,y)=(x^{2}+y^{2})-a-\gamma_{0}.$

It can be easily shown that $\dot{V}_{\ell}^{\prime}(x,y)\leq 0$ for any ${\rm
col}(x,y)\in\mathcal{A}_{\ell}^{\prime}$ and
$\dot{V}_{\ell}^{\prime\prime}(x,y)\leq 0$ for any ${\rm
col}(x,y)\in\mathcal{A}_{\ell}^{\prime\prime}$. As
$\mathcal{A}_{\ell}^{\prime}$ and $\mathcal{A}_{\ell}^{\prime\prime}$ can be
written as $\mathcal{A}_{\ell}^{\prime}=\\{x,y\in\mathbb{R}\colon 0\leq
V_{\ell}^{\prime}(x,y)<2\gamma_{0}\\}$ and
$\mathcal{A}_{\ell}^{\prime}=\\{x,y\in\mathbb{R}\colon
V_{\ell}^{\prime}(x,y)\geq 0\\}$, one can observe that both
$\mathcal{A}_{\ell}^{\prime}$ and $\mathcal{A}_{\ell}^{\prime\prime}$ are
invariant. Therefore, the set
$\mathcal{A}_{\ell}=\mathcal{A}_{\ell}^{\prime}\cup\mathcal{A}_{\ell}^{\prime\prime}$
is invariant, and $\dot{V}(x,y)=0$ only when ${\rm
col}(x,y)\in\mathcal{C}_{1}$. Starting from any initial condition in
$\mathcal{A}_{\ell}$, the solution converges to $\mathcal{C}_{1}$
asymptotically. $\square$

## References

  * [1] C. Begley, R. G. Wagner, A. Abraham, E. Beghi, C. Newton, C.-S. Kwon, D. Labiner, and A. S. Winkler, “The global cost of epilepsy: a systematic review and extrapolation,” _Epilepsia_ , vol. 63, no. 4, pp. 892–903, 2022\.
  * [2] P. Jiruska, M. De Curtis, J. G. Jefferys, C. A. Schevon, S. J. Schiff, and K. Schindler, “Synchronization and desynchronization in epilepsy: Controversies and hypotheses,” _The Journal of Physiology_ , vol. 591, no. 4, pp. 787–797, 2013.
  * [3] W. Löscher, H. Potschka, S. M. Sisodiya, and A. Vezzani, “Drug resistance in epilepsy: clinical impact, potential mechanisms, and new innovative treatment options,” _Pharmacological Reviews_ , vol. 72, no. 3, pp. 606–638, 2020.
  * [4] S. L. Moshé, E. Perucca, P. Ryvlin, and T. Tomson, “Epilepsy: New advances,” _The Lancet_ , vol. 385, no. 9971, pp. 884–898, 2015.
  * [5] E. M. Goldberg and D. A. Coulter, “Mechanisms of epileptogenesis: a convergence on neural circuit dysfunction,” _Nature Reviews Neuroscience_ , vol. 14, no. 5, pp. 337–349, 2013.
  * [6] R. A. Stefanescu, R. Shivakeshavan, and S. S. Talathi, “Computational models of epilepsy,” _Seizure_ , vol. 21, no. 10, pp. 748–759, 2012.
  * [7] F. Wendling, P. Benquet, F. Bartolomei, and V. Jirsa, “Computational models of epileptiform activity,” _Journal of Neuroscience Methods_ , vol. 260, pp. 233–251, 2016.
  * [8] F. L. Da Silva, W. Blanes, S. N. Kalitzin, J. Parra, P. Suffczynski, and D. N. Velis, “Epilepsies as dynamical diseases of brain systems: Basic models of the transition between normal and epileptic activity,” _Epilepsia_ , vol. 44, pp. 72–83, 2003.
  * [9] F. H. L. da Silva, W. Blanes, S. N. Kalitzin, J. Parra, P. Suffczynski, and D. N. Velis, “Dynamical diseases of brain systems: Different routes to epileptic seizures,” _IEEE Transactions on Biomedical Engineering_ , vol. 50, no. 5, pp. 540–548, 2003.
  * [10] E. Nozari and J. Cortés, “Oscillations and coupling in interconnections of two-dimensional brain networks,” in _American Control Conference_ , 2019, pp. 193–198.
  * [11] E. Nozari, R. Planas, and J. Cortés, “Structural characterization of oscillations in brain networks with rate dynamics,” _Automatica_ , vol. 146, p. 110653, 2022.
  * [12] M. McCreesh, T. Menara, and J. Cortés, “Sufficient conditions for oscillations in competitive linear-threshold brain networks,” _IEEE Control Systems Letters_ , 2023.
  * [13] F. Celi, A. Allibhoy, F. Pasqualetti, and J. Cortés, “Linear-threshold dynamics for the study of epileptic events,” _IEEE Control Systems Letters_ , vol. 5, no. 4, pp. 1405–1410, 2020.
  * [14] A. Allibhoy, F. Celi, F. Pasqualetti, and J. Cortés, “Optimal network interventions to control the spreading of oscillations,” _IEEE Open Journal of Control Systems_ , vol. 1, pp. 141–151, 2022.
  * [15] Y. Qin, D. S. Bassett, and F. Pasqualetti, “Vibrational control of cluster synchronization: Connections with deep brain stimulation,” in _61st IEEE Conference on Decision and Control_ , 2022, pp. 655–661.
  * [16] Y. Qin, A. M. Nobili, D. S. Bassett, and F. Pasqualetti, “Vibrational stabilization of cluster synchronization in oscillator networks,” _IEEE Open Journal of Control Systems_ , vol. 2, pp. 439–453, 2023.
  * [17] A. M. Nobili, Y. Qin, C. A. Avizzano, D. S. Bassett, and F. Pasqualetti, “Vibrational stabilization of complex network systems,” in _American Control Conference_ , 2023, pp. 1980–1985.
  * [18] E. A. Reed, G. Ramos, P. Bogdan, and S. Pequito, “Mitigating epilepsy by stabilizing linear fractional-order systems,” in _American Control Conference_ , 2023, pp. 2228–2233.
  * [19] F. Fröhlich, T. J. Sejnowski, and M. Bazhenov, “Network bistability mediates spontaneous transitions between normal and pathological brain states,” _Journal of Neuroscience_ , vol. 30, no. 32, pp. 10 734–10 743, 2010.
  * [20] G. Baier, M. Goodfellow, P. N. Taylor, Y. Wang, and D. J. Garry, “The importance of modeling epileptic seizure dynamics as spatio-temporal patterns,” _Frontiers in Physiology_ , vol. 3, p. 281, 2012.
  * [21] S. N. Kalitzin, D. N. Velis, and F. H. L. da Silva, “Stimulation-based anticipation and control of state transitions in the epileptic brain,” _Epilepsy & Behavior_, vol. 17, no. 3, pp. 310–323, 2010.
  * [22] O. Benjamin, T. H. Fitzgerald, P. Ashwin, and et al., “A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy,” _The Journal of Mathematical Neuroscience_ , vol. 2, pp. 1–30, 2012.
  * [23] S. H. Strogatz, _Nonlinear Dynamics and Chaos_. Boca Raton: CRC press, 2018.
  * [24] H. K. Khalil, _Nonlinear Systems_. Prentice Hall, 2002.
  * [25] M. A. Lopes, K. Hamandi, J. Zhang, and J. L. Creaser, “The role of additive and diffusive coupling on the dynamics of neural populations,” _Scientific Reports_ , vol. 13, no. 1, p. 4115, 2023.
  * [26] M. Stead, M. Bower, B. H. Brinkmann, and et al., “Microseizures and the spatiotemporal scales of human partial epilepsy,” _Brain_ , vol. 133, no. 9, pp. 2789–2797, 2010.