††thanks: These authors contributed equally.
$\dagger$<EMAIL_ADDRESS>
${\ddagger}$<EMAIL_ADDRESS>
$\sharp$<EMAIL_ADDRESS>
# Symmetry-protected higher-order exceptional points in staggered flatband
rhombic lattices

Yingying Zhang1∗, Shiqiang Xia1∗†, Xingdong Zhao1‡, Lu Qin1, Xuejing Feng1,
Wenrong Qi1, Yajing Jiang1,
Hai Lu1, Daohong Song2, Liqin Tang2♯, Zunlue Zhu1 and Yufang Liu1 1 School of
Physics, Henan Normal University, Xinxiang, Henan 453007, China
2 The MOE Key Laboratory of Weak-Light Nonlinear Photonics, TEDA Applied
Physics Institute and School of Physics, Nankai University, Tianjin 300457,
China

###### Abstract

Higher-order exceptional points (EPs), which appear as multifold degeneracies
in the spectra of non-Hermitian systems, are garnering extensive attention in
various multidisciplinary fields. However, constructing higher-order EPs still
remains as a challenge due to the strict requirement of the system symmetries.
Here we demonstrate that higher-order EPs can be judiciously fabricated in
$\mathcal{PT}$-symmetric staggered rhombic lattices by introducing not only
on-site gain/loss but also non-Hermitian couplings. Zero-energy flatbands
persist and symmetry-protected third-order EPs (EP3) arise in these systems
owing to the non-Hermitian chiral/sublattice symmetry, but distinct phase
transitions and propagation dynamics occur. Specifically, the EP3 arises at
the Brillouin zone (BZ) boundary in the presence of on-site gain/loss. The
single-site excitations display an exponential power increase in the
$\mathcal{PT}$-broken phase. Meanwhile, a nearly flatband sustains when a
small lattice perturbation is applied. For the lattices with non-Hermitian
couplings, however, the EP3 appears at the BZ center. Quite remarkably, our
analysis unveils a dynamical delocalization-localization transition for the
excitation of the dispersive bands and a quartic power increase beyond the
EP3. Our scheme provides a new platform towards the investigation of the
higher-order EPs, and can be further extended to the study of topological
phase transitions or nonlinear processes associated with higher-order EPs.

††preprint: APS/123-QED

##  I Introduction

Exceptional points (EPs), singularities at which eigenvalues and eigenvectors
simultaneously coalesce, have attracted considerable interest in the non-
Hermitian physics [1, 2, 3, 4, 5, 6, 7]. For instance, EPs in the parity–time
($\mathcal{PT}$)-symmetric optical arrangements [8, 9] have provided an
excellent platform for exploring intriguing phenomena such as the
$\mathcal{PT}$-symmetric lasers [10, 11] and topological phase transitions
[12, 13, 14]. In particular, recently, there has been a great deal of interest
in proposing various methods for constructing higher-order (greater than
second order) EPs. The eigenvalue shift near the higher-order EPs follows the
$\epsilon^{1/N}$ ($N$ is the order of the EPs) power law of the external
perturbation $\epsilon$, therefore leading to higher sensitivity of resonant
optical structures to small external disturbances [15, 16, 17, 18]. Another
interesting thing is that optical systems with $\mathcal{PT}$ symmetry can
show additional symmetries, such as time-reversal symmetry and chiral symmetry
(CS), making the higher-order EPs much more abundant [19, 20, 21]. In fact,
these symmetries ramify and show new features in non-Hermitian systems due to
the fact that the Hamiltonian $H\neq H^{\dagger}$. A crucial example is CS
which satisfies $CHC^{-1}=-H$. The CS coincides with sublattice symmetry (SLS)
($SH^{\dagger}S^{-1}=-H$) in the presence of Hermiticity. Obviously, they are
not equivalent for non-Hermitian Hamiltonians [22, 23, 24, 25].

The study of EPs has also been extended to flatband systems. A flatband is
dispersionless throughout the Brillouin zone (BZ) and it provides a unique
setting for exploring anomalous magnetic phases and strongly correlated states
of matter [26, 27, 28, 29]. In flatband systems, diffraction is suppressed due
to destructive interference of the Bloch wave functions, resulting in
typically compact localized states (CLSs) in real space [30, 31, 32, 33, 34,
35]. The most direct way to form $\mathcal{PT}$-symmetric flatband system is
to introduce on-site gain/loss to the lattices. Recently, the
$\mathcal{PT}$-symmetric flatband photonic waveguide arrays with tailored on-
site gain/loss distributions have been experimentally created using laser-
written technique [36, 37]. On the other hand, the non-Hermitian couplings
describe a situation where the mode amplitude undergoes gain/loss while
hopping between sites [38, 39, 40]. It has shown that the
$\mathcal{PT}$-symmetric frustrated lattices with fine-tuned non-Hermitian
couplings can also give rise to flatbands as well as EPs [41]. More
importantly, recent progress has demonstrated that even the higher-order EPs
can exist in the $\mathcal{PT}$-symmetric flatband lattices [20, 21, 42].
Nevertheless, due to the band touching (vanishing bandgap) between flat and
dispersive bands, so far, realizations of the non-zero valued higher-order EPs
have been limited to flatband systems with non-Hermitian couplings [43, 44,
45, 46]. Moreover, the distinct features of the two aspects of non-Hermiticity
(on-site gain/loss and non-Hermitian couplings) as well as the dynamical
behaviors near the phase transitions, also remain elusive.

In this work, we demonstrate that a higher-order EP (EP3) can occur in
staggered rhombic lattices by introducing the diagonal $\mathcal{PT}$ symmetry
(with on-site gain/loss) and the off-diagonal $\mathcal{PT}$ symmetry (with
non-Hermitian couplings) to the system. These $\mathcal{PT}$-symmetric
lattices possess both the non-Hermitian CS and SLS, thus leading to the
presence of zero-energy flatbands and non-zero EP3. However, the EP3 and the
propagation dynamics of the lattices, unexpectedly, show strikingly different
properties. For the diagonal $\mathcal{PT}$-symmetric lattices, a EP3 arises
at the BZ boundary. We find that the single-site excitations of different
sublattices display identical exponential power increases in the
$\mathcal{PT}$-broken phase. By contrast, a EP3 is generated in the BZ center
for the off-diagonal $\mathcal{PT}$-symmetric lattices. Intriguing propagation
dynamics, such as a dynamical delocalization-localization transition for the
excitation of dispersive bands and a quartic power increase beyond the EP3,
are observed. For both diagonal and off-diagonal $\mathcal{PT}$-symmetric
systems, the eigenvalue separations exhibit an $\sim\epsilon^{1/3}$ dependence
in the vicinity of the EP3. In particular, we show that a nearly flatband
sustains for the diagonal $\mathcal{PT}$-symmetric lattices in a small range
of perturbations.

##  II The diagonal $\mathcal{PT}$-symmetric flatband rhombic lattices

Figure 1: (a) Schematic of the diagonal $\mathcal{PT}$ -symmetric photonic
rhombic lattices formed by waveguide arrays without non-Hermitian coupling.
The unit cell consists of three sites: two edge sites with gain $\gamma$
($\bf{A}$, red) or loss -$\gamma$ ($\bf{C}$, blue) and one central site with a
neutral imaginary part ($\bf{B}$, gray). The staggered coupling coefficients
are $t_{L}=t(1+g)$ and $t_{R}=t(1-g)$, where $0\leq g<1$ and $t$ is the
coupling coefficient of uniform lattices ($g=0$). For simplicity, we set the
lattice period $l=1$. (b, c) Calculated real (b) and imaginary (c) parts of
the spectrum as a function of $\gamma/t$ for $g=0.2$. The red circles
represent the position of the EP3 and the corresponding critical value of
phase transition is $\gamma_{c}/t$=0.56. (d) Field distribution of a flatband
CLS, which occupies two unit cells. The characters show the corresponding
amplitudes.  Figure 2: Intensity profile $\left|\psi\right|^{2}$ of a single-
site $\bf{A}$ (a1-a3) or $\bf{B}$ (b1-b3) excitation. (a1) The $\bf{A}$
sublattice excitation for Hermitian case ($\gamma/t=0$) exhibits a suppression
of diffraction due to the excitation of flatband CLSs. (a2) The input remains
conserved below the EP3 threshold ($\gamma/t=0.5$). (a3) The input leads to an
exponentially growing total intensity in the $\mathcal{PT}$-broken phase
($\gamma/t=0.6$). (b1) The $\bf{B}$ site excitation for Hermitian case
($\gamma/t=0$), showing a typically discrete diffraction pattern. (b2) The
input excites both the flatband and dispersive bands, producing discrete
diffraction as well as a residual localized component ($\gamma/t=0.5$). (b3)
The input shows same exponential growth of total intensity as in (a3)
($\gamma/t=0.6$). The propagation length is $z$=20 $t^{-1}$. Note the
different scales of the vertical axis in these panels, though the initial
amplitudes of the excitation all equal 1.

Symmetry-protected flatband and EP3–We start from introducing the diagonal
$\mathcal{PT}$-symmetric staggered photonic rhombic lattices with on-site
gain/loss to elucidate the lattice symmetry, the formation of flatband and the
EP3. Figure 1(a) shows the lattice structure with each unit cell consisting of
an edge site $\bf{A}$ with gain ${\gamma}$, an edge site $\bf{C}$ with loss
${-\gamma}$, and a central neutral site $\bf{B}$. The sublattice $\bf{B}$ are
fourfold connected with the nearest neighbors, while the sublattices $\bf{A}$
and $\bf{C}$ have twofold connections with surrounding sites. Experimentally,
such non-Hermitian lattices can be obtained using the optical induced method
or the femtosecond laser-writing technique [36, 37, 14]. The couplings depend
on the spacing of the neighbor waveguides, thus a staggered arrangement of the
distance between the waveguides generates staggered couplings. Meanwhile, the
on-site gain/loss can be generated by inscribing the waveguides with tailored
loss distributions. The staggered couplings are $t_{L}=t(1+g)$ and
$t_{R}=t(1-g)$, where $t_{L}(t_{R})$ is the coupling on the left (right) side
of the site $\bf{B}$. Note that $t_{L}(t_{R})$ is real and represents
conserved (Hermitian) coupling. In the tight binding model, the Hamiltonian
can be written

$\displaystyle{H}=$
$\displaystyle{\sum_{n}}\big{(}t_{L}b_{n}^{\dagger}a_{n}+t_{R}b_{n}^{\dagger}a_{n+1}+t_{L}b_{n}^{\dagger}c_{n}+t_{R}b_{n}^{\dagger}c_{n+1}$
(1) $\displaystyle+{\rm
H.c.}\big{)}+i\gamma\left(a_{n}^{\dagger}a_{n}-c_{n}^{\dagger}c_{n}\right),$

where $a_{n}^{\dagger}$ ($a_{n}$), $b_{n}^{\dagger}$ ($b_{n}$) and
$c_{n}^{\dagger}$ ($c_{n}$) are the creation (annihilation) operators in the
$n$-th unit cell on the $\bf{A}$, $\bf{B}$ and $\bf{C}$ sites, respectively.
Transforming $H$ into momentum space as
$H=\sum_{k}\psi_{k}^{\dagger}H_{k}\psi_{k}$ with
$\psi_{k}=\left(a_{k},b_{k},c_{k}\right)^{T}$ , one can get the momentum space
Hamiltonian

$\displaystyle{H_{k}}=\begin{pmatrix}i\gamma&t_{L}+t_{R}{e^{ik}}&0\\\
t_{L}+t_{R}{e^{-ik}}&0&t_{L}+t_{R}{e^{-ik}}\\\
0&t_{L}+t_{R}{e^{ik}}&-i\gamma\end{pmatrix}.$ (2)

The diagonal terms of $H_{k}$ describe the propagation constant and on-site
gain/loss, whereas the the off-diagonal terms describe the lattice couplings.
Such a Hamiltonian possesses $\mathcal{PT}$ symmetry
$(PT){H_{k}}(PT)^{-1}={H_{-k}}$, and more importantly, a CS as well as a SLS
satisfying $CH_{k}C^{-1}=-H_{k}$ and $SH_{k}^{\dagger}S^{-1}=-H_{k}$, where
the unitary operator $C$ and $S$ are defined as

$\displaystyle{C}=\begin{pmatrix}0&0&1\\\ 0&-1&0\\\ 1&0&0\end{pmatrix},\
{S}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\ 0&0&1\end{pmatrix}.$ (3)

By calculating the eigenvalues of $H_{k}$, one can easily get the three
corresponding eigenvalues:

