# Coexistence of Multifold and Multidimensional Topological Phonons in KMgBO3

P. C. Sreeparvathy  Department of Physics, Indian Institute of Technology
Bombay, Powai, Mumbai 400076, India  Chiranjit Mondal  Department of Physics,
Indian Institute of Technology Bombay, Powai, Mumbai 400076, India  Center for
Correlated Electron Systems, Institute for Basic Science (IBS), Seoul 08826,
Korea  Chanchal K. Barman  Department of Physics, Indian Institute of
Technology Bombay, Powai, Mumbai 400076, India  Department of Physics,
Sungkyunkwan University, Suwon 16419, Republic of Korea  Aftab Alam
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Bombay,
Powai, Mumbai 400076, India

###### Abstract

Topological interpretations of phonons facilitate a new platform for novel
concepts in phonon physics. Though there are ubiquitous set of reports on
topological electronic excitations, the same for phonons are extremely
limited. Here, we propose a new candidate material, KMgBO3, which showcase the
co-existence of several multifold and multidimensional topological phonon
excitations, which are protected by spatial and non-spatial symmetries. This
includes zero dimensional double, triple and quadratic Weyl phonon nodes, one
dimensional nodal line/loop and two dimensional doubly degenerate nodal
surface states. Nodal line/loop emerges from the spin-$\frac{1}{2}$ phonon
nodes, while the two dimensional doubly degenerate nodal surface arises from a
combination of two fold screw rotational and time reversal symmetries.
Application of strain breaks the C3 rotational symmetry, which annihilates the
spin-1 double Weyl nodes, but preserves other topological features.
Interestingly, strain helps to create two extra single Weyl nodes, which in
turn preserve the total chirality. Alloying also breaks certain symmetries,
destroying most of the topological phonon features in the present case. Thus,
KMgBO3 is a promising candidate which hosts various Weyl points, large Fermi
arcs with a very clean phonon spectra and tunable topological phonon
excitations, and hence certainly worth for future theoretical/experimental
investigation of topological phononics.

Introduction: Topological properties of fermionic states enforced by spatial
and non-spatial symmetries have been ubiquitously explored in the
literature.Kaustuv ; YXia ; ShuoWang ; HaoZheng ; MasatoshiSato  The same for
bosons (photons, phonons), however, has started only recentlySebastian ; Roman
; Emil  and there are only limited studies along these lines. The concept of
chiral phonons is further extended to materials with non-symmorphic
crystalline symmetry, which unlocks the applications in THz frequency
range.Tiantian ; Jiangxu ; JiangxuLi  In addition, topological phonons also
find applications in unconventional heat transfer, electron phonon coupling,
and phonon diode.Jiangxu ; Emil ; Yizhou  The advantage of topological bosons
over fermions is the possibility of probing topological states in any
frequency range using experimental techniques unlike fermionic states, where
this is mostly possible only near the Fermi level (EF). In line with
electronic states, topological phonons also show anomalous transport
behaviour, such as anomalous phonon hall effect.Lifa ; KangtaiSun ; TaoQin
Over the years, diverse quasiparticle excitations have emerged as a
consequence of different band crossings which can host non-zero topological
Chern number. For instance, conventional spin-$\frac{1}{2}$ Weyl fermions
possess Chern number $\pm$1\. Moreover, unconventional chiral particles with
Chern number greater than unity have also been realised in crystalline
systems. Threefold spin-1 and fourfold charge-2 Dirac fermions are examples of
such higher Chern number quasiparticles.Zhicheng ; ShiXinZhang ; chanchal
These multifold band crossings yield a zero dimensional (0D) point degeneracy
in momentum space.MZahidHasan  There exists other types of topological band
crossings such as nodal line/loop and nodal surface which induce one
dimensional(1D) and two dimensional (2D) degeneracy respectively.XiaotianWang

A proper understanding and hunt for real materials which host multiple
topological phonons have emerged as a frontline area of research due to its
demand from both fundamental and application point of view. Several
3-dimensional and 2-dimensional materials are investigated for topological
phonons.RYWang ; ChengwuXie ; YuanjunJin  Coexistence of a few
multidimensional topological phonons have been identified in some
materials.Jianhua  However, to the best of our knowledge, coexistence of
multifold and multidimensional (such as nodal points(NP), nodal line(NL) and
nodal surfaces(NS)) phononic quasiparticles is never reported in a real
material.

Tuning the topological properties of a material via strain/pressure/doping is
a well accepted method which helps to modify the symmetry and hence the
electronic structure, facilitating the onset of trivial/nontrivial topological
excitations.CJEklund ; LeiJin  There has been several theoretical/experimental
studies which has showcased such symmetry mediated transition from topological
insulator to semimetal and vice versa, creation and annihilation of Weyl/Dirac
and other kinds of topological nodes.ChunLin ; SteveM ; YanSun  A similar band
engineering and the associated transition can be expected for phonon
quasiparticles as well.

In this letter, we present the co-existence of several topological phonon
excitations in a single compound, KMgBO3, which crystallizes in a cubic phase.
KMgBO3 is an experimentally synthesised chiral compound, which has been
explored for nonlinear optical properties.exp_paper ; JianghuiZheng  Recently,
this compound is also found to show interesting topological electronic
properties by us.chanchal  One of the main focus of the present study is to
investigate the topological phonon excitations and the effect of symmetry
breaking on topological properties of this material. Pristine KMgBO3 has a
rich phonon band dispersion, showing several topologically non-trivial
multifold/multidimensional band crossings in different frequency range. This
includes spin-1 Weyl, charge 2 Dirac, single Weyl, nodal line, nodal loop and
nodal surface states. These states are topologically protected by crystalline
and time reversal symmetries. Application of strain breaks the C3 rotational
symmetry, which annihilates the spin-1 double Weyl nodes but preserve the
charge 2 Dirac node and nodal surface states. We have further disturbed the
symmetry by doping with Rb, transforming it to a triclinic structure. Though
doping preserves the total chirality, it destroys all the topological phonon
features. Additionally, we have explored two other prototype compounds
RbMgBO3RVKurbatov  and CsCdBO3,HongweiYu  which show similar topological
features.

Figure 1: (Color online) For KMgBO3, (a) Crystal structure (space group P213(#
198)) (b) Brillouin zone(BZ) for bulk and (001) surface (c) Phonon dispersion
in a small frequency interval where a 3-fold spin -1 Weyl and 4-fold charge 2
Dirac point is observed. (d,f) Surface state around $\bar{\Gamma}$ and
$\bar{M}$ point on (001) surface, corresponding to the Spin-1 Weyl node and
charge 2 Dirac node respectively. (e,g) Surface arc at 18.76 THz and 18.73 THz
between $\bar{\Gamma}$ and $\bar{M}$ points.

Computational Details: First-principles calculations are done using the Vienna
Ab initio Simulation Package (VASP)Kresse1 ; Kresse2  and Phonopy.Togo  Other
computational details are provided in the supplementary material(SM).suppl

Crystal Structure : KMgBO3 crystallizes in cubic structure with space group
P213 (#198).exp_paper  The theoretically optimised lattice parameter (6.09
$\AA$) matches fairly well with the experiment (6.8345 $\AA$).exp_paper
Figure 1(a,b) show the real space crystal structure and bulk and (001) surface
Brillouin zone (BZ). The space group 198 holds the tetrahedron (T4) point
group symmetry, which provides two twofold screw rotations, S2z =
{C${}_{2z}|\frac{1}{2},0,\frac{1}{2}|$}, S2y
={C${}_{2z}|0,\frac{1}{2},\frac{1}{2}|$} and one threefold rotation, S3=
{C+3,111$|$0,0,0} as generators at $\Gamma$ point. At R point, the generators
are S2x = {C${}_{2x}|\frac{1}{2},\frac{3}{2},0|$}, S2y =
{C${}_{2y}|0,\frac{3}{2},\frac{1}{2}|$} and S3= {C-3,111$|$0,0,0}, while at X
point, they are S2y = {C${}_{2y}|0,\frac{1}{2},\frac{1}{2}|$} and S2z =
{C${}_{2z}|\frac{1}{2},0,\frac{1}{2}|$}.book  In addition, the space group
preserves time-reversal symmetry($T$).

Phonon spectra: KMgBO3 structure involves four formula units per unit cell
leading to 24 atoms in the cell, which generate 3 acoustic and 69 optical
phonon modes. The absence of imaginary modes confirms the dynamical stability
of the compound. Phonon dispersion over a wider frequency range (0-40 THz) is
presented in Fig. S1(a).suppl  A close inspection of bulk phonon dispersion
reveals the prosperity of different nature of band crossings in different
frequency range, yielding topological Weyl, nodal line/loop and nodal surface
states, as detailed below. We have also studied the surface phonon spectra.

