# APPLICATION OF ADAPTED-BUBBLES TO THE HELMHOLTZ EQUATION WITH LARGE WAVE
NUMBERS IN 2D

Adem Kaya<EMAIL_ADDRESS>Institut für Mathematik, Universität Potsdam
Karl-Liebknecht-Str. 24-25 14476 Potsdam/Golm Germany

###### Abstract

An adapted bubble approach which is a modifiation of the residual-free bubbles
(RFB) method, is proposed for the Helmhotz problem in 2D. A new two-level
finite element method is introduced for the approximations of the bubble
functions. Unlike the other equations such as the advection-diffusion
equation, RFB method when applied to the Helmholtz equation, does not depend
on another stabilized method to obtain approximations to the solutions of the
sub-problems. Adapted bubbles (AB) are obtained by a simple modification of
the sub-problems. This modification increases the accuracy of the numerical
solution impressively. The AB method is able to solve the Helmholtz equation
efficiently in 2D up to $ch=3.5$ where $c$ is the wave number and $h$ is the
mesh size. We provide analysis to show how the AB method mitigates the
pollution error.

###### keywords:

Helmholtz equation, adapted bubbles, residual-free bubbles, two-level finite
element

##  1  Introduction

Enriching linear finite element space with residual-free bubble functions is a
general framework for the discretizations of the problems[RFBhelmholtz,
choosingbuubles, nesliturkcdr, nesliturktez, unlocking, masslumping]. These
functions strongly satisfy the original differential equations and hence
obtaining the bubble functions is generally as difficult as solving the
original problem such as the convection-diffusion equation [choosingbuubles].
Unlike it was stated in [RFBhelmholtztwolevel], we will show that this is not
the case for the Helmholtz problem. Obtaining the bubble functions on
triangular elements is easier than solving the original problem. The standard
Galerkin finite element method can be used with a coarse mesh to obtain
efficient approximations to the bubble functions.

The residual-free bubbles method produces the exact solution of linear
dierential equations in the one-dimensional case. However, the method in
higher-dimensions is approximate and as we will show for the Helmholtz problem
in this article, its contribution to the stabilization of the standard
Galerkin method is very poor. We modify the residual-free bubbles (RFB) method
in 2D by multiplying the right-hand side of the bubble equations with a
constant. This operation impressively increase the accuracy. The new bubbles
are no more residual-free and we call them  adapted bubbles (AB). We provide
the optimal values of the constants for the triangular and rectangular
elements separately. We apply a two-level finite element method using the
standard Galerkin finite element method to get approximations to the bubble
functions.

We provide analysis to show how the AB method mitigates the pollution error.
To this end. we approximate the bubble functions with piece-wise defined
linear functions so-called pseudo-bubbles. The analysis give rise to a fourth
order finite difference scheme with seven-point stencil for plane waves. It is
perfectly applicable in polygonal and triangular domains. We use this method
to do comparison with the AB method.

Standard discretizations when applied to the Helmholtz problem suffer from the
pollution effect when the wave number is large [Babuska1997]. Moreover
standard iterative solvers are ineffective in obtaining the solutions of the
discrete Helmholtz equation [Ernst2012]. There is a great effort in literature
to overcome these difficulties. Among the discretization techniques, there are
finite difference [Singer1998, Feng2011], finite element [Ihlenburg1995,
Strouboulis2006, Alvarez2006], discontinuous Galerkin [paul, Haijun], virtual
element [Perugia2016], and boundary element methods [Kirkup2019]. At the same
time, there is a great effort to develop efficient preconditioners, such as
multigrid [Elman2002, Erlangga2006] and domain decomposition methods
[Gander2003, Stolk2013]. The AB method proposed in this article, does not
suffer from the pollution effect for very large wave numbers. It is by far
superior than the fourth order accurate scheme proposed here.

The rest of this paper is organized as follows. We review the RFB method for
the Helmholtz equation in Section 2. We explain how to implement two-level
finite element method in 1D and provide analysis to show the contribution of
the bubble functions in reducing the pollution error in Section 3. Section 4
is devoted for the analysis of the RFB method in 2D. We propose the AB method
for triangular elements in Section 5. Several numerical experiments are
provided in Section LABEL:sec:section5. The AB method is considered with
rectangular elements in Section LABEL:sec:section6. We finish with concluding
remarks in Section LABEL:sec:section7.

##  2 Residual-free bubbles method (RFB) for the Helmholtz equation

We start with considering the Helmholtz problem in 1D with Dirichlet boundary
conditions on unit interval.