$\displaystyle{{E_{0}(k)=0}},$ $\displaystyle{E_{\pm
1}\left(k\right)=\pm\sqrt{-\gamma^{2}+4t^{2}\left(1+{\rm
cos}k\right)+4t^{2}g^{2}\left(1-{\rm cos}k\right)}}.$ (4)

In the Hermitian case ($\gamma/t=0$), the CS/SLS symmetry requires all the
energy eigenvalues appear in pairs, namely, the three-band system must have a
zero-energy ﬂatband $E_{0}=0$ and two dispersive bands with opposite energies
$E_{1}=-E_{-1}$ [47, 48]. Figures 1(b) and 1(c) show the spectrum as a
function of $\gamma/t$ for $g$=0.2. It can be clearly seen that the Hermitian
spectrum in the limit $\gamma/t=0$ has three bands: a zero-energy flatband
located between two symmetric dispersive bands (dotted lines in Fig. 1(b)).
Unlike in uniform lattices where all the bands touch each other at $k=\pi$
[49, 50, 51], here the band degeneracy is lifted and two symmetric gaps with
gap width $E_{gap}=2\sqrt{-\gamma^{2}+8t^{2}g^{2}}$ emerge as a result of the
staggered couplings. When the on-site gain/loss is introduced ($\gamma/t\neq
0$), the energy eigenvalues become complex. As mentioned above, the CS
bifurcates into CS and SLS in the non-Hermitian case. However, we find both
the CS and SLS preserve. Meanwhile, they retain $E_{1}=-E_{-1}$ and non-
Hermitian zero mode occurs for the central band. As a consquence, the zero-
energy flatband survives and the spectrum stays symmetric about the imaginary
and real energy axis of the complex energy plane. The band gaps decrease by
increasing $\gamma/t$. More importantly, all the three eigenvalues coalesce at
$E_{0}=0$ for the critical value $\gamma_{c}/t=2\sqrt{2}g=0.56$ (Figs. 1(b)
and 1(c)). Note that $\gamma_{c}$ depends on the value of $g$, i.e., the
relative strength of the staggered couplings. The three eigenvectors also
coalesce at this specific degenerate point, and have the form
$\left[\psi_{n}^{A},\psi_{n}^{B},\psi_{n}^{C},\psi_{n+1}^{A},\psi_{n+1}^{C}\right]=\left[0.8,-0.2i,-0.8,1.2,-1.2\right]$.
Accordingly, the spectrum is completely real for $\gamma<\gamma_{c}$, and a
CS/SLS-protected EP3 emerges at $\gamma=\gamma_{c}$. For $\gamma>\gamma_{c}$,
the system turns into the $\mathcal{PT}$-broken phase. Figure 1(d) illustrates
the irreducible CLS of the flatband. As demonstrated in previous studies, CLSs
are localized in a single unit cell in Hermitian rhombic lattices [52, 53,
54]. Here, in contrast, the CLS is confined in the five-site-cross
configurations in two adjacent unit cells. Non-zero amplitudes appear in
$\bf{A}$ and $\bf{C}$ sites and the on-site gain/loss remodulates the
amplitude of the $\bf{B}$ site. The equal amplitude and opposite phase ensure
the destructive interference between the two legs formed by the $\bf{A}$ and
$\bf{C}$ sites.

Figure 3: (a, b) Analytical (dashed lines) and numerical (solid lines) results
for the real (a) and imaginary (b) parts of the spectrum as a function of the
detuning $\epsilon/t$ in the vicinity of EP3. Red lines, $E_{-1}$; green
lines, $E_{0}$; blue lines, $E_{1}$. (c) Corresponding results for
$E_{0}-E_{1}=\bigtriangleup E$, exhibiting a cube-root behavior. (d) The
results from (c) on a logarithmic scale. The slope of about 1/3 confirms the
cube-root response.

Propagation dynamics for single-site excitation–To verify our analysis, we
simulate the light propagation dynamics for initial excitations localized to
either an $\bf{A}$ or $\bf{B}$ site. A sufficiently large number of lattice
sites (101 unit cells) is considered. The input is $a_{n}(0)$ or
$b_{n}\left(0\right)=\left(\sqrt{\pi}/\sigma\right)e^{-\left(n^{2}/\sigma^{2}\right)}\left(\sigma=0.5\right)$,
i.e., only the central single site is excited. Figure 2 shows the evolution of
the beam intensity
$P\left(z\right)=\sum_{n}\left|\psi_{n}\left(z\right)\right|^{2}=\sum_{n}\left(\left|a_{n}\left(z\right)\right|^{2}+\left|b_{n}\left(z\right)\right|^{2}+\left|c_{n}\left(z\right)\right|^{2}\right)$,
where $z$ is an effective propagation distance. For comparison, we first
consider the result for the Hermitian case at $\gamma/t=0$. A suppression of
diffraction is observed for the excitation of $\bf{A}$ site (Fig. 2(a)). The
input excites a superposition of dispersive bands and flatband states. The
energy is mainly distributed in the flatband, therefore most of the power
remains confined to the initially excited site. When the system is operated
below the $\mathcal{PT}$ threshold value ($\gamma/t=0.5$), part of the energy
evolves into the dispersive bands though most of the energy still stays
localized (Fig. 2(a2)). It should be noted that, what is conserved here is the
quasi-energy
$Q\left(z\right)=\sum_{n}\psi\left(-n,z\right)\psi^{\ast}\left(n,z\right)$ as
opposed to the actual power $P(z)$. The presence of the gain/loss results in
exponential amplification $\left[P\left(z\right)\sim e^{z}\right]$ of the
total power when the non-Hermitian parameter is tuned to $\gamma/t$=0.6 (Fig.
2(a3)) (see Appendix A). A different scenario is found when a localized
$\bf{B}$ site is excited. In Fig. 2(b1), the input generates a conventional
discrete diffraction pattern with two ballistically expanding lobes for the
Hermitian lattices. In this case, only the dispersive bands are excited since
no flatband energy is distributed in the $\bf{B}$ sublattice (cf. Fig. 1(d)).
Nevertheless, the dispersive bands and flatband states, including the EP3
states are excited in the $\mathcal{PT}$ symmetry phase $\gamma/t$=0.5 (Fig.
2(b2)). The resulting diffraction consists of a discrete pattern and conserved
total power as required for a $\mathcal{PT}$-symmetric lattice in the unbroken
phase. Furthermore, there is a residual localized component due to the
excitation of the flatband. With the further increase of $\gamma/t$, the
energy is also obviously amplified when the system is operated in the
$\mathcal{PT}$ symmetry broken phase for $\gamma/t$=0.6 (Fig. 2(b3)).
Interestingly, we find that the power growth is the same as that in Fig.
2(a3).

Bifurcations of eigenvalues around the EP3–In practical implementations, the
propagation of photons is subject to inevitable perturbations. As mentioned
above, the eigenvalues bifurcation characteristics of the higher-order EPs can
enhance the response of the optical structure to small external disturbances.
In the following, we demonstrate the eigenvalues evolution of our ternary
$\mathcal{PT}$-symmetric photonic structure around the EP3 and show how a
small lattice perturbation affects the wave propagation dynamics.

Note that the perturbation could be introduced anywhere along the diagonal of
the matrix $H_{k}$. For simplicity and without loss of generality, we suppose
that the external perturbation $\epsilon$ is imposed on the sublattice
$\bf{A}$.

Figure 4: (a, b) Intensity profile for the $\bf{B}$ site excitation in the
presence of small perturbations $\epsilon/t=0.03$ and $\epsilon/t=0.06$,
respectively. For $\epsilon/t=0.03$, a similar pattern compared with Fig.
2(b2) is obtained. (c, d) Same as (a, b) but for the excitation of a flatband
CLS. The CLS can stay well localized despite the presence of the lattice
perturbations.

The energy bands can be analytically obtained through solving the cubic
characteristic equation (see Appendix B). Figures 3(a) and 3(b) illustrate
analytical and numerical solutions for the real and imaginary parts of the
eigenvalues as a function of $\epsilon/t$. The small perturbation ranges from
$10^{-10}$ to $10^{-1}$. The analytical solutions show a good agreement with
numerical results. The perturbation drives the system out of the degenerate
state and three complex eigenvalues are obtained. To have a quantitative
comparison, we calculate the numerical eigenvalue difference of the upper band
in Fig. 3(a) for $\epsilon/t=0.03$ and $\epsilon/t=0.06$, which are 0.1 and
0.15, respectively. The difference between the lower two eigenvalues ($E_{0}$
and $E_{1}$) is plotted in Fig. 3(c). Obviously, the eigenvalues splitting
increases along with the strength of the perturbation. By considering the
logarithmic behavior of this curve, we find that the slope is about 1/3, which
is exactly the characteristic of an EP3 (Fig. 3(d)). Another intriguing
feature is that both the real and imaginary parts of the central band change
slightly, especially for $\epsilon/t<0.03$. This indicates that a nearly
flatband can be sustained in such a range of perturbations.

The perturbation on the spectrum can also be reflected by light propagation
dynamics. For $\epsilon/t=0.03$, one can find a similar pattern for the
excitation of $\bf{B}$ sublattice in Fig. 4(a), which is compared with that in
Fig. 2(a2). An increase of the lattice perturbation ($\epsilon/t=0.06$)
results in an enhancement of wave intensity in Fig. 4(b). The input behaves
differently for the excitation of the flatband CLSs. An interesting property
of flatband is that, despite the non-zero couplings between lattice sites,
propagation of flatband eigenstates is completely suppressed. Therefore, if we
excite five sites according to Fig. 1(d), we can observe a completely
localized state with $z$-invariant intensity. In Fig. 4(c), a CLS excitation
can stay well localized in the presence of $\epsilon/t=0.03$, confirming that
the flatband can preserve in this case. Moreover, the input keeps nearly
localized even when the lattice detuning is increased up to $\epsilon/t=0.06$
as shown in Fig. 4(d).

##  III The off-diagonal $\mathcal{PT}$-symmetric flatband rhombic lattices

Figure 5:  (a) Schematic of the off-diagonal $\mathcal{PT}$-symmetric photonic
rhombic lattices formed by waveguide arrays with non-Hermitian coupling. The
complex-coupling coefficients are $T_{L}=t_{L}+i\gamma$,
$T_{L}^{*}=t_{L}-i\gamma$, $T_{R}=t_{R}+i\gamma$, and
$T_{R}^{*}=t_{R}-i\gamma$, respectively. Note that the real parts ($t_{L}$ and
$t_{R}$) are the same as in Fig. 1(a), which represent conserved staggered
couplings. (b, c) Calculated real (b) and imaginary (c) parts of the spectrum
as a function of $\gamma/t$ for $g=0.2$. The EP3 forms at $k=0$ and the
critical value of the phase transition is $\gamma_{c}=t$. (d) Field
distribution of the flatband CLS, which also occupies two unit cells but only
distributes in $\bf{A}$ and $\bf{C}$ sites.

Symmetry-protected flatband and EP3–The most intriguing peculiarity of our
scheme is that the non-Hermiticity can also be introduced by employing the
non-Hermitian couplings. Such a system also possesses an EP3 and a flatband
but exhibits different features. As illustrated in Fig. 5(a), we maintain the
geometric symmetry of the lattices in Fig. 1(a) and perform imaginary
processing on the coupling parameters. $T_{L}=t_{L}+i\gamma$,
$T_{L}^{*}=t_{L}-i\gamma$, $T_{R}=t_{R}+i\gamma$, and
$T_{R}^{*}=t_{R}-i\gamma$ are the complex-valued coupling coefficients on the
left (right) side of the site $\bf{B}$, respectively. These complex couplings
can be regarded as the phase of a hopping amplitude, which can be
experimentally realized by embedding amplifying/lossy media between adjacent
waveguides or periodically modulate the on-site gain/loss amplitude ratio [41,
55, 56, 57, 58]. Light propagates in engineered lattices with segmented
regions of alternating gain/loss and index contrast regions. The effective
gain/loss is sensitive to the relative phase between neighboring waveguides,
leading to stronger or weaker amplifications. In this case, the Hamiltonian
can be written as

$\displaystyle{H}=$
$\displaystyle{\textstyle\sum_{n}}\big{(}T_{L}b_{n}^{\dagger}a_{n}+T_{R}b_{n}^{\dagger}a_{n+1}+T^{*}_{L}b_{n}^{\dagger}c_{n}+T^{*}_{R}b_{n}^{\dagger}c_{n+1}$
(5) $\displaystyle+{\rm H.c.}\big{)}.$