Figure 2: (Color online) For KMgBO3, the phonon bands along (a) $\Gamma$-X and
(b) $\Gamma$-R, showing type-I and type II Weyl crossings respectively (c)
phonon surface spectra on (001) surface corresponding to type-I Weyl node (d)
phonon dispersion along R-$\Gamma$-X confirming a nodal loop behavior (e)
phonon dispersion along $\Gamma$-R confirming a nodal line behavior (f)
drumhead like surface state on (001) surface corresponding to the nodal loop
shown above (g) surface arc due to the nodal loop on kx-ky plane (h) phonon
dispersion along X-M, showing the doubly degenerate bands (i) surface spectra
due to the doubly degenerate bands.

Double Weyl Phonons : A symmetry enforced 3- or 4-fold band crossings can
induce unconventional topological excitations.BBradlyn  Let us first consider
the case at $\Gamma$ point. For phonons, the square of the time-reversal
symmetry operator($T$) is identity($I$). The generators S2z and S2y commute
each other at $\Gamma$ point, which can be understood by considering the
transformation of lattice coordinates as a function of the twofold screws C2z
and C2y. These transformations are shown in Ref. footnote0, .

To understand the degeneracy of bands at $\Gamma$, one need to consider the
eigenstates generated by these symmetry operators which commute with the
Hamiltonian of the system, S2z$|\psi>$ = $\lambda$1$|\psi>$, S2y$|\psi>$ =
$\lambda$2$|\psi>$, where $\lambda$1 =$\pm$1, $\lambda$2 =$\pm$1 according to
S22z = 1 and S22y = 1. S3 obeys S2zS3=S3S2y and S3S2zS2y = S2yS3, which impose
a threefold degeneracy at $\Gamma$ with eigenstates $|\psi>$, $|$S3$\psi$$>$
and $|$S23$\psi$$>$, when screws are nontrivial.chanchal  When screws are
trivial, the $\Gamma$ point can hold either nondgenerate or twofold
degeneracy.chanchal  The threefold degenerate bands at $\Gamma$ arises from a
combination of two highly linearised bands and one flat band, which can be
defined by spin-1 states with chirality +2,0 and -2.BBradlyn ; chanchal  Two
screws S2x and S2y anticommute each other at R point, and S22x = -1, S22y =
-1. This induces a twofold degenerate state (with eigenvalue $\pm$i) at R
point. Further, the combination of threefold rotation S3 and twofold screws
create two more distinct states, enabling a fourfold charge 2 Dirac point with
chirality $-2$,$-2$, $+2$ and $+2$. Further symmetry related details can be
found in our previous study performed for electronic excitation.chanchal

Figure 1(c) displays the phonon dispersion in a smaller frequency range where
a threefold spin-1 Weyl point (at $\Gamma$) and a fourfold charge 2 Dirac
point (at R) are observed. At $\Gamma$ point, near 18.76 THz, three bands
(band number 54, 55 and 56) cross degenerately, with Chern number -2 (for band
54), 0 (for 55) and +2 (for 56). At R point, near 18.73 THz, a fourfold
degenerate band crossing arises which involve 53, 54, 55 and 56th bands (Fig.
1(c)). This fourfold unconventional Dirac like node can be seen as a
combination of two identical spin- $\frac{1}{2}$ Weyl phonons, with a net
topological charge 2.BBradlyn  A similar spin-1 Weyl and charge $2$ Dirac node
also occur at other frequency range, one of which is shown in Fig. S1(b) of
SM.suppl  Simulated Berry curvature at $\Gamma$ point (for -2 chirality) in
the kx-ky plane is shown in Fig. S1(c),suppl  clearly confirming the flow of
Berry flux from $\Gamma$ point. Due to the large separation between bulk
double Weyl phonons, a clean surface spectra is expected as analysed below.

Surface states of the double Weyl phonons are examined on (001) and (111)
surfaces. Figure 1(d) and 1(f) show the (001) surface state corresponding to
spin-1 Weyl point and charge $2$ Dirac point projected on $\bar{\Gamma}$ and
$\bar{M}$ points respectively. Two linearly dispersed and one parabolic
surface state can be seen at $\bar{\Gamma}$ point, which are originated from
three bulk bands with chirality $-2$ ,$0$,$2$, while a highly linearised
surface state is observed at $\bar{M}$ point, which represent the charge 2
Dirac point. Since these two nodes at $\bar{\Gamma}$ and $\bar{M}$ possess
opposite chirality, one can expect a surface arc connecting these two points.
A surface arc simulated at 18.76 THz is shown in Fig. 1(e), which indirectly
connects $\bar{\Gamma}$ and $\bar{M}$ points. Figure 1(g) shows the surface
arc projected at a slightly lower frequency (18.73 THz) which directly
connects the two points, with a long diagonal. Such large surface arcs are
much easier to probe experimentally, as reported earlier for other topological
phonon materials.Tiantian  Further, the analysis of (111) surface shows the
presence of charge 2 Dirac node, see Fig. S1(d) of SM.suppl  The surface state
of spin-1 Weyl point located at higher phonon frequency is also presented in
SMsuppl  (see Fig. S1(e)).

Single Weyl Phonons : Pristine KMgBO3 also hosts two types of spin
$\frac{1}{2}$ topological states, (1) zero dimensional Weyl points and (2) one
dimensional Nodal line/Nodal ring. First let us discuss the single Weyl point.
Weyl nodes can be formed by breaking either the time-reversal symmetry or the
inversion symmetry. In the present case, it arises from the noncentrosymmetric
nature of KMgBO3. Weyl nodes can be further classified into type I and II,
according to the nature of band crossings.footnote1 ; Alexey ; BWXia  In
KMgBO3 phonon spectra, both types of Weyl nodes are observed, as shown in Fig.
2(a,b). Type I is observed on ’$xy$’ plane in the frequency range 18.71-18.74
THz (formed by the crossing of 54th and 55th phonon bands). While, type II is
observed along $\Gamma$-R line in the frequency range 6.76-6.84 THz (formed by
the crossing of 28th and 29th bands). We have chosen one of the Weyl points
(type-I) to present its surface state on (001) surface, as shown in Fig. 2(c).
Clearly, this is a highly linearised surface state with a pristine surface
arc.

Next, we will discuss two other types of spin-$\frac{1}{2}$ topological nodal
points. Presence of more than one nodal points at the same frequency but
different momentum forms one dimensional nodal line/loop like band crossing.
Figure 2(d) shows one such example where two tilted nodal points are observed
along two different high symmetry lines ($\Gamma$-X and $\Gamma$-R) at the
same frequency. This causes the formation of a nodal loop. Similarly, there
exists nodal loops along other directions such as $\Gamma$-Y and $\Gamma$-R,
$\Gamma$-Z and $\Gamma$-R. Additionally, the phonon dispersion along
$\Gamma$-R show two nodal points at same frequency which provides a nodal line
behaviour (see Fig. 2(e)). Since the bulk R point projects as $\bar{M}$ on
(001) surface, it gives a drumhead-like surface due to the nodal loop around
$\bar{\Gamma}$ point. This is shown in Fig. 2(f) and the corresponding iso-
surface arc in Fig. 2(g). A similar nature of iso-surface has been observed in
a carbon allotrope.JingYangYou

Nodal Surface : Higher dimensional band degeneracy has been theoretically
predictedQingBo  and experimentally realised in recent past.Yihao  There are
few symmetry mediated conditions which impose such high dimensional band
degeneracy.footnote2  The bands along X-M and X-R aligned in kx=$\pi$ plane
provide a minimum of double degeneracy due to the above condition.ChengwuXie ;
QingBo  Figure 2(h) shows the double degenerate band along X-M. Similarly, a
doubly degenerate phonon band can be seen in ky,z=$\pi$ planes. It is evident
from the figure that, these doubly degenerate bands are well separated and
hence can be easily realised experimentally. Figure 2(i) shows the
corresponding surface states, which is mostly flat with a low spread in
frequency range enabling KMbBO3 to be a promising candidate for ideal nodal
surface. In a similar way, one can observe a nodal surface on other two
surfaces. A combined visual of these nodal surface create and nodal cube
around the BZ centre. This is a unique surface feature of KMgBO3 which stands
out in comparison to other reported topological nodal surface
materials.Xiaotian

Figure 3: (Color online) For KMgBO3 under 3$\%$ compressive strain: phonon
bands (a) around $\Gamma$ point and (b) R point. (c) Surface state
corresponding to charge 2 Dirac node at R point. (d) Berry curvature at two
Weyl points on kx-kz plane. (e) Phonon dispersion in strained state (space
group P212121), representing the charge 2 Dirac point at R (with chirality -2)
and two Weyl nodes (with chirality +1) along -Z - $\Gamma$-Z. (f) Phonon
dispersion for Rb-substituted KMgBO3 (space group P1), representing the
annihilation of charge 2 Dirac point at R and the Weyl nodes along -Z -
$\Gamma$-Z direction.