$\displaystyle\Bigg{\\{}\begin{array}[]{ll}-u^{\prime\prime}-c^{2}u(x)=f(x),\qquad
x\in I\\\ u(0)=0,\qquad u(1)=1,\end{array}$ (3)

where we assume that the wave number $c$ is constant. RFB method can be
summarized as follows. Let’s start with recalling abstract variational
formulation of (3): Find $u\in H^{1}(I)$ such that

$a(u,v)=(f,v),\qquad\forall v\in H^{1}(I),$

where

$a(u,v)=\int_{I}u^{\prime}v^{\prime}dx-c^{2}\int_{I}uvdx$

and

$(f,v)=\int_{I}fvdx.$

Define $V_{h}\subset H^{1}(I)$ as a finite-dimensional space. Then the
Galerkin finite element method reads: Find $u_{h}\in V_{h}$ such that

$a(u_{h},v_{h})=(f,v_{h}),\qquad\forall v_{h}\in V_{h}.$

We now decompose the space $V_{h}$ as $V_{h}=V_{L}\bigoplus V_{B}$, where
$V_{L}$ is the space of continuous piecewise linear polynomials and
$V_{B}=\bigoplus_{K}B_{K}$ with $B_{K}=H_{0}^{1}(K)$. From this decomposition,
every $v_{h}\in V_{h}$ can be written in the form $v_{h}=v_{L}+v_{B}$, where
$v_{L}\in V_{L}$ and $v_{B}\in V_{B}$. Bubble component $u_{B}$ of $u_{h}$
satisfy the original differential equation in an element $K$ strongly, i.e.

$\mathcal{L}u_{B}=-\mathcal{L}u_{L}+f\qquad\text{in}\qquad K,$ (4)

subject to boundary condition,

$u_{B}=0\qquad\text{on}\qquad\partial K.$ (5)

Since the support of bubble $u_{B}$ is contained within the element $K$, we
can make a static condensation for the bubble part, getting directly the
$V_{L}$\- projection $u_{L}$ of the solution $u_{h}$ [RFBhelmholtz]. This can
be done as follows. Using $V_{h}=V_{L}\bigoplus V_{B}$, the finite element
approximation reads: Find $u_{h}=u_{L}+u_{B}$ in $V_{h}$ such that

$a(u_{L},v_{L})+a(u_{B},v_{L})=(f,v_{L}),\qquad\forall v_{L}\in V_{L}.$ (6)

##  3 A two level finite element method

In order to find bubble part $u_{B}$ of the solution, we need to solve
(4)-(5). The problems defined by equations (4)-(5) is addressed by solving
instead

$\displaystyle\left\\{\begin{array}[]{ll}-\varphi_{i}^{\prime\prime}-c^{2}\varphi_{i}=c^{2}\psi_{i}\quad\text{in}\quad
K,\quad(i=1,...,n_{en})\\\ \varphi_{i}=0\quad\text{on}\quad\partial
K,\end{array}\right.$ (9)

and

$\displaystyle\left\\{\begin{array}[]{ll}-\varphi_{f}^{\prime\prime}-c^{2}\varphi_{f}=f\quad\text{in}\quad
K,\\\ \varphi_{f}=0\quad\text{on}\quad\partial K.\end{array}\right.$ (12)

where $n_{en}$ is the number of element nodes. Thus if

$u_{L}=\sum\limits_{i}d_{i}^{K}\psi_{i},$ (13)

then

$u_{B}=\sum\limits_{i}d_{i}^{K}\varphi_{i}+\varphi_{f}.$ (14)

Substituting (13) and (14) into (6), we get the matrix formulation

$\sum_{K}\sum\limits_{i}^{n_{en}}d_{i}^{K}\left((\psi^{{}^{\prime}}_{i},\psi_{j}^{{}^{\prime}})-c^{2}(\psi_{i},\psi_{j})-c^{2}(\varphi_{i},\psi_{j}))=(f,\psi_{j})\right)$
(15)

at the global level where $d_{i}$ are the finite element approximations to the
solutions at the nodes.

Numerical solution of the bubble problems (9) and (12) generally requires
using a nonstandard method such as for the case of the advection-diffusion
equation. This makes the RFB method dependent on another stabilized method
when applying it as a two-level finite element method. In
[RFBhelmholtztwolevel], the Galerkin-least-squares method (GLS) was used to
get approximations to the bubble functions in solving the the Helmholtz
equation. Although this is true for the advection-diffusion equation, there is
no need to use a nonstandard method to get approximations to the bubble
functions when the Helmholtz problem is under consideration. We explain this
this fact in 1D. It is well known that standard discretizations are effective
up to $ch=0.6$. Suppose that we have a discretization of the domain such that
$ch=0.6$. Even if we use 3 nodes on the sub-domain (element), $ch_{e}$
($h_{e}$ is mesh size on the sub-domain) becomes less than $0.6$. If $ch=3$ on
the global mesh, then using 11 nodes for the sub-problems makes $ch_{e}=0.3$.
More precisely, it is always true that $ch_{e}<ch$.

It is true that GLS computation is known to incur at most marginal increase in
computational cost over the standard Galerkin method. However, GLS for the
sub-problems may lead to misinterpretations related to the bubble functions.

Another way of obtaining the bubble function is to use separation of variables
when rectangular elements are used [RFBhelmholtz]. However, this gives rise to
a series solution of the bubble function for which it must be truncated. For a
good accuracy, 200 terms are used in [RFBhelmholtz] which is computationally
not so effective. Another drawback is that this approach is limited to the
rectangular elements.

###  3.1 Analysis of the pollution effect of the sub-problems

The inequality $ch_{e}<ch$ is an indication that the sub-problems are easier
to solve; however, we must analyse the pollution effect for the sub-problems
for large wave numbers. It is well known that the condition $c^{2}h<1$ is
sufficient to guarantee that the error of the Galerkin solution is of the same
magnitude as the error of the best approximation [Babuska1997]. This condition
is necessary when the size the domain is fixed for increasing $c$. More
precisely, the exact solution is very oscillatory. In our case, the exact
solutions of the sub-problems are not oscillatory. When $ch<\pi$, the exact
solutions of the sub-problems are always in the form of a half wave as the
homogenous Dirichlet boundary condtion is applied everywhere on the boundary.
When $\pi<ch<2\pi$, the exact solutions of the sub-problems are always in the
form of a single wave. In this regime, the standard Galerkin method is
pollution free for the sub-problems for any wave number. Note that, in
simulations, 12 nodes per wave is generally chosen which correponds to
$ch\approx 0.62$.

We use the standard Galerkin finite element method with piecewise linear basis
functions to approximate the bubble functions. Note that the bubble problems
(9) and (12) can be solved independently and hence parallel processors can be
used to carry out these computations efficiently. When uniform meshes are used
and the right hand side function $f$ is constant, construction of the system
matrix is as cheap as construction of the system matrix of the standard
Galerkin finite element method.

###  3.2 Shape of the bubble functions and pseudo-bubbles

We have shown that the RFB method is not dependent on another stabilized
method to get approximations to the bubble functions when a two-level finite
element method is used due to the non-oscillatory behavior of the exact
solutions of the sub-problems. This non-oscillatory behavior of the exact
solutions opens a gateway to approximate these bubble functions with
piecewise-defined linear simple functions. These approximations are called
pseudo-bubbles and constructed considering the shape of the bubble functions.
Pseudo-bubbles were applied to the advection-diffusion-reaction equation in
[Asendur, Asendur2].

Here, we consider the case $ch<\pi$ for which the bubble functions are in the
form of a half wave. We present the bubble functions $\varphi_{1,2}$ in Figure
1 for $c=60,300$ when $h=0.01$. Efficient yet cheap approximations to these
bubble functions with piecewise-defined linear functions are given in Figure
2. While on the left in Figure 2, two pseudo-bubbles are used, it is possible
to approximate $\varphi_{1,2}$ with a single pseudo-bubble. The humps of the
bubble functions $\varphi_{1,2}$ come closer as $ch$ increases (see Figure 1.
Thus, for larger $ch$, we can derive more efficient approximations to the
bubble functions when a single pseudo-bubble is used. The advantage of using a
single pseudo-bubble is that the maximum of the pseudo-bubble occurs in the
middle of the element. Applying the minimization technique applied in
[Asendur], one can find the optimal heights and locations of the peaks of the
pseudo-bubbles. When a single pseudo-bubble is used, one can easily calculate
integrals in the finite element formulation. This will be important in
modifying the RFB method in 2D.

Figure 1: Bubble functions for $c=60$ (left) and $c=300$ (right).

Figure 2: Two different approaches to approximate the bubble functions.

###  3.3 Analysis of the pseudo-bubbles

In order to see how the residual-free bubbles method overcome the pollution
effect, we first consider linear finite element method for (3) when $f(x)=0$,
$u(0)=0$ and $u(1)=\sin(c)$. The exact solution is $u(x)=\sin(cx)$. We obtain
the truncation error and see how the error deteriorates as $c$ increases for
fixed $ch$ which is simply the pollution effect. To this end, we jump to the
finite difference equivalance of the linear finite method. Let $U_{j}$
represents the numerical solution and choose $n$ equally distributed nodes for
which $h=1/(n-1)$. Taking the integrals in linear finite element formulation
and scaling by $h$ gives

$-\frac{U_{j+1}-2U_{j}+U_{j-1}}{h^{2}}-c^{2}\frac{U_{j+1}+4U_{j}+U_{j-1}}{6}=0,\qquad
j=2,...,n-2.$ (16)

From the boundary conditions,

$U_{1}=0,\qquad U_{n}=\sin(c).$

Using the Taylor series expansion, we get the truncation error for (16).

$\tau(x)=-\frac{c^{2}h^{2}}{6}u^{\prime\prime}-\left(\frac{c^{2}h^{4}}{72}+\frac{h^{2}}{12}\right)u^{(4)}-\frac{h^{4}}{360}u^{(6)}+\mathcal{O}(h^{6}).$
(17)

The pollution effect can not be seen from (17). To see it, we substitute the
exact solution $u(x)=\sin(cx)$ into (17).

$\tau(x)=\frac{c^{4}h^{2}}{6}\sin(cx)-\left(\frac{c^{6}h^{4}}{72}+\frac{c^{4}h^{2}}{12}\right)\sin(cx)+\frac{c^{6}h^{4}}{360}\sin(cx)+\mathcal{O}(c^{8}h^{6}).$

Rearranging the above equation gives

$\tau(x)=\sin(cx)\left(\frac{c^{4}h^{2}}{12}-\frac{c^{6}h^{4}}{90}+\mathcal{O}(c^{8}h^{6})\right).$
(18)

When the exact solution is oscillatory, that is, $c$ is large, the term
$c^{4}h^{2}/12$ in (18) becomes large, that is, $\tau(x)$ is large, even if
$ch=\text{constant}$ is small. This is called the pollution effect. When $ch$
is sufficiently small, there is no phase error for the Dirichlet problem when
$c$ is large, however when a Neumann or Robin boundary condition is used,
phase error is also observed.

The simplest way to mitigate this pollustion effect is to choose $c^{2}h$
sufficiently small. However, this requires intractable matrices in higher
dimensions. The general idea in literature is to decrease the effect of the
first few terms in (18) so that tractable matrix sizes can be obtained. For
example, using higher order accurate methods of finite difference or higher
order polynomials finite element may allow to eliminate the first few terms.
If the first term can be eliminated, then the requirement to mitigate the
pollution effect reduces to $c^{3/2}h$ being sufficiently small. However,
higher order methods generally use more points and this increases the nonzero
entries of the matrices.

In order to get a deeper insight of the working principle of the residual-free
bubbles method in mitigating the pollution effect, we consider the pseudo-
bubbles in Figure 2 on the right. This choice allows us to take the integrals
containing the bubble functions, explicitly.

Figure 3: Basis functions employed in the approximation of bubble functions.

We can define the pseudo-bubbles $b_{1,2}$ using the basis functions
$b_{1,2}^{T}$ represented in Figure 3 and heights of $b_{1,2}$, i.e.,
$\alpha_{1,2}$. More precisely,

$b_{1,2}=\alpha_{1,2}b_{1,2}^{T},$ (19)

where $\alpha_{1}=\alpha_{2}$. Applying the technique proposed in [Asendur]
(set $\xi=h/2,\epsilon=1,\sigma=-c^{2}$ in equation (13) in [Asendur]) gives

$\alpha_{1}=\frac{3c^{2}h^{2}}{4(12-c^{2}h^{2})}.$ (20)

Taking the integrals in (6) making use of (19) gives the finite difference
formula

$-\frac{U_{j+1}-2U_{j}+U_{j-1}}{h^{2}}-c^{2}\frac{U_{j+1}+4U_{j}+U_{j-1}}{6}-\alpha_{1}c^{2}\frac{U_{j+1}+2U_{j}+U_{j-1}}{4}=0,\qquad
j=2,...,n-2.$ (21)

Using Taylor expansions of $u(x\pm h)$, definition of $\alpha_{1}$ given in
(20) and the exact solution $u(x)=\sin(cx)$, we obtain the truncation error

$\tau(x)=\sin(cx)\left(\frac{c^{4}h^{2}}{12}-\frac{c^{6}h^{4}}{90}+\mathcal{O}(c^{8}h^{6})\right)+\sin(cx)\left(-\frac{3c^{4}h^{2}}{4(12-c^{2}h^{2})}+\frac{3c^{6}h^{4}}{16(12-c^{2}h^{2})}+\mathcal{O}(c^{8}h^{6})\right).$

Rearranging the right-hand side of the above formula we end up with

$\tau(x)=\sin(cx)\left(c^{4}h^{2}\left(\frac{1}{12}-\frac{3}{4(12-c^{2}h^{2})}\right)+c^{6}h^{4}\left(-\frac{1}{90}+\frac{3}{16(12-c^{2}h^{2})}\right)+\mathcal{O}(c^{8}h^{6})\right).$
(22)

While the coefficient of $c^{4}h^{2}$ in (18) is $1/12$, the coefficient of
$c^{6}h^{4}$ is $1/90$ in magnitude. When $ch<\sqrt{15/2}\approx 2.73$, the
coefficient of $c^{4}h^{2}$ in (22) is smaller in magnitude and it is close to
zero when $ch\approx\sqrt{3}$. Moreover, when $ch<\sqrt{57/16}\approx 1.88$,
the coefficient of $c^{6}h^{4}$ in (22) becomes smaller in magnitude. The
approximate bubbles shows how the pollution effect is reduced. It is known
that RFB method for 1D linear equations is exact [RFBhelmholtztwolevel]. This
means that it automatically makes the coefficient of all powers
$c^{n+2}h^{n}$, $n=2,3,...$, zero. A good approximation to the residual-free
bubbles significantly reduce the pollution effect. The RFB method is
approximate in 2D. It is well known that the contributions of the residual-
free bubble functions to the stabilization of the Galerkin method is very
poor. The observations we made here will be helpful to further increase the
accuracy of the method in 2D. We will modify the sub-problems in 2D and use
adapted bubbles to further increase the accuracy of the bubble approach.

##  4 The RFB method in 2D

We have shown that the RFB method is able to solve the Helmholtz problem in 1D
cheaply and efficiently for very large wave numbers. As it was stated in
[RFBhelmholtztwolevel], RFB method is not as efficient in 2D as in 1D. To show
this fact, we consider the following problem on an equilateral triangular
shaped domain with vertices $(0,0)$, $(0,1)$ and $(0.5,\sqrt{3}/2)$.

$\displaystyle\left\\{\begin{array}[]{ll}-\Delta
u-c^{2}u=0,\quad\text{in}\quad\Omega,\\\
u(x,y)=\sin(cx\sin(\theta)+cy\cos(\theta)),\quad\text{on}\quad\partial\Omega_{D},\end{array}\right.$
(25)

where the exact solution is $u(x,y)=\sin(cx\sin(\theta)+cy\cos(\theta))$. We
use equilateral triangular elements with linear basis functions. We solve the
following equations on element level to get the RFB functions.

$\displaystyle\left\\{\begin{array}[]{ll}-\Delta\varphi_{i}-c^{2}\varphi_{i}=c^{2}\psi_{i}\quad\text{in}\quad
K,\quad(i=1,...,n_{en})\\\ \varphi_{i}=0\quad\text{on}\quad\partial
K,\end{array}\right.$ (28)

and

$\displaystyle\left\\{\begin{array}[]{ll}-\Delta\varphi_{f}-c^{2}\varphi_{f}=f\quad\text{in}\quad
K,\\\ \varphi_{f}=0\quad\text{on}\quad\partial K,\end{array}\right.$ (31)

where $n_{en}=3$. $\varphi_{i}$ and $\psi_{i}$, $(i=1,2,3)$ are the RFB and
the linear basis functions, respectively. Linear finite element method with a
coarse mesh can be used to obtain efficient approximations to the bubble
functions. To do some analyses, we approximate the RFB functions with
piecewise-defined linear functions with the maximum at the centroid of the
element. Let $b_{1,2,3}=\alpha_{2}b^{T}$ be the approximation to the bubble
functions where $b^{T}$ is the linear basis bubble function that assumes zero
at the vertices of the element and one at the centroid of the element.
Applying the same procedure we applied in 1D (see [Asendur2] for more
details), gives

$\alpha_{2}^{2}\int_{K}\nabla b^{T}\nabla
b^{T}dS-\alpha_{2}^{2}c^{2}\int_{K}b^{T}b^{T}dS=\alpha_{2}\int_{K}\psi_{i}b^{T}dS,\qquad
i=1,2,3.$ (32)

Solving the above equation for $\alpha_{2}$ and calculating the integrals for
$\psi_{1}$ gives

$\alpha_{2}=\frac{2c^{2}h^{2}}{3(72-c^{2}h^{2})}.$ (33)

Considering 6 adjacent elements as shown in Figure 4, the RFB method is
equivalent to the following finite difference scheme.

Figure 4: Equilateral triangles surrounding a node.

$\displaystyle\begin{array}[]{ll}\frac{6U_{i,j}-U_{i-1,j}-U_{i-1,j+1}-U_{i+1,j+1}-U_{i+1,j}-U_{i+1,j-1}-U_{i-1,j-1}}{\sqrt{3}h^{2}}-c^{2}\frac{6U_{i,j}+U_{i-1,j}+U_{i-1,j+1}+U_{i+1,j+1}+U_{i+1,j}+U_{i+1,j-1}+U_{i-1,j-1}}{8\sqrt{3}}\\\
-\alpha_{2}c^{2}\frac{3U_{i,j}+U_{i-1,j}+U_{i-1,j+1}+U_{i+1,j+1}+U_{i+1,j}+U_{i+1,j-1}+U_{i-1,j-1}}{6\sqrt{3}},\qquad
i,j=1,...,N_{\text{int}},\end{array}$ (36)

where $U_{i,j}\approx u(x,y)$, $U_{i-1,j}\approx u(x-h,y)$,
$U_{i-1,j+1}\approx u(x-\frac{h}{2},y+\frac{\sqrt{3}}{2}h)$,
$U_{i+1,j+1}\approx u(x+\frac{h}{2},y+\frac{\sqrt{3}}{2}h)$, $U_{i+1,j}\approx
u(x+h,y)$, $U_{i+1,j-1}\approx u(x+\frac{h}{2},y-\frac{\sqrt{3}}{2}h)$ and
$U_{i-1,j-1}\approx u(x-\frac{h}{2},y-\frac{\sqrt{3}}{2}h)$. Note that when
$\alpha_{2}=0$, (36) is equivalent to the linear finite element method with
equilateral triangular element.

To analyse the RFB method, we substitute the Taylor expansions of the exact
solution at the grid points. The derivation of the truncation error is given
in (42)-(51). In our analysis, we will examine the coefficients of
$c^{4}h^{2}$ and $c^{6}h^{4}$ in (51), that is, $C_{2}$ and $C_{1}$ in (52).

$\displaystyle\begin{array}[]{ll}\tau(x,y)=\frac{6u(x,y)-u(x-h,y)-u(x-\frac{h}{2},y+\frac{\sqrt{3}}{2}h)-u(x+\frac{h}{2},y+\frac{\sqrt{3}}{2}h)-u(x+h,y)-u(x+\frac{h}{2},y-\frac{\sqrt{3}}{2}h)-u(x-\frac{h}{2},y-\frac{\sqrt{3}}{2}h)}{\sqrt{3}h^{2}}-c^{2}\frac{6u(x,y)+u(x-h,y)+u(x-\frac{h}{2},y+\frac{\sqrt{3}}{2}h)+u(x+\frac{h}{2},y+\frac{\sqrt{3}}{2}h)+u(x+h,y)+u(x+\frac{h}{2},y-\frac{\sqrt{3}}{2}h)+u(x-\frac{h}{2},y-\frac{\sqrt{3}}{2}h)}{8\sqrt{3}}\\\
-c^{2}\alpha_{2}\frac{3u(x,y)+u(x-h,y)+u(x-\frac{h}{2},y+\frac{\sqrt{3}}{2}h)+u(x+\frac{h}{2},y+\frac{\sqrt{3}}{2}h)+u(x+h,y)+u(x+\frac{h}{2},y-\frac{\sqrt{3}}{2}h)+u(x-\frac{h}{2},y-\frac{\sqrt{3}}{2}h)}{6\sqrt{3}}=\\\
-\frac{189h^{4}u^{(0,6)}(x,y)+945h^{4}u^{(2,4)}(x,y)+315h^{4}u^{(4,2)}(x,y)+231h^{4}u^{(6,0)}(x,y)+7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)+120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)}{80640\sqrt{3}}\\\
-\frac{c^{2}\left(h^{2}\left(7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)+120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)\right)+967680u(x,y)\right)}{645120\sqrt{3}}\\\
-\frac{\text{$\alpha$2}c^{2}\left(h^{2}\left(7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)+120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)\right)+725760u(x,y)\right)}{483840\sqrt{3}}+\mathcal{O}(h^{5}).\end{array}$
(42)