The corresponding Hamiltonian in momentum space is

$\displaystyle{H_{k}}=\begin{pmatrix}0&T_{L}+T_{R}{e^{ik}}&0\\\
T_{L}+T_{R}{e^{-ik}}&0&T^{*}_{L}+T^{*}_{R}{e^{-ik}}\\\
0&T^{*}_{L}+T^{*}_{R}{e^{ik}}&0\end{pmatrix},$ (6)

which satisfies the $\mathcal{PT}$ symmetry
$\left(PT\right)H_{k}\left(PT\right)^{-1}=H_{-k}$. Interestingly, one can find
that such a arrangement also possesses CS and SLS. Here, the unitary operator
$C$ and $S$ can be set as

$\displaystyle{C}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\ 0&0&1\end{pmatrix},\
{S}=\begin{pmatrix}0&0&1\\\ 0&-1&0\\\ 1&0&0\end{pmatrix}.$ (7)

The energy eigenvalues of Hamiltonian (6) are

$\displaystyle{E_{0}(k)=0},$ $\displaystyle{E_{\pm 1}\left(k\right)=\pm
2\sqrt{\left(t^{2}-\gamma^{2}\right)\left(1+\cos
k\right)+g^{2}t^{2}\left(1-\cos k\right)}}.$ (8)

Figure 6: Same as Fig. 2, but for the off-diagonal $\mathcal{PT}$-symmetric
lattices. (a1) The $\bf{A}$ sublattice excitation oscillates periodically
below the EP3 threshold ($\gamma/t=0.45$). Most of the energy stay localized
due to the excitation of the flatband. (a2) The quasi-energy is conserved and
the oscillation period is doubled for $\gamma/t=0.9$. (a3) The input leads to
quartic total intensity grow once the system is operated in the $\mathcal{PT}$
symmetry broken phase ($\gamma/t=1.001$). (b1) The $\bf{B}$ input only excites
eigenstates of the dispersive bands and therefore experiences conventional
diffraction. (b2) The spreading velocity diminishes with the increase of
$\gamma/t$. (b3) The input becomes localized and the total intensity also show
a quartic increase in the $\mathcal{PT}$ symmetry broken phase.

For direct comparison of the two non-Hermitian arrangements, we also set
$g$=0.2 and plot the real and imaginary parts of the spectrum (as shown in
Figs. 5(b) and 5(c)). Due to the presence of CS and SLS, the central flatband
still survives and the three bands distribute symmetrically. With the increase
in $\gamma/t$, the bandgap tends to narrow and merge. In contrast, here, the
$\mathcal{PT}$ symmetry broken transition occurs at the center of the BZ
($k=0$) and the critical value is $\gamma_{c}=t$, i.e., is invariable and
independent of the value of $g$. The eigenmodes of the bands also coalesce
into a degenerated one, namely, a CS/SLS-protected EP3 is generated. Above the
EP3, the spectrum becomes complex. In addition, the irreducible flatband CLS
also occupies two unit cells but the energy only distribute in $\bf{A}$ and
$\bf{C}$ sites as shown in Fig. 5(d).

Figure 7: Same as Fig. 3, but for the off-diagonal $\mathcal{PT}$-symmetric
lattices.

Propagation dynamics for single-site excitation–In Fig. 6, we illustrate the
propagation dynamics of the system for the same initial excitations as in Fig.
2. Similar patterns as depicted in Figs. 2(a1) and 2(b1) can be obtained in
the Hermitian limit $\gamma/t=0$. However, diverse propagation dynamics can be
observed for $\gamma/t\neq 0$. The input experiences periodically oscillations
for the $\bf{A}$ site excitation due to the unfolding of the nonorthogonal
Floquet-Bloch modes in $\mathcal{PT}$ symmetry phase $\gamma/t<1$ (Figs. 6(a1)
and 6(a2)). The beat length, which is equal to ${L=2\pi/(\delta E)}$ (${\delta
E}$ is the difference between the propagation constants of two modes),
increases as $\gamma/t$ approaches the EP3 [9, 59, 60]. As can be found in
Fig. 6(a2), the oscillation is enhanced and has a double period for
$\gamma/t=0.9$. The beat length goes to infinity as the EP3 is approached.
Furthermore, we observe a sharp transition in Fig. 6(a3). The input becomes
localized and leads to a quartic power increase $\left[P(z)\sim z^{4}\right]$
in the $\mathcal{PT}$ symmetry broken phase (see Appendix A). The dynamics
exhibited in Figs. 6(b1-b3) is distinct from that in Figs. 6(a1-a3). Transport
is always ballistic in the $\mathcal{PT}$ symmetry phase ($\gamma/t<1$) for
the $\bf{B}$ site excitation. In this case, the input will never excite the
flatband CLSs (cf. Fig. 5(d)) and thus no localized component appears. More
saliently, the excitation displays a transition from ballistic wave-packet
spreading to dynamical localization, as the non-Hermitian parameter $\gamma/t$
is increased to the value at the EP3. In Figs. 6(b1) and 6(b2), the excitation
propagates bidirectionally along the lattice with the spreading speed
$\nu\sim\sigma\left(z\right)/z$, where the position operator
$\sigma^{2}={\textstyle\sum_{n}}n^{2}\psi_{n}^{2}/{\textstyle\sum_{n}}\psi_{n}^{2}$.
$\nu$ decreases as $\gamma/t$ is increased and the wave-packet spreading
decreases faster for large $\gamma/t$ (see Appendix C). When $\gamma/t$
approaches the EP3, the diffraction is prevented. In fact, the largest
velocity at which an excitation propagates along the lattice can also be
described by the group velocity
$\nu_{g}=Re\left\\{\left(dE/dk\right)_{k=\pi/2}\right\\}$. As the $\gamma/t$
is increased, the energy spectrum in Fig. 5(b) undergoes a deformation, which
changes the dispersion relation and results in a decrease of $\nu_{g}$ [61,
62]. This fact can be used for a direct observation of the EP: the propagation
constants approach when the non-Hermitian parameter is increased to the value
at the EP3. Though the total intensity in Fig. 6(b3) is much smaller than that
in Fig. 6(a3), we find that it also follows a quartic power increase in the
$\mathcal{PT}$ symmetry broken phase.

Bifurcations of eigenvalues around the EP3–Similarly, we calculate the
spectrum separation of the system in the presence of the external
perturbations. The eigenvalues difference depicted in Fig. 7 also has a 1/3
power law of the external perturbation (Figs. 7(c) and 7(d)). Note that the
changes of the actual value of both the dispersive bands and flatband are
greater compared with that in Fig. 3. For instance, the numerical eigenvalue
difference in Fig. 7(a) for $\epsilon/t=0.03$ and $\epsilon/t=0.06$ are 0.54
and 0.69, respectively. Besides, we find that the imaginary part of the
central band is identical to the upper band, indicating that the flatband
transforms into a dispersive band (Fig. 7(b)).

Figure 8 displays the corresponding propagation dynamics in the presence of
the lattice perturbations. For $\bf{B}$ site excitation, the energy becomes
obviously amplified for $\epsilon/t=0.06$ due to the complex spectrum (Fig.
8(b)). More importantly, one can clearly find that the CLSs no longer exist
even only a small perturbation is applied (Fig. 8(c)) and the pattern in Fig.
8(d) is similar to Fig. 8 (b). These results reveal that the flatband is
readily destroyed and evolves into a complex dispersive band as a result of
the lattice perturbations.

Figure 8: Same as Fig. 4, but for the off-diagonal $\mathcal{PT}$-symmetric
lattices.

##  IV Conclusion

In summary, we have demonstrated a feasible scheme to realize higher-order EPs
(EP3) in flatband rhombic lattices. By introducing on-site gain/loss and non-
Hermitian couplings, we have established a diagonal $\mathcal{PT}$-symmetric
and an oﬀ-diagonal $\mathcal{PT}$-symmetric rhombic lattices, respectively.
Furthermore, we have made careful consideration of the distributions of the
higher-order EPs and the dynamical behaviors around them. It has shown that
the systems possess both the non-Hermitian CS and SLS, thus leading to the
presence of zero-energy flatbands and non-zero EP3. We have also observed
varied propagation dynamics in these systems, such as sustained flatband CLSs,
exponential and quartic power increase as well as dynamical localization.
These results may be useful for realizing the high-sensitivity in a passive
wireless sensing system [63, 64, 65] and provide insightful information about
the underlying properties of the higher-order EPs in flatband systems. In
future works, it will also be interesting to investigate a wide range of
phenomena such as the topological transitions and nonlinear optical processes
in these systems [66, 67, 68, 69, 70].

Acknowledgments: This work was supported by the National Key R&D Program of
China (Grant No. 2017YFA0303800), the NSFC (Grant Nos. 11274096, 12134006,
11922408, and 12074105), the 111 Project (Grant No. B07013) in China.

##  Appendix A The energy amplification in the $\mathcal{PT}$ symmetry broken
phase

A localized single-site excitation in these $\mathcal{PT}$-symmetric lattices
can display exponential and quartic power increases in the $\mathcal{PT}$
symmetry broken phase. Specifically speaking, for the diagonal
$\mathcal{PT}$-symmetric lattices, the power increases for the excitation of
$\bf{A}$ and $\bf{B}$ sites are identical and follow an exponential
amplification $\left[P(z)=e^{0.39z}\right]$ (green line in Fig. 9(a)). For the
off-diagonal $\mathcal{PT}$-symmetric lattices, the power increase follows
quartic amplifications $\left[P(z)\sim z^{4}\right]$. However, the energy
amplification is greater for the excitation of $\bf{A}$ site (Fig. 9(b)). As
can be seen in Fig. 9(d), the corresponding intensities $P^{1/4}$ increase
linearly but the two lines have different slopes.

Figure 9: Energy evolution in the symmetry broken phase of the (a) diagonal
and (b-d) off-diagonal $\mathcal{PT}$-symmetric lattices. (a) The power
increase for the excitation of $\bf{A}$ and $\bf{B}$ sublattices are identical
and follow an exponential amplification $\left[P(z)=e^{0.39z}\right]$. (b-d)
The power increases for the excitation of (b) $\bf{A}$ and (c) $\bf{B}$
sublattices follow quartic amplifications $\left[P(z)\sim z^{4}\right]$, and
the corresponding intensities $P^{1/4}$ increase linearly. Note that the two
lines in (d) have different slops.

##  Appendix B  Eigenvalue evolutions of the $\mathcal{PT}$-symmetric lattices

Here, we provide an analytical study of the eigenvalue evolutions of our
$\mathcal{PT}$-symmetric structure in the vicinity of the EP3. For simplicity,
we only discuss the case of perturbing the gain sublattice $\bf{A}$. As a
result, the perturbation terms only appear along the diagonal element of
$H_{k}$. For the diagonal $\mathcal{PT}$-symmetric lattices, the Hamiltonian
$H_{k}$ can be written as

$\displaystyle{H_{k}}=\begin{pmatrix}i\gamma+\epsilon&t_{L}+t_{R}{e^{ik}}&0\\\
t_{L}+t_{R}{e^{-ik}}&0&t_{L}+t_{R}{e^{-ik}}\\\
0&t_{L}+t_{R}{e^{ik}}&-i\gamma\end{pmatrix}.$ (B.1)

The corresponding determinant of $H_{k}$ equated to zero is

$\displaystyle\begin{vmatrix}-E_{n}+i\gamma+\epsilon&t_{L}+t_{R}e^{ik}&0\\\
t_{L}+t_{R}e^{-ik}&-E_{n}&t_{L}+t_{R}e^{-ik}\\\
0&t_{L}+t_{R}e^{ik}&-E_{n}-i\gamma\end{vmatrix}=0.$ (B.2)

The energy bands can be analytically obtained through solving the cubic
characteristic equation, which is simplified to

$\displaystyle 0=E_{n}^{3}-E_{n}^{2}\epsilon-
E_{n}\left(i\gamma\epsilon-\gamma^{2}+2A\right)+\epsilon A,$ (B.3)

where $\gamma_{c}-2A=0$, because $A=t_{L}^{2}+t_{R}^{2}+2t_{L}t_{R}\cos
k=0.16t^{2}$ at the EP3 ($k=\pi,\gamma_{c}/t=2\sqrt{2}g$). By means of
Newton–Puiseux series $E_{n}\sim c_{1}\epsilon^{1/3}+c_{2}\epsilon^{2/3}$
($c_{1}$ and $c_{2}$ are complex constants), we can expand the above equation
by perturbation

