# Multi-orbital Flat Band Ferromagnetism with a Provable Percolation
Representation

Eric Bobrow  Junjia Zhang  Yi Li<EMAIL_ADDRESS>Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA

(November 20, 2020)

###### Abstract

We consider a two-layer multi-orbital system consisting of a $p_{x}$,
$p_{y}$-orbital honeycomb lattice layer and an $f$-orbital triangular lattice
layer with sites aligned with the centers of the honeycomb plaquettes. With an
appropriately tuned chemical potential difference between these two layers,
the system exhibits a flat band with provably ferromagnetic ground states at
half filling of the band in the presence of intra-orbital Hubbard interactions
and Hund’s coupling. Away from half filling, the interacting system admits a
percolation representation, where the ground state space is spanned by
maximum-spin clusters of localized single-particle states. A paramagnetic-
ferromagnetic transition occurs as the band approaches half filling and the
ground states become dominated by states with a large cluster.

##  I Introduction

Systems with flat bands can exhibit many interesting physical properties,
including ferromagnetism and exotic correlated states, where the flat band
enhances interaction effects, as well as topology in non-interacting flat
bands [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. In many models of flat band
systems, the flat band can be attributed to a large number of degenerate
states that are localized due to destructive interference. Tight-binding
models on line graphs [2] as well as decorated lattices such as the Tasaki
lattice [13] feature flat bands that can be understood in terms of a cell
construction [14] with destructively interfering hopping. In the presence of a
repulsive on-site Hubbard interaction, these models exhibit saturated
ferromagnetism when the flat band is half filled [15, 1]. For certain models
such as the Tasaki lattice, the intuition that overlapping single-particle
states maximize spin to avoid the repulsive interaction holds even below half
filling, where the interacting ground states are spanned by product states
formed from clusters of localized single-particle states, with the clusters
independently maximizing spin [16, 1]. Since the total spin of each cluster
depends only on its size, the result is a percolation representation that can
be efficiently simulated to find the transition from paramagnetic states with
small clusters at low filling to ferromagnetic states with large clusters at
high filling [17, 18].

Flat bands can also be realized in orbital-active systems, where orbital-
dependent hopping can prevent states from dispersing. A particular example of
this is the honeycomb lattice with $p_{x}$ and $p_{y}$ orbitals at each site
and orbital-dependent hopping along the bond directions, where it was shown in
Ref. [3] that both the lowest and highest bands are flat. The flat bands are
spanned by loop states around each honeycomb where each site on the loop
features a superposition of $p_{x}$ and $p_{y}$ orbitals to form a $p$ orbital
perpendicular to the outgoing bond that cannot hop out due to the bond-
projected hopping.

In order to find a percolation representation, Ref. [1] considers a flat band
with a basis of single-particle states that are quasilocal, meaning each basis
state is nonzero on a site where all other basis states vanish. When a
positive Hubbard $U$ interaction is turned on, the quasilocality condition
prevents double occupation of single-particle states. When the basis satisfies
the additional requirement that no more than two quasilocal states overlap at
any site, Hubbard interaction causes the overlapping states to maximize spin.
Breaking these conditions can introduce additional states into the interacting
ground state space that have lower spin within the cluster, making the
percolation picture non-rigorous. For the $p_{x}$\- and $p_{y}$-orbital model
of Ref. [3] together with an intra-orbital Hubbard interaction, the loop
states violate both conditions, with two or three loop states overlapping in
every orbital.

In this article, we study a two-layer system consisting of a $p_{x}$\- and
$p_{y}$-orbital honeycomb layer together with a triangular lattice $f$-orbital
layer. When the chemical potential of the $f$-orbital layer is appropriately
tuned, the system features a flat band spanned by localized states centered on
each $f$ orbital. We show that when this band is the highest energy band, the
system admits a percolation representation in the presence of Hund’s coupling
between the $p$ orbitals despite violating the condition in Ref. [1] that at
most two states overlap on site.

The remainder of this article is organized as follows. In Section II, we
introduce the multi-orbital model, find the condition under which the model
admits two bands, and discuss the localized single-particle states. In Section
III, we first review the percolation representation in Ref. [1] and then
extend it to our multi-orbital system with Hund’s coupling.

##  II Model Hamiltonian

The model Hamiltonian consists of spin-$1/2$ electrons in a two-layer system
with one layer a $p_{x}$, $p_{y}$-orbital honeycomb lattice and the other a
triangular $f$-orbital lattice with one $f_{y(3x^{2}-y^{2})}$ orbital per
site. The $f$-orbital lattice is arranged such that the $f$ orbitals are
aligned with the centers of the $p$-orbital honeycomb cells. The set of
p-orbital and f-orbital sites will respectively be labeled $\Lambda_{p}$ and
$\Lambda_{f}$, with these disjoint sets partitioning the overall lattice
$\Lambda=\Lambda_{p}\sqcup\Lambda_{f}$.

The kinetic Hamiltonian is $H_{K}=H_{K}^{p}+H_{K}^{f}+H_{K}^{fp}$ with
$p$-orbital part

$H_{K}^{p}=t_{p}\sum_{{\bf
r}\in\Lambda_{p}^{A}}\sum_{\sigma=\uparrow,\downarrow}\sum_{i=1}^{3}p^{\dagger}_{{\bf
r}+{\bf v}_{i},{\bf v}_{i},\sigma}p_{{\bf r},{\bf v}_{i},\sigma}+h.c.,$ (1)

where $\Lambda_{p}^{A}$ is the $A$ sublattice of the $p$-orbital layer and
$p_{{\bf r},{\bf v}_{i},\sigma}=\hat{v}_{i}\cdot\mathbf{p}_{{\bf r},\sigma}$
is the projection of the $p_{x}$, $p_{y}$ orbitals at site ${\bf r}$ in the
bond direction $\hat{{\bf v}}_{i}$, describing sigma bonding. The vectors
${\bf v}_{i}$ are the nearest neighbor vectors $v_{1}=(1,0)$,
$v_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$,
$v_{3}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$. Here we take the honeycomb lattice
bond length to be $1$. This $p_{x}$-$p_{y}$ orbital model has been previously
discussed in Refs. [3, 19, 20]. The $f$ orbital layer forms a triangular
lattice with sites ${\bf R}$ aligned with the centers of the $p$-orbital
honeycombs and are described by the kinetic Hamiltonian

$H_{K}^{f}=t_{f}\sum_{{\bf
R}\in\Lambda_{f}}\sum_{\sigma=\uparrow,\downarrow}\sum_{i=1}^{6}f^{\dagger}_{{\bf
R}+{\bf w}_{i},\sigma}f_{{\bf R},\sigma}+h.c.,$ (2)

where ${\bf w}_{i}$ are the six nearest-neighbor vectors on the $f$-orbital
triangular lattice, $w_{i}=\sqrt{3}(\cos\phi_{i},\sin\phi_{i})$ with
$\phi_{i}=\frac{\pi}{6}(2i-1)$ and $i=1,\cdots,6$.

Hopping between layers is described by

$H_{K}^{fp}=t_{fp}\sum_{{\bf
R}\in\Lambda_{f}}\sum_{\sigma=\uparrow,\downarrow}\sum_{i=1}^{6}(-1)^{i-1}f^{\dagger}_{{\bf
R},\sigma}p_{{\bf R}+{\bf u}_{i},{\bf u}^{\perp}_{i},\sigma}+h.c.$ (3)

where ${\bf u}_{i}=(\cos\theta_{i},\sin\theta_{i})$ with
$\theta_{i}=\frac{\pi}{3}(i-1)$, a reordering of the ${\bf v}$ vectors. In
particular, ${\bf R}+{\bf u}_{i}\in\Lambda_{p}$ when ${\bf R}\in\Lambda_{f}$.
Here ${\bf u}^{\perp}_{i}$ are defined to be unit vectors perpendicular to
${\bf u}^{\perp}_{i}$ with sign chosen so that ${\bf u}^{\perp}_{i+1}$ is a
$\pi/3$ rotation of ${\bf u}_{i}$. Thus, ${\bf
u}^{\perp}_{i}=(-\sin\theta_{i},\cos\theta_{i})$. The alternating sign in the
hopping amplitude is due to the fact that the $f_{y(3x^{2}-y^{2})}$ orbital
changes sign under $\frac{\pi}{3}$ rotation. To interpret the model as a two-
layer system, the ${\bf u}_{i}$ vectors can be thought of as the $xy$-plane
component of the vector between the $f$ orbitals on $\Lambda_{f}$ and the
adjacent $p$ orbitals on $\Lambda_{p}$, with the $z$ component being a small
inter-layer distance.

Typical hopping terms are depicted in Fig. 1, which shows that hopping between
$p$ orbitals with amplitude $t_{p}$ occurs between superpositions of $p_{x}$
and $p_{y}$ aligned along the bond direction and hopping between $p$ and $f$
orbitals involves superpositions of $p$ orbitals aligned perpendicular to the
in-plane hopping direction. Since the $p$ and $f$ orbitals are odd under
inversion, we take all hopping parameters to be positive,
$t_{p},t_{f},t_{fp}>0$.

In addition to the off-diagonal kinetic terms, we consider chemical potential
shifts for the $p$ and $f$ orbital layers,

$\displaystyle H_{\mu}$
$\displaystyle=-\mu_{p}\hat{N}_{p}-\mu_{f}\hat{N}_{f},$ (4)

where the total number operator for $p$-orbital electrons is
$\hat{N}_{p}=\sum_{{\bf
r}\in\Lambda_{p}}\sum_{p=p_{x},p_{y}}\sum_{\sigma=\uparrow\downarrow}n_{{\bf
r},p,\sigma}$ with $n_{{\bf r},p_{x},\sigma}=p^{\dagger}_{{\bf
r},\hat{x},\sigma}p_{{\bf r},\hat{x},\sigma}$ and a similar expression for
$n_{{\bf r},p_{y},\sigma}$. The total number operator for $f$-orbital
electrons is defined similarly, $\hat{N}_{f}=\sum_{{\bf
R}\in\Lambda_{f}}\sum_{\sigma=\uparrow\downarrow}n_{{\bf R},f,\sigma}$ with
$n_{{\bf R},f,\sigma}=f^{\dagger}_{{\bf R},\sigma}f_{{\bf R},\sigma}$.

The on-site intra-orbital Coulomb interaction is described by the Hubbard $U$
terms for $p$ and $f$ orbitals. Written in a particle-hole-symmetric form, the
interaction Hamiltonian $H_{U}$ is

$\displaystyle H_{U}$ $\displaystyle=U_{p}\sum_{{\bf
r}\in\Lambda_{p}}\sum_{p=p_{x},p_{y}}\left(n_{{\bf
r},p,\uparrow}-\frac{1}{2}\right)\left(n_{{\bf
r},p,\downarrow}-\frac{1}{2}\right)$ (5) $\displaystyle+U_{f}\sum_{{\bf
R}\in\Lambda_{f}}\left(n_{{\bf R},f,\uparrow}-\frac{1}{2}\right)\left(n_{{\bf
R},f,\downarrow}-\frac{1}{2}\right),$

which can also be written in terms of $H_{U}^{\prime}=U_{p}\sum_{{\bf
r}\in\Lambda_{p}}\sum_{p=p_{x},p_{y}}n_{{\bf r},p,\uparrow}n_{{\bf
r},p,\downarrow}+U_{f}\sum_{{\bf R}\in\Lambda_{f}}n_{{\bf
R},f,\uparrow}n_{{\bf R},f,\downarrow}$ as
$H_{U}=H_{U}^{\prime}-\frac{U_{p}}{2}\hat{N}_{p}-\frac{U_{f}}{2}\hat{N}_{f}+\frac{2U_{p}|\Lambda_{p}|+U_{f}|\Lambda_{f}|}{4}$.
We will consider only the repulsive case with $U_{p},U_{f}>0$.

As we will see, the intra-orbital Hubbard interactions alone are not
sufficient for the percolation representation we discuss in Section. III.2.
The ability of the Hund’s coupling in multi-orbital systems to polarize
electrons in degenerate orbitals will be essential. The on-site Hund’s
coupling between $p_{x}$ and $p_{y}$ orbitals is

$H_{J}=-J\sum_{{\bf r}\in\Lambda_{p}}\left(\mathbf{S}_{{\bf
r},p_{x}}\cdot\mathbf{S}_{{\bf r},p_{y}}-\frac{1}{4}n_{{\bf r},p_{x}}n_{{\bf
r},p_{y}}\right),$ (6)

where $S^{i}_{{\bf
r},p_{x/y}}=\frac{1}{2}\sum_{\mu,\nu=\uparrow,\downarrow}p^{\dagger}_{{\bf
r},\hat{x}/\hat{y},\mu}\sigma^{i}_{\mu\nu}p_{{\bf r},\hat{x}/\hat{y},\nu}$
with Pauli matrices $\sigma^{i}$, $i=x,y,z$, and the Hund’s coupling $J>0$.
The number operators without a spin index count both spin up and down, such as
$n_{{\bf r},p_{x}}=\sum_{\sigma=\uparrow,\downarrow}n_{{\bf r},p_{x},\sigma}$.
The Hund’s coupling energy from two electrons in an inter-orbital singlet
state at a site ${\bf r}$ is thus $J$, with an on-site triplet state having
$0$ Hund’s coupling energy.

Figure 1: Top view of the two-layer lattice. Red and blue respectively denote
positive and negative lobes for the $p$\- and $f$-orbital wavefunctions. (a)
Examples are shown of hopping $t_{f}$ between $f$ orbitals, $t_{p}$ between
$p$ orbitals along the honeycomb bond direction, and $t_{fp}$ between $f$ and
$p$ orbitals perpendicular to the hopping direction. (b) A localized state
formed by a superposition of an $f$ orbital and a surrounding loop of $p$
orbitals with alternating signs.

In the absence of the hoppings between $f$\- and $p$-orbitals, the decoupled
$p$-orbital Hamiltonian in Eq. (1), together with the chemical potential,
features two flat bands with energies $E_{p,\pm}=\pm\frac{3}{2}t_{p}-\mu_{p}$
as discussed in Refs. [3, 19, 20]. The lower flat band is spanned by localized
$p$-orbital loops with alternating sign $|\psi^{-(p)}_{{\bf
R},\sigma}\rangle=\frac{1}{\sqrt{6}}\sum_{i=1}^{6}(-1)^{i-1}p^{\dagger}_{{\bf
R}+{\bf u}_{i},{\bf u}_{i}^{\perp},\sigma}|0\rangle$ on each honeycomb, where
the $f$ orbital site ${\bf R}\in\Lambda_{f}$ is used to label the surrounding
honeycomb. The upper flat band is spanned by $|\psi^{+(p)}_{{\bf
R},\sigma}\rangle=\frac{1}{\sqrt{6}}\sum_{i=1}^{6}p^{\dagger}_{{\bf R}+{\bf
u}_{i},{\bf u}_{i}^{\perp},\sigma}|0\rangle$. These loop states feature a
superposition of $p_{x}$ and $p_{y}$ orbitals to form a $p$ orbital
perpendicular to the outgoing bond that cannot hop out due to the bond-
projected hopping. These states can also be thought of as localized due to
destructive interference in the $p_{x}$-$p_{y}$ basis.

In the presence of the $f$\- and $fp$-orbital hopping terms, these flat band
states can be extended to the full system. In general, $H_{fp}$ mixes the
$p$\- and $f$-orbital bands, potentially distorting the $p$-orbital flat
bands. We can find conditions under which the bands remain flat by first
assuming $|\psi_{{\bf R},\sigma}\rangle\equiv a|\psi^{-(p)}_{{\bf
R},\sigma}\rangle+b|f_{{\bf R},\sigma}\rangle$ and choosing coefficients $a,b$
such that $|\psi_{{\bf R},\sigma}\rangle$ is an eigenstate of $H_{K}$,

$\displaystyle H_{K}|\psi_{{\bf R},\sigma}\rangle$
$\displaystyle=(E_{p,-}-\mu_{p})a|\psi^{-(p)}_{{\bf
R},\sigma}\rangle+6t_{fp}a|f_{{\bf R},\sigma}\rangle$ (7)
$\displaystyle+t_{fp}b|\psi^{-(p)}_{{\bf R},\sigma}\rangle-\mu_{f}b|f_{{\bf
R},\sigma}\rangle$ $\displaystyle+(-at_{fp}+bt_{f})\sum_{i=1}^{6}|f_{{\bf
R}+{\bf w}_{i},\sigma}\rangle$ $\displaystyle=E_{fp}|\psi_{{\bf
R},\sigma}\rangle.$

The state $|\psi_{{\bf
R},\sigma}\rangle=t_{f}/{\mathcal{N}_{t}}|\psi^{-(p)}_{{\bf
R},\sigma}\rangle+t_{fp}/{\mathcal{N}_{t}}|f_{{\bf R},\sigma}\rangle$ with
normalization $\mathcal{N}_{t}=\sqrt{t_{f}^{2}+t_{fp}^{2}}$ is an eigenstate
with energy $E_{fp}=-(3/2)t_{p}-\mu_{p}+t_{fp}^{2}/t_{f}$ when the chemical
potential difference $\mu_{f}-\mu_{p}$ takes a critical value. As this can
also be thought of as tuning $\mu_{f}$ with $\mu_{p}$ fixed, we will consider
the $f$ orbital critical chemical potential

$\mu_{f}^{c}\equiv\mu_{p}+6t_{f}+\frac{3}{2}t_{p}-\frac{t_{fp}^{2}}{t_{f}},$
(8)

which can be satisfied in principle for any hopping amplitude values
$t_{p},t_{f},t_{fp}$. The localized states $|\psi_{{\bf R},\sigma}\rangle$ at
different ${\bf R}$ then span a flat band with energy $E_{fp}$.

One interesting feature of the model is that, for symmetry reasons, the upper
flat band of $H_{K}^{p}$ remains flat in the presence of $H_{K}^{f}$ and
$H_{K}^{fp}$ with energy $E_{p,+}=\frac{3}{2}t_{p}-\mu_{p}$, since the states
$|\psi^{+(p)}_{{\bf R},\sigma}\rangle$ centered at each honeycomb ${\bf R}$
destructively interfere at $|f_{{\bf R}}\rangle$ as well as at the neighboring
$f$ orbitals. This flat band requires no chemical potential tuning, though
whether it is the highest band of $H_{K}$ depends on the hopping amplitudes.

When $\mu_{f}=\mu_{f}^{c}$, there are five bands including the two flat bands
at energies $E_{fp}$ and $E_{p,+}$ and three dispersive bands with energies

$\displaystyle E_{2,\pm}({\bf k})$
$\displaystyle=-\mu_{p}\pm\frac{t_{p}}{2}\sqrt{4\cos{\frac{3k_{x}}{2}}\cos{\frac{\sqrt{3}k_{y}}{2}}+2\cos{\sqrt{3}k_{y}}+3}$
(9) $\displaystyle E_{1}({\bf k})$ $\displaystyle=-\frac{3t_{p}}{2}-\mu_{p}$
$\displaystyle+2t_{f}\left(2\cos{\frac{3k_{x}}{2}}\cos{\frac{\sqrt{3}k_{y}}{2}}+\cos{\sqrt{3}k_{y}}-3\right).$

In fact, $E_{2,\pm}$ are exactly the dispersive bands of the $p$-orbital model
$H^{p}_{K}$ [3]. When $\mu=\mu_{f}^{c}$, the addition of the $f$-orbital layer
leaves the dispersive bands and one flat band, $E_{p,+}$, of the $p$-orbital
layer unchanged, while introducing a dispersive band $E_{1}({\bf k})$ and
involving $f$ orbitals in the remaining flat band $E_{fp}$. These five bands
are shown in the first Brillouin zone in Fig. 2. The flat band of interest,
$E_{fp}$, is the highest energy band and does not feature any band touchings
when $t_{fp}^{2}>3t_{p}t_{f}$. We will assume that this is the case below. In
Fig. 2, the energy of the top band with a $\mu_{f}$ shift is plotted as a
function of $\mu_{f}-\mu_{f}^{c}$ and $k_{y}$, showing a flat band at
$\mu_{f}=\mu_{f}^{c}$.

Figure 2: For $\mu_{p}=0$, $t_{p}=5$, $t_{f}=1$, and $t_{fp}=4.5$, plots are
shown of (a) the band structure of $H_{K}$ in the first Brillouin zone with
$\mu_{f}=\mu_{f}^{c}$ and (b) the energy of the top band as a function of
$k_{y}$ and $\mu_{f}$. At $\mu_{f}=\mu_{f}^{c}$, the top band is flat.

##  III Ferromagnetic Percolation in Interacting Flat Bands

Before discussing the main result, we first briefly review the percolation
representation in the Hubbard model studied by Mielke and Tasaki [1]. We then
present and prove our results for the model in Sec. II, where the presence of
Hund’s coupling plays an important role.

###  III.1 Mielke-Tasaki Percolation in Flat Band Hubbard Models

For a system where the lowest-energy single-particle band is flat, a
percolation representation of ferromagnetism in the flat band was found in
Ref. [1] in terms of a linearly independent set of $N_{d}$ states
$\\{\varphi_{u}({\bf r})\\}$ spanning the space of the flat band. In this
notation, ${\bf r}$ is a lattice site and $u$ labels the state $\varphi_{u}$.
The index set of all state labels $u$ is denoted $\Lambda_{\varphi}$. The
percolation representation applies under two further conditions. First is
quasilocality, which requires each single-particle state to have a special
site in its support where every other state vanishes, i.e., for each $u$ there
is an ${\bf r}_{u}^{*}$ such that $\varphi_{v}({\bf r}_{u}^{*})\neq 0$ if and
only if $v=u$. Second is that no more than two states can overlap at any site,
i.e., for any site ${\bf r}$ there are at most two states $u$ and $v$ such
that $\varphi_{u}({\bf r})\neq 0$ and $\varphi_{v}({\bf r})\neq 0$. Two states
$u$ and $v$ are said to be directly connected or overlapping if there is such
a site ${\bf r}$ where both states are nonzero.

Under these conditions, Ref. [1] found that if the number of electrons
satisfies $N_{e}\leq N_{d}$, the ground states of the interacting system in
the presence of a repulsive single-orbital Hubbard interaction can be written
as a linear superposition of states formed from clusters of single-particle
states with maximum total spin.

The ground state space being spanned by the cluster states yields a
percolation representation where typical ground states are paramagnetic at low
filling and ferromagnetic at high filling. If filling is low, a typical ground
state consists of small, independent clusters that individually maximize spin
but can have low total spin. If filling is high, there is typically one
macroscopic cluster dominating the total spin. In particular, if
$N_{e}=N_{d}$, the ground state space is spanned by fully spin-polarized
states with trivial spin degeneracy. The above result from Ref. [1] thus
allows for a filling-dependent ferromagnetic transition to be studied by
treating the transition as a geometric percolation problem with additional
spin weights due to the cluster $S_{z}$ degeneracy [17].

This percolation representation was proven by constructing $N_{e}$-particle
interacting states from single-particle flat band states and requiring the
interacting states to have minimum Hubbard interaction energy, giving a ground
state of the interacting system since the flat band minimizes the kinetic
energy. In the proof of this result, the quasilocality condition is used to
prevent double occupation of single-particle states, and the requirement that
an overlap between single-particle states at a particular site involves only
two states $u$ and $v$ yields the fact that the overall state is symmetric
under swapping the spins of $u$ and $v$, indicating maximum total spin in the
cluster. Violating either of these conditions can lead to clusters with lower
total spin that still have zero interaction energy.

For the model in Sec. II, the localized states shown in Fig. 1 satisfy
quasilocality, as only one such state is nonzero at each $f$ orbital, but the
conditions of the above theorem do not hold since three single-particle states
have nonzero amplitude in the $p_{y}$ orbital component on each honeycomb
lattice site ${\bf r}\in\Lambda_{p}$. Thus, just the interaction $H_{U}$ is
insufficient to produce a percolation representation of maximum-spin clusters
in our multi-orbital model, and we will find that the combination of $H_{U}$
and $H_{J}$ allows for a percolation representation. We note as well that the
$f$ orbitals allowing quasilocality to be satisfied in the basis of localized
states $|\psi_{{\bf R},\sigma}\rangle$ is necessary to find a percolation
representation in terms of these localized states. Using features of the proof
of our main result in Sec. III.2, we demonstrate in Appendix A an explicit
example with a lower-than-maximum-spin cluster for the flat bands of the
p-orbital Hamiltonian $H^{p}_{K}$ even when intra-orbital Hubbard interactions
and Hund’s coupling are both considered.

###  III.2 Percolation representation with Hund’s coupling

We will now see, by studying the model in Sec. II that a multi-orbital
extension of the argument in Ref. [1] together with Hund’s coupling can lead
to a percolation representation that supports a ferromagnetic transition even
when there is an overlap between more than two localized states. The argument
then proceeds as follows. First, we perform a particle-hole transformation to
study a model with a lowest-energy flat band. Next, we write an arbitrary
ground state as a superposition of $N_{e}$-particle states in the localized
state basis. Then, we require the ground state to be a zero-energy eigenstate
of the intra-orbital Hubbard interaction in the $p_{x}$ and $p_{y}$ orbitals
as well as of the Hund’s coupling. This will show that when the particle-hole-
transformed flat band is at most half filled, the ground state space is
spanned by states with maximum-spin clusters. Finally, we invert the particle-
hole transformation to show that when the highest-energy flat band of the
original system is at least half filled, the ground state space is also
spanned by states with maximum-spin clusters.

Let us first define a creation operator $a^{\dagger}_{{\bf R},\sigma}$ for a
localized state centered at site ${\bf R}\in\Lambda_{f}$,

$\displaystyle a^{\dagger}_{{\bf R},\sigma}$ $\displaystyle\equiv\sum_{{\bf
R}\in\Lambda_{f}}\varphi_{\bf R}({\bf R},f)c^{\dagger}_{\sigma}({\bf
R},f)+\sum_{\begin{subarray}{c}{\bf r}\in\Lambda_{p}\\\
p=p_{x},p_{y}\end{subarray}}\varphi_{\bf R}({\bf
r},p)c^{\dagger}_{\sigma}({\bf r},p)$ (10)
$\displaystyle\equiv\bm{\varphi}_{{\bf R}}\cdot\mathbf{c}^{\dagger}_{\sigma},$

with $\mathbf{c}^{\dagger}_{\sigma}$ a compact notation for $f$-, $p_{x}$-,
and $p_{y}$-orbital electron creation operators on the entire lattice
$\Lambda$ as a $(|\Lambda_{f}|+2|\Lambda_{p}|)$-component vector. The
wavefunction $\bm{\varphi}_{{\bf R}}$ is compactly supported, taking nonzero
values only for the $f$ orbital at ${\bf R}$ and the $p_{x}$ and $p_{y}$
orbitals at ${\bf R}+{\bf u}_{i}$ with $i=1,2,\dots,6$. Due to the localized
state structure with $p$ orbitals perpendicular to the outgoing bonds, the
$p_{x}$ components vanish at ${\bf R}+{\bf u}_{1}$ and ${\bf R}+{\bf u}_{4}$
as well. In this notation, a component of $\mathbf{c}^{\dagger}_{\sigma}$ will
be labeled $c^{\dagger}_{\sigma}({\bf r},o_{\bf r})$ where ${\bf r}$ is a site
on the lattice and $o_{\bf r}$ is the orbital index at that site. Thus, if
${\bf r}={\bf R}\in\Lambda_{f}$, $o_{\bf r}=f$ and $c^{\dagger}_{\sigma}({\bf
R},f)=f^{\dagger}_{{\bf R},\sigma}$. Similarly, there are two operators
corresponding to $o_{\bf r}=p_{x}$ and $o_{\bf r}=p_{y}$ if ${\bf
r}\in\Lambda_{p}$. The dot product notation is a shorthand for summing over
all sites and orbitals. The subscript ${\bf R}\in\Lambda_{f}$ in
$\bm{\varphi}_{\bf R}$ identifies the localized state centered at ${\bf R}$.

The maximum-spin cluster states that will span the ground state space for
appropriate filling are defined by selecting a subset $A\subset\Lambda_{f}$
and placing one single-particle localized state $\bm{\varphi}_{\bf R}$ at each
${\bf R}\in A$. The set $A$ can then be partitioned into
$A=\bigsqcup_{k=1}^{n}C_{k}$, where $C_{k}$ are disjoint clusters, and two
states labeled by ${\bf R},{\bf R}^{\prime}\in A$ belong to the same cluster
if they are overlapping or connected by a path of overlapping states. Since
total spin is maximized within each cluster, these states can be constructed
from states $|\Phi_{A\uparrow}\rangle=\prod_{{\bf R}\in A}a^{\dagger}_{{\bf
R},\uparrow}|0\rangle$ or $|\Phi^{(h)}_{A\uparrow}\rangle=\prod_{{\bf R}\in
A}a_{{\bf R},\uparrow}|F\rangle$ with $|F\rangle=\prod_{{\bf
r}\in\Lambda}\prod_{o_{\bf r}}c_{\uparrow}^{\dagger}({\bf r},o_{\bf
r})c_{\downarrow}^{\dagger}({\bf r},o_{\bf r})|0\rangle$ the fully filled
state, depending on whether the flat band is the lowest- or highest-energy
band.

For the particle-hole transformed Hamiltonian, where the flat band is the
lowest-energy band, the $S_{z}$ spin of each cluster can be lowered by a
cluster spin-lowering operator $S^{-}_{C_{k}}=\sum_{({\bf r},o_{\bf r})\in
V_{k}}S^{-}_{{\bf r},o_{\bf r}}$ with $S^{-}_{{\bf r},o_{\bf
r}}=c^{\dagger}_{\downarrow}({\bf r},o_{\bf r})c_{\uparrow}({\bf r},o_{\bf
r})$ acting on $V_{k}=\\{({\bf r},o_{\bf r})|\varphi_{\bf R}({\bf r},o_{\bf
r})\neq 0\mathrm{\ for\ any\ }{\bf R}\in C_{k}\\}$, the set of orbitals where
at least one state in the cluster is nonzero. For the original Hamiltonian,
where the flat band is the highest-energy band, the analogous operator acting
on holes is $S^{(h)-}_{C_{k}}=\sum_{({\bf r},o_{\bf r})\in
V_{k}}c_{\downarrow}({\bf r},o_{\bf r})c^{\dagger}_{\uparrow}({\bf r},o_{\bf
r})$.

For the particle-hole transformed Hamiltonian, found by replacing creation and
annihilation operators $f_{{\bf R},\sigma}\leftrightarrow f^{\dagger}_{{\bf
R},\sigma}$, $p_{{\bf r},\hat{x},\sigma}\leftrightarrow p^{\dagger}_{{\bf
r},\hat{x},\sigma}$, and $p_{{\bf r},\hat{y},\sigma}\leftrightarrow
p^{\dagger}_{{\bf r},\hat{y},\sigma}$ in $H$, we find the following theorem.

###### Theorem 1

Consider the particle-hole transformed Hamiltonian
$H^{(ph)}=-H_{K}+H_{U}+H_{J}-\frac{J}{2}\hat{N}_{p}-H_{\mu}$ with
$t_{fp}^{2}>3t_{p}t_{f}$ and $t_{p},t_{f},t_{fp},U_{p},U_{f},J>0$. When
$N_{e}\leq|\Lambda_{f}|$ and
$\mu_{f}=\mu_{f}^{(ph),c}\equiv\mu_{f}^{c}+U_{f}/2-U_{p}/2+J/2$, with
$\mu_{f}^{c}$ defined in Eq. (8), the ground state space of $H^{(ph)}$ is
spanned by the states

$|\Phi^{(ph)}_{A,\\{m_{k}\\}}\rangle=\prod_{k=1}^{n}(S_{C_{k}}^{-})^{\frac{|C_{k}|}{2}-m_{k}}|\Phi_{A\uparrow}\rangle$
(11)

with $A\subset\Lambda_{f}$ and $|A|=N_{e}$.

One important detail is that the required $\mu_{f}$ in Theorem 1 is not the
$\mu_{f}^{c}$ for which $H_{K}$ has a highest-energy flat band. Instead, there
is a shift due to the additional chemical potential terms in the particle-
hole-symmetric $H_{U}$ and the particle-hole-transformed Hund’s coupling
$H_{J}-\frac{J}{2}\hat{N}_{p}$. When $\mu_{f}=\mu_{f}^{(ph),c}$, the single-
particle terms in $H^{(ph)}$ have a lowest-energy flat band spanned by the set
of states $\\{\bm{\varphi}_{\bf R}\\}_{{\bf R}\in\Lambda_{f}}$. We have
discarded additional constant terms in $H^{(ph)}$, as they will not affect the
spectrum or chemical potential condition.

For our original model, we find the following theorem.

###### Theorem 2

Consider $H$ defined in Sec. II with $t_{fp}^{2}>3t_{p}t_{f}$ and
$t_{p},t_{f},t_{fp},U_{p},U_{f},J>0$. When $N_{e}\geq
4|\Lambda_{p}|+|\Lambda_{f}|$ and $\mu_{f}=\mu_{f}^{(ph),c}$, defined in
Theorem. 1, the ground state space of $H$ is spanned by

$|\Phi_{A,\\{m_{k}\\}}\rangle=\prod_{k=1}^{n}(S^{(h)-}_{C_{k}})^{\frac{|C_{k}|}{2}-m_{k}}|\Phi^{(h)}_{A\uparrow}\rangle,$
(12)

where $A\subset\Lambda_{f}$, $|A|=N_{e}$, and $A=\bigsqcup_{k=1}^{n}C_{k}$
where $C_{k}$ are disjoint clusters.

We now proceed with the proof of Theorem 1, from which, our main result,
Theorem 2 immediately follows by a particle-hole transformation. Importantly,
the basis states $|\Phi_{A,\\{m_{k}\\}}\rangle$ feature maximum total spin
within each cluster, since $|\Phi_{A,\\{m_{k}\\}}\rangle$ is a particle-hole
transformation of $|\Phi^{(ph)}_{A,\\{m_{k}\\}}\rangle$ and total spin
commutes with particle-hole transformations, as reviewed in Appendix C.

Proof: Following Ref. [1], we construct operators canonically conjugate to
$a^{\dagger}_{{\bf R},\sigma}$ by defining

$\displaystyle b_{{\bf R},\sigma}$ $\displaystyle=\bm{\kappa}_{{\bf
R}}\cdot\mathbf{c}_{\sigma},$ (13) $\displaystyle\kappa_{{\bf R}}({\bf
r},o_{\bf r})$ $\displaystyle=\sum_{{\bf
R}^{\prime}\in\Lambda_{f}}(G^{-1})_{{\bf R},{\bf R}^{\prime}}\varphi_{{\bf
R}^{\prime}}({\bf r},o_{\bf r}),$

with $G_{{\bf R},{\bf R}^{\prime}}=\bm{\varphi}_{{\bf
R}}\cdot\bm{\varphi}_{{\bf R}^{\prime}}$ the Gram matrix for the states
$\bm{\varphi}_{{\bf R}}$. Thus, $\bm{\kappa}_{\bf R}\cdot\bm{\varphi}_{{\bf
R}^{\prime}}=\delta_{{\bf R},{\bf R}^{\prime}}$ implies the canonical
anticommutation relation $\\{b_{{\bf
R}^{\prime},\sigma^{\prime}},a^{\dagger}_{{\bf R},\sigma}\\}=\delta_{{\bf
R},{\bf R}^{\prime}}\delta_{\sigma,\sigma^{\prime}}$. The states
$\bm{\kappa}_{\bf R}$ serve as an alternate basis for the flat band, and this
Gram matrix procedure is essentially a method of constructing a dual basis
where each element of the dual basis is orthogonal to all but one of the
vectors in the original basis. This procedure is similar to the construction
of reciprocal lattice vectors from direct lattice vectors.

Now $c_{\sigma}({\bf r},o_{\bf r})$ can be expressed in terms of $b_{{\bf R}}$
by multiplying by $\varphi_{\bf R}({\bf r},o_{\bf r})$ and summing over ${\bf
R}$, an operator $d_{\sigma}({\bf r},o_{\bf r})$ can be defined to express the
electron annihilation operators as

$\displaystyle c_{\sigma}({\bf r}^{\prime},o^{\prime}_{{\bf
r}^{\prime}})=\sum_{{\bf R}\in\Lambda_{f}}\varphi_{{\bf R}}({\bf r},o_{\bf
r})b_{{\bf R},\sigma}-d_{\sigma}({\bf r},o_{\bf r}),$ (14)

where $d_{\sigma}({\bf r},o)=\sum_{{\bf r}^{\prime}\in\Lambda}\psi({\bf
r},o_{\bf r};{\bf r}^{\prime},o_{{\bf r}^{\prime}})c_{\sigma}({\bf
r}^{\prime},o^{\prime}_{{\bf r}^{\prime}})$ with $\psi({\bf r},o_{\bf r};{\bf
r}^{\prime},o^{\prime}_{{\bf r}^{\prime}})\equiv\delta_{{\bf r},{\bf
r}^{\prime}}\delta_{o_{\bf r},o^{\prime}_{{\bf r}^{\prime}}}-\sum_{{\bf
R}\in\Lambda_{f}}\varphi_{{\bf R}}({\bf r},o_{\bf r})\kappa_{{\bf R}}({\bf
r}^{\prime},o^{\prime}_{{\bf r}^{\prime}})$ a projection out of the flat band
spanned by $\bm{\varphi_{\bf R}}$ and $\bm{\kappa_{\bf R}}$. Thus,
$\\{d_{\sigma}({\bf r},o_{\bf r}),a^{\dagger}_{{\bf
R},\sigma}\\}=\\{d_{\sigma}({\bf r},o_{\bf r}),b^{\dagger}_{{\bf
R},\sigma}\\}=0$. Since the ground state will be expressed in terms of
$a^{\dagger}_{{\bf R},\sigma}$ operators, this construction allows the ground
state condition in the presence of interactions to be analyzed using only
states orthogonal or canonically conjugate to the single-particle localized
flat band states.

When it is possible to construct an $N_{e}$-electron state from single-
particle flat band states that simultaneously minimizes the interactions, such
a state will be a ground state, and the ground state space will be spanned by
the collection of these states. In particular, the
$n_{\uparrow}n_{\downarrow}$ interaction $H_{U}^{\prime}$ and $H_{J}$ are
positive semidefinite, so when there are states satisfying
$H_{U}^{\prime}|\Phi\rangle=H_{J}|\Phi\rangle=0$ with energy $-N_{e}E_{fp}$,
these states will be ground states of $H^{(ph)}$ and can be written

$|\Phi^{(ph)}\rangle=\sum_{A_{\uparrow},A_{\downarrow}\subset\Lambda_{f}}f(A_{\uparrow},A_{\downarrow})\prod_{{\bf
R}\in A_{\uparrow}\sqcup A_{\downarrow}}\prod_{\sigma_{\bf
R}}a^{\dagger}_{{\bf R},\sigma_{{\bf R}}}|0\rangle,$ (15)

with constraints on the coefficients $f(A_{\uparrow},A_{\downarrow})$ to be
determined by the zero-interaction-energy conditions. Here $A_{\uparrow}$ and
$A_{\downarrow}$ are subsets of $\Lambda_{f}$ and $\sigma_{\bf R}=\uparrow$ or
$\downarrow$ if ${\bf R}\in A_{\uparrow}$ or ${\bf R}\in A_{\downarrow}$. If
${\bf R}\in A_{\uparrow}\cap A_{\downarrow}$, the product over $\sigma_{\bf
R}$ includes both, with the spin-up operator to the left. The sum over
$A_{\uparrow},A_{\downarrow}$ is a sum over all possible such subsets
satisfying $|A_{\uparrow}|+|A_{\downarrow}|=N_{e}$. We will see that
$|\Phi^{(ph)}\rangle=0$ if $N_{e}>|\Lambda_{f}|$.

Since $H_{U}^{\prime}$ and $H_{J}$ are themselves sums of positive-
semidefinite operators at each site, we first consider the Hubbard interaction
on an $f$ orbital at ${\bf R}$, which gives the condition

$\displaystyle 0$ $\displaystyle=c_{\uparrow}({\bf R},f)c_{\downarrow}({\bf
R},f)|\Phi^{(ph)}\rangle$ (16) $\displaystyle=\left(d_{\uparrow}({\bf
R},f)+\sum_{{\bf R}^{\prime}\in\Lambda_{f}}\varphi_{{\bf R}^{\prime}}({\bf
R},f)b_{{\bf R}^{\prime},\uparrow}\right)$
$\displaystyle\times\left(d_{\downarrow}({\bf R},f)+\sum_{{\bf
R}^{\prime\prime}\in\Lambda_{f}}\varphi_{{\bf R}^{\prime\prime}}({\bf
R},f)b_{{\bf R}^{\prime\prime},\downarrow}\right)|\Phi^{(ph)}\rangle$
$\displaystyle\implies b_{{\bf R},\uparrow}b_{{\bf
R},\downarrow}|\Phi^{(ph)}\rangle=0,$

since $d_{\sigma}({\bf r},o_{\bf r})$ anticommutes with $b_{{\bf R},\sigma}$
and with the $a^{\dagger}_{{\bf R},\sigma}$ operators in $|\Phi\rangle$, and
$\varphi_{{\bf R}^{\prime}}({\bf R},f)\neq 0$ only when ${\bf R}^{\prime}={\bf
R}$. Thus, since this condition holds for any ${\bf R}\in\Lambda_{f}$, the
state $|\Phi^{(ph)}\rangle$ must satisfy $f(A_{\uparrow},A_{\downarrow})=0$ if
$A_{\uparrow}\cap A_{\downarrow}\neq\emptyset$, or, in other words, there must
be no double occupancy of localized single-particle states. If
$N_{e}>|\Lambda_{f}|$, this condition can only be met if
$|\Phi^{(ph)}\rangle=0$, meaning the ground state cannot be expressed solely
in terms of $N_{e}$ flat band states and must have energy higher than
$-N_{e}E_{fp}$.

Next, we examine the Hubbard interaction in the $p_{x}$ and $p_{y}$ orbitals.
Consider a honeycomb cell labeled by ${\bf R}_{0}\in\Lambda_{f}$ and examine
its rightmost vertex, ${\bf r}_{0}={\bf R}+{\bf u}_{1}$. Two additional
honeycomb cells, centered at ${\bf R}_{1}={\bf R}_{0}+{\bf w}_{1}$ and ${\bf
R}_{6}={\bf R}_{0}+{\bf w}_{6}$, share vertex ${\bf r}_{0}$. The corresponding
localized single-particle states $\bm{\varphi}_{{\bf R}_{i}}$ have nonzero
component in the $p_{y}$ orbital at ${\bf r}_{0}$ for all three of ${\bf
R}_{0}$, ${\bf R}_{1}$, and ${\bf R}_{6}$, but only ${\bf R}_{1}$ and ${\bf
R}_{6}$ have a nonzero $p_{x}$ component at ${\bf r}_{0}$. Excluding the
normalization factor $t_{f}/\sqrt{6(t_{f}^{2}+t_{fp}^{2})}$ on the $p$-orbital
components, the nonzero $p_{x}$ components are, $\varphi_{{\bf R}_{1}}({\bf
r}_{0},p_{x})=-\frac{\sqrt{3}}{2}$ and $\varphi_{{\bf R}_{6}}({\bf
r}_{0},p_{x})=\frac{\sqrt{3}}{2}$ while the nonzero $p_{y}$ components are
$\varphi_{{\bf R}_{0}}({\bf r}_{0},p_{y})=1$, $\varphi_{{\bf R}_{1}}({\bf
r}_{0},p_{y})=-\frac{1}{2}$, and $\varphi_{{\bf R}_{6}}({\bf
r}_{0},p_{y})=-\frac{1}{2}$. No other localized states $\bm{\varphi}_{\bf R}$
are nonzero at ${\bf r}_{0}$.

The zero-interaction-energy condition for the $p_{x}$-orbital Hubbard
interaction at site ${\bf r}_{0}$ is

$\displaystyle 0$ $\displaystyle=c_{\uparrow}({\bf
r}_{0},p_{x})c_{\downarrow}({\bf r}_{0},p_{x})|\Phi^{(ph)}\rangle$ (17)
$\displaystyle=\left(\sum_{i=1,6}\varphi_{{\bf R}_{i}}({\bf
r}_{0},p_{x})b_{{\bf R}_{i},\uparrow}\right)$
$\displaystyle\times\left(\sum_{i=1,6}\varphi_{{\bf R}_{i}}({\bf
r}_{0},p_{x})b_{{\bf R}_{i},\downarrow}\right)|\Phi^{(ph)}\rangle$
$\displaystyle\implies(b_{{\bf R}_{6},\uparrow}b_{{\bf
R}_{1},\downarrow}-b_{{\bf R}_{6},\downarrow}b_{{\bf
R}_{1},\uparrow})|\Phi^{(ph)}\rangle=0,$

where the last line follows from the no-double-occupancy condition Eq. (16)
and the fact that $\varphi_{{\bf R}_{0}}({\bf r}_{0},p_{x})=0$. The condition
in Eq. (17) essentially projects out states involving a spin-singlet component
between the states centered at ${\bf R}_{1}$ and ${\bf R}_{6}$. Explicitly,
this condition gives that for any configuration where
$A_{\uparrow}=B_{\uparrow}\sqcup\\{{\bf R}_{1}\\}$ and
$A_{\downarrow}=B_{\downarrow}\sqcup\\{{\bf R}_{6}\\}$ with ${\bf R}_{1},{\bf
R}_{6}\notin B_{\uparrow},B_{\downarrow}$, the coefficients are symmetric
under exchange of spins, $f(B_{\uparrow}\sqcup\\{{\bf
R}_{1}\\},B_{\downarrow}\sqcup\\{{\bf R}_{6}\\})=f(B_{\uparrow}\sqcup\\{{\bf
R}_{6}\\},B_{\downarrow}\sqcup\\{{\bf R}_{1}\\})$. Since the choice of ${\bf
R}_{0}$ is arbitrary, this gives the general condition that the state must
have spin exchange symmetry between localized single-particle states that are
nearest neighbors in the $y$ direction, states ${\bf R}$ and ${\bf
R}^{\prime}={\bf R}+{\bf w}_{2}$.

Eq. (16) and Eq. (17) are equivalent to the conditions resulting from
quasilocality and the requirement that no more than two single-particle states
overlap at any site in Ref. [1]. However, three states overlap in the $p_{y}$
orbitals in the model we consider. For the $p_{y}$ orbital at site ${\bf
r}_{0}$,

$\displaystyle 0$ $\displaystyle=c_{\uparrow}({\bf
r}_{0},p_{y})c_{\downarrow}({\bf r}_{0},p_{y})|\Phi^{(ph)}\rangle$ (18)
$\displaystyle=\varphi_{{\bf R}_{6}}({\bf r}_{0},p_{y})\varphi_{{\bf
R}_{0}}({\bf r}_{0},p_{y})(b_{{\bf R}_{6},\uparrow}b_{{\bf
R}_{0},\downarrow}-b_{{\bf R}_{6},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle+\varphi_{{\bf
R}_{1}}({\bf r}_{0},p_{y})\varphi_{{\bf R}_{0}}({\bf r}_{0},p_{y})(b_{{\bf
R}_{1},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{1},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle\implies(b_{{\bf
R}_{6},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{6},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle\hskip 22.76228pt+(b_{{\bf
R}_{1},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{1},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle=0,$

using Eq. $\eqref{eq:fcond}$ and Eq. $\eqref{eq:pxcond}$. While this condition
is satisfied for states that maximize spin (have spin exchange symmetry)
between single-particle states at ${\bf R}_{6}$ and ${\bf R}_{0}$ as well as
between those at ${\bf R}_{1}$ and ${\bf R}_{0}$, it is not the case that
every state $|\Phi\rangle$ satisfying Eq. (18) must have such spin exchange
symmetry. Thus, the percolation representation is not strictly valid when the
only interaction is $H_{U}^{\prime}$.

The final condition we consider is $H_{J}|\Phi\rangle=0$. This condition can
be written using only annihilation operators in $H_{J}$ by writing
$H_{J}=\frac{J}{2}\sum_{{\bf r}\in\Lambda_{p}}n_{{\bf r},S=0}$, where $n_{{\bf
r},S=0}=(p_{{\bf r},\hat{y},\downarrow}^{\dagger}p_{{\bf
r},\hat{x},\uparrow}^{\dagger}-p_{{\bf r},\hat{y},\uparrow}^{\dagger}p_{{\bf
r},\hat{x},\downarrow}^{\dagger})(p_{{\bf r},\hat{x},\uparrow}p_{{\bf
r},\hat{y},\downarrow}-p_{{\bf r},\hat{x},\downarrow}p_{{\bf
r},\hat{y},\uparrow})\equiv c^{\dagger}_{{\bf r},S=0}c_{{\bf r},S=0}$ is an
operator that counts whether there is a spin singlet between the $p_{x}$ and
$p_{y}$ orbitals at site ${\bf r}$. The derivation of this operator identity
is shown in Appendix B. In this form, $H_{J}$ is clearly a sum of positive
semidefinite operators, and $H_{J}|\Phi\rangle=0$ if and only if $c_{{\bf
r},S=0}|\Phi\rangle=0$ for every ${\bf r}\in\Lambda_{p}$.

In terms of localized state operators, the zero-interaction-energy condition
for Hund’s coupling at site ${\bf r}_{0}$ is

$\displaystyle 0$ $\displaystyle=\left[c_{\uparrow}({\bf
r}_{0},p_{x})c_{\downarrow}({\bf r}_{0},p_{y})-c_{\downarrow}({\bf
r}_{0},p_{x})c_{\uparrow}({\bf r}_{0},p_{y})\right]|\Phi^{(ph)}\rangle$ (19)
$\displaystyle=\varphi_{{\bf R}_{6}}({\bf r}_{0},p_{x})\varphi_{{\bf
R}_{0}}({\bf r}_{0},p_{y})(b_{{\bf R}_{6},\uparrow}b_{{\bf
R}_{0},\downarrow}-b_{{\bf R}_{6},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle+\varphi_{{\bf
R}_{1}}({\bf r}_{0},p_{x})\varphi_{{\bf R}_{0}}({\bf r}_{0},p_{y})(b_{{\bf
R}_{1},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{1},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle\implies(b_{{\bf
R}_{6},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{6},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle$ $\displaystyle\hskip 22.76228pt-(b_{{\bf
R}_{1},\uparrow}b_{{\bf R}_{0},\downarrow}-b_{{\bf R}_{1},\downarrow}b_{{\bf
R}_{0},\uparrow})|\Phi^{(ph)}\rangle,$

using the no-double-occupancy condition from Eq. (16) and the spin triplet
condition between ${\bf R}_{1}$ and ${\bf R}_{6}$ from Eq. (17). Taking the
sum and difference of the conditions in Eqs. (18) and (19) gives spin triplet
conditions between ${\bf R}_{0}$ and ${\bf R}_{1}$ and between ${\bf R}_{0}$
and ${\bf R}_{6}$. In other words, the spin degree of freedom must be fully
symmetrized among any overlapping states at ${\bf r}_{0}$. This argument holds
for any choice of ${\bf R}_{0}$, meaning it applies at any site in the same
honeycomb sublattice as ${\bf r}_{0}$. In fact, as can be seen by considering
the left-most site on the ${\bf R}_{0}$-centered honeycomb, ${\bf
r}_{0}^{\prime}={\bf R}_{0}+{\bf u}_{4}$, these spin symmetrization conditions
apply both sublattices of the $p$-orbital honeycomb lattice $\Lambda_{p}$. If
spin must be symmetrized between any two overlapping localized states, a
cluster of localized states will be fully spin symmetrized, since any two
localized states in a cluster can be connected by a path of overlapping
localized states each adjacent pair of which must have symmetrized spin. Thus,
clusters of localized states have maximum total spin
$S_{C_{k},tot}=\frac{|C_{k}|}{2}$ and ground states for
$N_{e}\leq|\Lambda_{f}|$ can be written in the form of Eq. (11).

Theorem 2 follows immediately by particle-hole transformation. In particular,
since Theorem 1 requires the lowest-energy flat band to be at most half
filled, Theorem 2 requires the highest-energy flat band to be at least half
filled. The clusters in Eq. (12) are then connected sets of singly-occupied
localized states surrounded by a doubly-occupied background.

In order to interpret Theorem 2 as a percolation representation, note that
when the highest-energy flat band is exactly half filled,
$|N_{e}|=4|\Lambda_{p}|+|\Lambda_{f}|$, there is a single cluster spanning the
system, and all ground states have total spin $\frac{|\Lambda_{f}|}{2}$. When
the system is close to fully filled, clusters are small and it is easy to find
combinations of basis states in Eq. (12) with low total spin. As the highest-
energy flat band approaches half filling, the ground state space becomes
dominated by states with a large cluster spanning the system and carrying
large spin. There is thus a paramagnetic-ferromagnetic transition as the
system approaches half filling of the top band from above in the sense that
sufficiently close to half filling, the ground state space is dominated by
states with macroscopic spin.

The percolation transition can be found through Monte Carlo simulation [17].
Since the localized single-particle states can be labeled by their $f$ orbital
sites, the percolation picture involves filled or empty sites on the
triangular lattice $\Lambda_{f}$, with adjacent filled sites on the triangular
lattice belonging to the same cluster. This can be thought of as a percolation
of hole states with $N_{h}$ holes in the otherwise fully-filled system. In
this percolation picture, physical quantities are computed by averaging over
the degenerate ground state space. In particular, the average expected value
of the spin in basis states with fixed $A$ is

$\displaystyle S_{A}^{2}$
$\displaystyle\equiv\frac{1}{W(A)}\sum_{\\{m_{k}\\}}\langle\Phi_{A,\\{m_{k}\\}}|S^{2}|\Phi_{A,\\{m_{k}\\}}\rangle$
(20)
$\displaystyle=\sum_{k=1}^{n}\frac{|C_{k}|}{2}\left(\frac{|C_{k}|}{2}+1\right),$

which depends only on the sizes of the $n$ clusters $C_{k}\subset A$. [1, 17].
This can be interpreted as a geometric weighted percolation picture on a
triangular lattice $\Lambda_{f}$ with filled sites specified by
$A\subset\Lambda_{f}$ and each geometric configuration having the weighting
factor $W(A)=\sum_{\\{m_{k}\\}}1=\prod_{k=1}^{n}(|C_{k}|+1)$. The percolation
transition is signalled by the filling at which $\langle
S^{2}\rangle/S^{2}_{max}$ becomes nonzero in the thermodynamic limit, where
the maximum value $S^{2}_{max}=\frac{N_{h}}{2}\left(\frac{N_{h}}{2}+1\right)$
occurs for configurations with only a single cluster.

##  IV Conclusions

We have studied a multi-orbital model with two layers, one a honeycomb lattice
layer with $p_{x}$ and $p_{y}$ orbitals and the other a triangular lattice
layer of $f_{y(3x^{2}-y^{2})}$ orbitals with triangular lattice vertices
aligned with the centers of the $p$-orbital honeycomb plaquettes. For an
appropriate chemical potential difference between the two layers, the system
admits a flat band spanned by single-particle states localized to each $f$
orbital and the surrounding $p$-orbital plaquette. We proved that when this
flat band is the highest-energy band, the ground state space of the system
with intra-orbital Hubbard interactions and Hund’s coupling is spanned by
states with maximum-spin clusters when the flat band is at least half filled.
The completeness of this cluster state basis allows for a percolation
representation, where the system is paramagnetic far from and ferromagnetic
near half filling of the flat band. In particular, when the flat band is
exactly half filled, the ground states are fully spin polarized.

## Acknowledgments

We thank Tyrel McQueen for helpful discussion. This work is supported by the
NSF CAREER grant DMR-1848349 and in part by the Alfred P. Sloan Research
Fellowships under grant FG-2018-10971.

##  Appendix A Lower-spin Clusters in $p$-orbital System

In this appendix, we consider just the $p$-orbital Hamiltonian
$H=H_{K}^{p}+H_{U}+H_{J}$, with the intra-orbital Hubbard interactions in the
$p_{x}$ and $p_{y}$ orbitals as well as Hund’s coupling. By use of a simple
example, we show that the loop state basis $|\psi_{{\bf
R},\sigma}^{-(p)}\rangle$ does not admit a simple percolation picture, as
there are states where clusters do not maximize total spin. For consistency,
we use the notation from Sec. III.2 with the modification that there are no
$f$ orbital sites and $\Lambda_{f}$ is simply the set of honeycomb plaquette
labels. Thus, $\varphi_{\bf R}({\bf r},o_{\bf r})$ describes the component at
orbital $o_{\bf r}=p_{x},p_{y}$ of site ${\bf r}$ of a loop state on the
plaquette ${\bf R}$. All operators $a_{{\bf R},\sigma}$ and $b_{{\bf
R},\sigma}$ are defined similarly to in Sec. III.2 using the wavefunctions
$\varphi_{\bf R}({\bf r},o_{\bf r})$ defined only at sites ${\bf
r}\in\Lambda_{p}$.

Consider an arbitrary plaquette centered at ${\bf R}_{0}$ together with the
six surrounding plaquettes centered at ${\bf R}_{i}={\bf R}_{0}+{\bf w}_{i}$
with $i=1,\dots,6$. The nonzero components at site ${\bf r}_{1}={\bf
R}_{0}+{\bf u}_{1}$ are $\varphi_{{\bf R}_{1}}({\bf
r}_{1},p_{x})=-\frac{\sqrt{3}}{2}$ and $\varphi_{{\bf R}_{6}}({\bf
r}_{1},p_{x})=\frac{\sqrt{3}}{2}$ while the nonzero $p_{y}$ components are
$\varphi_{{\bf R}_{0}}({\bf r}_{1},p_{y})=1$, $\varphi_{{\bf R}_{1}}({\bf
r}_{1},p_{y})=-\frac{1}{2}$, and $\varphi_{{\bf R}_{6}}({\bf
r}_{1},p_{y})=-\frac{1}{2}$, discarding the $1/\sqrt{6}$ normalization factor.

The equivalents of the interaction conditions in Eq. (17), (18), and (19) in
this case are

$\displaystyle 0=$ $\displaystyle\bigg{[}b_{{\bf R}_{1}\uparrow}b_{{\bf
R}_{1}\downarrow}+b_{{\bf R}_{6}\uparrow}b_{{\bf R}_{6}\downarrow}$ (21)
$\displaystyle-(b_{{\bf R}_{1}\uparrow}b_{{\bf R}_{6}\downarrow}-b_{{\bf
R}_{1}\downarrow}b_{{\bf R}_{6}\uparrow})\bigg{]}|\Phi\rangle$ $\displaystyle
0=$ $\displaystyle\bigg{[}4b_{{\bf R}_{0}\uparrow}b_{{\bf
R}_{0}\downarrow}+b_{{\bf R}_{1}\uparrow}b_{{\bf R}_{1}\downarrow}+b_{{\bf
R}_{6}\uparrow}b_{{\bf R}_{6}\downarrow}$ $\displaystyle+(b_{{\bf
R}_{1}\uparrow}b_{{\bf R}_{6}\downarrow}-b_{{\bf R}_{1}\downarrow}b_{{\bf
R}_{6}\uparrow})-2(b_{{\bf R}_{0}\uparrow}b_{{\bf R}_{6}\downarrow}-b_{{\bf
R}_{0}\downarrow}b_{{\bf R}_{6}\uparrow})$ $\displaystyle-2(b_{{\bf
R}_{0}\uparrow}b_{{\bf R}_{1}\downarrow}-b_{{\bf R}_{0}\downarrow}b_{{\bf
R}_{1}\uparrow})\bigg{]}|\Phi\rangle$ $\displaystyle 0=$
$\displaystyle\bigg{[}b_{{\bf R}_{1}\uparrow}b_{{\bf R}_{1}\downarrow}-b_{{\bf
R}_{6}\uparrow}b_{{\bf R}_{6}\downarrow}+(b_{{\bf R}_{0}\uparrow}b_{{\bf
R}_{6}\downarrow}-b_{{\bf R}_{0}\downarrow}b_{{\bf R}_{6}\uparrow})$
$\displaystyle-(b_{{\bf R}_{0}\uparrow}b_{{\bf R}_{1}\downarrow}-b_{{\bf
R}_{0}\downarrow}b_{{\bf R}_{1}\uparrow})\bigg{]}|\Phi\rangle.$

These conditions reduce to the conditions in the main text when there are no
doubly occupied states, meaning states where clusters of loop states maximize
spin remain ground states. There are, however, linearly independent ground
states that do not maximize spin due to having a doubly-occupied loop state.
One such example is

$|\Phi^{\prime}_{{\bf R}_{0}}\rangle=a^{\dagger}_{{\bf
R}_{0}\uparrow}\sum_{i=0}^{6}a^{\dagger}_{{\bf R}_{i}\downarrow}|0\rangle.$
(22)

It can be verified that the state $|\Phi^{\prime}_{{\bf R}_{0}}\rangle$
satisfies the three zero-interaction-energy conditions on every site. The
conditions in Eq. (21) are satisfied by $|\Phi^{\prime}_{{\bf R}_{0}}\rangle$,
and the conditions for the remaining sites can be shown to hold as well and
follow from simple mappings. For example, the conditions at site ${\bf
r}_{2}={\bf R}_{0}+{\bf u}_{2}$ follow from replacing ${\bf R}_{0}\to{\bf
R}_{1},$ ${\bf R}_{1}\to{\bf R}_{0}$, and ${\bf R}_{6}\to{\bf R}_{2}$ in Eq.
(21).

The two-electron state $|\Phi^{\prime}_{{\bf R}_{0}}\rangle$ includes a
component where the state $\bm{\varphi}_{{\bf R}_{0}}$ is doubly occupied and
has total spin $\langle\Phi^{\prime}_{{\bf
R}_{0}}|S^{2}_{tot}|\Phi^{\prime}_{{\bf
R}_{0}}\rangle/\langle\Phi^{\prime}_{{\bf R}_{0}}|\Phi^{\prime}_{{\bf
R}_{0}}\rangle=\frac{6}{7}$. The percolation representation with maximum-spin
clusters thus does not hold in the loop basis for $H_{K}^{p}$.

Figure 3: Two-electron state $|\Phi^{\prime}_{{\bf R}_{0}}\rangle$, which
consists of a central filled spin-up loop surrounded by a superposition of
spin-down loops on the central and six surrounding plaquettes. This state does
not maximize spin despite avoiding the interaction energy.

##  Appendix B Hund’s Coupling Singlet Operator

The Hund’s coupling term $H_{J}$ is a sum of positive semidefinite operators
at each honeycomb site,

$\displaystyle H_{J}$ $\displaystyle=\sum_{{\bf r}\in\Lambda_{p}}h_{J}({\bf
r}),$ (23) $\displaystyle h_{J}({\bf r})$
$\displaystyle\equiv-J\left(\mathbf{S}_{{\bf r},p_{x}}\cdot\mathbf{S}_{{\bf
r},p_{y}}-\frac{1}{4}n_{{\bf r},p_{x}}n_{{\bf r},p_{y}}\right)$

where $S^{i}_{{\bf
r},p_{x/y}}=\frac{1}{2}\sum_{\mu,\nu=\uparrow,\downarrow}p^{\dagger}_{{\bf
r},\hat{x}/\hat{y},\mu}\sigma^{i}_{\mu\nu}p_{{\bf r},\hat{x}/\hat{y},\nu}$ and
$n_{{\bf r},p_{x/y}}=\sum_{\mu=\uparrow,\downarrow}p^{\dagger}_{{\bf
r},\hat{x}/\hat{y},\mu}p_{{\bf r},\hat{x}/\hat{y},\mu}$. $h_{J}({\bf r})$
takes its minimum eigenvalue of $0$ when the $p_{x}$ and $p_{y}$ orbitals at
site ${\bf r}$ are singly occupied and form a spin triplet state, which can be
seen explicitly by using the Pauli matrix completeness identity
$\sum_{i=1}^{3}\sigma^{i}_{\alpha\beta}\sigma^{i}_{\mu\nu}=2\delta_{\alpha\nu}\delta_{\beta\mu}-\delta_{\alpha\beta}\delta_{\mu\nu}$
to write

$\displaystyle h_{J}({\bf r})$ $\displaystyle=-\frac{J}{4}\bigg{[}-n_{{\bf
r},p_{x}}n_{{\bf r},p_{y}}+$ (24)
$\displaystyle\sum_{\mu,\nu;\alpha,\beta=\uparrow,\downarrow}p^{\dagger}_{{\bf
r},\hat{x},\mu}p_{{\bf r},\hat{x},\nu}p^{\dagger}_{{\bf
r},\hat{y},\alpha}p_{{\bf
r},\hat{y},\beta}(2\delta_{\alpha\nu}\delta_{\beta\mu}-\delta_{\alpha\beta}\delta_{\mu\nu})\bigg{]}$
$\displaystyle=\frac{J}{2}\bigg{(}n_{{\bf r},p_{x},\uparrow}n_{{\bf
r},p_{y},\downarrow}+n_{{\bf r},p_{x},\downarrow}n_{{\bf r},p_{y},\uparrow}$
$\displaystyle-p^{\dagger}_{{\bf r},\hat{x},\uparrow}p_{{\bf
r},\hat{x},\downarrow}p^{\dagger}_{{\bf r},\hat{y},\downarrow}p_{{\bf
r},\hat{y},\uparrow}-p^{\dagger}_{{\bf r},\hat{x},\downarrow}p_{{\bf
r},\hat{x},\uparrow}p^{\dagger}_{{\bf r},\hat{y},\uparrow}p_{{\bf
r},\hat{y},\downarrow}\bigg{)}$ $\displaystyle=-\frac{J}{2}\bigg{(}p_{{\bf
r},\hat{x},\uparrow}^{\dagger}p_{{\bf r},\hat{y},\downarrow}^{\dagger}p_{{\bf
r},\hat{x},\uparrow}p_{{\bf r},\hat{y},\downarrow}+p_{{\bf
r},\hat{x},\downarrow}^{\dagger}p_{{\bf r},\hat{y},\uparrow}^{\dagger}p_{{\bf
r},\hat{x},\downarrow}p_{{\bf r},\hat{y},\uparrow}$ $\displaystyle-p_{{\bf
r},\hat{x},\uparrow}^{\dagger}p_{{\bf r},\hat{y},\downarrow}^{\dagger}p_{{\bf
r},\hat{x},\downarrow}p_{{\bf r},\hat{y},\uparrow}-p_{{\bf
r},\hat{x},\downarrow}^{\dagger}p_{{\bf r},\hat{y},\uparrow}^{\dagger}p_{{\bf
r},\hat{x},\uparrow}p_{{\bf r},\hat{y},\downarrow}\bigg{)}$
$\displaystyle=\frac{J}{2}(p_{{\bf r},\hat{y},\downarrow}^{\dagger}p_{{\bf
r},\hat{x},\uparrow}^{\dagger}-p_{{\bf r},\hat{y},\uparrow}^{\dagger}p_{{\bf
r},\hat{x},\downarrow}^{\dagger})(p_{{\bf r},\hat{x},\uparrow}p_{{\bf
r},\hat{y},\downarrow}-p_{{\bf r},\hat{x},\downarrow}p_{{\bf
r},\hat{y},\uparrow})$ $\displaystyle\equiv\frac{J}{2}n_{{\bf r},S=0}.$

The Hund’s coupling term can thus be written in terms of a sum of singlet
number operators.

##  Appendix C Particle-Hole Transformation and Total Spin

We show here for completeness that the particle-hole transformation and total
spin operators commute. For notational simplicity, consider electron operators
$c_{i\sigma}$, essentially absorbing both the site and orbital indices into
the single index $i$. The total spin operator can be written

$\displaystyle S^{2}_{tot}$ $\displaystyle=(\sum_{i}\mathbf{S}_{i})^{2}$ (25)
$\displaystyle=\sum_{i,j}\sum_{\mu\nu\alpha\beta=\uparrow,\downarrow}c^{\dagger}_{i\mu}c_{i\nu}c^{\dagger}_{j\alpha}c_{j\beta}\bm{\sigma}_{\mu\nu}\cdot\bm{\sigma}_{\alpha\beta}$
$\displaystyle=\sum_{i,j}\sum_{\mu\nu\alpha\beta=\uparrow,\downarrow}(-c_{i\nu}c^{\dagger}_{i\mu}+\delta_{\mu\nu})(-c_{j\beta}c^{\dagger}_{j\alpha}+\delta_{\alpha\beta})\bm{\sigma}_{\mu\nu}\cdot\bm{\sigma}_{\alpha\beta}$
$\displaystyle=\sum_{i,j}\sum_{\mu\nu\alpha\beta=\uparrow,\downarrow}c_{i\nu}c^{\dagger}_{i\mu}c_{j\beta}c^{\dagger}_{j\alpha}\bm{\sigma}_{\mu\nu}\cdot\bm{\sigma}_{\alpha\beta}$
$\displaystyle=U^{(ph)\dagger}S^{2}_{tot}U^{(ph)},$

which follows from the tracelessness of $\sigma^{i}$ and the fact that
$\bm{\sigma}_{\mu\nu}\cdot\bm{\sigma}_{\alpha\beta}=\bm{\sigma}_{\nu\mu}\cdot\bm{\sigma}_{\beta\alpha}$.
Thus, total spin is preserved by the particle-hole transformation.

## References

  * Mielke and Tasaki [1993] A. Mielke and H. Tasaki, Commun. Math. Phys. 158, 341 (1993).
  * Mielke [1991a] A. Mielke, J. Phys. A: Math. Gen. 24, L73 (1991a).
  * Wu  _et al._ [2007] C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Phys. Rev. Lett. 99, 070401 (2007).
  * Tang  _et al._ [2011] E. Tang, J.-W. Mei, and X.-G. Wen, Phys. Rev. Lett. 106, 236802 (2011).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 (2011).
  * Wang  _et al._ [2013] Z. F. Wang, N. Su, and F. Liu, Nano Letters 13, 2842 (2013), pMID: 23678979\.
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Bergman  _et al._ [2008] D. L. Bergman, C. Wu, and L. Balents, Phys. Rev. B 78, 125104 (2008).
  * Creutz [2001] M. Creutz, Rev. Mod. Phys. 73, 119 (2001).
  * Zurita  _et al._ [2020] J. Zurita, C. E. Creffield, and G. Platero, Advanced Quantum Technologies 3, 1900105 (2020).
  * Chiu  _et al._ [2020] C. S. Chiu, D.-S. Ma, Z.-D. Song, B. A. Bernevig, and A. A. Houck, (2020), arXiv:2010.11953 [cond-mat.mes-hall] .
  * Tasaki [1992] H. Tasaki, Phys. Rev. Lett. 69, 1608 (1992).
  * Tanaka [2020] A. Tanaka, Journal of Statistical Physics 181, 897 (2020).
  * Mielke [1991b] A. Mielke, J. Phys. A: Math. Gen. 24, 3311 (1991b).
  * Mielke [1992] A. Mielke, J. Phys. A: Math. Gen. 25, 4335 (1992).
  * Maksymenko  _et al._ [2012] M. Maksymenko, A. Honecker, R. Moessner, J. Richter, and O. Derzhko, Phys. Rev. Lett. 109, 096404 (2012).
  * Liu  _et al._ [2019] R. Liu, W. Nie, and W. Zhang, Science Bulletin 64, 1490 (2019).
  * Wu and Das Sarma [2008] C. Wu and S. Das Sarma, Phys. Rev. B 77, 235107 (2008).
  * Zhang  _et al._ [2010] S. Zhang, H.-h. Hung, and C. Wu, Phys. Rev. A 82, 053618 (2010).


11institutetext: Sberbank, Moscow 117997, Russia, 22institutetext: Artificial
Intelligence Research Institute, Moscow 105064, Russia, 33institutetext:
National University of Science and Technology “MISIS”, Moscow 119991, Russia
33email<EMAIL_ADDRESS>%**** ̵main.tex ̵Line ̵75 ̵****

# Logistics, Graphs, and Transformers:
Towards improving Travel Time Estimation

Natalia Semenova(✉) 1122 Vadim Porvatov 1133 Vladislav Tishin 1133
Artyom Sosedka 1133 Vladislav Zamkovoy 11

###### Abstract

The problem of travel time estimation is widely considered as the fundamental
challenge of modern logistics. The complex nature of interconnections between
spatial aspects of roads and temporal dynamics of ground transport still
preserves an area to experiment with. However, the total volume of currently
accumulated data encourages the construction of the learning models which have
the perspective to significantly outperform earlier solutions. In order to
address the problems of travel time estimation, we propose a new method based
on transformer architecture – TransTTE.

###### keywords:

Graph Embedding, Travel Time Estimation, Geospatial Linked Data.

##  1 Introduction

As long as ground transport dramatically increases its quantitative presence
in the cities, traffic management becomes more complex and hence less
predictable for drivers. In order to handle the escalation of such a negative
trend, it is important to effectively estimate the essential parameters
describing traffic dynamics. One of the most important values among such
quantities is the estimated time of arrival (ETA) which could be considered as
the expected time expenditure for a trip between two locations. Accurate
travel time estimation (TTE) is mostly challenging for cars due to the
presence of the extensive limitations induced by the road network structure.
These aspects of urban traffic ordinary require special spatio-temporal
methods implementation to be handled. The contributions of our work are the
following:

  1. 1.

We proposed TransTTE model that could utilize the spatio-temporal dependencies
and explored the capabilities of the transformer model application in the
domain of TTE via comparison with several baselines. Computational experiments
allow us to conclude that the TransTTE architecture achieves competitive
performance.

  2. 2.

We published a new dataset related to the city of Omsk.

  3. 3.

For the needs of demonstration, we developed a web service based on the
TransTTE model.

The application is available at http://transtte.online and the code could be
accessed from https://github.com/Vloods/TransTTE˙demo.

##  2 Related Work

Generally, the TTE methods could be divided into two categories related to the
different approaches of ETA computing. The first one is based on the
extraction of the total traveled time regarding each segment of the path [1].
Such models do not capture any global properties of the path which is
explicitly reflected in their performance. The second class of methods
utilizes the corresponding trip path as a whole [2]. The results achieved by
these approaches defined the mainstream of current research in the domain of
TTE and hence are widely established in the experiments section.

##  3 Framework design

In this section, we discuss the aspects of model design and deployment as a
part of a developed web service.

Task. Given an origin, destination, and departure time, our goal is to
estimate the duration using the set of historical trip dataset X and the
underlying road network G.

Model. The transformer architecture has recently become a prevalent approach
in many domains, such as natural language processing and computer vision. Yet,
it has not achieved competitive performance on popular leaderboards of graph-
level prediction compared to mainstream GNN variants. Therefore, it remains a
question of how transformers could perform well for graph representation
learning in the TTE task.

One of the graph-oriented aspects in the Graphormer architecture [8] is a
_centrality encoding_ which assigns each node two real-valued embedding
vectors according to its indegree and outdegree
$h_{i}^{(0)}=x_{i}+z^{-}+z^{+}$, where $z^{-},z^{+}\in\mathbb{R}^{d}$ are
learnable embedding vectors specified by the indegree $\text{deg}^{-}(v_{i})$
and outdegree $\text{deg}^{+}(v_{i})$ respectively. Along with centrality
encoding, _spatial encoding_ is used to capture the structural relation via
function $\phi(v_{i},v_{j}):V\times V\to\mathbb{R}$ which measures the spatial
relation between nodes $v_{i}$ and $v_{j}$ of road network $G$. Original
choice of such function $\phi(v_{i},v_{j})$ is the shortest distance between
$v_{i}$ and $v_{j}$ which further serves as a bias term in the self-attention
module
$A_{ij}=\frac{(h_{i}W_{Q})(h_{j}W_{K})^{T}}{\sqrt{d}}+b_{\phi(v_{i},v_{j})}$,
where $b_{\phi(v_{i},v_{j})}$ is a learnable scalar indexed by
$\phi(v_{i},v_{j})$ which is shared across all layers.

Data. In this paper, we use the dataset related to the road networks of Abakan
and Omsk111The full data could be requested from<EMAIL_ADDRESS>which
have different scales and road topology [6]. The datasets are collected in a
monthly period starting from December 1, 2020. The datasets consist of a road
network with corresponding to its segments features and trip part. The
presence of the noisy data among the trip part of the datasets encourages us
to apply filtering regarding the rebuild count feature, maximum/minimum
length, and the total time of the trips.

Web application. For the needs of model demonstration, we deployed our TTE
service in the platform of the Yandex Maps project, Figure 1. The interface
allows choosing between two cities and three types of routes regarding demands
of a user.

Figure 1: Application interface.

##  4 Results

In addition to the computing TTE for the shortest route, the framework also
evaluates the routes based on the metrics of picturesqueness and historicity.
Using the OpenStreetMap API, we managed to parse information about the
location of historical, cultural and natural objects. The number of certain
objects is used as the road segments weights for Dijkstra’s algorithm
$W_{i}=\frac{1}{1+C_{r}}$, where $W_{i}$ is the weight for the $i$’th segment
of the road, $C_{r}$ is the number of objects within a radius $r$ of a road
segment.

We reimplemented Graphormer architecture to accelerate the training process
and consider the peculiar properties of road trips. Due to caching spatial
encoding values, we were able to speed up training by almost 10 times. Several
baselines for TTE task were also implemented to verify the effectiveness of
the proposed model, Table 1. We considered results made by gradient boosted
decision trees along with the more sophisticated pipelines. WDR [7] uses
generalized linear model and LSTM together to compute travel time. MURAT [3],
in its turn, produces unsupervised representations due to DeepWalk [5] and
applies residual feedforward blocks to predict travel time and distance.

The best result was achieved on Graphormer${}_{\small\textsc{Slim}}$
($L=12,d=80$) version with reduced size of dimension. As the optimizer we used
AdamW [4]. Experiments were done with 5 Tesla V100 GPUs and 460 Gb of RAM. The
training time of the different configurations of TransTTE lies in the interval
from 2.5 hours up to 5 hours which is smaller than in case of WDR (7 hours)
and MURAT (5.5 hours).

##  5 Conclusion and Outlook

Dataset | Omsk | Abakan
---|---|---
Split | Train | Test | Train | Test
| MAE | RMSE | MAE | RMSE | MAE | RMSE | MAE | RMSE
GBDT | 403.921 | 582.011 | 408.644 | 573.559 | 244.119 | 449.250 | 248.862 | 399.534
MURAT | 279.616 | 438.228 | 286.491 | 443.397 | 179.037 | 285.003 | 185.153 | 286.934
WDR | 311.581 | 440.511 | 336.756 | 487.876 | 173.684 | 285.132 | 182.296 | 293.551
TransTTE | 101.381 | 387.241 | 105.464 | 261.103 | 81.048 | 285.032 | 83.616 | 168.421
Table 1: Evaluation of different pipelines and comparison with proposed method

In this paper, we proposed the new transformer-based approach to the computing
of ETA and explored its performance. The experiments revealed the perspective
of graph transformer utilization in the travel time estimation. In the
upcoming studies, we want to extend the current transformer architecture by
virtue of extra road network features and more precise work with the temporal
aspect of road trips. Future research should be devoted to the development of
a joint TransTEE model which could compute travel time for rides on different
city networks indeed.

## References

  * [1] Asghari, M., Emrich, T., Demiryurek, U., Shahabi, C.: Probabilistic estimation of link travel times in dynamic road networks. pp. 1–10 (11 2015)
  * [2] Jin, G., Wang, M., Zhang, J., Sha, H., Huang, J.: Stgnn-tte: Travel time estimation via spatial–temporal graph neural network. Future Generation Computer Systems 126, 70–81 (2022), https://www.sciencedirect.com/science/article/pii/S0167739X21002740
  * [3] Li, Y., Fu, K., Wang, Z., Shahabi, C., Ye, J., Liu, Y.: Multi-task representation learning for travel time estimation. In: International Conference on Knowledge Discovery and Data Mining (KDD ’18) (2018)
  * [4] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization (2019)
  * [5] Perozzi, B., Al-Rfou, R., Skiena, S.: Deepwalk: Online learning of social representations. In: Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 701–710 (2014)
  * [6] Porvatov, V., Semenova, N., Chertok, A.: Hybrid graph embedding techniques in estimated time of arrival task. In: Benito, R.M., Cherifi, C., Cherifi, H., Moro, E., Rocha, L.M., Sales-Pardo, M. (eds.) Complex Networks & Their Applications X. pp. 575–586. Springer International Publishing, Cham (2022)
  * [7] Wang, Z., Fu, K., Ye, J.: Learning to estimate the travel time. In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 858–866. KDD ’18, Association for Computing Machinery, New York, NY, USA (2018), https://doi.org/10.1145/3219819.3219900
  * [8] Ying, C., Cai, T., Luo, S., Zheng, S., Ke, G., He, D., Shen, Y., Liu, T.Y.: Do transformers really perform bad for graph representation? arXiv preprint arXiv:2106.05234 (2021)


# Ancillary Services Provision by Cross-Voltage-Level Power Flow Control using
Flexibility Regions ††thanks: The authors gratefully acknowledge funding by
the German Federal Ministry of Education and Research (BMBF) within the
Kopernikus project ENSURE ‘New ENergy grid StructURes for the German
Energiewende’.

Christian Holger Nerowski , Zongjun Li, Christian Rehtanz  0009-0001-1297-9936
0000-0002-8134-6841 Institute of Energy Systems, Energy Efficiency and Energy
Economics
TU Dortmund University - Dortmund, Germany
E-Mail<EMAIL_ADDRESS>
###### Abstract

The large-scale integration of distributed renewable energy sources into the
electricity grid requires the investigation of new methods to ensure
stability. For example, Active Distribution Networks (ADNs) can be used at
(sub-) transmission levels for emergency operation, provided robust and
efficient control is available. This paper investigates the use of Feasible
Operating Regions (FORs) and Flexibility Regions (FRs) for Cross-Voltage-Level
Power Flow Control (CPFC). The enhancement of network stability due to the
provision of ancillary services is illustrated, as is the need for
strengthened cooperation between Transmission (TSOs) and Distribution System
Operators (DSOs). Optimal power flow methods are considered, focusing on
computational advances through PieceWise Linearization (PWL) and convex
relaxation techniques aiming to speed up runtime while keeping high accuracy.
To illustrate the algorithms’ benefits and drawbacks, they are analyzed using
exemplary medium voltage grids.

###### Index Terms:

Active Distribution Network, Feasible Operating Region, Flexibility Region,
Optimal Power Flow, Piecewise Linearization, Second Order Cone Programming

##  I Introduction

With the energy transition, the restructuring of the electrical energy system
continues to progress. The shift from centralized to decentralized units in
electricity generation will lead to fewer controllable assets in transmission
grids, but in turn to more at distribution levels [1]. In future, the
installed capacity of renewable energies will exceed peak load demand in
Europe notably [2]. For this reason, situations will become frequent when
decentralized controllable assets have to ensure stability, particularly of
voltages and frequency. Therefore, network operation must adapt to these new
conditions [3].

To counteract the lack of control options in transmission grids, one approach
is to control renewable energy sources in clusters. Such groups of assets can
be allocated inside of Active Distribution Networks (ADNs) [4]. By controlling
inverter-based resources and on-load-tap changers, ADNs become able to provide
aforementioned support to higher levels. In principle, these ADNs depict
radial grids whose operating state is customizable using their renewable
energy sources. They rely heavily on information and communication technology
to guide their assets’ behavior. Said features enable the support of
superimposed levels through the Point of Common Coupling (PCC) as discussed in
numerous publications. Some approaches such as [5] are based on Model
Predictive Control or Optimal Power Flow (OPF) [6]. Others use PI-based
concepts such as [7], characterized by Cross-Voltage-Level Power Flow Control
(CPFC) enabling real-time management.

Due to increased prevalence of fast-acting inverter-based resources, providing
time-critical services will become possible. In principle, an emergency
requires the electrical system to be kept in balance within seconds to
minutes. The response time of a CPFC must therefore comply with hard limits,
as failure is not tolerated. As discussed in [3], reaction periods will change
significantly in the future, requiring rapid actions. According to [8], this
results in two possible CPFC-applications:

  * •

Frequency control (by providing active power control through subordinate ADNs
due to fast-acting assets)

  * •

Voltage regulation for transmission grids (by providing reactive power
regulation through subordinate ADNs)

This paper is structured as follows: In Section II, we define what
characterizes Flexibility Providing Units (FPUs), Feasible Operating Regions
(FORs) and Flexibility Regions (FRs). Henceforth, a conceptual framework for
TSO/DSO interaction is presented in Section III, as well as the ac-OPF for
benchmarking. In Section IV, various accelerated computation methods are
described and used in Section V to present test case results. Finally,
conclusions are drawn in Section VI.

##  II Modeling Flexibility Providing Units

As indicated, each asset needs to be managed in real-time to enable frequency
control and voltage regulation. These are often referred to as FPUs if they
are geographically distributed and remotely controllable via external signals.
As described in [9], there is a wide range of FPUs with various properties.
Depending on the ability to achieve certain active ($P$) and reactive ($Q$)
power combinations, five types are differentiated by their FORs. Here, any FOR
results from asset’s operating limits, which in turn are based on physical and
technical constraints. In short, FORs describe theoretically achievable
operational states of individual assets or entire network areas, reachable
only if no external time-dependent factors restrict them. However, if these
non-influenceable given factors, e.g. solar irradiation or wind strength, come
into play, they reduce available $PQ$-spaces resulting in FRs. FRs are thereby
subsets of FORs. As proposed in [10], we assume FPUs comprehend batteries (1),
controllable loads (1) with constant $\cos(\phi)$ (2), static synchronous
compensators (1) and synchronous generators (5) as shown in Fig. 1. Moreover,
photovoltaic systems (4) and wind turbines (3,4) are considered to be FPUs.

Figure 1: Classification of Flexibility Providing Units according to [10]
Figure 2: A Network’s Feasible Operating Regions and Flexibility Regions

Due to relevance for under- and overvoltage incidents [11], recent studies
have focused on calculating $P$\- and $Q$-flexibility. Furthermore, various
methods have been proposed to obtain FORs including those using computational
geometry, random sampling or OPFs. For example, Minkowski Sums handle generic
polygonal shapes to provide FORs as in [12]. However, crucial grid constraints
are not considered often causing problems. Applying random sampling is simple
but rather time consuming due to the vast number of calculations needed for
high-quality results. One such approach are Monte-Carlo simulations like in
[13], which define random setpoints for controllable assets before checking
solutions for the violation of limitations. Third, OPF methods have been
examined providing techniques for solving the problem efficiently. Advantages
are the reduction of runtime and improved effectiveness.

##  III Determination of Feasible Operating Regions and Flexibility Regions

###  III-A Setup

Within this paper, we consider different types of renewable energy sources to
be connected to a radial network as shown in Fig. 2. Each asset can be
represented by a characteristic FOR for which a time-dependent FR exists,
following [14]. Furthermore, such architecture is able to provide an
Interconnected Power Flow (IPF) at the PCC caused by aggregating assets’
current setpoints. Consequentially, Transmission (TSOs) and Distribution
System Operators (DSOs) must be coordinated to ensure ancillary services being
activated securely. Hence, we propose a conceptual framework where
responsibilities are managed upfront as illustrated in Fig. 3.

Figure 3: Conceptual Framework of Coordination Processes

Similarly to [15], we acknowledge the need for an ADN operator. Initially, it
must calculate the FOR once for a specific network topology and pool of
assets. If nothing changes in either aspect, the FOR remains the same
throughout upcoming processes. Following the FOR, FR calculations must be
performed with a pre-defined granularity, e.g. once per minute. The ADN
operator must then transmit obtained information to its TSO, where data of
other connected ADNs can be summarized. Depending on needs, the TSO optimizes
setpoints of individual ADNs, which again are to be addressed at their PCC as
a new IPF like proposed in [4]. By transmitting $P^{\text{ref}}$ and
$Q^{\text{ref}}$ back to each ADN operator, the TSO’s optimization and control
process concludes. ADN operators are ultimately responsible for breaking
target values down to each individual system handled. In the end, by utilizing
assets’ flexibility, new operating points are achieved without violating
constraints.

###  III-B Non-linear OPF exploration

Using ac-OPF to calculate FOR boundary points to determine the capabilities of
given networks requires formulating the problem with a set of linear/non-
linear equality and inequality constraints. In the course of which we use an
objective function (1) that includes $\alpha,\beta\in\\{-1,0,1\\}$ as weighing
factors like proposed in [16] to define the search direction.

$\text{minimize}\ f(x)=\alpha\cdot P_{\text{IPF}}(x)+\beta\cdot
Q_{\text{IPF}}(x)$ (1)

The power flow equations and operating constraints are represented by (2)-(9).
Let $\mathcal{N}$ be the set of buses, $\mathcal{E}$ the set of lines and
$\mathcal{G}$ the subset of buses with renewable energy sources connected. The
complex power $s_{i}=p_{i}+j\cdot q_{i}$, which is fed in at each bus
$i\in\mathcal{N}$, can then be expressed by

$\displaystyle
p_{i}(V,\theta)=|V_{i}|\sum_{j=1}^{n}y_{ij}|V_{j}|\cos(\theta_{ij}-\delta_{ij}),$
(2) $\displaystyle
q_{i}(V,\theta)=|V_{i}|\sum_{j=1}^{n}y_{ij}|V_{j}|\sin(\theta_{ij}-\delta_{ij}),$
(3)

where $n$ is the number of buses. Both equations depend on complex nodal
voltages $V_{i}=|V_{i}|\phase[2]{\theta_{i}}$, voltage angle differences
$\theta_{ij}$ and complex admittances $y_{ij}=|y_{ij}|\phase[2]{\delta_{ij}}$
of each line $(i,j)\in\mathcal{E}$. The nodal power balance in regular steady-
state operation is furthermore described by

$\displaystyle\Delta
p_{i}=-p_{d,i}+\sum_{g\in\mathcal{G}}p_{g,i}+p_{i}(V,\theta),$ (4)
$\displaystyle\Delta
q_{i}=-q_{d,i}+\sum_{g\in\mathcal{G}}q_{g,i}+q_{i}(V,\theta),$ (5)

where $p_{g,i},q_{g,i}$ represent generation and $p_{d,i},q_{d,i}$ consumption
for each bus $i\in\mathcal{N}$. These equations show the difference between
expected and calculated complex power at bus $i$.

The load flow problem aims to solve non-linear power flow equations, meaning
complex bus voltages $V_{i}$ and angles ${\theta_{i}}$ that solve $\Delta
p_{i}\approx 0$ and $\Delta q_{i}\approx 0$. For each line
$(i,j)\in\mathcal{E}$, we apply $V_{i}$ and ${\theta_{i}}$ to

$\displaystyle P_{ij}=$
$\displaystyle+g_{ij}\left(|V_{i}|^{2}-|V_{i}||V_{j}|\cos(\theta_{ij})\right)$
(6) $\displaystyle-b_{ij}|V_{i}||V_{j}|\sin(\theta_{ij}),$ $\displaystyle
Q_{ij}=$ $\displaystyle-
b_{ij}\left(|V_{i}|^{2}-|V_{i}||V_{j}|\cos(\theta_{ij})\right)$ (7)
$\displaystyle-g_{ij}|V_{i}||V_{j}|\sin(\theta_{ij}),$

where $g_{ij},b_{ij}$ are line series admittances and susceptances. This
permits the computation of line flows through all lines. The IPF, which is the
power flow through the slack bus, is included in those equations. Here, for
each line $(i,j)\in\mathcal{E}$, line capacity limits $S_{ij,\text{max}}$ are
enforced by

$P_{ij}^{2}+Q_{ij}^{2}\leq S_{ij,\text{max}}^{2}$ (8)

and voltage boundaries for each bus $i\in\mathcal{N}$ by

$V_{i,\text{min}}\leq V_{i}\leq V_{i,\text{max}}.$ (9)

##  IV Accelerated calculation methods

###  IV-A Piecewise linearization

With the strict PieceWise Linearization (PWL) as introduced in [17], the non-
linear power flow equations are discretized into linear segments. By
introducing auxiliary variables, the original non-linear problem is
transformed into a mixed-integer linear programming problem. The linear
approximation of the cosine function is realized by dividing it into multiple
segments. In (10), (11), $v_{i}=|V_{i}|^{2}$ depicts the squared nodal voltage
magnitudes. The linear combination of non-linear $\cos(\theta_{ij})$, which is
divided by breakpoints $\cos(\theta_{ij,s})$ with $N_{s}$ segments for all
$(i,j)\in\mathcal{E}$, is displayed by
$\sum_{s=0}^{N_{s}}w_{s}\cos(\theta_{ij,s})$ with coefficients
$w_{ij,s}\in[0,1]$ of breakpoints.

$\displaystyle P_{ij}=$
$\displaystyle+g_{ij}\left(\frac{v_{i}}{2}-\frac{v_{j}}{2}-\sum_{s=0}^{N_{s}}w_{ij,s}\cos(\theta_{ij,s})+1\right)$
(10) $\displaystyle-b_{ij}\sum_{s=0}^{N_{s}}w_{ij,s}\theta_{ij,s}$
$\displaystyle Q_{ij}=$ $\displaystyle-
b_{ij}\left(\frac{v_{i}}{2}-\frac{v_{j}}{2}-\sum_{s=0}^{N_{s}}w_{ij,s}\cos(\theta_{ij,s})+1\right)$
(11) $\displaystyle-g_{ij}\sum_{s=0}^{N_{s}}w_{ij,s}\theta_{ij,s}$

Furthermore, (12)-(15) represent the piecewise linearized quadratic line
capacity limits as in [18] for all $(i,j)\in\mathcal{E}$ with

$\displaystyle\sqrt{2}S_{ij,\text{max}}$ $\displaystyle\geq
P_{ij}+Q_{ij}\geq-\sqrt{2}S_{ij,\text{max}},$ (12)
$\displaystyle\sqrt{2}S_{ij,\text{max}}$ $\displaystyle\geq
P_{ij}-Q_{ij}\geq-\sqrt{2}S_{ij,\text{max}},$ (13) $\displaystyle
S_{ij,\text{max}}$ $\displaystyle\geq P_{ij}\geq-S_{ij,\text{max}},$ (14)
$\displaystyle S_{ij,\text{max}}$ $\displaystyle\geq Q_{ij}\geq-
S_{ij,\text{max}}.$ (15)

Further information on the introduction of auxiliary variables can be found in
[17]. In addition to the strict PWL, a relaxed PWL has been implemented
according to [19] and is analyzed hereafter. Compared to strict PWL, the
relaxed PWL transforms the mixed-integer linear programming problem into a
linear programming problem. Details are discussed in [19].

###  IV-B The DistFlow model

In the DistFlow model, the non-convexity of the ac-OPF is addressed by convex
relaxation. Since ADNs are mostly radial networks, we ignore phase angles of
voltages and currents according to [20]. Moreover, relaxing quadratic equality
constraints to inequality constraints yields a second order cone program as
shown in [20]. $P_{jk}$ depicts the power transfer going from $j$ to $k$,
whereas the power transfer coming from $i$ to $j$ is described by $P_{ij}$.
The power balance equations at each bus $j\in\mathcal{N}$ after angle
relaxation can then be described by

$\sum_{k}P_{jk}-\sum_{i}(P_{ij}-r_{ij}l_{ij})=p_{j},$ (16)
$\sum_{k}Q_{jk}-\sum_{i}(Q_{ij}-x_{ij}l_{ij})=q_{j},$ (17)

where $l_{ij}=|I_{ij}|^{2}$ are squared line current magnitudes and
$r_{ij},x_{ij}$ respective line series resistances and reactances describing
losses. In addition, Ohm’s law after angle relaxation for each line
$(i,j)\in\mathcal{E}$ is described by

$\displaystyle v_{j}=$ $\displaystyle
v_{i}-2(r_{ij}P_{ij}+x_{ij}Q_{ij})+(r_{ij}^{2}+x_{ij}^{2})l_{ij}.$ (18)

Furthermore, the line flow equation for each line $(i,j)\in\mathcal{E}$ after
angle and convex relaxation is described by

$\left\lVert\begin{array}[]{c}2P_{ij}\\\ 2Q_{ij}\\\
l_{ij}-v_{i}\end{array}\right\rVert_{2}\leq l_{ij}+v_{i},$ (19)

where $\left\lVert x\right\rVert_{2}$ is the second-order norm, or Euclidean
distance. Since (19) depicts a second order constraint, the feasible set is
convex. Furthermore, because the objective function is linear, the DistFlow
model represents a conic optimization.

##  V Test System Results

(a) CIGRE MV

(b) IEEE 33-bus

Figure 4: Test grids

In this section, we present FORs calculated by the algorithms from Sections
III-B, IV-A and IV-B using two test grids. The CIGRE benchmark network is
taken from pandapower [21] and adapted with photovoltaic, wind generation and
batteries as seen in Fig. 4. Also, the IEEE 33-bus network is taken from
MATPOWER [22] and supplied with photovoltaic and batteries accordingly. All
results are compared against each other regarding precision and runtime. Here,
the non-linear OPF is solved using CasADi [23]. Furthermore, both PWLs and the
DistFlow model are solved using Gurobi Optimizer [24]. Each test is conducted
on a PC with an Intel Core i5-1135G7 @ 2.4 GHz 8 GB RAM.

Figure 5: CIGRE’s Feasible Operating Regions using various Algorithms

The resulting FORs presented in Fig. 5 are obtained after five iterations.
Each iteration creates new boundary points, that improve the results’
accuracy. Starting with four boundary points initially, the number increases
exponentially with the number of iterations. The algorithms stop as soon as
the distance between neighboring boundary points crosses a
$5\leavevmode\nobreak\ \%$ threshold compared to the maximum values of $P$ and
$Q$ respectively. Finally, we connect all neighboring points obtaining a
polygon, the convex hull of the set, that serves as the FOR.

Figure 6: Algorithms’ runtimes for CIGRE’s Feasible Operating Regions

We can see that all FORs result in deviations compared to the non-linear case.
Both PWL methods lead to considerably shifted FORs if they contain 11 segments
instead of 21. Interestingly, the DistFlow model leads to errors at higher
values for $Q$, while $P$ fits well. In Fig. 6, the algorithms’ runtimes are
shown. The results indicate that the DistFlow model is the fastest method
examined, followed by relaxed PWL and strict PWL. Depending on the number of
segments, the PWL methods’ runtimes change significantly.

Table I: Runtimes of CIGRE’s and IEEE’s Feasible Operating Regions Algorithm ADN’s FOR |  CIGRE [s] \+ (speed-up) |  IEEE [s] \+ (speed-up)
---|---|---
Non-Linear | 14.24 (-) | 39.98 (-)
PWLstrict_11 | 2.95 (4.83) | 7.28 (5.49)
PWLstrict_21 | 4.63 (3.08) | 9.09 (4.4)
PWLrelaxed_11 | 2.28 (6.25) | 6.78 (5.9)
PWLrelaxed_21 | 2.8 (5.09) | 6.72 (5.95)
DistFlow | 1.44 (9.89) | 1.78 (22.46)

In Fig. 6, the acceleration obtained by all models stands out compared to the
non-linear case. In specific, Table I shows the substantial time saving
capabilities of the DistFlow model. Concerning the calculation of FRs,
exemplary results are shown in Fig. 7, obtained by using time series data from
[21]. Here, a single day was split into 96 intervals of 15 minutes each,
displaying 10 exemplary FRs. The increasing output of photovoltaic systems
during the middle of the day can be seen.

Figure 7: CIGRE’s Flexibility Regions from Time Series Data

##  VI Conclusion and Outlook

In this paper, the applicability of various approaches to calculate FORs/FRs
within an existing adaptive control system was discussed. The goal was to
improve grid stability by providing ancillary services. For this purpose,
several calculation methods, from linearization to relaxation, were compared
against each other. The simulation results show that the calculation of
FORs/FRs to be used in a CPFC is possible within a reasonable time frame.
Therefore, it would be advantageous if these methods were utilized more
frequently. The performance of ADNs thereby essentially depends on available
resources, while the computation methods of FORs/FRs can be accelerated
significantly by linearization and relaxation.

In future work, the increasing complexity of greater networks should be
thematized, accompanied by the challenges of communication. In addition,
laboratory experiments and field tests are planned to extend the analysis
proposed in this paper. In particular, laboratory tests based on Power-
Hardware-in-the-Loop will be conducted to characterize the actual FORs/FRs on
real components under different scenarios.

## Acknowledgment

Special thanks to Prof. Dr. Costas Vournas from the National Technical
University of Athens for valuable discussions.

## References

  * [1] ENTSO-E  “Ten Year Network Development Plan 2020: Completing the Map”, 2020  URL: https://eepublicdownloads.entsoe.eu/tyndp-documents/IoSN2020/200810_IoSN2020mainreport_beforeconsultation.pdf
  * [2] ENTSO-E  “Ten Year Network Development Plan 2022: Scenario Report”, 2022  URL: https://2022.entsos-tyndp-scenarios.eu/wp-content/uploads/2022/04/TYNDP2022_Joint_Scenario_Full-Report-April-2022.pdf
  * [3] Federal Ministry for Economic Affairs and Climate Action  “System Stability Roadmap”, 2023  URL: https://www.bmwk.de/Redaktion/EN/Publikationen/Energie/system-stability-roadmap.pdf?__blob=publicationFile&v=5
  * [4] Daniel Mayorga Gonzalez et al.  “Determination of the Time-Dependent Flexibility of Active Distribution Networks to Control their TSO-DSO Interconnection Power Flow”  In _20th Power Systems Computation Conference (PSCC)_ IEEE, 2018  DOI: 10.23919/PSCC.2018.8442865
  * [5] Gustavo Valverde and Thierry van Cutsem  “Model Predictive Control of Voltages in Active Distribution Networks”  In _IEEE Transactions on Smart Grid_ 4.4, 2013  DOI: 10.1109/TSG.2013.2246199
  * [6] Adrien Saint-Pierre and Pierluigi Mancarella  “Active Distribution System Management: A Dual-Horizon Scheduling Framework for DSO/TSO Interface under Uncertainty”  In _IEEE Transactions on Smart Grid_ 8.5, 2017  DOI: 10.1109/TSG.2016.2518084
  * [7] Lena Robitzky et al.  “Impact of Active Distribution Networks on Voltage Stability of Electric Power Systems”  In _2017 IREP Bulk Power Systems Dynamics and Control Symposium_ , 2017  URL: https://www.researchgate.net/profile/lena-robitzky/publication/319527382_impact_of_active_distribution_networks_on_voltage_stability_of_electric_power_systems
  * [8] Jannik Zwartscholten and Christian Rehtanz  “Adaptive and Robust Cross-Voltage-Level Power Flow Control of Active Distribution Networks”  In _2022 IREP Bulk Power Systems Dynamics and Control Symposium_ , 2022  URL: http://arxiv.org/pdf/2207.05462
  * [9] Daniel Contreras, Ouafa Laribi, Manswet Banka and Krzysztof Rudion  “Assessing the Flexibility Provision of Microgrids in MV Distribution Grids”  AIM, 2018  DOI: 10.34890/153
  * [10] Daniel A. Contreras and Krzysztof Rudion  “Improved Assessment of the Flexibility Range of Distribution Grids using Linear Optimization”  In _20th Power Systems Computation Conference (PSCC)_ IEEE, 2018  DOI: 10.23919/pscc.2018.8442858
  * [11] Thierry van Cutsem and Costas Vournas  “Voltage Stability of Electric Power Systems”  Boston, MA: Springer US, 2008  DOI: 10.1007/978-0-387-75536-6
  * [12] Soumya Kundu, Karanjit Kalsi and Scott Backhaus  “Approximating Flexibility in Distributed Energy Resources: A Geometric Approach”  In _20th Power Systems Computation Conference (PSCC)_ IEEE, 2018  DOI: 10.23919/PSCC.2018.8442600
  * [13] Miguel Heleno et al.  “Estimation of the Flexibility Range in the Transmission-Distribution Boundary”  In _2015 IEEE Eindhoven PowerTech_ IEEE, 2015  DOI: 10.1109/PTC.2015.7232524
  * [14] Daniel A. Contreras and Krzysztof Rudion  “Impact of Grid Topology and Tap Position Changes on the Flexibility Provision from Distribution Grids”  In _New Businesses for Energy Transition_ IEEE, 2019  DOI: 10.1109/ISGTEurope.2019.8905479
  * [15] Giorgos Prionistis, Costas Vournas and Maria Vrakopoulou  “A Fast Method to Approximate the Flexibility Region of an Active Distribution Network in PQ Space”  In _2023 IEEE Belgrade PowerTech_ IEEE, 2023  DOI: 10.1109/PowerTech55446.2023.10202983
  * [16] Joao Silva et al.  “Estimating the Active and Reactive Power Flexibility Area at the TSO-DSO Interface”  In _IEEE Transactions on Power Systems_ 33.5, 2018  DOI: 10.1109/TPWRS.2018.2805765
  * [17] Mang Jiang, Qinglai Guo, Hongbin Sun and Huaichang Ge  “Decoupled Piecewise Linear Power Flow and its Application to Under Voltage Load Shedding”  In _CSEE Journal of Power and Energy Systems_ , 2020  DOI: 10.17775/CSEEJPES.2019.02900
  * [18] Xin Chen, Wenchuan Wu and Boming Zhang  “Robust Restoration Method for Active Distribution Networks”  In _IEEE Transactions on Power Systems_ 31.5, 2016  DOI: 10.1109/TPWRS.2015.2503426
  * [19] Guanghui Sun et al.  “Decoupled Piecewise Linear Relaxation for the AC Optimal Power Flow Problem”  In _2022 5th International Conference on Renewable Energy and Power Engineering (REPE)_ IEEE, 2022  DOI: 10.1109/REPE55559.2022.9949533
  * [20] Masoud Farivar and Steven H. Low  “Branch Flow Model: Relaxations and Convexification—Part I”  In _IEEE Transactions on Power Systems_ 28.3, 2013  DOI: 10.1109/TPWRS.2013.2255317
  * [21] Leon Thurner et al.  “Pandapower—An Open-Source Python Tool for Convenient Modeling, Analysis, and Optimization of Electric Power Systems”  In _IEEE Transactions on Power Systems_ 33.6, 2018  DOI: 10.1109/TPWRS.2018.2829021
  * [22] Ray Daniel Zimmerman, Carlos Edmundo Murillo-Sanchez and Robert John Thomas  “MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education”  In _IEEE Transactions on Power Systems_ 26.1, 2011  DOI: 10.1109/TPWRS.2010.2051168
  * [23] Joel A.. Andersson et al.  “CasADi: a Software Framework for nonlinear Optimization and Optimal Control”  In _Mathematical Programming Computation_ 11.1, 2019  DOI: 10.1007/s12532-018-0139-4
  * [24] Gurobi Optimization, L. L. C.  “Gurobi Optimizer”  URL: https://www.gurobi.com