Since

$-\frac{120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)}{80640\sqrt{3}}-\frac{967680u(x,y)}{645120\sqrt{3}}=0,$

(42) becomes

$\displaystyle\begin{array}[]{ll}\tau(x,y)=-\frac{189h^{4}u^{(0,6)}(x,y)+945h^{4}u^{(2,4)}(x,y)+315h^{4}u^{(4,2)}(x,y)+231h^{4}u^{(6,0)}(x,y)+7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)}{80640\sqrt{3}}\\\
-\frac{c^{2}\left(h^{2}\left(7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)+120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)\right)\right)}{645120\sqrt{3}}\\\
-\frac{\text{$\alpha$2}c^{2}\left(h^{2}\left(7560h^{2}u^{(0,4)}(x,y)+15120h^{2}u^{(2,2)}(x,y)+7560h^{2}u^{(4,0)}(x,y)+120960u^{(0,2)}(x,y)+120960u^{(2,0)}(x,y)\right)+725760u(x,y)\right)}{483840\sqrt{3}}+\mathcal{O}(h^{5}).\end{array}$
(46)

Substituting the exact solution $u(x,y)=\sin(c\cos(\theta)x+c\sin(\theta)y)$
(note that more general solutions can be chosen) into (46) gives

$\displaystyle\begin{array}[]{ll}\tau(x,y)=-\frac{1}{322560\sqrt{3}}(c^{2}(35c^{5}h^{5}\cos(cx\cos(\theta)+cy\sin(\theta)+\theta)+21c^{5}h^{5}\cos(cx\cos(\theta)+cy\sin(\theta)+3\theta)+7c^{5}h^{5}\cos(cx\cos(\theta)+cy\sin(\theta)+5\theta)+\\\
c^{5}h^{5}\cos(cx\cos(\theta)+cy\sin(\theta)+7\theta)+35c^{5}h^{5}\cos(-cx\cos(\theta)-cy\sin(\theta)+\theta)+21c^{5}h^{5}\cos(-cx\cos(\theta)-cy\sin(\theta)+3\theta)+7c^{5}h^{5}\cos(-cx\cos(\theta)-cy\sin(\theta)+5\theta)+\\\
c^{5}h^{5}\cos(-cx\cos(\theta)-cy\sin(\theta)+7\theta)+5040\text{$\alpha$2}c^{4}h^{4}\sin(cx\cos(\theta)+cy\sin(\theta))+2940c^{4}h^{4}\sin(cx\cos(\theta)+cy\sin(\theta))-42c^{4}h^{4}\sin(cx\cos(\theta)+cy\sin(\theta)+6\theta)+\\\
42c^{4}h^{4}\sin(-cx\cos(\theta)-cy\sin(\theta)+6\theta)-80640\text{$\alpha$2}c^{2}h^{2}\sin(cx\cos(\theta)+cy\sin(\theta))-30240c^{2}h^{2}\sin(cx\cos(\theta)+cy\sin(\theta))+483840\text{$\alpha$2}\sin(cx\cos(\theta)+cy\sin(\theta))+\mathcal{O}(h^{5}).\end{array}$
(51)

The last 6 terms in (51) can be written in the form

$\sin(cx\cos(\theta)+cy\sin(\theta))\left(c^{6}h^{4}\underbrace{\frac{-2940+\frac{80640\alpha
2}{c^{2}}+84\cos(6\theta)}{322560\sqrt{3}}}_{C_{2}}+c^{4}h^{2}\underbrace{\frac{30240-\frac{483840\alpha
2}{c^{2}h^{2}}}{322560\sqrt{3}}}_{C_{1}}\right).$ (52)

Figure 5 shows the graph of $C_{1}$ and $C_{2}$ for $\alpha_{2}=0$ (standard
Galerkin) and for $\alpha_{2}$ in (33) (pseudo-RFB) . We set $\theta=\pi/3$ to
plot the graph of $C_{2}$. The slight decreases in $C_{1}$ and $C_{2}$ in
magnitude for $0<ch<3$, explains why the RFB method is not effective in 2D.

Figure 5: Comparison of the linear Galekin method and the pseudo-RFB method
for the coefficients $C_{1}$ (left) and $C_{2}$ (right).

One way to improve the accuracy of the bubble approach in 2D is to modify the
right-hand side of the bubble equations in (28) by multiplying with a
constant, say $\mu$. Then, $\alpha_{2}$ becomes

$\alpha_{2}=\frac{2\mu c^{2}h^{2}}{3(72-c^{2}h^{2})}.$

After this modification, the bubble functions are no more residual-free. We
call these modified functions as adaptive bubble functions. We call the
piecewise-defined linear approximations to these adaptive bubble functions as
pseudo-adaptive bubble functions.

We give two examples here to validate the approach. Figure 6 shows the graph
of $C_{1}$ and $C_{2}$ when $\mu=6.8$. It is clear that $C_{1}$ is decreased
in magnitude substantially. It is almost zero when $ch$ is close to zero.
There is not much change in $C_{2}$ in magnitude.

Figure 6: Comparison of the linear Galekin method and the pseudo-adaptive
bubbles method ($\mu=6.8$) for the coefficients $C_{1}$ (left) and $C_{2}$
(right).

For the second example, consider $\alpha_{2}=0.0625c^{2}h^{2}$ which makes
$C_{1}$ zero for all values of $ch$. In this case, the finite difference
scheme in (36) is a fourth order scheme with seven points for plane waves.
Figure 7 shows the graph of $C_{1}$ and $C_{2}$ when
$\alpha_{2}=0.0625c^{2}h^{2}$. While $C_{1}$ is zero for all values of $ch$,
there is only a slight change in $C_{2}$ in magnitude. This fourth order
accurate finite difference scheme can be easily applied in triangular,
trapezoidal and polygonal domains. Our main aim in this article is to propose
adaptive bubbles approximated by standard Galerkin method. However, the above
two method will be used to compare the success of the AB method.

Figure 7: Comparison of the linear Galekin method and the fourth order
accurate method for the coefficients $C_{1}$ (left) and $C_{2}$ (right).

##  5 Adaptive bubbles (AB) in 2D with triangular elements

We have shown using the pseudo-bubbles and the truncation error that the RFB
method is not effective in 2D. However, it is possible to increase its
accuracy with a simple modification, that is, multiplying the right-hand side
of the bubble problems with a constant. We proposed two methods using this
approach; a pseudo-adaptive bubbles method and a fourth order accurate finite
difference scheme that uses seven points. However, our main aim is to obtain
more accurate solutions by approximating the adaptive bubble functions with
linear finite element method on a coarse mesh.

We follow an empirical way to determine the optimal values of $\mu_{i}$ for
varying $cm_{i}$ where $m_{i}$ is the median of the global triangular element
(see Figure 8, (left)). This is actually a necessity because we have to use a
coarse mesh for the sub-problems and shapes of the bubble functions may change
significantly when a small change occurs in the number of mesh used for the
sub-problems.

Figure 8: A global mesh (left) and a decomposition of a global mesh with
triangular elements when $N_{s}=10$ (right).

Considering the problem in (25), we report the optimal values of $\mu_{i}$
($i=1,2,3$) for different values of $cm_{i}$ in Table 1 for an equilateral
triangular element. The optimality criterion in determining these values is
minimization of the error in inifinity norm by doing many tests. The values of
$\mu_{i}$ are optimal for $\theta=0,\pi/3,3\pi/2$. We decompose each global
mesh into triangular elements by choosing $N_{s}$ uniformly distributed nodes
on all edges of a triangular element (see Figure 8). We set $N_{s}=10$ when
$cm_{i}\leq 2.577$ and set $N_{s}=15$ when $cm_{i}>2.577$. While choosing
$N_{s}=10$ amounts to $28\times 28$ matrices, choosing $N_{s}=15$ amounts to
$78\times 78$ matrices on element level. It is possible to choose smaller
values of $N_{s}$, especially for smaller $ch$, but one has to report more
optimal values of $\mu$ in this case.

Table 1: Optimal values of $\mu_{i}$ for varying $cm_{i}$ ($i=1,2,3$) when $\theta=0,\pi/3,3\pi/2$. | $cm_{i}$ | $\mu_{i}$ | $N_{s}$
---|---|---
$\leq$ 0.57 | 5.4 | 10
0.583 | 5.43 | 10
0.644 | 5.45 | 10
0.71 | 5.5 | 10
0.7876 | 5.5 | 10
0.859 | 5.51 | 10
0.893 | 5.58 | 10
0.930 | 5.6 | 10
1.002 | 5.65 | 10
1.074 | 5.7 | 10
1.145 | 5.75 | 10
1.217 | 5.8 | 10
1.288 | 5.87 | 10
1.360 | 5.95 | 10
1.431 | 6.05 | 10
1.503 | 6.1 | 10
1.575 | 6.17 | 10
1.646 | 6.3 | 10
1.718 | 6.4 | 10
1.789 | 6.5 | 10
1.861 | 6.6 | 10
| $cm_{i}$ | $\mu_{i}$ | $N_{s}$
---|---|---
1.933 | 6.75 | 10
2.004 | 6.9 | 10
2.076 | 7.05 | 10
2.147 | 7.2 | 10
2.219 | 7.35 | 10
2.291 | 7.52 | 10
2.362 | 7.7 | 10
2.434 | 7.9 | 10
2.505 | 8.05 | 10
2.577 | 8.3 | 10
2.577 | 7.8 | 15
2.649 | 7.95 | 15
2.72 | 8.1 | 15
2.75 | 8.21 | 15
2.79 | 8.25 | 15
2.863 | 8.4 | 15
2.93 | 8.5 | 15
3.007 | 8.55 | 15
3.078 | 8.6 | 15
3.15 | 8.65 | 15

When $cm_{i}$ is between any of the successive two values in Table 1, we use
linear interpolation to get $\mu_{i}$. When shape of the global mesh changes,
the bubble functions behave differently, and hence it becomes more difficult
to find the optimal values of $\mu_{i}$ for each bubble functions. We
therefore expect deterioration of the AB method when nonuniform mesh is used,
especially for large wave numbers. Rectangular elements require solving 5
different bubble problems on each element and hence it becomes more
complicated to determine the optimal values of $\mu_{i}$ on nonuniform mesh.
Hence, we expect the triangular elements to be more efficient than the
rectangular elements on nonuniform mesh.

###### Remark 1

The optimal values in Table 1 were determined when $\theta=0$ and hence they
are true values for $\theta=\pi/3,2\pi/3$. The optimal values can be find for
other values of $\theta$. We will show by numerical test that the values of
$\mu_{i}$ in Table 1 can be used in any direction when $ch<1$. It is possible
to obtain good approximations up to $ch=2$ in any direction when $c<200$.


# Analysis of Load-Altering Attacks Against Power Grids: A Rare-Event Sampling
Approach

Maldon Patrice Goodridge1, Subhash Lakshminarayana2 and Christopher Few3
1Global Development Initiatives, Queen Mary University of London
2School of Engineering, University of Warwick, UK
3National Grid, UK
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

By manipulating tens of thousands of internet-of-things (IoT) enabled high-
wattage electrical appliances (e.g., WiFi-controlled air-conditioners), large-
scale load-altering attacks (LAAs) can cause severe disruptions to power grid
operations. In this work, we present a _rare-event sampling_ approach to
identify LAAs that lead to critical network failure events (defined by the
activation of a power grid emergency response (ER)). The proposed sampler is
designed to ‘skip’ over LAA instances that are of little interest (i.e., those
that do not trigger network failure), thus significantly reducing the
computational complexity in identifying the impactful LAAs. We perform
extensive simulations of LAAs using the Kundur two-area system (KTAS) power
network while employing the rare-event sampler. The results help us identify
the victim nodes from which the attacker can launch the most impactful attacks
and provide insights into how the spatial distribution of LAAs triggers the
activation of ERs.

##  I Introduction

Cyber attacks against power system grids can have significant social and
economic consequences. The threats can be broadly divided into two categories
– (i) attacks that directly target the power grid’s supervisory control and
data acquisition (SCADA) system, and (ii) attacks that indirectly target the
power grid’s control loops by manipulating end-user internet-of-things (IoT)
enabled electrical appliances. Direct attacks against the SCADA system, such
as false data injection attacks and/or coordinated cyber-physical attacks
against power grid state estimation have received significant attention [1, 2,
3]. In contrast, indirect attacks that target a large number of demand-side
appliances in a Botnet-type attack have been studied only recently [4, 5].
Unlike the SCADA assets, these devices cannot be monitored continuously due to
their large numbers.

The focus of this work is on _load-altering attacks_ (LAAs), which refer to a
sudden and abrupt change in the power grid demand by synchronously turning
on/off a large number of IoT-enabled high-wattage appliances [4, 5, 6, 7, 8].
LAAs pose a major threat to power grid operations since they can potentially
disrupt the balance between supply and demand. It has been shown [4, 5] that
such attacks can lead to unsafe frequency excursions, line outages, and/or
increase the grid’s operational costs. Moreover, dynamic LAAs, in which the
attacker injects a series of load perturbations over time, can also
destabilize the power grid’s frequency control loop [7]. Subsequent work has
also focused on detecting LAAs using data-driven approaches based on the data
gathered from phasor measurement units [9, 10].

Understanding the impact of LAAs is an important component in risk analysis.
Existing work on quantifying the impact of LAAs can be categorized into two
approaches – (i) simulation-based approach [4, 5, 6, 7], and (ii) analytical
approach [8]. Under the former, the attack impact is assessed by simulating
the power grid’s control loops (e.g., frequency dynamics) [4, 5, 6] or
computing the system’s eigenvalues under LAAs [7]. However, the results
presented in these works correspond to only a few specific LAA scenarios
(i.e., specific load perturbations injected at the victim nodes). Evaluating
the attack impact under all possible spatial distributions of LAAs (over the
victim nodes) requires performing extensive simulations considering different
combinations of the victim nodes and attack magnitudes, which can be
computationally prohibitive. To overcome these issues, an analytical approach
based on the theory of second-order dynamical systems was proposed in [8]. The
closed-form analytical functions to evaluate the attack impact (in terms of
the dynamic response and eigensolutions of the power grid’s frequency control
loop) only need to be computed once, thus avoiding the requirement for
repeated simulations. However, the analytical approach is restricted to a
second-order model with a direct current (DC) power flow model. Extending
these results to higher-order models (e.g., one that considers both frequency
and voltage dynamics) and involving the non-linear alternating current (AC)
power flow model is non-trivial. This is important to obtain a realistic
assessment of the LAA impact.

To overcome the aforementioned limitations, we apply a sampling approach to
map LAA magnitudes and their spatial distributions (across the different
victim nodes) to their attack impact. In particular, the focus of this work is
on LAAs that lead to the activation of power system emergency responses (ERs)
which disconnect critical power system components (e.g.,
generators/load/transmission lines). A common sampling approach is _Monte
Carlo_ simulations, which would apply randomly generated realisations of LAAs
from an underlying distribution to a simulated model of the physical system.
However, power grid design philosophies, such as $N-1$ scheduling, make the
grid resilient to various contingencies (including cyber attacks). This
implies component disconnections induced by LAAs can be extremely rare.
Consequently, Monte Carlo sampling can be computationally expensive, as the
majority of sampled LAAs will not result in component disconnections,
requiring a large number of realisations applied to the power system model to
generate a sample of the rare event (i.e., the activation of an ER).

To avoid an exhaustive search to locate potential LAAs that lead to network
failures, we instead employ a novel methodology in the context of LAAs based
on a Markov chain Monte Carlo (MCMC) approach for rare-event sampling, known
as the _skipping sampler_. The proposed approach is designed to reduce the
time and computational effort spent on evaluating LAAs that are of little
interest (i.e. those that do not result in an emergency response), allowing it
to more efficiently construct a sample of a rare event. The skipping sampler
has been applied to draw samples of low probability, high impact events in
power networks in the literature (see [12] and [13]). Section III-B provides a
detailed discussion on the skipping sampler.

We evaluate the framework by performing extensive simulations using the Kundur
two-area system (KTAS) [11]. We simulate the power grid’s transient dynamics
by a third-order model, which accounts for both the frequency as well as the
voltage dynamics [12]. The LAAs that perturb these dynamics are modelled
according to the _Log-normal_ distribution (since real-world cyber attack
magnitudes are well modelled by this distribution [13]). If the local
frequency metrics at any node exceed pre-set tolerances, appropriate ERs, such
as generation/load shedding, are activated. Our results show that in the KTAS
network, most instances of ERs are activated for two specific spatial
distributions of LAAs, namely (i) when the attacker increases the load
throughout the system or (ii) when the attacker decreases the load in the
over-provisioned area (with excess generation) and increases the load in the
under-provisioned area (with excess load). In particular, attacks that
exacerbate the power imbalance in the system (i.e., type (ii)) can trigger
inter-connector line disconnections, and lead to other network failures.

The rest of the paper is organised as follows. Section II introduces the
system model; Section III presents the statistical model and details of the
rare event sampling approach. Section IV describes the simulation results and
Section V concludes. The simulation parameters are provided in an online
Appendix found in [14].

##  II System model

###  II-A Power Grid Model

Using the third-order model for the generator, the power system model for
rare-event sampling simulations will also include a model for governor action,
automatic voltage regulation and a model of protection system ERs if the
frequency or RoCoF exceeds pre-defined thresholds. We consider a power grid
represented by $\mathcal{G}=\\{\mathcal{N},\mathcal{W}\\},$ where
$\mathcal{N}$ is the set of buses and $\mathcal{W}$ is the set of transmission
lines. The set of buses $\mathcal{N}$ consists of $N$ generation buses and $L$
load buses, with $|\mathcal{N}|=N+L$. At each generation bus $i=1,\dots N$,
the dynamics for the phase angle $\delta_{i}$, voltage magnitude $E_{i}$ and
governor action $\rho_{i}$ are given respectively by:

$\displaystyle
M(\psi)\ddot{\delta}_{i}+D\dot{\delta}_{i}=\psi_{i}\chi^{G}_{i}-\chi^{L}_{i}(\mathcal{R}_{i})-$
$\displaystyle\hskip
85.35826ptE_{i}\sum_{j=1}^{N+L}B_{ij}(\Omega_{ij})E_{j}\sin(\delta_{ij})$ (1a)
$\displaystyle
S_{i}\dot{E}_{i}=\psi_{i}(E_{f,i}-\text{v}_{i})-E_{i}+X_{i}\sum_{j=1}^{N+L}B_{ij}(\Omega_{ij})E_{j}\cos(\delta_{ij})$
(1b)
$\displaystyle\dot{\rho}_{i}=-A_{i}\dot{\delta}_{i}(1-1_{\mathcal{\mathcal{\mathcal{W}}}}[\dot{\delta}_{i}]).$
(1c)

In a similar manner, the dynamics for $\delta_{i}$ and $E_{i}$ at each load
bus $i=N+1,\dots,N+L$ are given by:

$\displaystyle
M(\psi)\ddot{\delta}_{i}+D\dot{\delta}_{i}=-\chi^{L}_{i}(\mathcal{R}_{i})-E_{i}\sum_{j=1}^{N+L}B_{ij}(\Omega_{ij})E_{j}\sin(\delta_{ij})$
(2a) $\displaystyle
S_{i}\dot{E}_{i}=\psi_{i}E_{f,i}-E_{i}+X_{i}\sum_{j=1}^{N+L}B_{ij}(\Omega_{ij})E_{j}\cos(\delta_{ij}.)$
(2b)

In equations (1) and (2), $\psi_{i}$, $\Omega_{ij}$ and $\mathcal{R}_{i}$ are
indicator variables associated with generator, line and load disconnections
respectively, explained in Section II-C. The system angular momentum,
$M(\psi)=\sum_{j=1}^{N}\psi_{j}H_{j}$ is given by the sum of each generator’s
inertia constant, $H_{i}$ (see online Appendix for parameter values).

TABLE I: Variables used in (1) and (2).  Symbol | Meaning | Units
---|---|---
$A_{i}$ | Governor’s droop response | MW/rad
$B_{ij}(\Omega_{ij})$ | Susceptance matrix | p.u.
$\chi_{i}^{G}$ | Net generation at node $i$ | p.u.
$\chi_{i}^{L}(\mathcal{R})$ | Net loads at node $i$ | p.u.
$D$ | System damping | %
$\delta_{i}$ | Phase angle | p.u
$\delta_{ij}$ | $\delta_{i}-\delta_{j}$ | p.u.
$\dot{\delta}_{i}$ | Frequency | p.u
$\ddot{\delta}_{i}$ | Rate of change of frequency (RoCoF) | p.u.
$E_{i}$ | Voltage | p.u.
$E_{f,i}$ | Machine $i$ rotor field voltage | p.u.
$M(\psi)$ | System angular momentum | Ws2
$\Omega_{ij}$ | Line disconnection indicator | -
$\psi_{i}$ | Generator shed indicator | -
$\mathcal{R}_{i}$ | UFLS counter | -
$S_{i}$ | Machine $i$ transient time constant | s
$X_{i}$ | Machine $i$ equivalent reactance | ohms
$\mathcal{W}$ | Governor’s deadband frequency range | Hz

As we assume the network to be lossless, the elements of $B_{ij}(\Omega_{ij})$
correspond to the imaginary part of the elements of the network’s admittance
matrix [11]. The net generation at node $i$ is given by
$\chi^{G}_{i}=\min\\{P_{i}^{\text{max}},P^{G}_{i}+\rho_{i}\\}$, where
$P_{i}^{\text{max}}$ is the nominal maximum power output of generator $i$,
$P_{i}^{G}$ is the equilibrium power of the generator and $\rho_{i}$ is the
power contributed by a governor unit, whose dynamics are given in (1c). The
variable vi accounts for the action of automatic voltage regulation (see
online Appendix). The net load at node $i$, $\chi^{L}_{i}$, is inclusive of
the LAA and a load disconnection scheme, and is discussed in Sections II-B and
II-C. The remaining parameters are given in Table I.

###  II-B Load-Altering Attack Model

Several security vulnerabilities have been identified in IoT-enabled high-
wattage consumer appliances (see, e.g., [5]). These vulnerabilities can be
exploited by a strategic attacker to cause security incidents such as
information disclosure and privilege escalation, leading to a change in the
device’s operational settings (e.g., switch ON/OFF or change the mode of
operation). Considering a $2~{}$kW power rating for the ACs, and a Botnet-
scale attack that potentially compromises tens of thousands of such devices,
LAAs can lead to a sudden load change of several MWs of power [5].

Not all loads are expected to be susceptible to LAAs, thus we decompose
$P^{L}_{i}$, the equilibrium load at bus $i$, into a vulnerable part, given by
$\nu P^{L}_{i}$, where $\nu\in[0,1]$ denotes the proportion of equilibrium
loads in the network vulnerable to an LAA; and a secure part (i.e., protected
or non-smart loads) $(1-\nu)P^{L}_{i}$. LAAs at node $i$, denoted $u_{i}$, are
modelled as $u_{i}\coloneqq\eta_{i}\nu P^{L}_{i}$, where $\eta_{i}\in[-1,1]$
is the proportional change to equilibrium vulnerable load. Thus, the load at
node $i$, inclusive of the LAA, is given by

$\displaystyle\chi^{L}_{i}$ $\displaystyle=(1-\nu)P^{L}_{i}+\nu
P^{L}_{i}+u_{i}$ $\displaystyle=(1-\nu)P^{L}_{i}+(1+\eta_{i})\nu P^{L}_{i}.$
(3)

This constrains the authority of the attacker to alter vulnerable loads at
node $i$ to a minimum of $0$ MW (i.e. $\min(u_{i})=-P^{L}_{i}\nu$), or, at
maximum, double vulnerable load demand (i.e. $\max(u_{i})=P^{L}_{i}\nu$). This
restriction to the maximum LAA reflects the finite capacity of inactive loads
the attacker can activate during an LAA.

###  II-C Emergency Responses

ERs refer to systems designed to protect sensitive power system components
from excessive frequency deviations following a change in the active power
balance. In this section we provide a brief discussion of ER employed, and
refer the reader to [15] and [16] for a detailed mathematical description.

_Generation shedding:_ To protect synchronous generators, we model two
independent schemes intended to disconnect the generator from the network: (i)
RoCoF-induced generation shedding (RIGS) - the generator is disconnected when
nodal RoCoF $|\ddot{\delta_{i}}|$ exceeds an upper threshold; (ii) over
frequency generation shedding (OFGS) - generation is shed when nodal frequency
$\dot{\delta_{i}}$ exceeds a pre-set upper limit. The binary variable
$\psi_{i}$ models the activation of generation shedding at node $i$ during the
simulation: $\psi_{i}=1$ initially under normal operation; however, when
either threshold is met and generator $i$ is disconnected, $\psi_{i}$ is set
to 0 until the end of the simulation.

_Under-frequency load shedding (UFLS):_ We model UFLS schemes as a progressive
disconnection of loads when the frequency $\dot{\delta}_{i}$ falls below a
strictly decreasing sequence of four frequency thresholds
$F^{U}\coloneqq\\{F^{U}_{1},\dots,F^{U}_{4}\\}$ where $F^{U}_{j-1}>F^{U}_{j}$.
In our model, at each frequency threshold, 10% of equilibrium loads
$P^{L}_{i}$ is automatically disconnected to arrest the decline in nodal
frequency. Letting $\mathcal{R}_{i}\in\\{0,1,2,3,4\\}$ count the total number
of UFLS activations at node $i$ at each time step $t$ in the simulation, the
net load is a dynamic variable in the power system model:

$\chi^{L}_{i}(\mathcal{R}_{i})=\bigg{(}1-0.1\mathcal{R}_{i}\bigg{)}\bigg{(}(1-\nu)P^{L}_{i}+(1+\eta_{i})\nu
P^{L}_{i}\bigg{)}$ (4)

_Line disconnection:_ In the KTAS network, we model the disconnection of the
line connecting Areas 1 and 2 when the power flow through the line, given by
$\phi_{ij}\coloneqq B_{ij}E_{i}E_{j}sin(\delta_{i}-\delta_{j})$, exceeds a
pre-set power threshold $P^{\phi}$. When excess power flow is detected through
the inter-connector line, the indicator $\Omega_{ij}$ switches from 1 to 0 for
the remainder of the simulation, setting the $ij^{th}$ element of $B_{ij}$ to
0 [16].

The ER model inspects the continuous time variables $\dot{\delta}_{i}(t)$,
$\ddot{\delta}_{i}(t)$ and $\phi_{ij}(t)$ from the power system model (1) at
regular time intervals. Once a criteria for activation is observed, the
corresponding ER is activated. This is represented in (1) as a discontinuity,
where changes to the relevant input variables (power injection, load or
network topology) are applied. Subsequently, the simulation is resumed with
the new network parameters.

##  III Statistical Model for LAAs

In this section, we present the statistical model for the distribution of LAAs
and describe the proposed rare-event sampling approach to identify the
impactful LAAs.

###  III-A Modeling the Unconditional Distribution of LAAs

There are several studies that document the frequency and magnitude of cyber
breaches in enterprise networks [13]. For instance, [13] demonstrates the size
of data breaches can be well-modelled by the log-normal family of
distributions. Assuming nodal LAAs magnitudes are independent and follow a
similar distribution, we model
$U\in\mathbb{R}^{N+L}\coloneqq[|u_{1}|,\dots,|u_{N+L}|]$ as

$U\sim\prod_{i=1}^{N+L}\text{Lognormal}(\mu_{i},\sigma_{i}^{2}).$ (5)

We also note that our analysis is not restricted to the log-normal family of
distributions, and can be extended to any underlying distribution in a
straightforward manner.

###  III-B Rare-Event Sampler for LAAs

If given an unconditional density $\rho$ over $\mathbb{R}^{N+L}$ and a rare
event of interest $C\subset\mathbb{R}^{N+L}$, rare-event sampling involves
drawing elements from $\rho$ conditioned on the occurrence $C$. The density of
this conditional distribution for the element $U\in\mathbb{R}^{N+L}$ is:

$\pi(U)=\frac{\rho(U)\mathbbm{1}_{C}(U)}{\rho(C)},$ (6)

where $\rho(C)$ is the probability of the event $C$ occurring, and

$\mathbbm{1}_{C}(U)=\begin{cases}1&U\in C\\\ 0&U\notin C.\end{cases}$ (7)

In the context of our research, $C$ is the set of LAA magnitudes which result
in the activation of at least one ER and $\rho$ is the distribution in (5). As
$C$ is expected to be rare, we employ the skipping sampler MCMC algorithm to
efficiently draw samples of $U\in C$. The skipping sampler is formalised in
Algorithm 1. We provide an intuitive explanation of the algorithm in the
following.

1 Input: initial state $U_{1}$;

2 for i = 1 to n do: Generate an initial proposal $Z_{1}$ distributed
according to the density $q(y-U_{i})dy$;

3 Calculate the direction
$\Phi=\left(Z_{1}-U_{i}\right)/\left\|Z_{1}-U_{i}\right\|$;

4 Generate a halting index $K\sim K_{\varphi}$;

5 Set $k=1$ and:

6 while  _$Z_{k}\notin C\textbf{{ and }}$ $k<K$_ do

7 Generate a distance increment $R$ distributed according to
$q_{r|\Phi}\left(r|\Phi\right)$;

8 Set $Z_{k+1}=Z_{k}+\Phi R$;

9 k=k+1;

10

11 end while

12Set $Z\coloneqq Z_{k}$;

13 Evaluate the acceptance probability:
$\alpha(U_{i},Z)=\begin{cases}\min\left(1,\frac{\pi(Z)}{\pi(U_{i})}\right)&\text{
if }\pi(U_{i})\neq 0,\\\ 1,&\text{ otherwise, }\end{cases}$ (8) Generate a
uniform random variable $V$ on $(0,1)$;

14 if  _$V\leq\alpha(U_{i},Z)$_ then

15 $U_{i+1}=Z$;

16 else

17 $U_{i+1}=U_{i}$;

18 end if

19return _$U_{i+1}$_.

Output: final sample $[U_{1},\dots,U_{h}]$

Algorithm 1 Skipping sampler algorithm

As a Metropolis-class algorithm [17], the skipping sampler can be understood
as a two step procedure: (1) a proposal step where, starting from a state
$U_{n}\in\mathbb{R}^{N+L}$, a potential new state $Z$ for the final sample is
generated; (2) an acceptance/rejection step, which determines whether the
proposed state $Z$ is included in the final sample, according to a specified
acceptance probability. This ensures the distribution of the sample follows
the desired target distribution $\pi$. If it is accepted, the proposal is
included in the final sample and becomes the starting state for the next
proposal step. This procedure is repeated a desired number of times, after
which the final sample is returned.

The skipping sampler improves the sampling of $C$ by using a specialised
proposal function which ‘skips’ over $C^{c}$ \- the set of LAAs which do not
lead to the activation of an ER (which are not of interest); until the rare
event is sampled or the skipping process is terminated. Thus, the proposal of
the skipping sampler was designed to enable efficient transitions between
connected components of $C$. Denoting the current state $U_{n}$, if the
initial proposal $Z_{1}\notin C$, we update the initial proposal (or ‘skip’)
by an adding an independent random distance increment $R_{2}$ in the direction
$\Phi=\frac{Z_{1}-U_{n}}{|Z_{1}-U_{n}|}$, where $R_{k}$ has the conditional
distribution of $\left\|Z_{1}-U_{n}\right\|$ conditioned on $\Phi$. The
proposal function continues this linear update procedure until either $C$ is
entered or the budget for skipping is exhausted [17].

##  IV Simulations

###  IV-A Simulation Settings

Our case study is based on the KTAS power grid [11]. We take a Kron reduced
version consisting of $N=4$ generation buses and $L=2$ load buses as shown in
Figure 1. At $t=0^{-}$, the system is modelled in equilibrium, with power
flows from Area 1 to Area 2 through the line connecting nodes $5$ and $6$ (tie
line). The system parameters are such that the system is $N-1$ secure, in the
sense that the loss of a generator (in the absence of any other disturbance)
does not trigger an ER. The initial conditions of the above system of
equations, denoted $\delta_{i}\left(0\right)$, $E_{i}\left(0\right)$,
$\rho_{i}(0)$, $P_{i}^{G}$ and $P^{L}_{i}$ are set equal to equilibrium states
which can be determined numerically, such that $\ddot{\delta}_{i}\approx 0$.
These values along the those of the parameters of (1) and (2) can be found in
the online Appendix [14].

Figure 1: Schematic drawing of the Kundur two-area 4 node network after Kron
reduction. Generator buses (green circles) correspond to nodes $i=[1,\dots,4]$
and load buses (brown circles) correspond to nodes $i=5,6$. Line lengths are
indicated.

To generate LAA instances, we implement the skipping sampler proposed in
Section III-B. During each proposal step, we sample an $N+L-$dimension LAA
vector $U$ from Lognormal distribution as in (5). We investigated various
values for $\sigma_{i}\in[1,8]$ which controls the rareness of large LAAs. For
this study, we present the results for $\mu_{i}=0$ and $\sigma_{i}=4$ for
$i=1,\dots,N+L$, and reserve a detailed sensitivity analysis for a dedicated
manuscript. We apply $U$ as an input to the power system model (1) at $t=0$,
with frequency dynamics simulated for $15$ seconds following the LAA using
MATLAB. We conduct $n=60,000$ proposals, which generated a final sample of
$h\approx 10,500$ LAAs conditioned on the activation of at least one ER. This
is an acceptance rate of 17.5%, within the 15 - 48% rate considered optimal
for exploring of a sample [18]. Following [17], no burn-in nor thinning was
required, thus all samples collected were available for analysis. Since these
responses are an undesirable event (from a system operator’s point of view),
we label such instances as a “network failure”.

Figure 2: Distribution of the absolute value LAAs at each node of the KTAS
network conditioned on a network failure event, for different levels of
network vulnerability $\nu$ and the unconditional distribution ($\rho$). Note
- bins for successive values of $\nu$ are offset slightly to the right to
improve readability.

###  IV-B Simulation Results

We evaluate the susceptibility of the network to LAAs at both a local and
global scope for three regimes of network vulnerability to LAAs- a ‘secure
network’ ($\nu\leq 45\%)$, a ‘moderately vulnerable network’ ($45\%<\nu\leq
65\%)$ and a ‘highly vulnerable network’ ($\nu>65\%)$.

Local Analysis – Identifying the vulnerable nodes: Figure 2 plots the
distribution of LAA magnitudes at each node, conditioned on the occurrence of
a network failure, for different degrees of network vulnerability $\nu$. We
observe the following : (i) only at nodes $5$ and $6$ does the conditional
distribution of LAA magnitudes differ significantly from $\rho$, with network
failures associated with larger magnitude LAAs, located in the low density
region of $\rho$. This is driven by the design of the KTAS network, where most
loads are concentrated at nodes $5$ and $6$, giving the attacker sufficient
leverage over system frequency to trigger a network failure. Thus, network
failures are primarily driven by LAAs at these nodes, mostly independent of
the LAA magnitude at nodes $1-4$; (ii) as network vulnerability $\nu$
increases, the distribution of LAA magnitudes at node 6 shifts rightwards,
implying larger magnitude LAAs at node $6$ become more prevalent in the
sample; however (iii) we note that large magnitude LAAs are rare under their
assumed Log-normal distribution. Thus, Fig. 2 reveals network failures, driven
by large LAA magnitudes at node $6$, are indeed low probability events.
Summarily, these results suggest nodes $5$ and $6$ are the critical nodes in
the KTAS network from which the attacker can launch the most impactful
attacks, regardless of the degree of network vulnerability.

Figure 3: Average number of each ER per sample $\epsilon_{i},i=1,\dots,4$
(defined in (9)) in the KTAS network for different degrees of network
vulnerability to LAAs.

Figure 3 illustrates the average number of activations of each ER per sample
for different levels of network vulnerability, defined as

$\displaystyle\epsilon_{f}=\frac{\text{Total ER${}_{f}$ activations in the
sample}}{\text{Sample size}}.$ (9)

Herein, the indices $f=1,\dots,4$ correspond to RoCoF-induced generation
shedding (RIGS), over-frequency generator shedding (OFGS), UFLS activations,
and inter-connector line trips respectively. Note that $\epsilon_{f}\in[0,1]$
for $f=1,2,4$. However, $\epsilon_{3}\in[0,4]$, as each node can experience a
maximum of 4 UFLS events ($f=3$) during the simulation (see Section II-C).

For secure networks, i.e, when $\nu\leq 45\%$, the average number of RIGS is
comparatively small, as low-magnitude LAAs are unable to induce sufficiently
large RoCoF deviations to trigger generator disconnection. Additionally, we
observe all generator nodes experience OFGS and UFLS events at similar rates
of $\epsilon_{2}\approx 0.1$ and $\epsilon_{3}\approx 1$ respectively.
Together, these imply for a secure KTAS network, each class of ER is similarly
likely at any node in the network. Thus, to counter the potential threat of
LAAs in a secure KTAS network will require a system-wide solution, e.g. - a
coordinated automatic generation control (AGC) system.

As $\nu$ increases, it grants the attacker greater authority over network
loads and the ability to significantly disrupt the active power balance of the
network. This is associated with increased rates of RIGS, UFLS and line
disconnections, as these are triggered by large power deviations. However,
when $\nu>65\%$, the average number of OFGS responses declines, as large
changes in loads result in RIGS dominating generation shedding events. In
contrast to secure networks, where nodes are similarly vulnerable to UFLS
events under LAAs, Figure 3 reveals a differential in the nodal risk of an
UFLS event in moderately and highly vulnerable networks. For example, when
$\nu>65\%$, nodes $3,4$ and $6$ experience approximately twice as many UFLS
events as nodes $1,2$ and $5$. This can be best understood through a global
analysis of the KTAS network, which is discussed next.

Global Analysis - Spatial Distribution of LAAs: Recall that the KTAS network
is comprised of two areas- Area 1 (nodes 1, 2 and 5) with excess generation,
and Area 2 (nodes 3, 4 and 6) with excess demand. With sufficient authority
over network loads, an attacker can exploit the global pre-LAA power
imbalances of each area to trigger an ER by decreasing loads in Area 1 to
exacerbate the excess generation, and increasing loads in Area 2 to exacerbate
the generation deficit. Thus, as $\nu$ increases, we observe fewer UFLS
responses in Area 1, more UFLS responses in Area 2, and a substantial increase
in power transferred across the inter-connector, increasing the rates of line
disconnections (Figure 3).

The global susceptibility of the KTAS network to LAAs is also illustrated by
Figure 4, which plots the probability of the relative changes in loads in each
Area of the KTAS network, conditioned on the occurrence of a network failure.
First, we note the KTAS network is secure against global reduction of loads in
both Areas, as it is likely governors are able to adjust generator output to
handle the loss of loads in all vulnerability regimes. Instead, we observe
that failure events occur primarily for large increases in loads, dominated by
two scenarios: (i) when the LAA increases the demand in both Areas- this
scenario is the primary driver of failure events when KTAS is secure, as the
attacker increases loads at both nodes 5 and 6 beyond the maximum power
capabilities of generators, triggering a failure anywhere in the network; (ii)
the attacker decreases loads in Area 1 and increases them in Area 2- this
scenario is rare when the KTAS is secure, as the attacker generally lacks
sufficient authority to exploit the imbalances of the KTAS network. However,
as $\nu$ (and the attacker’s authority) increases, so too does the conditional
probability of scenario (ii) (see Figure 3).

Figure 4: Probability of network failure for increase/decrease of system load
in each area for the KTAS network.

##  V Conclusions and Future Research

In this work, we present a framework to evaluate the impact of LAAs on power
grid operations using a rare-event sampling approach. The proposed approach
provides a comprehensive framework to examine the impact of LAAs under
different potential spatial distributions of LAAs across the power network.
Our results identify the nodes from which the attacker can launch the most
impactful LAA and further illustrate how the attacker can exploit the inter-
area power imbalances in the network to trigger ER events. Future work
includes (i) an extension of the analysis to dynamic LAAs against power grids
[7], (ii) considering correlations in LAA injections, (iii) showing
scalability of the proposed method to large-dimensional systems.

## References

  * [1] Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against state estimation in electric power grids,” in _Proc. ACM Conference on Computer and Communications Security (CCS)_ , 2009, pp. 21–32.
  * [2] S. Lakshminarayana, E. V. Belmega, and H. V. Poor, “Moving-target defense against cyber-physical attacks in power grids via game theory,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 6, pp. 5244–5257, 2021.
  * [3] S. Lakshminarayana, A. Kammoun, M. Debbah, and H. V. Poor, “Data-driven false data injection attacks against power grids: A random matrix approach,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 1, pp. 635–646, 2021.
  * [4] A. Dabrowski, J. Ullrich, and E. R. Weippl, “Grid shock: Coordinated load-changing attacks on power grids: The non-smart power grid is vulnerable to cyber attacks as well,” in _Proc. ACSAC_ , 2017, pp. 303–314.
  * [5] S. Soltan, P. Mittal, and H. V. Poor, “BlackIoT: IoT botnet of high wattage devices can disrupt the power grid,” in _Proc. USENIX Security Symposium_ , Baltimore, MD, Aug. 2018, pp. 15–32.
  * [6] B. Huang, A. A. Cardenas, and R. Baldick, “Not everything is dark and gloomy: Power grid protections against iot demand attacks,” in _Proc. USENIX Security Symposium_ , Aug. 2019, pp. 1115–1132.
  * [7] S. Amini, F. Pasqualetti, and H. Mohsenian-Rad, “Dynamic load altering attacks against power system stability: Attack models and protection schemes,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 4, pp. 2862–2872, July 2018.
  * [8] S. Lakshminarayana, S. Adhikari, and C. Maple, “Analysis of IoT-based load altering attacks against power grids using the theory of second-order dynamical systems,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 4415–4425, 2021.
  * [9] S. Amini, F. Pasqualetti, M. Abbaszadeh, and H. Mohsenian-Rad, “Hierarchical location identification of destabilizing faults and attacks in power systems: A frequency-domain approach,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 2, pp. 2036–2045, 2019.
  * [10] S. Lakshminarayana, S. Sthapit, H. Jahangir, C. Maple, and H. V. Poor, “Data-driven detection and identification of iot-enabled load-altering attacks in power grids,” _IET Smart Grid_ , pp. 1–16, 2022.
  * [11] P. Kundur and N. Balu, _Power System Stability and Control_. McGraw-Hill, 1994.
  * [12] K. Schmietendorf, J. Peinke, R. Friedrich, and O. Kamps, “Self-organized synchronization and voltage stability in networks of synchronous machines,” _The European Physical Journal Special Topics_ , vol. 223, no. 12, pp. 2577–2592, 2014.
  * [13] B. Edwards, S. Hofmeyr, and S. Forrest, “Hype and heavy tails: A closer look at data breaches,” _Journal of Cybersecurity_ , vol. 2, no. 1, pp. 3–14, 12 2016.
  * [14] “Online appendix,” https://tinyurl.com/ycknrnzv.
  * [15] M. Patrice Goodridge, J. Moriarty, and A. Pizzoferrato, “Distributions of cascade sizes in power system emergency response,” pp. 1–6, 2020.
  * [16] ——, “A rare-event study of frequency regulation and contingency services from grid-scale batteries,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 379, no. 2202, p. 20190433, 2021.
  * [17] J. Moriarty, J. Vogrinc, and A. Zocca, “The Skipping Sampler: A new approach to sample from complex conditional densities,” pp. 1–20, 2019. [Online]. Available: http://arxiv.org/abs/1905.09964
  * [18] D. Gamerman and H. Lopes, _Markov Chain monte Carlo, Stochastic Simulatoin for Bayesian Inference_ , 2nd ed. Boca Raton: Chapman & Hall/CRC, 2006.


[toc] language=isabelle, escapeinside=&&, columns=fixed, extendedchars,
basewidth=0.5em,0.45em, basicstyle=, mathescape, tabsize=4

SCHOOL OF COMPUTATION, INFORMATION AND TECHNOLOGY — INFORMATICS & MATHEMATICS

TECHNISCHE UNIVERSITÄT MÜNCHEN

Bachelor’s Thesis in Informatics

Bachelor’s Thesis in Mathematics

On the Formalization of Martingales

Eine Formalisierung von Martingalen

Author: Ata Keskin Supervisor: Prof. Dr. Tobias Nipkow Advisor: M. Sc.
Katharina Kreuzer Submission Date: 15 September 2023

##  Chapter 0 Abstract

This thesis presents a formalization of martingales in arbitrary Banach spaces
using Isabelle/HOL. We begin by examining formalizations in prominent proof
repositories and extend the definition of the conditional expectation operator
from the real numbers to general Banach spaces, drawing inspiration from prior
work. We define filtered measure spaces, adapted, progressively measurable and
predictable processes and rigorously formalize martingales, submartingales,
and supermartingales. Additionally, our contributions expand the scope of
Bochner integration techniques to general Banach spaces, and introduces
additional lemmas and induction schemes for integrable functions. Our
formalization provides a robust framework for future formalizations within the
theory of stochastic processes.

###### Contents

  1. 0 Abstract


* [95] D. Vollhardt, From Gutzwiller Wave Function to Dynamical Mean-Field Theory, in: Lecture Notes of the Autumn School on Correlated Electrons 2014, Vol. 4, Eds. E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein (Verlag des Forschungszentrum Jülich, 2014).
  * [96] D. Vollhardt, From Infinite Dimensions to Real Materials, in: Lecture Notes of the Autumn School on Correlated Electrons 2018, Vol. 8, Eds. E. Pavarini, E. Koch, A. Lichtenstein, and D. Vollhardt (Verlag des Forschungszentrum Jülich, 2018).
  * [97] G. Kotliar, Electronic Structure of Correlated Materials: Slave-Boson Methods and Dynamical Mean-Field Theory, in: Lecture Notes of the Autumn School on Correlated Electrons 2018, Vol. 8, Eds. E. Pavarini, E. Koch, A. Lichtenstein, and D. Vollhardt (Verlag des Forschungszentrum Jülich, 2018).
  * [98] J. K. Freericks and V. Zlatic, Exact dynamical mean-field theory of the Falicov-Kimball model, Rev. Mod. Phys. 75, 1333 (2003).
  * [99] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * [100] L. de’Medici, A. Georges, and S. Biermann, Orbital-selective Mott transition in multiband systems: Slave-spin representation and dynamical mean-field theory, Phys. Rev. B 72, 205124 (2005).
  * [101] E. Jakobi, N. Blümer, and P. van Dongen, Orbital-selective Mott transitions in a doped two-band Hubbard model with crystal field splitting, Phys. Rev. B 87, 205135 (2013).
  * [102] M. Kollar and D. Vollhardt, Exact analytic results for the Gutzwiller wave function with finite magnetization, Phys. Rev. B 65, 155121 (2002).
  * [103] R. Frésard and G. Kotliar, Interplay of Mott transition and ferromagnetism in the orbitally degenerate Hubbard model, Phys. Rev. B 56, 12909 (1997).
  * [104] T. Obermeier, T. Pruschke, and J. Keller, Ferromagnetism in the large-$U$ Hubbard model, Phys. Rev. B 56, R8479 (1997).
  * [105] P. Fazekas, Lecture Notes on Electron Correlation and Magnetism (World Scientific, Singapore, 1999).
  * [106] R. Bulla, Zero temperature metal-insulator transition in the infinite-dimensional Hubbard model, Phys. Rev. Lett. 83, 136 (1999).
  * [107] M. J. Rozenberg, R. Chitra, and G. Kotliar, Finite Temperature Mott Transition in the Hubbard Model in Infinite Dimensions, Phys. Rev. Lett. 83, 3498 (1999).
  * [108] R. Bulla, T. A. Costi, and D. Vollhardt, Finite-temperature numerical renormalization group study of the Mott transition, Phys. Rev. B 64, 045103 (2001).
  * [109] E. G. C. P. van Loon, F. Krien, and A. A. Katanin, Bethe-Salpeter Equation at the Critical End Point of the Mott Transition, Phys. Rev. Lett. 125, 136402 (2020).
  * [110] K. I. Kugel, D. I. Khomskii, A. O. Sboychakov, and S. V. Streltsov, Spin-orbital interaction for face-sharing octahedra: Realization of a highly symmetric SU(4) model, Phys. Rev. B 91, 155125 (2015).
  * [111] V. E. Valiulin, A. V. Mikheyenkov, K. I. Kugel, and A. V. Barabanov, Thermodynamics of Symmetric Spin-Orbital Model: One- and Two-Dimensional Cases, JETP Letters 109, 546-551 (2019).
  * [112] N. J. Laurita, J. Deisenhofer, L.-D. Pan, C. M. Morris, M. Schmidt, M. Johnsson, V. Tsurkan, A. Loidl, and N. P. Armitage, Singlet-Triplet Excitations and Long-Range Entanglement in the Spin-Orbital Liquid Candidate FeSc2S4, Phys. Rev. Lett. 114, 207201 (2015).
  * [113] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Kitaev quantum spin liquid—concept and materialization, Nature Rev. Phys. 1, 264-280 (2019).
  * [114] J. Bünemann, K. Jávorne-Radnóczi, P. Fazekas, and F. Gebhard, Orbital order in degenerate Hubbard models: A variational study, J. Phys.: Condensed Matter 19, 326217 (2007).
  * [115] Y. Nagaoka, Ferromagnetism in a narrow, almost half-filled $s$ band, Phys. Rev. 147, 392 (1966).
  * [116] J. van den Brink and D. Khomskii, Orbital ordering of complex orbitals in doped Mott insulators, Phys. Rev. B 63, 140416(R) (2001).
  * [117] https://www.humboldt-foundation.de/web/humboldt-preis.html