$\displaystyle 0=(c_{1}^{3}+A)\epsilon^{3/3}+(3c_{1}^{2}c_{2}-i\gamma
c_{1})\epsilon^{4/3}+(3c_{1}c_{2}^{2}-c_{1}^{2}-c_{2}i\gamma)$
$\displaystyle\epsilon^{5/3}+(c_{2}^{3}-2c_{1}c_{2})\epsilon^{6/3}-c_{2}^{2}\epsilon^{7/3}.$
(B.4)

Forcing the coefficients of the first two terms to be zero, we obtain three
sets of values for the coefficients $c_{1}$ and $c_{2}$, corresponding to the
three eigenvalues

$\displaystyle
E_{-1}\sim\sqrt[3]{\frac{4}{25}}e^{i\pi/3}\epsilon^{1/3}+\sqrt[3]{\frac{4}{25}}\frac{i\sqrt{2}}{3}e^{-i\pi/3}\epsilon^{2/3},$
$\displaystyle
E_{0}\sim\sqrt[3]{\frac{4}{25}}e^{-i\pi/3}\epsilon^{1/3}+\sqrt[3]{\frac{4}{25}}\frac{i\sqrt{2}}{3}e^{i\pi/3}\epsilon^{2/3},$
$\displaystyle
E_{1}\sim-\sqrt[3]{\frac{4}{25}}e^{i\pi/3}\epsilon^{1/3}-\sqrt[3]{\frac{4}{25}}\frac{i\sqrt{2}}{3}e^{-i\pi/3}\epsilon^{2/3}.$
(B.5)

Similarly, we can get the results of the off-diagonal $\mathcal{PT}$-symmetric
system. The Hamiltonian $H_{k}$ can be written as

$\displaystyle{H_{k}}=\begin{pmatrix}\epsilon&T_{L}+T_{R}{e^{ik}}&0\\\
T_{L}+T_{R}{e^{-ik}}&0&T^{*}_{L}+T^{*}_{R}{e^{-ik}}\\\
0&T^{*}_{L}+T^{*}_{R}{e^{ik}}&0\end{pmatrix}.$ (B.6)

The determinant of $H_{k}$ equated to zero is

$\displaystyle\begin{vmatrix}-E_{n}+\epsilon&T_{L}+T_{R}e^{ik}&0\\\
T_{L}+T_{R}e^{-ik}&-E_{n}&T_{L}+T_{R}e^{-ik}\\\
0&T_{L}+T_{R}e^{ik}&-E_{n}\end{vmatrix}=0.$ (B.7)

Here, $A=T_{L}^{2}+T_{R}^{2}+2T_{L}T_{R}\cos k=8i$ and
$A^{*}=T_{L}^{*2}+T_{R}^{*2}+2T_{L}^{*}T_{R}^{*}\cos k=-8i$ at the EP3
($k=0,\gamma_{c}/t=1$). So, the determinant can be simplified to

$\displaystyle 0=E_{n}^{3}-E_{n}^{2}\epsilon+\epsilon A^{*}.$ (B.8)

Using the Newtonian series expansion for $E_{n}\sim
c_{1}\epsilon^{1/3}+c_{2}\epsilon^{2/3}$, one can get

$\displaystyle 0=$
$\displaystyle(c_{1}^{3}-8i)\epsilon^{3/3}+3c_{1}^{2}c_{2}\epsilon^{4/3}+(3c_{1}c_{2}^{2}-c_{1}^{2}-c_{2}i\gamma)\epsilon^{5/3}$
(B.9)
$\displaystyle+(c_{2}^{3}-2c_{1}c_{2})\epsilon^{6/3}-c_{2}^{2}\epsilon^{7/3}.$

The bifurcations in the eigenvalues have the form

$\displaystyle E_{-1}\sim 2e^{i5\pi/6}\epsilon^{1/3},$ $\displaystyle
E_{0}\sim 2e^{i\pi/6}\epsilon^{1/3},$ $\displaystyle
E_{1}\sim-2i\epsilon^{1/3}.$ (B.10)

##  Appendix C The spreading velocity $\nu$ for the $\bf{B}$ site excitation
of the off-diagonal $\mathcal{PT}$-symmetric lattices

For the off-diagonal $\mathcal{PT}$-symmetric lattices, the $\bf{B}$ site
excitation displays a transition from ballistic wave-packet spreading to
dynamical localization. As can be seen in Fig. 10(a)), the diffraction is
prevented when the non-Hermitian parameter $\gamma/t$ approaches the value at
the EP3. Meanwhile, the spreading speed $\nu$ decreases faster for large
$\gamma/t$ as shown in Fig. 10(b))

Figure 10: (a) Intensity profile $\left|\psi\right|^{2}$ of a localized
$\bf{B}$ sublattice excitation for $\gamma/t=0.98$. When $\gamma/t$ approaches
the EP3, the diffraction is prevented. (b) Behavior of the spreading velocity
$\nu$ as a function of $\gamma/t$, computed for the largest propagation length
$z$ = 20 $t^{-1}$. $\nu$ is expressed in units of $t$.

Disclosures: The authors declare no conflicts of interest.

Data Availability Statement: Data underlying the results presented in this
paper are not publicly available at this time but may be obtained from the
authors upon reasonable request.