Strain effects : Both uni-axial (along $c$ axis) and bi-axial (along $a$ and
$b$ axes) compressive and tensile strains are applied on KMgBO3 to illustrate
the effect of symmetry breaking on topological phonon excitations. Such strain
breaks the cubic symmetry and transforms the structure to an orthorhombic
space group P212121 (# 19), which lacks the threefold C3 rotational
symmetrybook  along with deltahedron type (D42) point group (see
footnotefootnote3  for more symmetry related details).

Figure 4: (Color online) Phonon dispersion for (a) RbMgBO3 and (b) CsCdBO3.

Let us now understand the emergence of different topological phonon
excitations as a function of uni-axial strain. Here we applied a maximum of
3$\%$ strain (both in compressive and tensile forms) in an interval of 1$\%$
along $c$ axis. At ambient state (cubic phase), KMgBO3 hosts two kinds of
double Weyl phonons i.e. spin-1 type at $\Gamma$ and charge 2 Dirac type at R
points. It is clear that, in strained system, the generators at $\Gamma$ will
enforce the bands to become non-degenerate (since the two screw operators
commute), which should destroy the spin-1 type double Weyl phonon nature. The
eigenvalues can be $\pm$1 at $\Gamma$ point.RMatthias  This is precisely shown
in Fig. 3(a) (under 3$\%$ compressive strain). Coming to R point, the possible
combinations of eigenvalues are ($1$,-$1$,-$1$), (-$1$,$1$,-$1$),
(-$1$,-$1$,$1$) and ($1$,$1$,$1$)RMatthias  which provide a fourfold
degeneracy and assure the possibility of preserving the charge 2 Dirac point.
The bulk phonon dispersion at R point is shown in Fig. 3(b)(under 3$\%$
compressive strain), confirming the persistence of charge 2 Dirac state at R.
To confirm the double Weyl phonon nature, the surface state of charge 2 Dirac
point is simulated, which shown in Fig. 3(c). A highly linearised surface
state along with a clear surface arc is clearly visible from this figure.

Further, let us analyse the effect of strain on spin-$\frac{1}{2}$ Weyl
phonons. The analysis of band topology between 54th and 55th bands shows the
presence of two Weyl points with +1 chirality around $\Gamma$-point (with
locations (0,0,+kz) and (0,0,-kz)). This balances the chirality generated due
the double Weyl phonons at R point. Figure 3(d) displays the computed Berry
curvature on kx-kz plane, which shows the outcome as expected. Importantly,
the combination of double and single Weyl phonons conserves the Ninomiya
theoremHBNielsen  in strained state. To further understand this point, we have
disturbed the symmetry by substituting one of the K element by Rb. This
transforms the system to a low symmetry triclinic space group (P1), which
lacks all the screw rotation symmetries. The phonon dispersion of strained
(P212121) and the doped phase (P1) are presented in Fig. 3(e) and 3(f)
respectively. The Dirac nature at R point is defaced in the P1 phase, which
arises purely due to the breaking of screw symmetry, with the simultaneous
defacing of +1 chiral Weyl points along -Z- $\Gamma$-Z direction.

Further, the ‘1D’ nodal ring/loop is destroyed under uni-axial strain, with
the opening of a gap along $\Gamma$-R. Regarding the ‘2D’ nodal surface, the
presence of screws and time-reversal symmetry induce a Kramers like doubly
degeneracy at ki=$\pi$ which cause the formation of nodal surface. The
application of biaxial strain alters the topological phonon excitations in a
similar way. The phonon dispersion around the $\Gamma$ and R point under 3$\%$
bi-axial strain are shown in SM.suppl  Changes in other topological features
remain in line with those of uni-axial strain.

Other materials : We have also analysed the topological phonon excitations of
two other experimentally synthesized prototype compounds RbMgBO3 and CsCdBO3.
Optimised lattice parameter for RbMgBO3 and CsCdBO3 are 7.02 $\AA$ and 7.629
$\AA$ respectively, which are in good agreement with experiment.RVKurbatov ;
HongweiYu  A similar nature of phonon dispersion is observed for these two
compounds, see Fig. 4(a,b). They show various topological phonons, including
double Weyl, single Weyl, nodal line/loop and nodal surface. Full phonon
dispersion in larger frequency interval and bulk/surface states of spin-1 type
double Weyl nodes are shown in SM.suppl

Summary : We propose three promising candidate materials (KMgBO3, RbMgBO3, and
CsCdBO3) for topological phononic applications in different THz frequency
range. At ambient condition, these materials provide a rich platform to
simultaneously host several multifold and multidimensional topological Weyl
phonons, such as spin-1 Weyl, charge 2 Dirac, single Weyl, nodal line/loop and
nodal surface. These topological states are protected by crystalline and time
reversal symmetries. Application of strain breaks the C3 rotational symmetry
which mediates the annihilation of spin-1 double Weyl nodes while preserves
the other topological features. Substitution at K/Rb/Cs-site also lowers the
symmetry to triclinic structure, and destroys most of the topological
features. The present work, for the first time, showcase a novel candidate
material where multiple topological phonon excitations coexist and hence can
be an interesting platform for future theoretical/experimental investigations.

Acknowledgements: AA acknowledges DST SERB, India (Grant No. CRG/2019/002050)
for funding to support this research. SPC thank IIT Bombay for institute
postdoctoral fellowship and computing facility.C.M. would like to thank
Institute for Basic Science in Korea (Grant No. IBS-R009-D1) for funding.

## References

  * (1) Kaustuv Manna, Yan Sun, Lukas Muechler, Jürgen Kübler and Claudia Felser, Nature Reviews Materials , 3, 244-256 (2018).
  * (2) Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava and M. Z. Hasan, Nature Physics, 5, 98-402 (2009).
  * (3) Shuo Wang, Ben-Chuan Lin, An-Qi Wang, Da-Peng Yu and Zhi-Min Liao, Advances in Physics: X, 2, 518-544 (2017).
  * (4) Hao Zheng and M. Zahid Hasan, Physics: X, 3, 466661 (2018).
  * (5) Masatoshi Sato and Yoichi Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * (6) Sebastian D. Huber, Nature Physics, 12, 621-623 (2016).
  * (7) Roman Süsstrunk and Sebastian D. Huber, PNAS, 113, E4767-E4775 (2016).
  * (8) Emil Prodan, Kyle Dobiszewski, Alokik Kanwal, John Palmieri, Camelia Prodan, Nature Communications 8, 14587 (2017).
  * (9) Tiantian Zhang, Zhida Song, A. Alexandradinata, Hongming Weng,Chen Fang, Ling Lu, and Zhong Fang, Phys. Rev. Lett. 120, 016401 (2018).
  * (10) Jiangxu Li, Jiaxi Liu, Stanley A. Baronett, Mingfeng Liu, Lei Wang, Ronghan Li, Yun Chen, Dianzhong Li, Qiang Zhu and Xing-Qiu Chen, Nature Communications, 12, 1204 (2021).
  * (11) Jiangxu Li, Qing Xie, Sami Ullah, Ronghan Li, Hui Ma, Dianzhong Li, Yiyi Li, and Xing-Qiu Chen, Phys. Rev. B 97, 054305 (2018).
  * (12) Yizhou Liu, Yong Xu, Shou-Cheng Zhang, and Wenhui Duan, Phys. Rev. B 96, 064106 (2017).
  * (13) Lifa Zhang, Jie Ren, Jian-Sheng Wang, and Baowen Li Phys. Rev. Lett. 105, 225901 (2010).
  * (14) Kangtai Sun, Zhibin Gao and Jian-Sheng Wang, Phys. Rev. B 103, 214301 (2021).
  * (15) Tao Qin, Jianhui Zhou,and Junren Shi, Phys. Rev. B 86, 104305 (2012).
  * (16) Zhicheng Rao, Hang Li, Tiantian Zhang, Shangjie Tian, Chenghe Li, Binbin Fu, Cenyao Tang, Le Wang, Zhilin Li, Wenhui Fan, Jiajun Li, Yaobo Huang, Zhehong Liu, Youwen Long, Chen Fang, Hongming Weng, Youguo Shi, Hechang Lei, Yujie Sun, Tian Qian and Hong Ding, Nature, 567, 496-499 (2019).
  * (17) Shi-Xin Zhang, Shao-Kai Jian, and Hong Yao, Phys. Rev. B 96, 241111(R) (2017).
  * (18) Chanchal K. Barman, Chiranjit Mondal, Sumiran Pujari, Biswarup Pathak, and Aftab Alam, Phys. Rev. B 102, 155147 (2020).
  * (19) M. Zahid Hasan, Guoqing Chang, Ilya Belopolski, Guang Bian, Su-Yang Xu, Jia-Xin Yin, Nature Reviews Materials, 6, 784-803 (2021)
  * (20) Xiaotian Wang, a Zhenxiang Cheng, Gang Zhang, Biao Wang, Xiao-Lin Wang and Hong Chen, Nanoscale, 12, 8314 (2020)
  * (21) R. Y. Wang, Z. J. Chen, Z. Q. Huang, B. W. Xia, and H. Xu, Phys. Rev. Materials 5, 084202 (2021).
  * (22) Chengwu Xie, Ying Liu, Zeying Zhang, Feng Zhou, Tie Yang, Minquan Kuang, Xiaotian Wang, and Gang Zhang, Phys. Rev. B 104, 045148 (2021).
  * (23) Yuanjun Jin, Rui Wang and Hu Xu, Nano Lett. 18, 7755-7760 (2018).
  * (24) Jianhua Wang, Hongkuan Yuan,Minquan Kuang , Tie Yang, Zhi-Ming Yu , Zeying Zhang, and Xiaotian Wang, Phys. Rev. B . 104, L041107 (2021).
  * (25) C.-J. Eklund, C. J. Fennie, and K. M. Rabe, Phys. Rev. B 79, 220101 (R) (2009).
  * (26) Lei Jin, Xiaoming Zhang, Ying Liu, Xuefang Dai, Liying Wang, and Guodong Liu, Phys. Rev. B 102, 195104 (2020).
  * (27) Chun Lin, Masayuki Ochi, Ryo Noguchi, Kenta Kuroda, Masahito Sakoda, Atsushi Nomura, Masakatsu Tsubota, Peng Zhang, Cedric Bareille, Kifu Kurokawa, Yosuke Arai, Kaishu Kawaguchi, Hiroaki Tanaka, Koichiro Yaji, Ayumi Harasawa, Makoto Hashimoto, Donghui Lu, Shik Shin, Ryotaro Arita, Satoshi Tanda and Takeshi Kondo, Nature Materials, 20, 1093-1099 (2021).
  * (28) Steve M. Young, Sugata Chowdhury, Eric J. Walter, Eugene J. Mele, Charles L. Kane, and Andrew M. Rappe, Phys. Rev. B 84, 085106 (2011).
  * (29) Yan Sun, Shu-Chun Wu, Mazhar N. Ali, Claudia Felser, and Binghai Yan, Phys. Rev. B 92, 161107(R) (2015).
  * (30) L.Wu,J.C.Sun,Y.Zhang,S.F.Jin,Y.F.Kong,andJ.J.Xu, Inorg. Chem. 49, 2715 (2010).
  * (31) JianghuiZheng, Qijin Cheng, Shunqing Wu, Yixi Zhuang, Ziquan Guo, Yijun Lu, Chao Chen, Materials Chemistry and Physics ,165 168-176 (2015).
  * (32) R.V.Kurbatov, Leonid A. Solovyov, B.G. Bazarov, A.K. Subanakov, J.G. Bazarova, Solid State Communications, 172 33-36 (2013).
  * (33) Hongwei Yu, Hongping Wu, Shilie Pan, Bingbing Zhang, Ming Wen, Zhihua Yang, Hongyi Li, and Xiangzhan Jiang, Eur. J. Inorg. Chem. 2013, 5528–5533 (2013).
  * (34) G. Kresse and J. Hafner, Phys. Rev. B 47, 558(R) (1993).
  * (35) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * (36) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * (37) Supplemental Material at URL[]. It contains computational details, explanation about high symmetry points in BZ, few essential figures which supports the main text for KMgBO3 and other prototype compounds RbMbgBO3 and CsCdBO3.
  * (38) P. Cracknell and Christopher J. Bradley, The Mathematical Theory of Symmetry in Solids: Representation Theory for Point Groups and Space Groups, (Oxford: Clarendon Press, 1972).
  * (39) B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser,R. J. Cava, B. Andrei Bernevig, Science 353, aaf5037 (2016).
  * (40) The transformations are C2z($x$,$y$,$z$) $\rightarrow$ ($-x$,$-y$,$z$) and C2y($x$,$y$,$z$) $\rightarrow$ ($-x$,$y$,$-z$), which define S22z = 1, S22y = 1 and S2zS2y = S2yS2z.book  The lattice coordinate transformation due to C3 is C3,111($x$,$y$,$z$) = ($z$,$x$,$y$), which shows S32 = 1.
  * (41) The general definition of the type I and type II Weyl nodes in a fermionic systems can be understand the following way: The Weyl point, which possess a zero like Fermi surface is defined as type I and the type II shows a tilted dispersion along the momentum space and make an overlap between the electron and hole pockets.Alexey  A phonon counterpart for both type I and type II can be understand in the similar way.BWXia
  * (42) Alexey A. Soluyanov, Dominik Gresch, Zhijun Wang, , QuanSheng Wu, Matthias Troyer, Xi Dai, B. Andrei Bernevig, Nature 527, 495-498 (2015).
  * (43) B. W. Xia, R. Wang, Z. J. Chen, Y. J. Zhao, and H. Xu, Phys. Rev. Lett. 123, 065501 (2019).
  * (44) Jing-Yang You, Xian-Lei Sheng, Gang Su, Phys. Rev. B 103, 165143 (2021).
  * (45) Qing-Bo Liu, Zhe-Qi Wang, and Hua-Hua Fu, Phys. Rev. B 104, L041405 (2021).
  * (46) Yihao Yang, Jian-ping Xia, Hong-xiang Sun, Yong Ge, Ding Jia, Shou-qi Yuan, Shengyuan A. Yang, Yidong Chong and Baile Zhang, Nature Communications, 10, 5185 (2019).
  * (47) The combination of twofold screw rotational symmetry and the time-reversal symmetry operators leads to the formation of an anti-unitary operator ($T$$C$2i)2 = $-1$, which opens up the possibility of a higher dimensional (2D) band degeneracy at ki=$\pi$ plane. As a result of ($T$$C$2i)2 = $-1$, at $k$i=$\pi$ high symmetry plane the bands act like Kramers degenerate with eigenvalue $\pm$i.
  * (48) Xiaotian Wang, Feng Zhou, Tie Yang, Minquan Kuang, Zhi-Ming Yu, Gang Zhang, Phys. Rev. B 104, 041104 (2021)
  * (49) Space group P2${}_{1}2_{1}2_{1}$ has Deltahedron type (D${}^{4}_{2}$) point group involving three twofold screw rotation symmetries. At $\Gamma$ point, the generators are S2z = {C${}_{2z}|\frac{1}{2}|,0,\frac{1}{2}|$} and S2y = {C${}_{2z}|\frac{1}{2}|,\frac{1}{2},0|$}. These two satisfy the commutation relation at this point [S2y ,S2z ] =0. At R point, the symmetries are S2x = {C${}_{2z}|0,\frac{1}{2},\frac{1}{2}|$}, S2x = {C${}_{2z}|\frac{1}{2},\frac{1}{2},0|$}, and identity operator. The two screws are anticommute at R point. At X point the symmetries are S2x = {C${}_{2x}|0,\frac{1}{2},\frac{1}{2}|$} and S2z = {C${}_{2x}|\frac{1}{2},0,\frac{1}{2}|$}
  * (50) R. Matthias Geilhufe1, Stanislav S. Borysov 1, Adrien Bouhon2 and Alexander V. Balatsky, Sci Rep, 7, 7298 (2017).
  * (51) H. B. Nielsen and M. Ninomiya, Phys. Lett. B 105, 219 (1981).


# Predicting the power grid frequency of European islands

Thorbjørn Lund Onsaker  Faculty of Science and Technology, Norwegian
University of Life Sciences, 1432 Ås, Norway  Heidi S. Nygård  Faculty of
Science and Technology, Norwegian University of Life Sciences, 1432 Ås, Norway
Damià Gomila  Instituto de Física Interdisciplinar y Sistemas Complejos, IFISC
(CSIC-UIB), Campus Universitat Illes Balears, E-07122 Palma de Mallorca, Spain
Pere Colet  Instituto de Física Interdisciplinar y Sistemas Complejos, IFISC
(CSIC-UIB), Campus Universitat Illes Balears, E-07122 Palma de Mallorca, Spain
Ralf Mikut  Institute for Automation and Applied Informatics, Karlsruhe
Institute of Technology, 76344 Eggenstein-Leopoldshafen, Germany  Richard
Jumar  Institute for Automation and Applied Informatics, Karlsruhe Institute
of Technology, 76344 Eggenstein-Leopoldshafen, Germany  Heiko Maass  Institute
for Automation and Applied Informatics, Karlsruhe Institute of Technology,
76344 Eggenstein-Leopoldshafen, Germany  Uwe Kühnapfel  Institute for
Automation and Applied Informatics, Karlsruhe Institute of Technology, 76344
Eggenstein-Leopoldshafen, Germany  Veit Hagenmeyer  Institute for Automation
and Applied Informatics, Karlsruhe Institute of Technology, 76344 Eggenstein-
Leopoldshafen, Germany  Dirk Witthaut  Forschungszentrum Jülich, Institute for
Energy and Climate Research - Systems Analysis and Technology Evaluation (IEK-
STE), 52428 Jülich, Germany  Institute for Theoretical Physics, University of
Cologne, 50937 Köln, Germany  Benjamin Schäfer  Institute for Automation and
Applied Informatics, Karlsruhe Institute of Technology, 76344 Eggenstein-
Leopoldshafen, Germany

###### Abstract

Modelling, forecasting and overall understanding of the dynamics of the power
grid and its frequency is essential for the safe operation of existing and
future power grids. Much previous research was focused on large continental
areas, while small systems, such as islands are less well-studied. These
natural island systems are ideal testing environments for microgrid proposals
and artificially islanded grid operation. In the present paper, we utilize
measurements of the power grid frequency obtained in European islands: the
Faroe Islands, Ireland, the Balearic Islands and Iceland and investigate how
their frequency can be predicted, compared to the Nordic power system, acting
as a reference. The Balearic islands are found to be particularly
deterministic and easy to predict in contrast to hard-to-predict Iceland.
Furthermore, we show that typically 2-4 weeks of data are needed to improve
prediction performance beyond simple benchmarks.

##  1 Introduction

An increasingly interdisciplinary community is contributing to the
understanding of power grid dynamics and stability, including dynamical
modelling, simulations and machine-learning-based forecasts, see e.g. [1] for
a recent review. The energy system is a very important and challenging
research area as it constitutes a complex system with many interacting
entities, e.g. its energy markets, generators, consumers, different control
areas and control mechanisms.

In this article, we investigate the dynamics and predictability of the power
grid frequency in five grids of very different sizes. The grid frequency is
critical to ensure a safe and reliable supply of electricity due to its
central role in power system control [2]. The electric frequency of the
alternating current power grids mirrors the current balance in supply and
demand: An excess of demand effectively draws from the kinetic energy stored
in rotating turbines and slows these down. Vice versa, an excess of generation
cannot be stored directly in the grid and thereby leads to an increased
frequency (unless batteries, pumped hydropower or other storage options are
set in operation to restore the power balance). Critically, the frequency has
to be controlled tightly using numerous control mechanisms [2, 3, 4], as
otherwise consumers or generators have to be disconnected. Desired frequency
ranges depend on the synchronous area, i.e. the connected power grid sharing
the same frequency, and range from $f=50\pm 0.2$ Hz [5] to $f=50\pm 0.05$ Hz
[6] for the areas considered in the present paper. With an increasing share of
volatile renewable generation, in particular wind and solar power, the total
inertia of the power system decreases, while fluctuations overall increase [7,
8]. This further amplifies the need for accurate forecasts and a solid
understanding of power systems.

Islands and their power grids are particularly interesting for two reasons:
First, there exist many natural islands which are unconnected to any other
power grid or will only be connected via HVDC (high-voltage-directed-current)
links, still ensuring a decoupled frequency dynamics [9]. Small island grids
typically have a lower inertia and display larger frequency deviations than
large integrated networks [10]. Hence they may provide valuable insights for
variable grid operation. Second, microgrid proposals [11, 12, 13], which
suggest establishing a local cell of the power grid, are receiving much
attention. Depending on the specific scenario, cells might be mostly
autonomous, have their own market [14] or only decouple during emergencies to
stabilize the grid and avoid cascading failures [15, 16].

Fully modelling the short-term frequency dynamics is very challenging due to
its combined stochastic-deterministic nature and the numerous external
influences [17, 18, 10], suggesting the use of data-driven approaches [19] and
the usage of machine learning to predict time series as well as to understand
and control energy systems [20, 21, 22, 23]. These data-driven approaches
complement modelling-based approaches [10, 24] as they do not make any
assumptions about the governing dynamical equations but still provide
forecasts and explanations of the system.

Any data-based forecasting method depends on the availability of data with a
sufficient quality covering a sufficiently long period. Unfortunately, such
data are still not easily available. Even research projects, such as the
GridEye/FNET-initiative [25], are not always open, which limits their value
for the research community. In recent work, we have pushed for more open data
of power grid frequency measurements [26, 27, 28], including coverage of power
grids on islands. Still, with only limited measurement devices available, some
questions arise: Are previously developed methods which forecast the power
grid frequency [20] applicable to islands? How much data are necessary to
develop data-driven approaches that can outperform simple benchmarks? How can
additional information, e.g. from transparency platforms [29], be integrated
into frequency forecasts?

To answer these questions, the remaining paper is structured as follows. We
start with an overview of the data obtained from the Nordic system, the Faroe
Islands, Ireland, the Balearic Islands and Iceland in Sec. 2. Next, we provide
an introduction to the nearest-neighbour predictor used for forecasting in
Sec. 3.1. Then, we investigate the performance of the predictor, compared to
benchmarks and as a function of available data in Sec. 3.2 to then discuss the
possibility of extending the predictor by including additional data in Sec.
3.3. We close with a discussion in Sec. 4.

##  2 Data overview

In the present article, we analyse the power grid frequency time series from
four different European islands in addition to the Nordic synchronous area.
These islands are the Faroe Islands, Ireland, the Balearic Islands and
Iceland, see also Fig 1, showing good geographical spread and a difference in
size. We consider time series from 6 days (Faroe Islands) up to 450 days
(Balearic Islands) for our analysis. The five regions (Nordic and the islands)
are also very different in terms of their energy mix, regulations and
interconnections:

The Nordic power system includes Norway, Sweden and Finland as well as the
Eastern part of Denmark, a total of approximately 24 million people. The
energy mix for the Nordic electricity supply consists of a significant share
of renewable energy, in particular wind and hydropower, with a steady increase
in recent years [30].

The Faroe Islands are an archipelago of 18 major islands in the Atlantic
Ocean, with approximately 50 000 inhabitants. The power generation relies
mostly on fossil fuels [31].

The Irish power system covers both the Republic of Ireland and Northern
Ireland, with about 6.8 million inhabitants total. Ireland’s generation is
tailed towards a high share of wind energy and in 2020 the renewable
electricity generation reached 43% of annual generation [32]. Ireland is
connected to Great Britain via HVDC cables [33].

The Balearic Islands, belonging to Spain, are located in the Mediterranean Sea
with a total population of about 1.2 million. Generation relies mostly on
fossil fuels such as gas (combined cycle gas turbines) and diesel engines
leading to a share in generation of more than 70% [34]. Furthermore, an HVDC
cable connecting the islands to mainland Spain is an important factor in the
electric system.

Finally, Iceland is located in the North Atlantic Ocean and has roughly 300
000 inhabitants. Its main source of electricity generation is hydropower,
making up 75% of the total generation. In addition, Iceland’s volcanic
activities allow extensive geothermal power generation. These cheap energy
sources have attracted energy-intensive industries such as ferrosilicon,
aluminium smelters and data centres [35].

Figure 1: Frequency data are taken from several European islands. We display
the synchronous areas analysed within the present article: Ireland, Balearic
Islands, Iceland and Faroe Islands as the island systems with the Nordic
synchronous area as a continental reference system. Map created via Basemap
and Python 3.

The data for the Nordic grid are based on measurements provided by the Finish
operator Fingrid [36], see also [26, 27, 28]. Meanwhile, the time series of
the islands are obtained via our measurement device, the Electrical Data
Recorder (EDR), which has been developed at KIT. The EDR records the voltage
with a sampling frequency of $25\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ and
assigns a standard time to each data point based on a GPS recorder. Then, the
frequency is extracted based on zero-crossings of the voltage signal, see [37,
38, 27, 39] for details. Finally, the frequency data are pre-processed to
identify acquisition interruptions and to remove outliers. Unfortunately, some
single devices had some longer measurement outages, leading to up to 25% of
missing data related to the installed time period in the Balearic Islands and
Ireland, respectively, but typically less than 1% for the other locations, see
also code for details [40].

To obtain an initial impression of the frequency dynamics in the different
regions, we compute the daily profiles in Fig. 2, noting stark differences.
The daily profile is the average frequency for each second of a day, averaged
over all days in the data set, defined in detail in the next section. The
daily profile in the Faroe Islands displays the largest deviations and seems
to be somewhat unpredictable. This is in part due to the short time series
available (6 days) and might further be linked to the small population size.
Meanwhile, Ireland, the Balearic Islands and the Nordic grid all display clear
patterns with the frequency jumping up and down at specific hours. These jumps
are associated with typical power dispatch actions and load and generation
ramps at full hours [41, 42]. Finally, the Icelandic daily profile is mostly
flat and does not display clear patterns, likely because the dominant
hydropower and geothermal power plants are easily controllable on the
generation side.

While the daily profile (mean daily frequency values) provides insights into
what typical frequency values to expect, the daily standard deviation reveals
the spread of frequency values over the day, see Fig. 3. Technically, instead
of computing the mean value, we compute the standard deviation for each second
of a day. Analysing the standard deviation confirms a previous observation
from the daily profile: The Balearic Islands, the Nordic and to a slightly
lesser extent the Irish systems show clean curves with regular distinct peaks
at the hourly dispatch intervals. This is consistent with earlier
observations: Larger deterministic frequency excursions at the start of an
hour are complemented with smaller but random fluctuations in between [42, 17,
10]. Meanwhile, the Faroe Islands show little structure and the largest
absolute standard deviation, consistent with the large frequency deviations
observed earlier. Iceland is special in that the standard deviation is on
average higher than in Ireland, the Balearic Islands or the Nordic grid and
occasionally, very large standard deviations are observed due to extreme
frequency deviation of the order of Hertz, instead of deviations of tens or
maybe one hundred millihertz as would be observed in the other synchronous
areas.

Complementing the daily profile and standard deviation, we analyse how much
history is preserved in the trajectories by computing the autocorrelation
decay [17]. For a generic stochastic time series, we might expect an
exponential decay of the autocorrelation [43] and we do indeed observe such a
(quick) decay in all time series, see Fig. 4. For all grids except Iceland, we
also observe regular peaks in the autocorrelation on a daily basis,
reinforcing the idea that the frequency follows daily patterns. Iceland is
special in that the frequency displays almost no correlation and seems to be
much more stochastic than the frequency in the other regions. Note that the
Faroe Islands are omitted due to their short time series.

Figure 2: The daily profile acts as the benchmark for all regions. We define
the daily profile as the average of all recordings for every second of the
day, averaged over all days in the data set, see also eq. (7). We note that
the daily profiles in the Nordic and Balearic grids are not significantly
affected if the time series is truncated to shorter periods, see also the code
[40].  Figure 3:  The daily standard deviation reflects the variability within
the frequency samples at each second of the day. Faroe Islands and Iceland
show particularly large deviations, in contrast to the deterministic profiles
in Ireland, the Balearic Islands and the Nordic grid.  Figure 4: The
autocorrelation reveals daily and other regular patterns in three out of the
four island grids. We plot the autocorrelation function with time lags up to
20 days and notice very pronounced patterns in the Balearic Islands but almost
no patterns in Iceland.

##  3 Forecasting frequency on islands

###  3.1 Introducing the WNN predictor

Figure 5: The WNN predictor predicts the next hour trajectory based on similar
patterns from the past. The past hour is compared with patterns from the
training set, for which the next hour is available, and patterns most similar
to the past hour are identified. Here, as an example, the three nearest
neighbours are shown. A weighted average of the nearest neighbours’ next hour
trajectory is computed as the prediction and can thus be compared to the
actual trajectory from the test set.

To forecast the power grid frequency, we expand upon an earlier paper [20],
which investigated numerous aspects of frequency predictability with a
weighted nearest neighbour (WNN) predictor. This earlier paper [20] focused on
the three largest synchronous areas in Europe, for which several years of data
is available. In contrast, measurements form islands, investigated in the
present article, cover few weeks or months up to at most 450 days in the case
of the Balearic islands. The WNN algorithm identifies similar patterns in the
past to forecast future dynamics, see Fig. 5 for an illustration.

Let us specify the details of the prediction process: For each synchronous
area, the time series is split into chunks of equal length which are then used
to find similar chunks (neighbors) from the past. More formally, we define the
pattern vectors

$\mathbf{F}_{n}=\begin{pmatrix}f(t_{0}-(n+1)\gamma\tau)\\\
f(t_{0}-(n+1)\gamma\tau+\tau)\\\ f(t_{0}-(n+1)\gamma\tau+2\tau)\\\ {\dots}\\\
f(t_{0}-n\gamma\tau-\tau)\end{pmatrix},$ (1)

where $t_{0}$ is the last time point for which we have data in our initial
pattern $\mathbf{F}_{0}$. Furthermore, the time delay $\tau$ is identical with
the time resolution $\tau=1$s and we choose $\gamma=3600$ data points so that
$\gamma\tau=3600$s covers one hour, a central time scale for market dynamics
and regular grid patterns [42]. Hence, each pattern $\mathbf{F}_{n}$ is a
vector with 3600 frequency measurements and the entire time series is covered
is covered by combining non-overlapping sets as

$\mathcal{F}=\\{\mathbf{F}_{n}|\exists i\in\mathbb{N}:n\gamma\tau=i\cdot
24\textrm{h}\\}$ (2)

and the index $n$ is used to distinguish between different patterns. Given the
past hour pattern as $\mathbf{F}_{0}$, the next hour of the frequency time
series is predicted by searching the set $\mathcal{F}$ for the most similar
patterns to this initial pattern. We implement this search by computing the
Euclidean distance between the initial pattern and all patterns in
$\mathcal{F}$:

$\textrm{d}(\mathbf{F}_{n})=\|\mathbf{F}_{n}-\mathbf{F}_{0}\|.$ (3)

The patterns in $\mathcal{F}$ are then sorted based on these distances in
increasing order, thus resulting in a new set of patterns. This set only
contains the $k$ most similar patterns to $\mathbf{F}_{0}$ given as

$\mathcal{N}_{k}=\\{n_{1},n_{2},\ldots,n_{k}|\mathbf{F}_{n_{j}}\in\mathcal{F}\\}.$
(4)

The final prediction $f_{wnn}(t_{0}+\Delta t)$ with time steps $\Delta t\in$
{1s, 2s, 3s, …, 3600s} is performed by a weighted average of the $k$ nearest
neighbours’ next hour trajectories:

$f_{wnn}(t_{0}+\Delta
t)=\frac{1}{\sum_{j=1}^{k}\alpha_{j}}\sum_{j=1}^{k}\alpha_{j}f(t_{0}-n_{j}\gamma\tau+\Delta
t).$ (5)

Note that by construction of $\mathcal{F}$, we restrict our search for
frequency patterns to the same hour, e.g. we predict the frequency at 9pm only
using training data from 9pm as well, see also [20, 40]. The weighting
$\alpha$ decreases with greater distance between the patterns, since close
neighbours should have a stronger influence on the prediction than more
distant ones:

$\alpha_{j}=\frac{\textrm{d}(\mathbf{F}_{n_{k}})-\textrm{d}(\mathbf{F}_{n_{j}})}{\textrm{d}(\mathbf{F}_{n_{k}})-\textrm{d}(\mathbf{F}_{n_{1}})}.$
(6)

The number of nearest neighbours $k$ is a hyperparameter of the WNN predictor,
which has to be optimised. Within this paper, we utilize the $adaptive$-$k$
approach [20], which optimises $k$ for every time step as follows: The optimal
$k$-value is given by the lowest $\textrm{MSE}_{\Delta t}(f_{wnn})$ of the
predictor for each $\Delta t$, i.e. we obtain an $adaptive$-$k$ vector with
3600 values. The $adaptive$-$k$ is smoothed with a sliding window of one
minute to avoid sudden jumps and noise for any time steps, see [40] for
details.

As a benchmark, we compare the WNN performance with the daily profile, which
effectively is an unweighted $k=\infty$ predictor (including all patterns)
defined as

$f_{dp}(t_{0}+\Delta
t)=\frac{1}{|\mathcal{F}|}\sum_{n\in\mathcal{F}}f(t_{0}-n\gamma\tau+\Delta
t).$ (7)

###  3.2 Performance of the WNN predictor

Let us now investigate whether the WNN predictor can outperform common
benchmarks on islands, as it did in continental regions [20]. We consider
three prediction models: first: a naive and constant 50Hz predictor, as this
is the typical mean value for each hour, second: the daily profile, which
takes all patterns into account equally and third: the WNN, which selects
fewer patterns and weights them to obtain a forecast. For all predictors, we
use about 70% of the data for training (i.e. data to select patterns from or
data to compute the daily profile), 15% as validation (to determine the
optimal number of nearest neighbours $k$) and 15% as test (to evaluate the
performance), see the code [40] for details.

The WNN predictor typically outperforms these benchmarks even when little data
are available, see Fig. 6. We compute the root-mean-squared-error (RMSE) for
each predictor for each minute, averaged over all hours. Consistent with the
daily profiles, we note that both benchmarks display a large error at the
start of an hour, as dispatch actions such as generation or load ramps lead to
a large deviation in the frequency. Especially at the beginning of an hour,
the WNN predictor is more precise than the benchmarks, leveraging its specific
information from the previous hour. The further away the prediction is
evaluated from the start of the hour, the less specific information is
contained in the previous patterns and the WNN approaches the daily profile
behaviour, as also seen for continental regions [20]. There are two more
interesting observations: First, even with only 6 days of training and test
data available, the WNN outperforms the daily profile and 50Hz predictor in
the Faroe Islands. Second, all predictors fare approximately equal for
Iceland, where the grid frequency time series shows almost now characteristic
patterns but large stochastic deviations.

In many situations, we will not have access to years or months of data, due to
a large variety of potential reasons: Measurement time is limited, large parts
of the data are corrupted, the specific grid is very new, either because the
grid is newly built or major changes were made, such as the expansion of a
synchronous area [44] or a change in regulations. Regardless of what limits
the access to data, we need to understand how well the WNN predicts the
frequency with shorter periods of training data available.

Interestingly, there is no general threshold of how many weeks of data are
necessary for the WNN to outperform the daily profile, see Fig. 7. There, we
systematically investigate the prediction quality as a function of the length
of available data by computing the RMSE for the constant 50Hz prediction, the
daily profile and the WNN. The latter two predictors are based on data ranging
from one week up to 6 months (where available) and always utilised an 80-20
split for training and validation, hence reporting the validation RMSE. Again,
the Faroe Islands are excluded as less than one week of data are available. We
typically note that for very little data, the constant prediction is the best,
followed by the daily profile and the WNN coming in last. If measurement times
are too short, the unspecific, average prediction is superior and the WNN
restricts itself to finding patterns in its very limited training set, leading
to a bad performance. Vice versa, if sufficient data are available, the WNN
correctly identifies characteristic patterns and is capable to make the best
predictions.

Notably, this transition from "unspecific information being best" to "specific
historic patterns being best" depends on the power grid at hand: The Balearic
islands with their very deterministic frequency profile, only require 2 weeks
for the WNN to outperform the daily profile, compared to about 4 weeks
necessary in Ireland, while in Iceland even with 2 months of data, the
constant 50 Hertz prediction is still better than daily profile or the WNN
predictor.

Figure 6: For three out of four islands, the WNN predictor outperforms the two
benchmark models. We plot the Root Mean Square Error (RMSE) of the WNN
prediction model versus the daily profile and a constant 50Hz prediction. In
all regions except Iceland, the WNN predictor outperforms the null models,
especially the first 15 minutes of the hour. The Balearic Islands have the
lowest absolute error, while the Faroe Islands have the highest. The y-axis
varies between the subplots due to the varying error between the regions.
Figure 7: With sufficient data, the WNN predictor outperforms the daily
profile at three out of four islands. We plot the average root mean square
error for the WNN predictor and the two null models for different time
intervals. As the data size increases, the WNN predictor is the best model for
all regions except Iceland, where the constant 50 Hz prediction performs best.
The Faroe Islands are not included in the analysis as the region’s frequency
time series is shorter than the minimum interval of one week.

###  3.3 Additional features

So far, we have only used historic frequency information to forecast frequency
time series. However, there are an increasing amount of energy system data
available, so we ask: How can additional information, such as the generation
of coal power plants or volatile renewable generation (wind and solar power),
be integrated into the frequency prediction, and how do they influence the
performance of the WNN predictor?

To answer this question, we utilise additional information on the Balearic
energy system, obtained from the transmission system operator of Spain, Red
Eléctrica de España, and their real-time electricity demand and generation
tracker [45]. It seems intuitive that (frequency) patterns of the power system
will differ based on the total load or the mix of generation types.
Mathematically, we extend the WNN predictor as follows: Let
$a_{\text{origin}}$ be the original additional time series. To compensate for
the large deviations in the typical values of $a_{\text{origin}}$, we carry
out a min-max scaling:

$a=\frac{a_{\text{origin}}-a_{\text{min}}}{a_{\text{max}}-a_{\text{min}}},$
(8)

where $a$ is the transformed data point and $a_{\text{max}}$ and
$a_{\text{min}}$ are the maximum and minimum values of the original time
series. Next, the additional feature is split into non-overlapping patterns
$\mathbf{A}_{n}$, similar to $\mathbf{F}_{n}$ from eq. (1):

$\mathbf{A}_{n}=\begin{pmatrix}a(t_{0}-(n+1)\eta\lambda)\\\
a(t_{0}-(n+1)\eta\lambda+\lambda)\\\ a(t_{0}-(n+1)\eta\lambda+2\lambda)\\\
{\dots}\\\ a(t_{0}-n\eta\lambda-\lambda)\end{pmatrix},$ (9)

where $\eta$ is the number of additional data points in each pattern, and
$\lambda$ is the time resolution. $\eta\lambda$ is thus the prediction window
size and must be the same size as the frequency patterns to represent the same
specific time window of 3600 seconds, matching $\mathbf{F}_{n}$. Hence,
$\eta\lambda\stackrel{{\scriptstyle!}}{{=}}\gamma\tau$. When the vectors
$\mathbf{A}_{n}$ are concatenated with $\mathbf{F}_{n}$, the result is new
vectors given as:

$\mathbf{G}_{n}=\begin{pmatrix}\mathbf{F}_{n}\\\ \mathbf{A}_{n}\end{pmatrix}.$
(10)

As the extended model is also limited to searching through patterns from the
same time of the day, $\mathcal{G}$ is a new set similar to $\mathcal{F}$ from
(2). To now choose the nearest neighbours in terms of distance, the previous
distance computation eq. (3) is extended to

$\textrm{d}(\mathbf{G}_{n})=\|\mathbf{F}_{n}-\mathbf{F}_{0}\|+\beta\cdot\|\mathbf{A}_{n}-\mathbf{A}_{0}\|,$
(11)

where $\|\cdot\|$ denotes the Euclidean distance, and $\beta$ is a tunable
weight to adjust the influence of the additional feature. With new distances,
the weighting $\alpha$ from (6) is now computed as

$\alpha_{j}=\frac{\textrm{d}(\mathbf{G}_{n_{k}})-\textrm{d}(\mathbf{G}_{n_{j}})}{\textrm{d}(\mathbf{G}_{n_{k}})-\textrm{d}(\mathbf{G}_{n_{1}})}.$
(12)

The weights are then used as in eq. (5) to forecast the frequency time series.

Additional information might improve performance slightly for certain time
periods, see Fig. 8. To obtain these results, we optimized the hyperparameter
$\beta$ via a grid search $\beta\in{0.1,0.2,...,1.5}$ Including renewable
generation accounted for the best performance improvement, closely followed by
improvements by including data from combined cycle or coal [46] generation,
all of which lead to a reduction in RMSE of more than 1%. Interestingly, the
additional features decreased the performance during the first 15 minutes,
while it improved after that. Without additional features, the predicted
trajectory has its best performance at the start of an hour and then
deteriorates since the information of the previous hour is slowly being
invalidated by additional, unforeseen events, see also [20]. The additional
feature balances this effect: In the first few minutes, the additional feature
dilutes the specific information, leading to a larger RMSE, while in the later
parts of the hour, the additional information is more useful. Further research
will be necessary to clarify the exact effects of additional features included
in the WNN predictor, especially when applying it to other islands, such as
Iceland, for which unfortunately no public data on demand and generation were
available when writing this article. Moreover, in our current implementation,
we use the actual generation for the hour for which we forecast the frequency.
Thereby, we mix past and future information. Future research could rely on
forecasts or scheduled generation, which are not always publicly available but
might be available to system operators wishing to implement the WNN predictor.

Figure 8: With additional features included in the WNN predictor for the
Balearic Islands, the prediction performance is increased. Presented are the
three best performing features, all improving the performance by an average of
1-1.5%. Across the full hour, the first 15 min got worse, but the last 45 min
significantly improved.

##  4 Discussion and Conclusion

The power grid frequency in islands may be predicted via a weighted-nearest
neighbour (WNN) predictor [20]. Compared to the continental areas considered
[20] and the Nordic area analysed in the current article, the four
investigated islands displayed more extreme properties: The daily profile in
the Balearic islands is the most pronounced and regular behaviour we have
observed so far, while the frequency dynamics in Iceland shows no
characteristic patters but strong stochastic deviations, making it
particularly tricky to predict. We also note some dependency on the population
size and the composition of the power generation and consumption: the Faroe
Islands and Iceland are the hardest to predict. These two islands are the
least populated areas and have very special generation and demand patterns,
influenced e.g. by large industry in Iceland. Meanwhile, the very good
predictability of the Balearic Islands might at least be partially explained
by its strong interconnection with the Continental European synchronous area
and the predominance of combined cycle production with its fast control
capabilities [47]. Overall, we have shown that the WNN predictor outperforms
the daily profile even when only limited data are available in three out of
the four island grids. Furthermore, we have demonstrated how additional time
series, such as generation, could be included as additional information when
forecasting the frequency.

The WNN predictor discussed in the present article could be a useful tool for
the planning and operation of islanded grids, let these be natural islands or
islanded microgrids [12]. Notably, almost no prior knowledge about the system
at hand is necessary nor do we require many assumptions to apply the WNN. Even
the difference between weekends and weekdays is implicitly included in the
WNN: Given sufficient data, it will select suitable patterns. Furthermore, the
WNN could be used as a diagnostic tool, e.g. to analyse systematic offsets in
the frequency or external perturbations and thereby detect anomalies.

Further research is necessary to fully characterise to which extent power
grids are predictable and which external factors should be included. Also,
further ML methods, beyond nearest-neighbour techniques could be considered
[48], e.g. to handle any concept drift in long-term forecasting or to provide
further monitoring capabilities.

###### Acknowledgements.

This project has received funding from the Helmholtz Association under grant
no. VH-NG-1727, the Program “Energy System Design” and the Helmholtz
Association’s Initiative and Networking Fund through Helmholtz AI.

## References

  * Witthaut et al. [2022] D. Witthaut, F. Hellmann, J. Kurths, S. Kettemann, H. Meyer-Ortmanns, and M. Timme, Reviews of Modern Physics 94, 015005 (2022).
  * Machowski et al. [2020] J. Machowski, Z. Lubosny, J. W. Bialek, and j. R. Bumby, _Power system dynamics stability and control_ (John Wiley & Sons Ltd, 2020).
  * Modig et al. [2022] N. Modig, R. Eriksson, P. Ruokolainen, J. N. Ødegård, S. Weizenegger, and T. D. Fechtenburg, _Overview of Frequency Control in the Nordic Power System_ , https://www.epressi.com/media/userfiles/107305/1648196866/overview-of-frequency-control-in-the-nordic-power-system-1.pdf (2022), [Online; accessed 25-April-2022].
  * Tosatto et al. [2020] A. Tosatto, M. Dijokas, T. Weckesser, S. Chatzivasileiadis, and R. Eriksson, _Sharing Reserves through HVDC: Potential Cost Savings in the Nordic Countries_ (2020).
  * Landsnet [2020] Landsnet, _Tíðnistýring_ , https://www.landsnet.is/arsskyrslur/arsskyrsla-2021/frammistoduskyrsla/tidnistyring-og-spennugaedi/ (2020), [Online; accessed June-2022].
  * ENTSO-E [2021] ENTSO-E, _Operational Limits and Conditions for Mutual Frequency support over HVDC_ (2021).
  * Ulbig et al. [2014] A. Ulbig, T. S. Borsche, and G. Andersson, IFAC Proceedings Volumes 47, 7290 (2014).
  * Milano et al. [2018] F. Milano, F. Dörfler, G. Hug, D. J. Hill, and G. Verbič, in _2018 power systems computation conference (PSCC)_ (IEEE, 2018), pp. 1–25.
  * Pillai and Bak-Jensen [2010] J. R. Pillai and B. Bak-Jensen, in _IEEE PES General Meeting_ (IEEE, 2010), pp. 1–8.
  * Gorjão et al. [2020] L. R. Gorjão, M. Anvari, H. Kantz, C. Beck, D. Witthaut, M. Timme, and B. Schäfer, IEEE Access 8, 43082 (2020).
  * Parhizi et al. [2015] S. Parhizi, H. Lotfi, A. Khodaei, and S. Bahramirad, IEEE Access 3, 890 (2015).
  * Lidula and Rajapakse [2011] N. Lidula and A. Rajapakse, Renewable and Sustainable Energy Reviews 15, 186 (2011).
  * Lagrange et al. [2020] A. Lagrange, M. de Simón-Martín, A. González-Martínez, S. Bracco, and E. Rosales-Asensio, International Journal of Electrical Power & Energy Systems 119, 105865 (2020).
  * Palizban et al. [2014] O. Palizban, K. Kauhaniemi, and J. M. Guerrero, Renewable and Sustainable Energy Reviews 36, 428 (2014).
  * Teymouri et al. [2017] F. Teymouri, T. Amraee, H. Saberi, and F. Capitanescu, IEEE Transactions on Smart Grid 10, 1735 (2017).
  * Schäfer et al. [2018a] B. Schäfer, D. Witthaut, M. Timme, and V. Latora, Nature Communications 9, 1 (2018a).
  * Schäfer et al. [2018b] B. Schäfer, C. Beck, K. Aihara, D. Witthaut, and M. Timme, Nature Energy 3, 119 (2018b).
  * Schafer et al. [2018] B. Schafer, M. Timme, and D. Witthaut, in _2018 IEEE PES Innovative Smart Grid Technologies Conference Europe (ISGT-Europe)_ (IEEE, 2018), pp. 1–5.
  * Vorobev et al. [2019] P. Vorobev, D. M. Greenwood, J. H. Bell, J. W. Bialek, P. C. Taylor, and K. Turitsyn, IEEE Transactions on Power Systems 34, 3098 (2019).
  * Kruse et al. [2020] J. Kruse, B. Schäfer, and D. Witthaut, IEEE Access 8, 149435 (2020).
  * Hossain et al. [2019] E. Hossain, I. Khan, F. Un-Noor, S. S. Sikander, and M. S. H. Sunny, IEEE Access 7, 13960 (2019).
  * Yang et al. [2021] S.-G. Yang, B. J. Kim, S.-W. Son, and H. Kim, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 123127 (2021).
  * Kruse et al. [2021a] J. Kruse, B. Schäfer, and D. Witthaut, Patterns 2, 100365 (2021a).
  * Kyesswa et al. [2020] M. Kyesswa, H. Çakmak, U. Kühnapfel, and V. Hagenmeyer, IET Renewable Power Generation 14, 4010 (2020).
  * Chai et al. [2016] J. Chai, Y. Liu, J. Guo, L. Wu, D. Zhou, W. Yao, Y. Liu, T. King, J. R. Gracia, and M. Patel, in _2016 Power Systems Computation Conference (PSCC)_ (2016).
  * Rydin Gorjão et al. [2020] L. Rydin Gorjão, R. Jumar, H. Maass, V. Hagenmeyer, G. C. Yalcin, J. Kruse, M. Timme, C. Beck, D. Witthaut, and B. Schäfer, Nature Communications 11, 6362 (2020).
  * Jumar et al. [2021] R. Jumar, H. Maaß, B. Schäfer, L. R. Gorjão, and V. Hagenmeyer, arXiv preprint arXiv:2006.01771 (2021).
  * [28] L. R. Gorjão, B. Schäfer, and G. Hassan, _Power-grid frequency databasepermalink_ , https://power-grid-frequency.org/.
  * ENT [2022] _ENTSO-E Transparency Platform_ , https://transparency.entsoe.eu/ (2022).
  * Research [2020] N. E. Research, _The Nordics: Share of renewables has increased_ , https://www.nordicenergy.org/figure/share-of-renewables-has-increased/ (2020), [Online; accessed 30-March-2022].
  * Tróndheim et al. [2021] H. M. Tróndheim, B. A. Niclasen, T. Nielsen, F. F. D. Silva, and C. L. Bak, IEEE Open Access Journal of Power and Energy 8, 23 (2021).
  * SONI [2022] SONI, _What is the DS3 Programme?_ , https://www.soni.ltd.uk/__uuid/a6f0ce76-c5a9-4120-8b53-cac2daa99840/ds3-programme/ (2022), [Online; accessed 07-April-2022].
  * Pütz et al. [2022] S. Pütz, B. Schäfer, D. Witthaut, and J. Kruse, arXiv preprint arXiv:2204.10727 (2022).
  * Curto and Trapanese [2018] D. Curto and M. Trapanese, in _2018 IEEE International Conference on Environment and Electrical Engineering and 2018 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I CPS Europe)_ (2018), pp. 1–6.
  * [35] O. N. E. Authority, _Electricity_ , https://nea.is/hydro-power/electric-power/nr/69, [Online; accessed 10-March-2022].
  * Fingrid [2022] Fingrid, _Nordic power system frequency measurement data_ (2022), URL https://data.fingrid.fi/en/dataset/frequency-historical-data.
  * Maaß et al. [2013] H. Maaß, H. K. Çakmak, W. Suess, A. Quinte, W. Jakob, K.-U. Stucky, and U. G. Kuehnapfel, IEEE Transactions on Instrumentation and Measurement 62, 2384 (2013).
  * Maaß et al. [2015] H. Maaß, H. K. Cakmak, F. Bach, R. Mikut, A. Harrabi, W. Süß, W. Jakob, K.-U. Stucky, U. G. Kühnapfel, and V. Hagenmeyer, EURASIP Journal on Advances in Signal Processing 2015, 14 (2015).
  * Förstner et al. [2022] E. Förstner, R. Jumar, O. D. Tanrikulu, H. Maaß, U. Kühnapfel, and V. Hagenmeyer, in _The 2022 International Conference on Smart Grid Synchronized Measurements and Analytics – SGSMA_ (2022).
  * [40] T. L. Onsaker, _Code for predicting islanded power grid frequency_ , https://github.com/thlonsaker/predictability-of-island-power-grids.git.
  * Kruse et al. [2021b] J. Kruse, B. Schäfer, and D. Witthaut, in _2021 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_ (IEEE, 2021b), pp. 133–139.
  * Weißbach and Welfonder [2009] T. Weißbach and E. Welfonder, in _2009 IEEE/PES Power Systems Conference and Exposition_ (IEEE, 2009), pp. 1–6.
  * Gardiner et al. [1985] C. W. Gardiner et al., _Handbook of stochastic methods_ , vol. 3 (springer Berlin, 1985).
  * Böttcher et al. [2022] P. C. Böttcher, L. R. Gorjão, C. Beck, R. Jumar, H. Maass, V. Hagenmeyer, D. Witthaut, and B. Schäfer, arXiv preprint arXiv:2204.07508 (2022).
  * de España. [2022] R. E. de España., _Islas Baleares - Electricity demand tracking in real time_ , https://demanda.ree.es/visiona/baleares/baleares/total (2022), [Online; accessed 08-February-2022].
  * Eléctrica [2022] R. Eléctrica, _Singular features of the Balearic system_ , https://www.ree.es/en/about-us/ree-2-minutes (2022), [Online; accessed 08-February-2022].
  * Martínez-Barbeito et al. [2021] M. Martínez-Barbeito, D. Gomila, and P. Colet, in _ENERGY 2021, The Eleventh International Conference on Smart Grids, Green Communications and IT Energy-aware Technologies, 2021, p. 13-18_ (2021).
  * Lim and Zohren [2021] B. Lim and S. Zohren, Philosophical Transactions of the Royal Society A 379, 20200209 (2021).