## References

  * Feng  _et al._ [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian photonics based on parity–time symmetry, Nature Photonics 11, 752 (2017).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and $\mathcal{P}\mathcal{T}$ symmetry, Nature Physics 14, 11 (2018).
  * Özdemir  _et al._ [2019] Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Parity–time symmetry and exceptional points in photonics, Nature Materials 18, 783 (2019).
  * Miri and Alu [2019] M.-A. Miri and A. Alu, Exceptional points in optics and photonics, Science 363, eaar7709 (2019).
  * Parto  _et al._ [2021] M. Parto, Y. G. Liu, B. Bahari, M. Khajavikhan, and D. N. Christodoulides, Non-Hermitian and topological photonics: optics at an exceptional point, Nanophotonics 10, 403 (2021).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Reviews of Modern Physics 93, 015005 (2021).
  * Xiao  _et al._ [2021] L. Xiao, T. Deng, K. Wang, Z. Wang, W. Yi, and P. Xue, Observation of Non-Bloch parity-time symmetry and exceptional points, Phys. Rev. Lett. 126, 230402 (2021).
  * Makris  _et al._ [2008] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Beam dynamics in $\mathcal{P}\mathcal{T}$ symmetric optical lattices, Phys. Rev. Lett. 100, 103904 (2008).
  * Klaiman  _et al._ [2008] S. Klaiman, U. Günther, and N. Moiseyev, Visualization of branch points in $\mathcal{P}\mathcal{T}$-symmetric waveguides, Phys. Rev. Lett. 101, 080402 (2008).
  * Hodaei  _et al._ [2014] H. Hodaei, M.-A. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, Parity-time–symmetric microring lasers, Science 346, 975 (2014).
  * Feng  _et al._ [2014] L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, Single-mode laser by parity-time symmetry breaking, Science 346, 972 (2014).
  * Malzard  _et al._ [2015] S. Malzard, C. Poli, and H. Schomerus, Topologically protected defect states in open photonic systems with non-Hermitian charge-conjugation and parity-time symmetry, Phys. Rev. Lett. 115, 200402 (2015).
  * Shen  _et al._ [2018] H. Shen, B. Zhen, and L. Fu, Topological band theory for non-Hermitian Hamiltonians, Phys. Rev. Lett. 120, 146402 (2018).
  * Xia  _et al._ [2021a] S. Xia, D. Kaltsas, D. Song, I. Komis, J. Xu, A. Szameit, H. Buljan, K. G. Makris, and Z. Chen, Nonlinear tuning of $\mathcal{P}\mathcal{T}$ symmetry and non-Hermitian topological states, Science 372, 72 (2021a).
  * Chen  _et al._ [2017] W. Chen, Ş. Kaya Özdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature 548, 192 (2017).
  * Hodaei  _et al._ [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature 548, 187 (2017).
  * Hokmabadi  _et al._ [2019] M. P. Hokmabadi, A. Schumer, D. N. Christodoulides, and M. Khajavikhan, Non-Hermitian ring laser gyroscopes with enhanced sagnac sensitivity, Nature 576, 70 (2019).
  * Duggan  _et al._ [2022] R. Duggan, S. A. Mann, and A. Al$\grave{\rm u}$, Limitations of sensing at an exceptional point, ACS Photonics 9, 1554 (2022).
  * Longhi [2018] S. Longhi, $\mathcal{P}\mathcal{T}$ symmetry and antisymmetry by anti-Hermitian wave coupling and nonlinear optical interactions, Optics Letters 43, 4025 (2018).
  * Mandal and Bergholtz [2021] I. Mandal and E. J. Bergholtz, Symmetry and higher-order exceptional points, Phys. Rev. Lett. 127, 186601 (2021).
  * Delplace  _et al._ [2021] P. Delplace, T. Yoshida, and Y. Hatsugai, Symmetry-protected multifold exceptional points and their topological characterization, Phys. Rev. Lett. 127, 186602 (2021).
  * Kawabata  _et al._ [2019a] K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Topological unification of time-reversal and particle-hole symmetries in non-Hermitian physics, Nature Communications 10, 1 (2019a).
  * Kawabata  _et al._ [2019b] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019b).
  * Yoshida  _et al._ [2021] T. Yoshida, T. Mizoguchi, Y. Kuno, and Y. Hatsugai, Square-root topological phase with time-reversal and particle-hole symmetry, Phys. Rev. B 103, 235130 (2021).
  * Rivero and Ge [2021] J. D. H. Rivero and L. Ge, Chiral symmetry in non-Hermitian systems: Product rule and clifford algebra, Phys. Rev. B 103, 014111 (2021).
  * Sun  _et al._ [2011] K. Sun, Z. Gu, H. Katsura, and S. Das Sarma, Nearly flatbands with nontrivial topology, Phys. Rev. Lett. 106, 236803 (2011).
  * Kopnin  _et al._ [2011] N. B. Kopnin, T. T. Heikkilä, and G. E. Volovik, High-temperature surface superconductivity in topological flat-band systems, Phys. Rev. B 83, 220503 (2011).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Wang  _et al._ [2020] P. Wang, Y. Zheng, X. Chen, C. Huang, Y. V. Kartashov, L. Torner, V. V. Konotop, and F. Ye, Localization and delocalization of light in photonic moiré lattices, Nature 577, 42 (2020).
  * Vicencio  _et al._ [2015] R. A. Vicencio, C. Cantillano, L. Morales-Inostroza, B. Real, C. Mejía-Cortés, S. Weimann, A. Szameit, and M. I. Molina, Observation of localized states in Lieb photonic lattices, Phys. Rev. Lett. 114, 245503 (2015).
  * Mukherjee  _et al._ [2015] S. Mukherjee, A. Spracklen, D. Choudhury, N. Goldman, P. Öhberg, E. Andersson, and R. R. Thomson, Observation of a localized flat-band state in a photonic Lieb lattice, Phys. Rev. Lett. 114, 245504 (2015).
  * Whittaker  _et al._ [2018] C. E. Whittaker, E. Cancellieri, P. M. Walker, D. R. Gulevich, H. Schomerus, D. Vaitiekus, B. Royall, D. M. Whittaker, E. Clarke, I. V. Iorsh, I. A. Shelykh, M. S. Skolnick, and D. N. Krizhanovskii, Exciton polaritons in a two-dimensional Lieb lattice with spin-orbit coupling, Phys. Rev. Lett. 120, 097401 (2018).
  * Leykam and Flach [2018] D. Leykam and S. Flach, Perspective: photonic flatbands, APL Photonics 3, 070901 (2018).
  * Tang  _et al._ [2020] L. Tang, D. Song, S. Xia, S. Xia, J. Ma, W. Yan, Y. Hu, J. Xu, D. Leykam, and Z. Chen, Photonic flat-band lattices and unconventional light localization, Nanophotonics 9, 1161 (2020).
  * Poblete [2021] R. A. V. Poblete, Photonic flat band dynamics, Advances in Physics: X 6 (2021).
  * Weimann  _et al._ [2017] S. Weimann, M. Kremer, Y. Plotnik, Y. Lumer, S. Nolte, K. G. Makris, M. Segev, M. C. Rechtsman, and A. Szameit, Topologically protected bound states in photonic parity–time-symmetric crystals, Nature Materials 16, 433 (2017).
  * Biesenthal  _et al._ [2019] T. Biesenthal, M. Kremer, M. Heinrich, and A. Szameit, Experimental realization of $\mathcal{P}\mathcal{T}$-symmetric flat bands, Phys. Rev. Lett. 123, 183601 (2019).
  * Jin [2019] L. Jin, Flat band induced by the interplay of synthetic magnetic flux and non-Hermiticity, Phys. Rev. A 99, 033810 (2019).
  * Qi  _et al._ [2018] B. Qi, L. Zhang, and L. Ge, Defect states emerging from a non-Hermitian flatband of photonic zero modes, Phys. Rev. Lett. 120, 093901 (2018).
  * Longhi and Feng [2018] S. Longhi and L. Feng, Invited article: Mitigation of dynamical instabilities in laser arrays via non-Hermitian coupling, APL Photonics 3, 060802 (2018).
  * Leykam  _et al._ [2017] D. Leykam, S. Flach, and Y. D. Chong, Flat bands in lattices with non-Hermitian coupling, Phys. Rev. B 96, 064305 (2017).
  * Xia  _et al._ [2021b] S. Xia, C. Danieli, Y. Zhang, X. Zhao, H. Lu, L. Tang, D. Li, D. Song, and Z. Chen, Higher-order exceptional point and Landau–Zener Bloch oscillations in driven non-Hermitian photonic Lieb lattices, APL Photonics 6, 126106 (2021b).
  * Xia  _et al._ [2018] S. Xia, A. Ramachandran, S. Xia, D. Li, X. Liu, L. Tang, Y. Hu, D. Song, J. Xu, D. Leykam, S. Flach, and Z. Chen, Unconventional flatband line states in photonic Lieb lattices, Phys. Rev. Lett. 121, 263902 (2018).
  * Rhim  _et al._ [2020] J.-W. Rhim, K. Kim, and B.-J. Yang, Quantum distance and anomalous landau levels of flat bands, Nature 584, 59 (2020).
  * Ge [2018] L. Ge, Non-Hermitian lattices with a flat band and polynomial power increase, Photonics Research 6, A10 (2018).
  * Xiao  _et al._ [2020] Y.-X. Xiao, K. Ding, R.-Y. Zhang, Z. H. Hang, and C. T. Chan, Exceptional points make an astroid in non-Hermitian Lieb lattice: Evolution and topological protection, Phys. Rev. B 102, 245144 (2020).
  * Ramachandran  _et al._ [2017] A. Ramachandran, A. Andreanov, and S. Flach, Chiral flat bands: Existence, engineering, and stability, Phys. Rev. B 96, 161104 (2017).
  * Kolovsky  _et al._ [2018] A. R. Kolovsky, A. Ramachandran, and S. Flach, Topological flat wannier-stark bands, Phys. Rev. B 97, 045120 (2018).
  * Xia  _et al._ [2020] S. Xia, C. Danieli, W. Yan, D. Li, S. Xia, J. Ma, H. Lu, D. Song, L. Tang, S. Flach, _et al._ , Observation of quincunx-shaped and dipole-like flatband states in photonic rhombic lattices without band-touching, APL Photonics 5, 016107 (2020).
  * Kremer  _et al._ [2020] M. Kremer, I. Petrides, E. Meyer, M. Heinrich, O. Zilberberg, and A. Szameit, A square-root topological insulator with non-quantized indices realized with photonic Aharonov-Bohm cages, Nature Communications 11, 1 (2020).
  * Zhang and Jin [2020] S. M. Zhang and L. Jin, Localization in non-Hermitian asymmetric rhombic lattice, Phys. Rev. Research 2, 033127 (2020).
  * Xia  _et al._ [2021c] S. Xia, Y. Zhang, Z. Li, L. Qin, C. Yang, H. Lu, J. Zhang, X. Zhao, Z. Zhu, _et al._ , Band evolution and Landau-Zener Bloch oscillations in strained photonic rhombic lattices, Optics Express 29, 37503 (2021c).
  * Mukherjee  _et al._ [2018] S. Mukherjee, M. Di Liberto, P. Öhberg, R. R. Thomson, and N. Goldman, Experimental observation of Aharonov-Bohm Cages in photonic lattices, Phys. Rev. Lett. 121, 075502 (2018).
  * Mukherjee and Thomson [2015] S. Mukherjee and R. R. Thomson, Observation of localized flat-band modes in a quasi-one-dimensional photonic rhombic lattice, Optics Letters 40, 5443 (2015).
  * Alexeeva  _et al._ [2014] N. V. Alexeeva, I. V. Barashenkov, K. Rayanov, and S. Flach, Actively coupled optical waveguides, Phys. Rev. A 89, 013848 (2014).
  * Longhi [2016] S. Longhi, Non-Hermitian tight-binding network engineering, Phys. Rev. A 93, 022102 (2016).
  * Longhi  _et al._ [2015a] S. Longhi, D. Gatti, and G. Della Valle, Non-Hermitian transparency and one-way transport in low-dimensional lattices by an imaginary gauge field, Phys. Rev. B 92, 094204 (2015a).
  * Ding  _et al._ [2021] L. Ding, Z. Lin, S. Ke, B. Wang, and P. Lu, Non-Hermitian flat bands in rhombic microring resonator arrays, Optics Express 29, 24373 (2021).
  * Makris  _et al._ [2010] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, $\mathcal{P}\mathcal{T}$-symmetric optical lattices, Physical Review A 81, 063807 (2010).
  * Schnabel  _et al._ [2017] J. Schnabel, H. Cartarius, J. Main, G. Wunner, and W. D. Heiss, $\mathcal{P}\mathcal{T}$-symmetric waveguide system with evidence of a third-order exceptional point, Phys. Rev. A 95, 053868 (2017).
  * Longhi  _et al._ [2015b] S. Longhi, D. Gatti, and G. D. Valle, Robust light transport in non-Hermitian photonic lattices, Scientific Rep. 5, 1 (2015b).
  * Longhi [2021] S. Longhi, Phase transitions in a non-Hermitian Aubry-André-Harper model, Phys. Rev. B 103, 054203 (2021).
  * Zeng  _et al._ [2019] C. Zeng, Y. Sun, G. Li, Y. Li, H. Jiang, Y. Yang, and H. Chen, Enhanced sensitivity at high-order exceptional points in a passive wireless sensing system, Optics Express 27, 27562 (2019).
  * Zhang  _et al._ [2021] L. Zhang, Y. Yang, Z. Jiang, Q. Chen, Q. Yan, Z. Wu, B. Zhang, J. Huangfu, and H. Chen, Demonstration of topological wireless power transfer, Science Bulletin 66, 974 (2021).
  * Yu  _et al._ [2020] S. Yu, Y. Meng, J.-S. Tang, X.-Y. Xu, Y.-T. Wang, P. Yin, Z.-J. Ke, W. Liu, Z.-P. Li, Y.-Z. Yang, G. Chen, Y.-J. Han, C.-F. Li, and G.-C. Guo, Experimental investigation of quantum $\mathcal{P}\mathcal{T}$-enhanced sensor, Phys. Rev. Lett. 125, 240506 (2020).
  * Ezawa [2022] M. Ezawa, Nonlinear non-Hermitian higher-order topological laser, Phys. Rev. Research 4, 013195 (2022).
  * Hang  _et al._ [2021] C. Hang, D. A. Zezyulin, G. Huang, and V. V. Konotop, Nonlinear topological edge states in a non-Hermitian array of optical waveguides embedded in an atomic gas, Phys. Rev. A 103, L040202 (2021).
  * Pyrialakos  _et al._ [2022] G. G. Pyrialakos, H. Ren, P. S. Jung, M. Khajavikhan, and D. N. Christodoulides, Thermalization dynamics of nonlinear non-Hermitian optical lattices, Phys. Rev. Lett. 128, 213901 (2022).
  * Komis  _et al._ [2022] I. Komis, D. Kaltsas, S. Xia, H. Buljan, Z. Chen, and K. Makris, Robustness versus sensitivity in non-Hermitian topological lattices probed by pseudospectra, arXiv:2203.02519 (2022).
  * Ding  _et al._ [2016] K. Ding, G. Ma, M. Xiao, Z. Q. Zhang, and C. T. Chan, Emergence, coalescence, and topological properties of multiple exceptional points and their experimental realization, Phys. Rev. X 6, 021007 (2016).


# Long Term Drought Prediction using Deep Neural Networks based on Geospatial
Weather Data

Vsevolod Grabar  Alexander Marusov  Yury Maximov  Nazar Sotiriadi  Alexander
Bulkin  Alexey Zaytsev

###### Abstract

The accurate prediction of drought probability in specific regions is crucial
for informed decision-making in agricultural practices. In particular, for
long-term decisions it is important to make predictions for one year ahead.
However, forecasting this probability presents challenges due to the complex
interplay of various factors within the region of interest and neighboring
areas. In this study, we propose an end-to-end solution based on various
spatio-temporal neural networks to address this issue. Considered models focus
on predicting the drought intensity based on Palmer Drought Severity Index
(PDSI) for subregions of interest, leveraging intrinsic factors and insights
from climate models to enhance drought predictions.

Comparative evaluations demonstrate the superior accuracy of Convolutional
LSTM and Transformer models compared to baseline Gradient Boosting and
Logistic Regression solutions. The two former models achieved impressive
$ROC\,AUC$ scores from $0.90$ to $0.70$ for forecast horizons from 1 to 6
months, outperforming baseline models. The transformer showed superiority for
shorter horizons, while ConvLSTM took the lead for longer horizons. Thus, we
recommend selecting the model in this way for long-term drought forecasting.

To ensure the broad applicability of considered models, we conduct extensive
validation across regions worldwide, considering different environmental
conditions. Moreover, we also run several ablation and sensitivity studies to
challenge our findings and provide additional information on an appropriate
way to solve the problem at hand.

###### keywords:

weather , climate , drought forecasting , deep learning , long-term
forecasting

††journal: Computers and Electronics in Agriculture

[inst1]organization=Skolkovo Institute of Science and Technology

[inst2]organization=Los Alamos National Laboratory, Theoretical Division,
[inst3]organization=Sberbank of Russia PJSC, [inst4]organization=Moscow State
University,

##  1 Introduction

Drought forecasting is one of the crucial problems nowadays [1]. Indeed,
monitoring and forecasting droughts hold significant importance due to their
high occurrence across diverse terrains [2]. These natural climate events are
costly and have far-reaching impacts on populations and various economic
sectors [3]. Additionally, it is worth noting that global climate change will
likely increase the probability of droughts [4]. So, drought forecasting is an
essential but complicated task since the drought’s beginning, ending, and
duration are often hard to predict [5].

A quantitative measure of the drought is needed to solve the task described
above. However, it is well-known that drought depends on many other climatic
factors such as temperature and precipitation [6]. Hence, the drought measure
(indices) can be calculated in many ways. Ones of the most fundamental and
mostly known indices [7] are the Standardized Precipitation Index (SPI) [8]
and Palmer Drought Severity Index (PDSI) [9]. In our studies, we use monthly
PDSI for several reasons. First of all, our forecasting horizon is quite long
- 12 months. According to recent work [10], PDSI is a good choice for long-
term forecasting. Also, PDSI has a long recording history and allows us to
consider the effect of global warming, which is quite a significant problem in
our days [11].

More and more machine learning algorithms are involved in solving drought
forecasting tasks. Regarding classic approaches, in [5] prominent stochastic
models, AutoRegressive Integrated Moving Average (ARIMA) and multiplicative
Seasonal AutoRegressive Integrated Moving Average (SARIMA), were exploited to
predict drought via SPI index. Multivariate regression having a history of
PDSI and some other global climate indexes on input was used to predict the
PDSI index in South-Central Oklahoma [10]. Gradient boosting algorithm shown
effectiveness for geospatial data [12, 13] and can handle imbalanced problems
[14] encountered in drought prediction. Besides, this algorithm has proved its
power in drought classification in Turkish regions [15]. Authors [16] used a
logistic regression model for binary drought classification task via SPI
index.

Moving on to deep learning methods, the authors in [17] proposed to use a
recursive multistep neural network and direct multistep neural network for
drought forecasting via the SPI index. One of the most prominent approaches
for time series data is Recurrent Neural Network (RNN) [18]. It is known that
RNNs are often better than classical approaches [19]. Indeed, research [20]
shows that Long Short-Term Memory (LSTM) [21] outperforms ARIMA for long-term
forecasting of SPI index, but in short-term prediction, ARIMA has quite well
results. To exploit the advantages of ARIMA and LSTM, authors of [22] proposed
to use a hybrid ARIMA-LSTM model. The approaches described above used only
historical (temporal) information about the data. Nevertheless, besides
temporal information, the data has spatial dependencies. One of the most
famous methods to account for both temporal and spatial information is
ConvLSTM [23]. This method applies to many tasks, such as earthquake
prediction [24]. Authors of [7] used ConvLSTM for short-term drought
forecasting (8 days) of satellite-based drought indices - Scaled Drought
Condition Index (SDCI) and SPI. Another idea is to use Transformer
architectures. Initially developed for Natural Language Processing (NLP) tasks
[25] now these models are widely used in a wide variety of domains, e.g.,
processing video frames in Computer Vision (CV) [26]. Consequently, adopting
attention-based architectures to spatio-temporal modeling is a natural idea.
One of the most known and modern weather transformers are EarthFormer [27] and
FourCastNet [28]. These models solve different regression tasks (e.g.,
precipitation nowcasting) and show superior performance for various spatio-
temporal problems. We adopted these architectures to our task.

The applications of deep learning models in conjunction with geospatial data
extend beyond drought forecasting and can be effectively employed in diverse
weather scenarios. Unlike computationally intensive numerical weather
prediction methods based on solving partial differential equations, recent
research by [29] demonstrates the enhanced capabilities of transformer-based
models in delivering accurate medium-range predictions for various upper-air
and surface variables. Furthermore, [30] successfully combined multiple deep-
learning models to rectify short-term wind predictions derived from
conventional numerical algorithms. In the realm of meteorological forecasting
complexity, [31] harnessed the power of Generative Adversarial Networks (GANs)
and Variational Autoencoder GANs (VAE-GANs) to generate spatial forecasts of
precipitation in the UK. Similarly, in the context of the contiguous United
States, [32] showcased the superior performance of a 3D convolutional network
in short-term precipitation predictions. The existing literature highlights
the strengths of deep learning models for generating short-term and medium-
term forecasts, spanning from durations of several hours to several days or
even weeks. This diverse range of applications underscores the adaptability
and potential of deep learning in enhancing meteorological prediction
methodologies.

Besides weather forecasting problems, deep learning models for geospatial data
have found many applications. To begin with, this type of data could be
transformed into images and then used to forecast traffic with CNNs, as in
[33]. Naturally, one would exploit spatial dependencies to predict epidemic
propagation, as in the case of the Madrid region and COVID-19 in the paper
[34]. More recently, [35] reviewed various machine learning methods for
forecasting dengue risk, including CNNs and RNNs. Next, in the paper by [36],
geospatial LSTMs were exploited to predict invasive species outbreaks.
Finally, various convolutional and recurrent neural network architectures were
successfully exploited in crop yield forecasts, thus immediately providing
economic value for the agricultural sector.

Our paper considers different classic (Gradient boosting, logistic regression)
and deep-learning (ConvLSTM, EarthFormer, and FourCastNet) methods for long-
term drought forecasting using PDSI values. The main contribution of our
research are:

  1. 1.

Deep learning models solve the problem of medium (up to 12 months) drought
forecasting. We formulate the problem as the prediction of the probability of
a drought.

  2. 2.

Medium-term forecasting (up to 6 months ahead) is better to provide via
EarthFormer meanwhile Long-term forecasting is quite suitable for ConvLSTM.
So, we recommend to use them accordingly.

  3. 3.

Surprisingly, we observed pretty good performance of Logistic regression and
Gradient Boosting models, especially in Short-term forecasting, if correct
features are used.

Detailed explanations and intuition on achieved results are in 4.1 section.

The list of key novelties of our work is the following:

  1. 1.

We adopted and explored traditional spatio-temporal methods (ConvLSTM,
gradient boosting, logistic regression) as modern approaches (EarthFormer,
FourCastNet). According to our knowledge, we are the first to use these models
for long-term forecasting of PDSI values.

  2. 2.

Unlike many previous works, which tried to predict PDSI value directly
(regression task), we treated the drought forecasting task as a binary
classification task. We chose $-2$ as threshold according to [10]. All PDSI
values below $-2$ are interpreted as a drought. So, now a model report a
probability of a drought, making it more suitable for estimation of financial
risk for particular

  3. 3.

We provided extensive experiments using freely accessible data from diverse
regions and continents.

##  2 Data

We have used publicly available geospatial data from Google Earth Engine [37].
Specifically, for obtaining the PDSI data, we employed the TerraClimate
Monthly dataset [38]. Our PDSI data encompasses a comprehensive range of
climatic values spanning from 1958 to 2022, covering the entirety of Earth. To
test the consistency of considered models, we look at regions across
continents and various climate zones - from the state of Missouri and Eastern
Europe to India. The considered regions are depicted in Figure 1. Their
characteristics are given in Table 1.

Figure 1: Regions chosen for PDSI forecast

The input data is downloaded as a tif file and later restructured as a 3D
tensor, with one dimension representing the time scale (monthly intervals). In
contrast, the other two dimensions correspond to spatial coordinates (x and y)
of the grid. The resolution of the regional dataset varies, where grid
dimensions range from 30 by 60 up to 128 by 192 (spatial resolution of a
single cell of TerraClimate data is roughly 5 km), accommodating different
levels of granularity and spatial detail.

Region | Span, |  $\%$ of normal |  $\%$ of drought | Spatial
---|---|---|---|---
| months |  PDSI $\geq-2$ |  PDSI $\leq-2$ | dimensions, km
Missouri, USA | 754 | 74.91 | 25.09 | 416x544
Madhya Pradesh, India | 754 | 70.66 | 29.34 | 512x768
Goias, Brazil | 754 | 68.97 | 31.03 | 640x640
Northern Kazakhstan | 742 | 68.70 | 31.30 | 256x480
Eastern Europe | 742 | 66.28 | 33.72 | 352x672
Table 1: Regions’ summary statistics

##  3 Methods

We compare deep learning approaches, including Convolutional LSTM and novel
transformer models, such as FourCastNet from Nvidia and EarthFormer from
Amazon, with classic methods, including baseline model, gradient boosting, and
logistic regression.

###  3.1 Baseline

As a global baseline and a coherence check, we took the most prevalent class
from train data and compared it with actual targets from the test subset. We
also checked a rolling baseline - i.e., the most frequent class from recent
history (from 6 to 24 months), but the results were almost indistinguishable
from the global baseline, so we kept them out of our tables and graphs.

###  3.2 Basic methods: Logistic regression and Gradient boosting

Both logistic regression and gradient boosting cannot work with raw data
itself as it is. Hence, we created a data module that treated each grid cell
as an individual value and transformed our task into a typical time series
forecasting problem. To benefit from spatial correlations, we incorporate
values from neighboring cells. For example, if we consider a 3x3 neighborhood,
this includes eight additional time series. It is important to note that for
”edge” cells, some of the neighboring cells may contain all zeros due to data
limitations.

##### Logistic Regression

Logistic regression is usually the natural choice for tasks with linear
dependence or as a baseline model. The novel research [39] shows that time
series linear models are a good choice.

##### Gradient boosting

We adopted the gradient boosting of decision trees, implemented using the
well-established library XGBoost [40]. XGBoost, renowned for its speed and
efficiency across a wide range of predictive modeling tasks, has consistently
been favored by data science competition winners. It operates as an ensemble
of decision trees, with new trees aimed at rectifying errors made by existing
trees in the model. Trees are successively added until no further improvements
can be made.

###  3.3 Transformer-based methods

##### FourCastNet

FourCastNet (short for Fourier ForeCasting Neural Network) is a weather
forecasting model developed by authors from Nvidia [28]. It constitutes part
of their Modulus Sym deep learning framework for solving various applied
physics tasks. It combines Adaptive Fourier Neural Operator (AFNO) from [41]
with Vision Transformer, and the model is both computationally, reducing
spatial mixing complexity to $O(N\log N)$, and memory efficient, where $N$ is
the sequence length. In the original papers, authors produced high-resolution
short-term forecasts of wind speed and precipitation. We modified the last
layer of the model to switch it from regression to classification task and
evaluate it for the long-term forecasting problem.

##### EarthFormer

Vanilla transformers have $O(N^{2})$ complexity, and consequently, it is hard
to apply them to weather spatio-temporal data because of their large
dimensionality. Authors [27] propose to use the ”divide and conquer” method:
they divide original data into non-intersecting parts (called cuboids) and use
a self-attention mechanism on each cuboid separately in parallel. Such an
approach allows significantly reduced complexity. The authors introduced
EarthFormer for regression tasks, so we adopted it for a classification
problem.

###  3.4 ConvLSTM

This model draws inspiration from the work presented in the paper by Kail et
al. (2021) [24], and it represents a modified version of the Convolutional
LSTM architecture proposed by Shi et al. (2015) [23]. To capture temporal
dependencies, we adopt Recurrent Neural Networks (RNNs). Specifically, we
employ Long Short-Term Memory (LSTM), a type of RNN that utilizes an
additional state cell to facilitate long-term memory retention. We extend the
traditional one-dimensional hidden states of LSTM to two-dimensional feature
maps, enabling grid-to-grid transformations, which are crucial for our task.
We process grids of PDSI values with various dimensions, ranging from 9x16 to
40x200.

To exploit spatial dependencies among drought severity values, we incorporate
Convolutional Neural Networks (CNNs), well-suited for image and two-
dimensional signal processing, including geospatial data. This approach
combines the strengths of both RNN and CNN architectures by passing
information through the RNN component as a feature map obtained using CNN
operations. This integration allows us to capture and utilize temporal and
spatial information within the following architecture.

##### Details of architecture

Convolutional LSTM follows the pipeline:

  1. 1.

We represent data as a sequence of grids: for each cell, we specify a value of
PDSI for a particular month; the input grid at each time moment has dimension
$grid_{h}\times grid_{w}$ (varying from 50x50 to 200x200 for different regions
of interest).

  2. 2.

We pass the input grid through a convolutional network to create an embedding
of grid dimensionality size with 16 channels. As an output of LSTM at each
time moment, we have a hidden representation (short-term memory) of size
$h\times grid_{h}\times grid_{w}$, cell (long-term memory) representation of a
similar size, and the output of size $h\times grid_{h}\times grid_{w}$.

  3. 3.

We transform the output to the size $1\times grid_{h}\times grid_{w}$ using
convolution $1\times 1$ to receive probabilities of the drought for each cell
as a final prediction (or to $k\times grid_{h}\times grid_{w}$, where $k>2$ is
a number of classes of drought conditions that we are trying to predict). As
an additional hyperparameter, we vary the forecasting horizon - i.e., we can
forecast PDSI for the next month or $k$-th month.

##  4 Results

##### Description of experiment

For a binary classification problem, we arbitrarily set a drought threshold as
$-2$, which can be justified by PDSI bins from Table 2 and later adjusted for
further experiments. For each cell in a considered region, the problem to
solve is a binary classification problem with the true value is the indicator
that PDSI exceeds $-2$. By construction, the model outputs the probability of
a serious drought from the model output.

For validation purposes, we have divided our dataset into the train (70%) and
test (30%) subsets for all of the five world regions (Missouri, Northern
Kazakhstan, Madhya Pradesh, Eastern Europe, and Goias). The model is trained
via a train subset of the data, the quality metrics are calculated via an
independent test subset. The validation is out-of-time: the test subset
follows the training subset in time.

##### Evaluation procedure

We use $ROC\,AUC$, $PR\,AUC$ and $F1$ score to evaluate the model. During
validation for early stopping and hyperparameter optimization, we chose
$ROC\,AUC$ score. All these scores are the medians - i.e., concatenating
spatial prediction, we receive a temporal vector at every cell of spatial
prediction. Next, we receive a single score for each cell, so we end up with a
grid full of metrics, and finally, we compute a median. Higher values of all
scores correspond to better models. $ROC\,AUC$ scores has the perfect value
being $1$ and the value for a random prediction being $0.5$.

##### Compared methods

For our main experiments, we have exploited the performance of baseline (most
frequent class from historical data), XGBoost (gradient boosting of random
forests), LogReg (logistic regression), Convolutional LSTM (ConvLSTM),
FourCastNet (transformer model from Nvidia) and Eartformer (transformer model
from Amazon).

XGBoost and LogReg are two basic algorithms often used as strong baselines.
The last three are variants of neural networks that showed strong performance
in various geospatial problems. They represent two dominating architectures in
geospatial modeling: ConvLSTM is a combination of recurrent and convolutional
neural networks, FourCastNet and Eartformer

###  4.1 Main results

PDSI value | Drought severity class
---|---
4.00 and above | Extreme wet spell
3.00-3.99 | Severe wet spell
2.00-2.99 | Moderate wet spell
1.00-1.99 | Mild wet spell
-1.00 to 0.99 | Normal
-1.00 to -1.99 | Mild dry spell
-2.00 to -2.99 | Moderate dry spell
-3.00 to -3.99 | Severe dry spell
-4.00 and below | Extreme dry spell
Table 2: Classification of various PDSI values, source [42]

##### Analysis of results

Our findings indicate that ConvLSTM model outperformed the competitors,
including the gradient boosting tree algorithm, achieving an impressive
$ROC\,AUC$ score of $0.9$ for a one-month prediction. Notably, ConvLSTM
exhibits a gradual decline in performance, reaching $0.6-0.65$ for longer
forecasting horizons (ranging from 1 month to 12 months). The standard
gradient boosting approach initially yielded a similar $ROCAUC$ score of $0.9$
but sharply dropped to $0.5$ as the forecasting horizon extended. EarthFormer
outperforms other approaches for shorter horizons and falling short of
ConvLSTM at longer horizons of 9-12 months. These results are visually
depicted in Figure 2. Additionally, we present the results for six-month
prediction by regions in Figure 4, where the variation in scores for different
geographies across models is visible.

Why does transformer fail in long-term prediction? Our assumption of such
behavior is the permutation-invariance of the attention mechanism. Despite
positional encoding, transformers can’t effectively extract temporal
information from long input sequences. Since a long input sequence is
essential for long-term forecasting, transformers don’t show good results.
Similar results for different time-series tasks were observed in [39].
ConvLSTM naturally extracts temporal information via LSTM.

Why is logistic regression impressive? Since Gradient boosting results are
almost identical to Logistic regression, we want to discuss only Logistic
regression performance. First of all, the power of linear models was already
shown in [39], where they beat modern Transformer architectures on almost all
datasets. In our experiments, linear models are worse than other models on
Long-term prediction, but on Short-term scale, we can see quite comparable
results. We tried different history lengths, but our results show that it is
enough to take only the nearest to the future horizons element. Our intuition
is that the nearest future predictor variables are closely (particularly
linearly) related to the history element. For example, PDSI in July is close
to PDSI in August but far away from December. Hence, linear models are good at
Short-term predictions but poor at Long-term forecasting.

Figure 2: Median metrics: ROC-AUC, PR-AUC, F1 from left to right for different forecast horizons averaged over five considered regions, binary drought classification Horizon, months | 1 | 3 | 6 | 9 | 12
---|---|---|---|---|---
Median ROC AUC: |  |  |  |  |
Baseline | 0.5 | 0.5 | 0.5 | 0.5 | 0.5
LogReg | 0.886 | 0.774 | 0.640 | 0.546 | 0.518
XGBoost | 0.878 | 0.754 | 0.628 | 0.568 | 0.542
FourCastNet | 0.881 | 0.711 | 0.624 | 0.561 | 0.536
EarthFormer | 0.948 | 0.840 | 0.690 | 0.556 | 0.480
ConvLSTM | 0.887 | 0.802 | 0.693 | 0.650 | 0.617
Median PR AUC: |  |  |  |  |
Baseline | 0.355 | 0.354 | 0.354 | 0.355 | 0.356
LogReg | 0.766 | 0.598 | 0.470 | 0.394 | 0.360
XGBoost | 0.752 | 0.574 | 0.44 | 0.39 | 0.372
FourCastNet | 0.776 | 0.546 | 0.455 | 0.402 | 0.382
EarthFormer | 0.880 | 0.650 | 0.514 | 0.438 | 0.362
ConvLSTM | 0.772 | 0.689 | 0.565 | 0.505 | 0.452
Median F1: |  |  |  |  |
Baseline | 0.645 | 0.646 | 0.646 | 0.645 | 0.644
LogReg | 0.846 | 0.698 | 0.480 | 0.226 | 0.094
XGBoost | 0.836 | 0.674 | 0.480 | 0.366 | 0.292
FourCastNet | 0.831 | 0.603 | 0.460 | 0.366 | 0.314
EarthFormer | 0.816 | 0.604 | 0.448 | 0.226 | 0.178
ConvLSTM | 0.784 | 0.698 | 0.600 | 0.558 | 0.543
Table 3: Median Metrics vs Forecast Horizon, binary classification; best values are in bold, second best are underlined Region | Northern | Eastern | Madhya | Goias | Missouri
---|---|---|---|---|---
| Kazakhstan | Europe | Pradesh |  |
median ROC AUC: |  |  |  |  |
Baseline | 0.5 | 0.5 | 0.5 | 0.5 | 0.5
LogReg | 0.60 | 0.63 | 0.61 | 0.61 | 0.75
XGBoost | 0.59 | 0.64 | 0.59 | 0.61 | 0.71
FourCastNet | 0.60 | 0.55 | 0.65 | 0.61 | 0.69
EarthFormer | 0.54 | 0.69 | 0.84 | 0.65 | 0.73
ConvLSTM | 0.71 | 0.68 | 0.71 | 0.60 | 0.75
median PR AUC: |  |  |  |  |
Baseline | 0.37 | 0.43 | 0.35 | 0.39 | 0.23
LogReg | 0.46 | 0.67 | 0.24 | 0.55 | 0.43
XGBoost | 0.46 | 0.63 | 0.20 | 0.54 | 0.37
FourCastNet | 0.46 | 0.55 | 0.36 | 0.54 | 0.37
EarthFormer | 0.47 | 0.22 | 0.77 | 0.52 | 0.59
ConvLSTM | 0.55 | 0.65 | 0.57 | 0.54 | 0.50
median F1: |  |  |  |  |
Baseline | 0.63 | 0.57 | 0.65 | 0.61 | 0.77
LogReg | 0.42 | 0.62 | 0.35 | 0.41 | 0.60
XGBoost | 0.45 | 0.58 | 0.30 | 0.54 | 0.53
FourCastNet | 0.42 | 0.39 | 0.46 | 0.51 | 0.52
EarthFormer | 0.03 | 0.18 | 0.80 | 0.55 | 0.65
ConvLSTM | 0.65 | 0.64 | 0.64 | 0.52 | 0.56
Table 4: Median Metrics vs Region, binary classification, 6 months horizon;
best values are in bold, second best are underlined

###  4.2 Predictions and errors for a particular region

To assess the performance of the Convolutional LSTM algorithm (which proved to
be the most stable and promising for drought forecasting), we focused on the
region of Missouri state, where we will run several ablation studies. As an
illustration, the spatial distribution of $ROCAUC$ scores is depicted in
Figure 3. Notably, we observed a non-uniform distribution of $ROCAUC$ values
across the cells within the region. The standard deviation of the scores is
substantial, and individual values range from those close to random predictors
($ROC\,AUC=0.6$) to near-perfect scores approaching $1.0$. This variability
highlights the diverse predictive capabilities of our algorithm across
different spatial locations within Missouri.

Figure 3: Spatial distribution of $\mathrm{ROC\,AUC}$ for 6 month forecast,
Missouri, ConvLSTM

####  4.2.1 Performance Evaluation for Cropped Region

##### Description of experiment

As is typical with ROC AUC maps, the worst predictions are found on the edges
and some corners. We have observed that this behavior is consistent regardless
of the region being studied. Consequently, making predictions for a larger
region and cropping the desired region of interest may be advantageous. We
have conducted a study to test this hypothesis for the attached map, and the
results are shown in Table 5.

Percent of map cropped | 0 | 10 | 20 | 30 | 40
---|---|---|---|---|---
median ROC AUC | 0.7525 | 0.7592 | 0.7665 | 0.7749 | 0.7834
Percent of map cropped | 50 | 60 | 70 | 80 | 90
median ROC AUC | 0.7886 | 0.7899 | 0.7880 | 0.7838 | 0.7825
Table 5: ROC AUC score vs crop percentage, 6 month forecast, ConvLSTM model
for Missouri

##### Analysis of results

Based on the findings of this specific experiment, we deduce that cropping
approximately 40-50% of the initially selected region maximizes our score. In
other words, choosing a region that is initially 1.6-2 times larger than our
target region is advisable. However, the precise amount of zoom required for
optimal results has to be determined through further experiments.
Next, we conducted several similar experiments to investigate how the model
predictions change with the decrease in the total squared area. For this, we
took the same geographic region, Missouri, and examined various combinations
of history length and forecast horizon. We trained a new checkpoint for each
variant of history length, forecast horizon, and region area (varying from
entire state to about a quarter of it).

The results are summarized in Table 6, and in more detail, they are presented
in Tables 8, 9, 10, 11. We observed that predicting at a smaller area with
pre-trained model from a larger area generally works better. However, the
degree of improvement is marginal, usually not exceeding 0.5-1%. Figure 4
presents the evolution of spatial maps.

Region area | 100% | 75% | 53% | 27%
---|---|---|---|---
$h=6$, $f=1$ | 0.9292 | 0.9466 | 0.9432 | 0.9326
$h=12$, $f=3$ | 0.8746 | 0.8525 | 0.8595 | 0.8710
$h=9$, $f=6$ | 0.7525 | 0.6901 | 0.7597 | 0.7544
$h=24$, $f=12$ | 0.6117 | 0.5890 | 0.6126 | 0.6264
Table 6: ROC AUC score vs zoomed region area (h - length of input history, f -
forecast horizon)

####  4.2.2 Standard deviations and ensembling

##### Description of experiment

To assess the consistency of our findings, we repeated the experiment of PDSI
binary classification for Missouri state with five random seeds for each
configuration of history length and forecast horizon. In addition, we tested
an averaged ensemble of these five trained models.

##### Analysis of results

The results for horizons of 1, 3, 5, 6, 9, and 12 months can be found in
Figure 5. Ensemble scores are denoted by the red crosses on these plots.
Overall, the numbers are varied, and there is no definitive optimal choice of
history length for any horizon. However, we observed that extreme values, such
as the shortest or longest history lengths (e.g., 6, 9, 21, 24 months), yield
better performance. Notably, in most cases, the averaged ensemble of models
outperforms the individual underlying models. We note, that this effect is
more sound for neural networks, as they provide more diverse predictions, even
if there are no difference in the architecture and the only difference is a
starting point for training [43]. Contrary, logistic regression ensemble is
similar to a single model, as the optimization problem is convex [44], and
Gradient boosting is an ensemble per se [15]. Obtained ensembles can also be
used to access the uncertainty of predictions by machine learning models
improving the decision-making process [45].

Figure 4: Evolution of spatial maps with zooming in: top to bottom area
decreases from 100% of region to 27%, left to right forecast horizons are 1,
3, 6 and 12 months correspondingly. Figure 5: Variation of ROC-AUC scores for
five different random seeds; ConvLSTM binary classification for Missouri
state. Metrics for an ensemble of five models are marked as red crosses.

##  5 Conclusion

Droughts have emerged as severe natural disasters that significantly impact
the economy, influencing sectors such as agriculture and the well-being of
populations. The frequency and intensity of droughts have been amplified
worldwide due to climate change, as exemplified by the summer of 2022 in the
Northern Hemisphere.

Accurate prediction of future drought occurrences and adequate preparation for
their consequences is crucial in mitigating the adverse impacts of climate
change. Our research study has shown substantial improvements in forecasting
capabilities for PDSI indicator, including the Palmer Drought Severity Index
(PDSI), by combining convolutional neural networks and recurrent neural
networks. This Convolutional LSTM model surpasses the performance of standard
algorithms, such as gradient boosting, and even more advanced transformer
models for long-term predictions.

By enhancing our understanding and predictive abilities related to droughts,
we can proactively address the challenges of climate change. We hope that the
insights gained from this study will contribute to better preparedness and
mitigation strategies to alleviate the adverse effects of droughts on the
economy and society.

## References

  * [1] S. Mohammed, et al., A comparative analysis of data mining techniques for agricultural and hydrological drought prediction in the eastern mediterranean, Computers and Electronics in Agriculture 197 (2022) 1–19.
  * [2] A. Ghozat, et al., A novel intelligent approach for predicting meteorological drought based on satellite-based precipitation product: Application of an emd-dfa-dbn hybrid model, Computers and Electronics in Agriculture 211 (2023).
  * [3] M. Ali, et al., Multi-stage committee based extreme learning machine model incorporating the influence of climate parameters and seasonality on drought forecasting, Computers and Electronics in Agriculture 152 (2018) 149–165.
  * [4] C. Xiujia, et al., Application of wnn-pso model in drought prediction at crop growth stages: A case study of spring maize in semi-arid regions of northern china, Computers and Electronics in Agriculture 199 (2022).
  * [5] A. K. Mishra, V. R. Desai, Drought forecasting using stochastic models, Stochastic environmental research and risk assessment 19 (2005) 202–216.
  * [6] M. Özger, et al., Comparison of wavelet and empirical mode decomposition hybrid models in drought prediction, Computers and Electronics in Agriculture 179 (2020).
  * [7] S. e. a. Park, Short-term forecasting of satellite-based drought indices using their temporal patterns and numerical model output, Remote Sensing 12 (21) (2020) 3499.
  * [8] T. B. McKee, N. J. Doesken, J. Kleist, The relation of drought frequency and duration to time scales, Proceedings of the Eighth Conference on Applied Climatology (1993) 179–184.
  * [9] W. M. Alley, The Palmer drought severity index: limitations and assumptions, Journal of Applied Meteorology and Climatology 23 (7) (1984) 1100–1109.
  * [10] I. L. C. McPherson, Renee A., M. B. Richman, A place‐based approach to drought forecasting in south‐central oklahoma, Earth and Space Science 9 (10) (2022).
  * [11] A. Dai, Characteristics and trends in various forms of the palmer drought severity index during 1900–2008, Journal of Geophysical Research: Atmospheres 116 (2011).
  * [12] P. Proskura, A. Zaytsev, I. Braslavsky, E. Egorov, E. Burnaev, Usage of multiple rtl features for earthquakes prediction, in: International Conference on Computational Science and Its Applications, Springer, 2019, pp. 556–565.
  * [13] D. Koldasbayeva, P. Tregubova, D. Shadrin, M. Gasanov, M. Pukalchik, Large-scale forecasting of heracleum sosnowskyi habitat suitability under the climate change on publicly available data, Scientific reports 12 (1) (2022) 6128\.
  * [14] N. Kozlovskaia, A. Zaytsev, Deep ensembles for imbalanced classification, in: 2017 16th IEEE International Conference on Machine Learning and Applications (ICMLA), IEEE, 2017, pp. 908–913.
  * [15] A. Danandeh Mehr, Drought classification using gradient boosting decision tree, Acta Geophysica 69 (2021) 909–918.
  * [16] R. Niaz, et al, Logistic regression analysis for spatial patterns of drought persistence, Complexity (2021) 1–13.
  * [17] A. K. Mishra, V. R. Desai, Drought forecasting using feed-forward recursive neural network, Ecological Modelling 198 (2006) 127–138.
  * [18] D. Rumelhart, et al., Learning internal representations by error propagation, California Univ San Diego La Jolla Inst for Cognitive Science (1985).
  * [19] H. Hewamalage, et al., Recurrent neural networks for time series forecasting: current status and future directions, IJF (2021).
  * [20] S. Poornima, M. Pushpalatha, Drought prediction based on spi and spei with varying timescales using lstm recurrent neural network, Soft Computing (2019) 8399–8412.
  * [21] J. Schmidhuber, S. Hochreiter, Long short-term memory, Neural Comput (1997) 1735–1780.
  * [22] D. Xu, et al., Application of a hybrid arima-lstm model based on the spei for drought forecasting, Environmental Science and Pollution Research (2022) 4128–4144.
  * [23] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, W.-c. Woo, Convolutional lstm network: A machine learning approach for precipitation nowcasting, Advances in neural information processing systems 28 (2015).
  * [24] R. Kail, E. Burnaev, A. Zaytsev, Recurrent convolutional neural networks help to predict the location of earthquakes, IEEE Geoscience and Remote Sensing Letters 19 (2021) 1–5.
  * [25] A. Vaswani, et al., Attention is all you need, Advances in Neural Information Processing Systems 30 (2017).
  * [26] C. Zheng, et al., 3d human pose estimation with spatial and temporal transformers, Proceedings of the IEEE/CVF International Conference on Computer Vision (2021) 11656–11665.
  * [27] Z. Gao, et al., Earthformer: Exploring space-time transformers for earth system forecasting, Advances in Neural Information Processing Systems 35 (2022).
  * [28] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, P. Hassanzadeh, K. Kashinath, A. Anandkumar, Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators (2022).  arXiv:2202.11214.
  * [29] Accurate medium-range global weather forecasting with 3d neural networks, Nature 619 (2023) 533–538.
  * [30] Y. Han, L. Mi, L. Shen, C. Cai, Y. Liu, K. Li, G. Xu, A short-term wind speed prediction method utilizing novel hybrid deep learning algorithms to correct numerical weather forecasting, Applied Energy 312 (2022) 118777.
  * [31] L. Harris, A. T. McRae, M. Chantry, P. D. Dueben, T. N. Palmer, A generative deep learning approach to stochastic downscaling of precipitation forecasts, Journal of Advances in Modeling Earth Systems 14 (10) (2022) e2022MS003120.
  * [32] G. Chen, W.-C. Wang, Short-term precipitation prediction for contiguous united states using deep learning, Geophysical Research Letters 49 (8) (2022) e2022GL097904.
  * [33] W. Jiang, L. Zhang, Geospatial data to images: A deep-learning framework for traffic forecasting, Tsinghua Science and Technology 24 (1) (2019) 52–64.  doi:10.26599/TST.2018.9010033.
  * [34] M. Muñoz-Organero, P. Queipo-Álvarez, Deep spatiotemporal model for covid-19 forecasting, Sensors 22 (9) (2022).  doi:10.3390/s22093519.
URL https://www.mdpi.com/1424-8220/22/9/3519

  * [35] Z. Li, J. Dong, Big geospatial data and data-driven methods for urban dengue risk forecasting: A review, Remote Sensing 14 (19) (2022) 5052.
  * [36] N. Elias, Deep learning methodology for early detection and outbreak prediction of invasive species growth, in: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, 2023, pp. 6335–6343.
  * [37] N. Gorelick, M. Hancher, M. Dixon, S. Ilyushchenko, D. Thau, R. Moore, Google earth engine: Planetary-scale geospatial analysis for everyone, Remote sensing of Environment 202 (2017) 18–27.
  * [38] J. T. Abatzoglou, S. Z. Dobrowski, S. A. Parks, K. C. Hegewisch, Terraclimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015, Scientific data 5 (1) (2018) 1–12.
  * [39] A. Zeng, et al., Are transformers effective for time series forecasting?, Proceedings of the AAAI conference on artificial intelligence (2023) 11121–11128.
  * [40] T. Chen, C. Guestrin, Xgboost: A scalable tree boosting system, in: Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, 2016, pp. 785–794.
  * [41] J. Guibas, M. Mardani, Z. Li, A. Tao, A. Anandkumar, B. Catanzaro, Adaptive fourier neural operators: Efficient token mixers for transformers (2022).  arXiv:2111.13587.
  * [42] Y. Liu, L. Ren, M. Ma, X. Yang, F. Yuan, S. Jiang, An insight into the palmer drought mechanism based indices: comprehensive comparison of their strengths and limitations, Stochastic Environmental Research and Risk Assessment 30 (02 2015).  doi:10.1007/s00477-015-1042-4.
  * [43] S. Fort, H. Hu, B. Lakshminarayanan, Deep ensembles: A loss landscape perspective, arXiv preprint arXiv:1912.02757 (2019).
  * [44] C. M. Bishop, N. M. Nasrabadi, Pattern recognition and machine learning, Vol. 4, Springer, 2006.
  * [45] S. Jain, G. Liu, J. Mueller, D. Gifford, Maximizing overall diversity for improved uncertainty estimates in deep ensembles, in: Proceedings of the AAAI conference on artificial intelligence, Vol. 34, 2020, pp. 4264–4271.
  * [46] Y. Maximov, D. Reshetova, Tight risk bounds for multi-class margin classifiers, Pattern Recognition and Image Analysis 26 (2016) 673–680.
  * [47] Y. Maximov, M.-R. Amini, Z. Harchaoui, Rademacher complexity bounds for a penalized multi-class semi-supervised algorithm, Journal of Artificial Intelligence Research 61 (2018) 761–786.
  * [48] M.-R. Amini, V. Feofanov, L. Pauletto, E. Devijver, Y. Maximov, Self-training: A survey, arXiv preprint arXiv:2202.12040 (2022).

##  6 Appendix

###  6.1 Multiclass problem

##### Description of experiment

As an additional experiment and study of models’ limits, we looked at the
multiclass classification problem for Missouri state. For the 3 class problem,
we arbitrarily set up thresholds of PDSI values as $-1,\,\,1$ and for the 5
class as $-3,\,\,-1,\,\,1,\,\,3$. As an evaluation metric, we use median
accuracy over all celled predictions. For gradient boosting and logistic
regression we use default implementations for multiclass predictions, while
for neural networks-based model we replace two possible output cells with the
number of cells equal to the number of classes.

##### Analysis of results

Results, provided in Figure 6 and Table 7, are similar to a binary problem
while multi-class classification typically has worse risk estimates [46, 47,
48]. Logistic regression and gradient boosting hold their superiority longer,
Convolutional LSTM score is relatively stable, and only transformers’
prediction disappoints, falling to the level of historical baseline.

Figure 6: Accuracy vs forecast horizon, 3 and 5 class task, Missouri Horizon, months | 1 | 3 | 6 | 9 | 12
---|---|---|---|---|---
median Accuracy, 3 class: |  |  |  |  |
Baseline | 0.372 | 0.366 | 0.358 | 0.358 | 0.358
LogReg | 0.86 | 0.71 | 0.59 | 0.52 | 0.46
XGBoost | 0.86 | 0.72 | 0.61 | 0.53 | 0.49
FourCastNet | 0.451 | 0.409 | 0.509 | 0.444 | 0.431
EarthFormer | 0.37 | 0.36 | 0.37 | 0.33 | 0.36
ConvLSTM | 0.637 | 0.598 | 0.561 | 0.500 | 0.451
median Accuracy, 5 class: |  |  |  |  |
Baseline | 0.221 | 0.223 | 0.226 | 0.225 | 0.223
LogReg | 0.73 | 0.56 | 0.41 | 0.32 | 0.26
XGBoost | 0.75 | 0.58 | 0.42 | 0.32 | 0.29
FourCastNet | 0.367 | 0.305 | 0.319 | 0.266 | 0.279
EarthFormer | 0.17 | 0.18 | 0.26 | 0.19 | 0.20
ConvLSTM | 0.504 | 0.455 | 0.389 | 0.344 | 0.312
Table 7: Median Accuracy vs Forecast Horizon, 3 and 5 possible class problems, Missouri; best values are in bold, second best are underlined Region | 104x136 | 88x120 | 72x104 | 48x80
---|---|---|---|---
area | (100%) | (75%) | (53%) | (27%)
104x136 (100%) | 0.9292 | 0.9302 | 0.9327 | 0.9356
88x120 (75%) | - | 0.9466 | 0.9487 | 0.9515
72x104 (53%) | - | - | 0.9432 | 0.9438
48x80 (27%) | - | - | - | 0.9326
Table 8: ROC AUC score trained on subset (left) and evaluated on different region (top) (length of input history = 6, forecast horizon=1) Region | 104x136 | 88x120 | 72x104 | 48x80
---|---|---|---|---
area | (100%) | (75%) | (53%) | (27%)
104x136 (100%) | 0.8746 | 0.8775 | 0.8816 | 0.8829
88x120 (75%) | - | 0.8525 | 0.8536 | 0.8508
72x104 (53%) | - | - | 0.8595 | 0.8655
48x80 (27%) | - | - | - | 0.8710
Table 9: ROC AUC score trained on subset (left) and evaluated on different region (top) (length of input history = 12, forecast horizon=3) Region | 104x136 | 88x120 | 72x104 | 48x80
---|---|---|---|---
area | (100%) | (75%) | (53%) | (27%)
104x136 (100%) | 0.7525 | 0.7624 | 0.7722 | 0.7897
88x120 (75%) | - | 0.6901 | 0.6974 | 0.7081
72x104 (53%) | - | - | 0.7597 | 0.7794
48x80 (27%) | - | - | - | 0.7544
Table 10: ROC AUC score trained on subset (left) and evaluated on different region (top) (length of input history = 9, forecast horizon=6) Region | 104x136 | 88x120 | 72x104 | 48x80
---|---|---|---|---
area | (100%) | (75%) | (53%) | (27%)
104x136 (100%) | 0.6117 | 0.6154 | 0.6142 | 0.6210
88x120 (75%) | - | 0.5890 | 0.5980 | 0.6033
72x104 (53%) | - | - | 0.6126 | 0.6184
48x80 (27%) | - | - | - | 0.6264
Table 11: ROC AUC score trained on subset (left) and evaluated on different
region (top) (length of input history = 24, forecast horizon=12)