11institutetext:  1 University of New South Wales, Australia
<EMAIL_ADDRESS>22institutetext: 2 Division of Fire Safety
Engineering, Lund University, Sweden<EMAIL_ADDRESS>
# Revisiting the paper ”Simulating dynamical features of escape panic”: What
have we learnt since then?

Milad Haghani1  Enrico Ronchi2

###### Abstract

The paper ”Simulating dynamical features of escape panic” by Helbing, Farkas,
and Vicsek, published over two decades ago in Nature, has left an indelible
mark on the field of crowd dynamics. With nearly 3,000 citations to date,
according to the Web of Science records, and significant influence, it has
shaped the crowd dynamics field. This analysis investigates the overall
influence of this paper through a variety of indicators, mapping its reach
across research areas. The intellectual foundation of the paper is traced,
examining the references cited. The terminological impact is also explored,
showing how the paper made use of terms like ”panic” and ”herding”. Moreover,
the alignment of the assumptions of the paper with empirical evidence is
discussed, finding discrepancies in key assertions about panic behaviour. The
numerical simulations of the paper and observations have significantly
influenced the field, such as for the case of the ”faster-is-slower”
phenomenon. The paper remains a key pillar in crowd dynamics, nevertheless, we
advocate for a new course of the field shifting away from the terminology
adopted in the paper and focusing more on empirical evidence.

###### keywords:

Social force model pedestrian dynamics crowd dynamics evacuation simulation

##  1 Introduction

More than two decades have passed since the publication of one of the most
influential papers in the field of crowd dynamics authored by Helbing, Farkas,
and Vicsek, which was featured in Nature (Helbing et al. [1]). Remarkably,
this paper has obtained nearly 3,000 citations according to Web of Science and
6,000 citations according to Google Scholar up to now. Its influence extends
beyond mere numerical counts, shaping the trajectory of crowd dynamics
research in profound ways since its inception. The paper is not just the most
cited paper in the field; it has also fundamentally shaped the discourse and
development of crowd dynamics as a standalone discipline. This prompts us to
consider how differently our understanding of crowd dynamics would have
evolved had this groundbreaking paper not appeared on the pages of Nature. Has
this influence predominantly steered the field in a positive direction, or
have there been negative consequences? If the impact has been mixed, what
aspects have been beneficial and what aspects have not? Our investigation
operates on two distinct levels: a bibliometric analysis and a content-focused
examination. We undertake several key inquiries:
(1) Overall influence: Through analysis of bibliometric indicators, we
quantify the influence of the paper based on objective metrics, while also
mapping the extent of influence that this paper has had and the areas of
research to which its influence has permeated.
(2) Intellectual Foundation: We trace the intellectual lineage that led to the
genesis of this paper by dissecting its references, and in some cases, the
references within those references. By doing so, we aim to unveil the origins
and foundations of concepts like escape panic as referred by the authors.
(3) Terminological Influence: We explore how the language and terminology
employed in this paper have shaped the broader lexicon of crowd dynamics. To
gauge this influence, we analyse the language employed in the field prior to
and subsequent to the publication of this paper.
(4) Assumptions vs. Empirical Evidence: Helbing et al. [1] outlined a set of
assumptions characterising crowd (escape) panic. Since then, numerous
empirical studies have emerged, testing these assumptions through various
experimental methods. We critically assess whether this accumulated body of
evidence aligns with the conceptual characteristics laid out in the original
paper.
(5) Numerical Simulations and Experimental Replication: Following the
introduction of the social force modelling concept, Helbing et al. [1]
conducted a series of numerical simulations that yielded noteworthy
conclusions. These findings have permeated the crowd dynamics literature. We
evaluate the extent to which subsequent experimental studies have investigated
these numerical observations.
Overall, this study endeavours to provide an evidence-based revisiting of the
influence, intellectual underpinnings, terminological evolution, and
validation of assumptions and numerical findings within the referenced work.
Such an assessment can help guide the trajectory of the field, offering
insights into whether a conscious paradigm shift is necessary to address the
issues associated with its influence. From here on, we refer to the paper as
the HFV paper (after the initials of the authors’ surnames).

##  2 Overall Influence

At the time of this analysis, we identified nearly 3,000 papers that have
referenced the HFV paper, according to the Web of Science, making it the
second most cited paper ever published in the crowd dynamics literature. It
follows an earlier work by Helbing and Molnar [2] where they laid the
foundation of the Social Force (SF) model of pedestrians and introduced the
model for the first time. The rate of annual citations for the HFV paper
exhibited a consistent linear increase since its publication and that trend
continued until 2015, when there was a noticeable change in this growing trend
(See fig:1). This may coincide with a paradigm shift in crowd dynamics
research around that year. Researchers began to shift their focus
significantly towards experimental work, diverging from the predominant trend
of numerical work up until that point.

Figure 1: The number of citing articles of HFV paper every year since its
publication (part (a)), the countries of origin of the citing articles and the
clusters of citing articles and their relative prominence (part (c)).

One noteworthy aspect is the substantial portion of these citations
originating from papers authored by researchers based in China, indicating the
exceptional popularity and influence of this paper within this community. The
majority of citing articles belong to the fields of Physics and Computer
Science, with a notable concentration of these articles in the journal Physica
A. Within the field of crowd dynamics, the HFV paper has also left its mark on
influential papers such as [3] and [4] which introduced the concept of
cellular automata modelling to the field in 2001. To gain a deeper
understanding of the areas where the HFV paper has exerted its influence, we
conducted a document co-citation analysis among the reference lists of the
nearly 3,000 citing papers. This analysis identifies clusters of references
that are co-cited by these papers. In other words, it reveals groups of
references that tend to appear together in the citations to the HFV paper by
its citing articles. We identified ten such clusters, and by analysing the key
phrases in the titles of these citing articles, we identified and categorised
the themes of studies that have referenced the HFV paper. These themes are
listed in app:1 in order of significance, showing how the influence of this
paper has extended beyond the field of crowd dynamics. Additionally, app:1
provides a list of references that have been most frequently co-cited with the
HFV paper.

##  3 Intellectual Foundation

The HFV paper cites a total of twenty references. In our analysis, we do not
intend to scrutinise every single reference, as not all appear to be equally
fundamental in shaping its intellectual foundation. Notably, the initial
references cited in the HFV paper have played a pivotal role in setting the
stage and justifying the study and modelling of escape panic. An intriguing
observation concerning these initial references is that not all of them appear
to support the argument; in fact, some seem to directly contradict the
statements they accompany. For instance, the paper by Keating [5] is cited
following a statement asserting ”Sometimes this behaviour [panic] is triggered
in life-threatening situations such as fires”. The reference title [5], The
myth of panic, and its content strongly suggests that it does not align with
the statement, but rather points in the opposite direction. Similarly, the
subsequent statement reads, ”At other times, stampedes can arise during the
rush for seats”. One of the two references cited following this statement is
the 1987 study by Johnson [6]. However, the content and conclusions of
Johnson’s paper appear to contradict the statement. This study [6] analyses
empirical evidence related to a tragic incident before a rock concert and
reports “evidence showing that panic did not cause the death and injury of
numerous young people”. Furthermore, it highlights that post-disaster
interviews and event reconstructions revealed no signs of stampede; instead,
most competition stemmed from people’s attempts to escape a crowd crush. This
contradicts the notion that panic can arise from people’s rush for seats, as
stated by HFV. Interestingly, in an editorial written by David Low about the
HFV paper in the same issue of Nature, Johnson’s paper is inaccurately cited
following this statement, ”The consequences of crushing, trampling, and panic
in crowds are well known”. References [5] and [6] are repeatedly cited by HFV
throughout the paper without them offering support in favour of the arguments
that they accompany.

Another example arises in the statement, ”Panicking individuals tend to show
maladaptive and relentless mass behaviour like jamming and life-threatening
overcrowding, which has often been attributed to social contagion”. The 1957
work of Quarantelli [7] is frequently cited following many of these
statements. Published in 1957, it has been hard to trace this reference in
contemporary databases. However, we managed to locate a more recent working
paper by Quarantelli from 2001 [8], which we believe would naturally reflect
the author’s observations and reports from earlier stages of their career. The
working paper is titled The Sociology of Panic and it clearly indicates that
the Quarantelli’s findings could not have possibly been consistent with the
portrayal presented in the HFV paper. Taken together, these observations cast
doubt on some of the foundational ideas upon which the HFV paper appears to
rely. In fact, certain aspects of the work cited seem to fundamentally
contradict key premises of the paper or at least reflect a possible inaccurate
use of terminology.

##  4 Terminological Influence

An examination of the text within the article reveals the prominent usage of
two terms: panic (mentioned 37 times) and herding (mentioned 10 times). The
notable prevalence of these terminologies within the field of pedestrian
dynamics can be directly attributed to the influential impact of this article.
As shown in [9], prior to the publication of this article, there was minimal
mention of these two terms in the context of studies on crowds and
evacuations. Similarly, the terminology of crowd dynamics (mentioned 3 times),
which has been detected in the titles, abstracts, or keywords of more than 700
articles since 2000, was scarcely identifiable in any publications within the
aforementioned domains before the publication of this article. It is
reasonable to conclude that the nomenclature by which this field is recognised
owes its existence to the terminological influence exerted by this article.

##  5 Assumptions vs. Empirical Evidence

The foundation of the HFV paper hinges on the concept of panic, which it
delineates through a set of defining characteristics. In this analysis, we
prospectively reevaluate some of these key assumptions to determine whether
they align with the empirical evidence that has emerged in the field since
that publication. People move or try to move considerably faster than normal.
This behaviour does not necessarily indicate panic per se. Individuals start
pushing, and interactions among people become physical in nature. Empirical
evidence does not support the notion that people engage in competitive
behaviour when responding to a life-threatening situation [10, 11]. At exits,
arching and clogging are observed. Arching is not a particular characteristic
of panic; rather, it is observable to some extent in laboratory experiments
when a crowd of pedestrians is instructed to pass through a narrow bottleneck
[12]. Jams build up. Similar to the previous point, the formation of a
pedestrian traffic jam is linked to a situation where the inflow of
pedestrians exceeds the capacity of a passage. People show a tendency toward
mass behaviour, that is, to do what other people do. This assumption is
particularly important because it is readily operationalisable in evacuation
models. In many cases, in fact, crowd models have been assessed on their
ability to replicate such phenomenon. However, empirical evidence accumulated
since the HFV paper suggests that the matter of imitation is complex. Social
influence has been observed in several experimental studies, but there are
several factors and variables affecting its extent, e.g., crowd size,
familiarity with the environment, and perceived urgency [13]. While there is
evidence of some level of imitation regarding decision adaptation, it does not
lead to blind herd-like behaviour where all individuals follow a single action
(e.g., as visually depicted by the Figure 1 of the Editorial piece by David
Low [14] on the HFV paper, marking one of the most influential illustrations
produced in this field). In addition, empirical evidence has overwhelmingly
suggested that humans take into consideration a range of factors in flights
situations and not purely the decisions of others [15].

##  6 Numerical Simulations and Experimental Replication

The HFV article subjected its proposed crowd escape panic model to numerical
testing, leading to recommendations that have significantly influenced both
the field and public perception. Among these observations, the most
consequential is the so-called faster-is-slower phenomenon. The study
considered the evacuation of a crowd of 200 agents through a 1m-wide exit
while gradually increasing a parameter known as desired velocity. The outcome
revealed that as the desired velocity increased, system efficiency (in this
context, the inverse of total evacuation time) initially improved. However,
further increases in this parameter led to reduced system efficiency, implying
longer evacuation times. This finding was translated into the recommendation
that the most efficient way to evacuate and survive a crowded space with
limited exit capacity, relative to occupancy, is for individuals to remain
patient and keep lower speeds when passing through bottlenecks. Another
suggestion was that asymmetrical placement of columns in front of exits can
improve outflow. Several dimensions warrant consideration: (1) Empirical
testing: Empirical studies have in some instances confirmed and in other
contradicted the faster-is-slow and column phenomena [16, 17, 18]. In a
controlled experiment, participants were instructed to exit forcefully
(without resorting to aggression), resulting in significantly higher flow
efficiency compared to a calm and patient manner [16]. (2) Model results need
context: Numerical findings could be very parameter specific. According to
Figure 1 (c) of the original paper, the minimum time required for 200 people
to evacuate a room with a single 1m-wide exit is approximately 120s. Empirical
observations indicate that this number of people may complete this evacuation
in different (even faster) times depending on several factors. (3) Desired
velocity in non-free flow conditions: The concept of desired velocity is hard
to grasp physically, except in free-flow conditions. In scenarios with
bottlenecks or obstacles, it remains unclear how different levels of desired
flow correspond to various levels of escape motivation. It is possible that in
reality, assuming that the model components can manifest physically in the
real world (at least in terms of desired velocity), we are always within the
monotonically decreasing limb of the graph in Figure 1(c) of the HFV article.
This would imply that an increase in evacuation time does not manifest for
typical (non-aggressive) forces applied by people at a bottleneck. This
observation has had profound implications and has influenced both research
studies that employed ants and mice as models for human crowds [18] as well as
public perceptions regarding efficient evacuations. The main argument is that
the panic concept assumed by the model may be misinterpreted as an
encouragement to stay on the beginning of the left-hand side of the
flow/density relationship, whereas we know that flow/density reaches a maximum
with a speed which is at an optimal density value, and not at slow speed
level.

##  7 Conclusions

Our analysis suggests that the HFV paper stands as one of the most fundamental
and influential contributions to the field of crowd dynamics. It has not only
profoundly shaped the terminologies employed within this field but has also
significantly influenced the trajectory of research undertaken by scholars. It
is evident that the introduction of the social force model of pedestrian
dynamics marks a major paradigm shift in crowd research. This model introduced
a paradigm that has spurred numerous developments in agent-based modelling of
pedestrian traffic and beyond, establishing itself as a standard benchmark.
Patterns of referencing to the article demonstrate its broad-reaching impact,
reflecting its use as an intellectual foundation and source of inspiration in
a multitude of domains. Our work highlights aspects that facilitate the
interpretation of the HFV paper: (1) Model limitations: There are a set of
misconceptions regarding the model applicability. The model excels mostly when
integrated with other layers of modelling, addressing various aspects of
pedestrian behaviour and decision-making beyond the operational level and
step-taking behaviour [19]. (2) Calibration procedure: Presenting the model
along with a detailed sensitivity analysis to identify critical model
parameters would have helped understanding their impact on simulation output.
A standardised calibration procedure could have facilitated future model
applications. In the social force model paradigm, the forces lack physical
manifestations, making them metaphorical or imaginary concepts. This renders
them unmeasurable (though some have attempted to measure these forces in
experimental settings [20]), making it exceedingly difficult to link
experimental observations to the model components and calibrate parameters
using established procedures. This led to partial calibration attempts [19].
(3) Model validation: The numerical case studies presented in the paper,
including the faster-is-slower phenomenon and the obstacle effect, have been
later scrutinised showing they may not hold for all scenarios. These findings
have directed research efforts toward replicating them, in some instances
relying on non-valid methods, such as animal experiments [18]. Relying on
assumptions about panic, unsupported by references, has significantly
influenced the direction of research in this field.

In summary, the HFV paper has played a pivotal role in the establishment of
crowd dynamics as a research field. While we recognise the methodological and
modelling benefits derived from this paper, we advocate for a paradigm shift
in crowd dynamics research. It is essential to recognise that the mere
publication of this paper in Nature should not serve as an unequivocal
endorsement of the quest to delineate and explore the concept of crowd panic.
Thus, we call for a reimagining of the field, focused on empirical evidence to
chart a new course for research in crowd dynamics.

Milad Haghani: Conceptualization, Methodology, Software, Validation, Formal
Analysis, Investigation, Writing – Original Draft, Visualization Enrico
Ronchi: Conceptualization, Methodology, Validation, Formal Analysis,
Investigation, Writing – Review and Editing

## References

  * [1] Helbing, D., Farkas, I., Vicsek, T.: Simulating dynamical features of escape panic.  Nature 407(6803), 487–490 (2000)
  * [2] Helbing, D., Molnar, P.: Social force model for pedestrian dynamics.  Physical review E 51(5), 4282 (1995)
  * [3] Helbing, D.: Traffic and related self-driven many-particle systems.  Reviews of modern physics 73(4), 1067 (2001)
  * [4] Burstedde, C., Klauck, K., Schadschneider, A., Zittartz, J.: Simulation of pedestrian dynamics using a two-dimensional cellular automaton.  Physica A: Statistical Mechanics and its Applications 295(3-4), 507–525 (2001)
  * [5] Keating, J.P.: The myth of panic.  Fire Journal 76(3), 57–61 (1982)
  * [6] Johnson, N.R.: Panic at “the who concert stampede”: an empirical assessment.  Social Problems 34(4), 362–373 (1987)
  * [7] Quarantelli, E.L.: The behavior of panic participants.  Sociology and Social Research 41, 187–194 (1957)
  * [8] Quarantelli, E.L.: Panic, sociology of (2001)
  * [9] Feliciani, C., Corbetta, A., Haghani, M., Nishinari, K.: How crowd accidents are reported in the media: Lexical and sentiment analyses.  arXiv preprint arXiv:2309.14633 (2023)
  * [10] Drury, J., Cocking, C., Reicher, S., Burton, A., Schofield, D., Hardwick, A., Graham, D., Langston, P.: Cooperation versus competition in a mass emergency evacuation: A new laboratory simulation and a new theoretical model.  Behavior research methods 41(3), 957–970 (2009)
  * [11] Bartolucci, A., Casareale, C., Drury, J.: Cooperative and competitive behaviour among passengers during the costa concordia disaster.  Safety science 134, 105055 (2021)
  * [12] Uesten, E., Schumann, J., Sieben, A.: Exploring the dynamic relationship between pushing behavior and crowd dynamics.  Collective Dynamics 8, 1–29 (2023)
  * [13] Haghani, M., Cristiani, E., Bode, N.W., Boltes, M., Corbetta, A.: Panic, irrationality, and herding: three ambiguous terms in crowd dynamics research.  Journal of advanced transportation 2019 (2019)
  * [14] Low, D.J.: Following the crowd.  Nature 407(6803), 465–466 (2000)
  * [15] Kinateder, M., Comunale, B., Warren, W.H.: Exit choice in an emergency evacuation scenario is influenced by exit familiarity and neighbor behavior.  Safety science 106, 170–175 (2018)
  * [16] Haghani, M., Sarvi, M., Shahhoseini, Z.: When ‘push’does not come to ‘shove’: Revisiting ‘faster is slower’in collective egress of human crowds.  Transportation research part A: policy and practice 122, 51–69 (2019)
  * [17] Zuriguel, I., Echeverría, I., Maza, D., Hidalgo, R.C., Martín-Gómez, C., Garcimartín, A.: Contact forces and dynamics of pedestrians evacuating a room: the column effect.  Safety science 121, 394–402 (2020)
  * [18] Haghani, M.: Empirical methods in pedestrian, crowd and evacuation dynamics: Part i. experimental methods and emerging topics.  Safety science 129, 104743 (2020)
  * [19] Haghani, M., Sarvi, M.: Crowd model calibration at strategic, tactical, and operational levels: Full-spectrum sensitivity analyses show bottleneck parameters are most critical, followed by exit-choice-changing parameters.  Transportation Letters pp. 1–28 (2023)
  * [20] Zhao, Y., Lu, T., Su, W., Wu, P., Fu, L., Li, M.: Quantitative measurement of social repulsive force in pedestrian movements based on physiological responses.  Transportation research part B: methodological 130, 1–20 (2019)

##  Appendix A Appendix 1

Clusters of topics and research areas where the HFV paper has been cited,
based on a document co-citation analysis of the citing articles of the HFV
paper:
Cluster 1: pedestrian evacuation; social force model; bidirectional pedestrian
flow; new lattice; visual field.
Cluster 2: crowd space; predictive crowd analysis technique; dense crowd;
emotional contagion; personality trait.
Cluster 3: jamming transition; pedestrian evacuation; fire evacuation model;
pedestrian dynamics; self-driven particle.
Cluster 4: collective escape; pedestrian counter flow; emergency escape;
pedestrian contact force; pedestrian dynamics.
Cluster 5: collective motion; self-propelled particle; collective decision;
nonlinear dynamics.
Cluster 6: pedestrian trajectory prediction; evacuation crowd dynamics; human
behaviour; social interaction.
Cluster 7: dynamic decision behaviour; information service; discrete opinion
dynamics; critical market crash.
Cluster 8: modelling pedestrian crowd; panic condition; crowd guidance;
disaster evacuation; aircraft disembarking.
Cluster 9: switching phenomena; financial market; correlated randomness; trend
switching processes; trend switching.
Cluster 10: granular material; clogging transition; horizontal hopper; hopper
angle; large reflective obstacle.

References most frequently co-cited by the HFV paper:
1104 Co-citations: Helbing D, 1995, PHYS REV E, V51, P4282,
DOI 10.1103/PhysRevE.51.4282
565 Co-citations: Burstedde C, 2001, PHYSICA A, V295, P507,
DOI 10.1016/S0378-4371(01)00141-8
392 Co-citations: Kirchner A, 2002, PHYSICA A, V312, P260,
DOI 10.1016/S0378-4371(02)00857-9
380 Co-citations: Helbing D, 2005, TRANSPORT SCI, V39, P1,
DOI 10.1287/trsc.1040.0108
306 Co-citations: Helbing D, 2007, PHYS REV E, V75, P0,
DOI 10.1103/PhysRevE.75.046109
291 Co-citations: Helbing D, 2001, REV MOD PHYS, V73, P1067,
DOI 10.1103/RevModPhys.73.1067
285 Co-citations: Hughes RL, 2002, TRANSPORT RES B-METH, V36, P507,
DOI 10.1016/S0191-2615(01)00015-7
270 Co-citations: Reynolds C.W., 1987, ACM SIGGRAPH COMPUTE, V21, P25,
DOI 10.1145/37402.37406
268 Co-citations: Moussaid M, 2011, P NATL ACAD SCI USA, V108, P6884,
DOI 10.1073/pnas.1016507108
261 Co-citations: Helbing D, 2002, PEDESTRIAN AND EVACUATION DYNAMICS, V0, P21
229 Co-citations: Muramatsu M, 1999, PHYSICA A, V267, P487,
DOI 10.1016/S0378-4371(99)00018-7
207 Co-citations: Vicsek T, 1995, PHYS REV LETT, V75, P1226,
DOI 10.1103/PhysRevLett.75.1226
201 Co-citations: Helbing D, 2000, PHYS REV LETT, V84, P1240,
DOI 10.1103/PhysRevLett.84.1240
193 Co-citations: Kirchner A, 2003, PHYS REV E, V67, P0,
DOI 10.1103/PhysRevE.67.056122
181 Co-citations: Henderson LF, 1971, NATURE, V229, P381,
DOI 10.1038/229381a0
169 Co-citations: Treuille A, 2006, ACM T GRAPHIC, V25, P1160,
DOI 10.1145/1141911.1142008
167 Co-citations: Helbing D, 2003, PHYS REV E, V67, P0,
DOI 10.1103/PhysRevE.67.067101
166 Co-citations: Hughes RL, 2003, ANNU REV FLUID MECH, V35, P169,
DOI 10.1146/annurev.fluid.35.101101.161136
163 Co-citations: Zheng XP, 2009, BUILD ENVIRON, V44, P437,
DOI 10.1016/j.buildenv.2008.04.002


#  Compact localised states in magnonic Lieb lattices

Grzegorz Centała  Jarosław W. Kłos<EMAIL_ADDRESS>Institute of Spintronics
and Quantum Information, Faculty of Physics, Adam Mickiewicz University,
Poznań, Uniwersytetu Poznańskiego 2, Poznań 61-614, Poland

(March 26, 2023)

###### Abstract

Lieb lattice is one of the simplest bipartite lattices where compact localized
states (CLS) are observed. This type of localisation is induced by the
peculiar topology of the unit cell, where the modes are localized only on one
sublattice due to the destructive interference of partial waves. The CLS exist
in the absence of defects and are associated with the flat bands in the
dispersion relation. The Lieb lattices were successfully implemented as
optical lattices or photonic crystals. This work demonstrates the possibility
of magnonic Lieb lattice realization where the flat bands and CLS can be
observed in the planar structure of sub-micron in-plane sizes. Using forward
volume configuration, we investigated numerically (using the finite element
method) the Ga-dopped YIG layer with cylindrical inclusions (without Ga
content) arranged in a Lieb lattice of the period 250 nm. We tailored the
structure to observe, for the few lowest magnonic bands, the oscillatory and
evanescent spin waves in inclusions and matrix, respectively. Such a design
reproduces the Lieb lattice of nodes (inclusions) coupled to each other by the
matrix with the CLS in flat bands.

Keywords

flat bands, compact localized states, Lieb lattice, spin waves, finite element
method

††preprint: APS/123-QED

##  I  Introduction

There are many mechanisms leading to wave localization in systems with long-
range order, i.e. in crystals or quasicrystals. The most typical of these
require (i) the local introduction of defects, including the defects in the
form of surfaces or interfaces [1] (ii) the presence of global disorder [2],
(iii) the presence of external fields [3] or (iv) the existence of many-body
phenomena [4]. However, since at least the late 1980s, it has been known that
localization can occur in unperturbed periodic systems in the absence of
fields and many-body effects, and is manifested by the presence of flat, i.e.,
dispersion-free bands in the dispersion relation. The pioneering works are
often considered to be the publications of B. Sutherland [5] and E. H. Lieb
[6], who found the flat bands of zero energy [7] for bipartite lattices with
use of the tight-binding model Hamiltonians, where the hoppings occur only
between sites of different sublattices. The simplest realization of this type
of system is regarded as the Lieb lattice [6, 8], where the nodes of one
square sublattice, of coordination number $z=4$, connect to each other only
via nodes with a coordination number $z=2$ from other two square sublattices
(Fig. 1). In the case of extended Lieb lattices [9, 10], the nodes of $z=2$
form chains: dimmers, trimmers, etc.(Fig. 2). An intuitive explanation for the
presence of the flat bands is the internal isolation of excitations located in
one of the sublattices. The cancelling of excitations at one sublattice is the
result of forming destructive interference and local symmetry within the
complex unit cell [11]. When only one of the sublattices is excited, the other
sublattice does not mediate the coupling between neighbouring elementary
cells, and the phase difference between the cells is irrelevant to the energy
(or the frequency) of the eigenmode on the whole lattice - i.e. the Bloch
function. Modes of this type are therefore degenerated for different wave
vector values in infinite lattices. We are dealing here with the localization
on specific arrangements of structure elements, which are isolated from each
other. Such kinds of modes are called compact localized states (CLS) [12, 13,
14, 15, 16] and show a certain resistance to the introduction of defects [17,
18]. The flat band systems with CLS are the platform for the studies of
Anderson localization [19], and unusual properties of electric conductivity
[20]. A similar localization is observed in the quasicrystals, where the
arrangements of the elements composing the structure are replicated
aperiodically and self-similarly throughout the system [21, 22] and the
excitation can be localized on such patterns. The CLS in finite Lieb lattices
have a form of loops (plaquettes) occupying the majority nodes ($z=2$). These
states are linearly dependent and do not form a complete basis for the flat
band. Therefore occupancy gaps need to be filled (for infinite lattice) by
states occupying only one sublattice of majority nodes, localizes at lines,
called noncontractible loop states (NLS) [14, 15, 23].

The topic of Lieb lattices and other periodic structures with compact
localization and flat bands was renewed [8] about 10 years ago when physical
realizations of synthetic Lieb lattices began to be considered for electronic
systems [24, 25], optical lattices [26, 27], superconducting systems [28, 29],
in phononics [30] and photonics [31, 14]. In a real system, where the
interaction cannot be strictly limited to the nearest elements of the
structure, the bands are not perfectly flat. Therefore, some authors use the
extended definition of the flat band to consider the bands that are flat only
along particular directions or in the proximity of high-symmetry Brillouin
zone points [32]. In tight-binding models, this effect can be included by
considering the hopping to at least next-nearest-neighbours [33, 34].
Similarly, the crossing of the flat band by Dirac cones can be transformed
into anti-crossing and lead to opening gaps, separating the flat band from
dispersive bands. This effect can in induced by the introduction of spin-orbit
term to tight-binding Hamiltonian (manifested by the introduction of Peierls
phase factor to the hopping) or by dimerization of the lattice (by alternative
changes of hoppings or site energies) [33, 34, 35, 36, 37]. The later scenario
can be easily observed in real systems where the position of rods/wells
(mimicking the sites of Lieb lattice) and contrast between them can be easily
altered [38]. Opening the narrow gap between flat band and dispersive bands
for Lieb lattice is also fundamentally interesting because it leads to the
appearance of so-called Landau-Zener Bloch oscillations [39].

The isolated and perfectly flat bands for Lieb lattices are topologically
trivial – their Chern number is equal to zero [40]. For weakly dispersive
(i.e. almost flat) bands the Chern numbers can be non-zero [41]. However, when
the flat band is intersected by dispersive bands then it can exhibit the
discontinuity of Hilbert–Schmidt distance between eigenmodes corresponding to
the wave vectors just before and just after the crossing. Such an effect is
called singular band touching [16]. This limiting value of Hilbert–Schmidt
distance is bulk invariant, different from the Chern number.

One of the motivations for the photonic implementation of systems with flat,
or actually nearly flat bands [42], was the desire to reduce the group
velocity of light in order to compress light in space, which leads to the
concentration of the optical signal and an increase in the light-matter
interaction, or the enhancement of non-linear effects. Another, more obvious
application is the possibility of realizing delay lines that can buffer the
signal to adjust the timing of optical signals [43]. A promising alternative
to photonic circuits are magnonic systems, which allow signals of much shorter
wavelengths to be processed in devices several orders of magnitude smaller
[44, 45]. For this reason, it seems natural to seek a magnonic realization of
Lieb lattices.

In this paper, we propose the realization of such lattices based on a magnonic
structure in the form of a perpendicularly magnetized magnetic layer with
spatially modulated material parameters or spatially varying static internal
field. Lieb lattices have been studied also in the context of magnetic
properties, mainly due to the possibility of enhancing ferromagnetism in
systems of correlated electrons [46], where the occurrence of flat bands with
zero kinetic energy was used to expose the interactions. There are also known
single works where the spin waves have been studied in the Heisenberg model in
an atomic Lieb lattice, such as the work on the magnon Hall effect [47]. But
the comprehensive studies of spin waves in nanostructures that realize
magnonic Lieb lattices and focus on wave effects in a continuous model have
not been carried out so far. In this work, we demonstrate the possibility of
realization of magnonic lattices in planar structure based on low spin wave
damping material: yttrium iron garnet (YIG) where the iron is partially
substituted by gallium (Ga). We present the dispersion relation with a weakly
depressive (flat) band exhibiting the compact localized spin waves. The flat
is almost intersected at the $M$ point of the 1st Brillouin zone by highly
dispersive bands, similar to Dirac cones. We discuss the spin wave spectra and
compact localized modes both for simple and extended Lieb lattices.

The introduction is followed by the section describing the model and numerical
method we used, which precedes the main section where the results are
presented and discussed. The paper is summarized by conclusion and
supplemented with additional materials where we showed: (A) the results for
extended Lieb-7 lattice, (B) an alternative magnonic Lieb lattice design via
shaping the demagnetizing field, and (C) an attempt of formation magnonic Lieb
lattice by dipolarly coupled magnetic nanoelements, (D) discussion of small
differences in the demagnetizing field of majority and minority nodes
responsible for opening a small gap in the Lieb lattice spectrum.

##  II  Structure

Magnonic crystals (MCs) are regarded as promising structures for magnonic-
based device applications [48, 45]. In our studies, we consider planar MCs to
design the magnonic Lieb lattice, owing to the relative ease of fabrication of
such structures and their experimental characterization [49, 50, 51]. We
proposed realistic systems that mimic the main features of the tight-binding
model of Lieb lattice [33, 16].

Investigated MCs consist of yttrium iron garnet doped with gallium (Ga:YIG)
matrix and yttrium iron garnet (YIG) cylindrical inclusions arranged in Lieb
lattice Fig. 1. Doping YIG with Gallium is a procedure where magnetic ${\rm
Fe}^{3+}$ ions are replaced by non-magnetic ${\rm Ga}^{3+}$ ions. This method
not only decreases saturation magnetization $M_{S}$ but, simultaneously,
arises uni-axial out-of-plane anisotropy, that ensures the out-of-plane
orientation of static magnetization in Ga:YIG layer at a relatively low
external field applied perpendicularly to the layer. Discussed geometry, i.e.
forward volume magnetostatic spin wave configuration, does not introduce an
additional anisotropy in the propagation of spin waves, related to the
orientation of static magnetization.

Figure 1: Basic magnonic Lieb lattice. The planar magnonic structure consists
of YIG cylindrical nanoelements embedded within Ga:YIG. Dimensions of the
ferromagnetic unit cell are equal to 250x250x59 nm and the unit cell contains
three inclusions of 50 nm diameter. (a) The structure of basic Lieb lattice,
and (b) top view of the Lieb lattice unit cell where the node (inclusion) from
minority sublattice $A$ and two nodes (inclusions) from two majority
sublattices $B$ are marked.

The design of the Lieb lattice requires the partial localization of spin wave
in inclusions, which can be treated as an approximation of the nodes from the
tight-binding model. Furthermore, the neighbouring inclusions in the lattice
have to be coupled strongly enough to sustain the collective spin wave
dynamics, and weakly enough to minimize the coupling between further
neighbours. Therefore, the geometrical and material parameters were selected
to ensure the occurrence of oscillatory excitations in the (YIG) inclusions
and exponentially evanescent spin waves in the (Ga:YIG) matrix. The size of
inclusions was chosen small enough to separate three lowest magnonic bands
with almost uniform magnetization precession inside the inclusion from the
bands of higher frequency, where the spin waves are quantised inside the
inclusions. Also, the thickness of the matrix and inclusion was chosen in a
way that there are no nodal lines inside the inclusion. The condition which
guarantee the focussing magnetization dynamics inside the inclusions is
fulfilled in the frequency range below the ferromagnetic resonance (FMR)
frequency of the out-of-plane magnetized layer made of Ga:YIG (matrix
material): $f_{\rm FMR,Ga:YIG}=$4.95 GHz and above the FMR frequency of out-
of-plane magnetized layer made of YIG (inclusions material): $f_{\rm
FMR,YIG}=2.42$ GHz. These limiting values were obtained using the Kittel
formula for out-of-plane magnetised film: $f_{\rm
FMR}=\frac{\gamma}{2\pi}\text{\textbar}{\mu_{0}H_{\rm 0}+\mu_{0}H_{\rm
ani}-\mu_{0}M_{S}\text{\textbar}}$, where we used the following values of
material parameters [52] for YIG: gyromagnetic ratio $\gamma=177~{}\frac{\rm
GHz}{\rm T}$, magnetization saturation $\mu_{0}M_{S}=182.4$ mT, exchange
stiffness constant $A=3.68~{}\frac{\rm pJ}{\rm m}$, (first order) uni-axial
anisotropy field $\mu_{0}H_{\rm ani}=-3.5~{}$mT, and for Ga:YiG:
$\gamma=179~{}\frac{\rm GHz}{\rm T}$, $\mu_{0}M_{S}=20.2$ mT,
$A=1.37~{}\frac{\rm pJ}{\rm m}$, $\mu_{0}H_{\rm ani}=94.1$ mT. Since the
greatest impact of the first order uniaxial anisotropy field ($\mu_{0}H_{\rm
ani}$), we decided to neglect higher order terms of uni-axial anisotropy and
cubic anisotropy of (Ga:)YIG. Due to the presence of out-of-plane anisotropy
and relatively low saturation magnetization, we could consider a small
external magnetic field $\mu_{0}H_{0}=100$ mT to reach saturation state.

Figure 2: Extended magnonic Lieb lattice – Lieb - 5. Dimensions of the unit
cell are 375x375x59 nm and contain 5 inclusions of size 50 nm in diameter.
Also, we maintain the same separation (distance between centres of
neighbouring sites is 125 nm) as for considered basic Lieb lattice – Fig. 1.
(a) The structure of Lieb-5 lattice, and (b) top view on Lieb-5 lattice unit
cell where the node (inclusion) from minority sublattice $A$ and four nodes
(inclusions) from two majority sublattices $B$ are marked.

It is worth noticing that without the evanescent spin waves in the
ferromagnetic matrix, the appropriate coupling between inclusions would not be
possible. Therefore the realization of the Lieb lattice in form of the array
of ferromagnetic nanoelemets embedded in air/vacuum seems to be very
challenging – see the exemplary results in Supplementary Information C.

We also tested the possibility of other realizations of magnonic Lieb
lattices. One solution seemed to be the design of a structure in which the
concentration of the spin wave amplitude in the Lieb lattice nodes would be
achieved through an appropriately shaped profile of the static demagnetizing
field – Supplementary Information B. However, the obtained results were not as
promising as for YIG/Ga:YIG system.

In the main part of the manuscript, we present the results for the basic Lieb
lattice (showed in Fig. 1) and extended Lieb-5 lattice (showed in Fig. 2),
based on YIG/Ga:YIG structures. The further extension of the Lieb lattice may
be realized by increasing the number of $B$ nodes between neighbouring $A$
nodes. Supplementary Information A presents the results for Lieb-7, where for
each site (inclusion) from minority sublattice $A$, we have six nodes
(inclusions), grouped in three-element chains, from majority sublattices $B$.

##  III  Methods

The spin waves spectra and the spatial profiles of their eigenmodes were
obtained numerically in a semi-classical model, where the dynamics of
magnetization vector $\textbf{M}(\textbf{r},t)$ is described by the Landau-
Lifshitz equation [53]:

$\frac{d\textbf{M}}{dt}=-\gamma\mu_{0}[\textbf{M}\times\textbf{H}_{\rm
eff}+\frac{\alpha}{M_{S}}\textbf{M}\times(\textbf{M}\times\textbf{H}_{\rm
eff})].$ (1)

The symbol $\textbf{H}_{\rm eff}(\textbf{r},t)$ denotes effective magnetic
field.

In numerical calculations, we neglected the damping term since $\alpha$ is
small both for YIG and for YIG with Fe substituted partially by Ga (for
$\alpha_{\rm Ga:YIG}=$6.1\text{\times}{10}^{-4}$$ and $\alpha_{\rm
YIG}=$1.3\text{\times}{10}^{-4}$$ [52]). The effective magnetic field $H_{\rm
eff}$ contains the following components: the external field $H_{0}$, exchange
field $H_{\rm ex}$, bulk uniaxial anisotropy field $H_{\rm ani}$ and dipolar
field $H_{\rm d}$:

$\textbf{H}_{\rm
eff}(\textbf{r},t)=H_{0}\hat{\mathbf{z}}+\frac{2A}{\mu_{0}M^{2}_{S}}\laplacian\textbf{M}(\textbf{r},t)+H_{\rm
ani}(\textbf{r})\hat{\mathbf{z}}-\gradient\varphi(\textbf{r},t),$ (2)

where the $z-$direction is normal to the plane of the magnonic crystal. We
assume that the sample is saturated in $z-$direction and magnetization vector
precesses around this direction. The material parameters ($M_{S}$, $A$,
$\alpha$ and $\gamma$) are constant within matrix and inclusions.

Using the magnetostatic approximation the dipolar term of the effective
magnetic field may be expressed as a gradient of magnetic scalar potential:

$\textbf{H}_{\rm d}(\textbf{r},t)=-\gradient\varphi(\textbf{r},t)$ (3)

By using the Gauss equation magnetic scalar potential may be associated with
magnetisation as follows:

$\laplacian\varphi(\textbf{r},t)=\divergence\textbf{M}(\textbf{r},t)$ (4)

Spin-wave dynamics is calculated numerically using the finite-element method
(FEM). We used the COMSOL Multiphysics [54] to implement the Landau-Lifshitz
equation (Eq. 1) and performed FEM computation for the defined geometry of
magnonic Lieb lattices. The COMSOL Multiphysics is the software used for
solving a number of physical problems, since many implemented modules it
becomes more and more convenient. Nevertheless, all the equations were
implemented in the Mathematics module which contains different forms of
partial differential equations. Eq. 1 was solved by using eigenfrequency
study, on the other hand, to solve Eq. 4 we used stationary study. To obtain
free decay of scalar magnetic potential in the model we applied $5~{}\mu$m of
a vacuum above and underneath the structure. At the bottom and top surface of
the model with vacuum, we applied the Dirichlet boundary condition. We use the
Bloch theorem for each variable (magnetostatic potential and components of
magnetization vector) at the lateral surfaces of a unit cell. We selected the
symmetric unit cell with minority node $A$ in the centre to generate a
symmetric mesh which does not perturb the four-fold symmetry of the system –
this approach is of particular importance for the reproduction of the
eigenmodes profiles in high-symmetry points. In our numerical studies, we used
2D wave vector $\textbf{k}=k_{x}\hat{\textbf{x}}+k_{y}\hat{\textbf{y}}$ as a
parameter for eigenvalue problem which was selected along the high symmetry
path $\Gamma-X-M-\Gamma$ to plot the dispersion relation. We considered the
lowest 3, 5 and 7 bands for basic Lieb lattice, Lieb-5 lattice and Lieb-7
lattice, respectively.

##  IV  Results

The tight-biding model of the basic Lieb lattice, with hopping restricted to
next-neighbours gives three bands in the dispersion relation. The top and
bottom bands are symmetric with respect to the second, perfectly flat band,
and intersect with this dispersionless band at $M$ point of 1st Brillouin
zone, with constant slope forming two Dirac cones[27, 8]. In a realistic
magnonic system, the spin wave spectrum showing the particle-hole symmetry
with a zero energy flat band is difficult to reproduce because (i) the
dipolarly dominated spin waves, propagating in magnetic film, experience a
significant reduction of the group velocity with an increase of the wave
vector (this tendency is reversed for much larger wave vectors were the
exchange interaction starts to dominate) [53], (ii) the dipolar interaction is
long-range. The first effect makes the lowest band wider than the third band,
and the latter one – induces the finite width of the second band [33]. We are
going to show, that this weakly dispersive band supports the existence of CLS.
Therefore, we will still refer to it as flat band, which is a common practice
for different kinds of realization of Lieb lattices in photonics or optical
lattices.

The results obtained for the basic magnonic Lieb lattice, (Fig. 1), are shown
in Fig. 3. As we predicted, three lowest bands form a band structure which is
similar to the dispersion relation known from the tight-binding model [10].
However, in a considered realistic system there is an infinite number of
higher bands, not shown in Fig. 3(a). For higher bands, spin waves can
propagate in an oscillatory manner in the matrix hence the system does not
mimic the Lieb lattice where the excitations should be associated with the
nodes (inclusions) of the lattice.

Figure 3: Dispersion relation for the basic magnonic Lieb lattice, containing
three inclusions in the unit cell: one inclusion $A$ from minority sublattice
and two inclusions $B$ from majority sublattices (see Fig. 1). (a) The
dispersion relation is plotted along the high symmetry path
$\Gamma$-X-M-$\Gamma$ (see the inset). The lowest band (blue) and the highest
band (red) create Dirac cones almost touching (b) in the M point. The middle
band (green) is relatively flat in the vicinity of the M point.

Due to the fourfold symmetry of the system, the dispersion relation could be
inspected along the high symmetry path $\Gamma-X-M-\Gamma$. Frequencies of the
first three bands are in the range $f_{\rm FMR,YIG}~{}-~{}f_{\rm FMR,Ga:YIG}$.
Their total width is about $\approx 0.78$ GHz. The first and third band form
Dirac cones at $M$ point, separated by a tiny gap $\approx 15$ MHz. The
possible mechanism responsible for opening the gap is a small difference in
the demagnetizing field in the areas of inclusions $A$ (from the minority
lattice) and inclusions $B$ (from two majority sublattices) – see
Supplementary Information D. Inclusions $A$ ($B$) have four (two) neighbours
of type $B$ ($A$). Although inclusions $A$ and $B$ have the same size and are
made of the same material, the static field of demagnetization inside them
differs slightly due to the different neighbourhoods. This effect is
equivalent to the dimerization of the Lieb lattice by varying the energy of
the nodes in the tight-binding model, which leads to the opening of a gap
between Dirac cones and parabolic flattening of them in very close proximity
to the $M$-point. It is worth noting that in the investigated system, the gap
opens between the first and second bands, while the second and third bands
remain degenerated at point $M$, with numerical accuracy.

The middle band can be described as weakly dispersive. The band is more flat
on the $X-M$ path and, in particular, in the vicinity of $M$ point – see Fig.
3(b). The small width of the second band can be attributed to long-range
dipolar interactions which govern the magnetization dynamics in a considered
range of sizes and wave vectors. It is known that even the extension of the
range of interactions to next-nearest-neighbours in the tight-binding model
induces the finite width of the flat band for the Lieb lattice.

Figure 4: The spin wave profiles obtained for the basic magnonic Lieb lattice,
composed of three inclusions in the unit cell (see Fig. 1). The modes are
presented for each band exactly at $M$ (left column) and in its proximity
($M^{\leftarrow}$) on the path $M-\Gamma$ (right column). In the presented
profiles, the saturation and the colour denote the amplitude and phase of the
dynamic, in-plane component of magnetization. The compact localized states
(CLS) are presented at the point $M^{\leftarrow}$ for the second band – right
column. The CLS do not occupy minority sublattice $A$. The inclusions $B$, in
which the magnetization dynamics is focused, are quite well isolated from each
other. One can easily notice that the lattice is decorated by loops (marked by
grey patches) where the phase of the precessing magnetization flips between
inclusions ($+$ and $-$ signs). Exactly at point $M$ – left column, we observe
the degeneracy of the second and third bands. The spin waves occupy $B$
inclusions only in one majority sublattice, i.e. along vertical or horizontal
lines, filliping the phase from inclusion to inclusion which gives the pattern
characteristic to noncontractible loop states (NLS) - marked by grey stripes.

To prove that the second band supports the CLS regardless of its finite width,
we plotted the profiles of spin wave eigenmodes at $M$ point and in its close
vicinity. The results are presented in Fig. 4. The profiles were shown for
infinite lattice and are presented in the form of square arrays containing 3x3
unit cells, where the dashed lines mark their edges. It is visible that the
spin waves are concentrated in the cylindrical inclusions, where the amplitude
and phase of precession is quite homogeneous. In calculations, we used the
Bloch boundary conditions applied for a single unit cell, which means that at
$M$ point the Bloch function is flipped after translation by lattice period,
in both principal directions of the lattice and we will not see the single
closed loops of CLS or lines of NLS. Exactly at $M$ point, all three bands
have zero group velocity. Therefore, the corresponding modes (left column) are
not propagating. The lowest band ($M_{1}$) occupy only inclusions $A$ from the
minority sublattice where the static demagnetizing field is slightly lower
than inside inclusions $B$ (see Supplementary Information D), which justifies
its lower frequency and lifting the degeneracy with two higher modes $M_{2}$
and $M_{3}$ of the same frequency. Each of the modes $M_{2}$ and $M_{3}$
occupy only one of two sublattices $B$, therefore they can be interpreted as
NLS. To observe the pattern typical for CLS, we need to move slightly away
from $M$ point. The first and third modes have then the linear dispersion with
high group velocity and the second band remains flat. We selected the point
$M^{\leftarrow}$ shifted from $M$ point toward $\Gamma$ point by 5% of
$M-\Gamma$ distance (right column). We can see that the first and third modes
$M^{\leftarrow}_{1}$, $M^{\leftarrow}_{3}$ occupy now all inclusions and the
mode $M^{\leftarrow}_{2}$ from the flat band has a profile typical for CLS,
predicted by tight-binding models [55, 8, 56, 9, 10]:

$\text{\textbar}m_{\textbf{k}}>=[\underbrace{-e^{i\frac{k_{y}}{2}}}_{B},\underbrace{0}_{A},\underbrace{e^{i\frac{k_{x}}{2}}}_{B}]$
(5)

where $m_{\textbf{k}}$ is the complex amplitude of the Bloch function in the
base of unit cell (i.e. on two inclusion $B$ from majority sublattices and one
inclusion $A$ from minority sublattice), $\textbf{k}=[k_{x},k_{y}]$ is
dimensionless wave vector. From (Eq. 5), we can see that (i) CLS do not occupy
the minority nodes $A$ and (ii) close to $M$ point the phases at two nodes
$B$, from different majority sublattices, are opposite. These two features are
reproduced for $M^{\leftarrow}_{2}$ mode in investigated magnonic Lieb
lattice. In the profile of this mode, we marked (by a grey patch) the
elementary loop of CLS which is easily identified in finite systems. Here, in
an infinite lattice with Bloch boundary conditions, the loops are infinitely
replicated with $\pi$ phase shift after each translation $x-$ and
$y-$direction. The localization at the inclusions $B$ and the absence of the
spin wave dynamics in inclusions $A$ is observed regardless of the wave
vector. Therefore, the coupling can take place only between the next
neighbours (inclusions $B$), i.e. on larger distances and mostly due to
dipolar interactions, that makes the second band not perfectly flat.

Figure 5: Dispersion relation for the extended magnonic Lieb lattice Lieb-5,
containing five inclusions in the unit cell: one inclusion $A$ from minority
sublattice and four inclusions $B$ from majority sublattices (see Fig. 2). (a)
The dispersion relation is plotted along the high symmetry path
$\Gamma$-X-M-$\Gamma$ (see the inset). The first, third and fifth bands (dark
blue, red and cyan) are strongly dispersive bands, while the second and fourth
bands (green and magenta) are less dispersive and related to the presence of
CLS. The system does not support the appearance of Dirac cones, even in case
when the interaction is fictitiously limited only to inclusions, according to
tight-binding model. (b) The zoomed regions in the vicinity of $\Gamma$ (in
dark green frame) and $M$ points show the essential gaps with relatively low,
parabolic-like curvatures for top and bottom bands.

Let’s discuss now the presence of flat bands and CLS in an extended magnonic
Lieb lattice (Lieb-5), containing five inclusions in the unit cell: one
inclusion $A$ form minority sublattice and four inclusions $B$ from majority
sublattices, as it is presented in Fig. 2. In the considered structure, we add
two additional inclusions $B$ into the unit cell in such a way that
neighbouring inclusions $A$ are linked by the doublets of inclusions $B$. The
sizes of inclusions, distances between them, the thickness of the layer and
the material composition of the structure remained the same as for the basic
Lieb lattice, discussed earlier (Fig. 1).

The dispersion relation obtained for the magnonic Lieb-5 lattice can be found
in Fig. 5(a). The properties of the extended Lieb lattices are well described
in the literature [10, 57, 58, 59]. The tight-binding model description of
Lieb-5 lattices, with information about their dispersion relation and the
profiles of the eigenmodes, are presented in numerous papers[10, 55, 16].
Therefore, it is possible to compare the obtained results with the theoretical
predictions of the tight-binding model.

Figure 6: The profiles obtained for the extended Lieb lattice consisted of 5
inclusions in the unit cell. The modes are presented for bands No. 3-5 in
$\Gamma$ point and its proximity $\Gamma^{\leftarrow}$ (the first and second
column). In the third and fourth columns, we presented the profiles for bands
No. 1-3 at $M$ point and its vicinity $M^{\leftarrow}$. Each profile of
eigenmode is presented on a grid composed of 3x3 unit cells - dashed lines
mark the edges of unit cells. The scheme of the unit cell is presented in top-
left corner. Exactly at $\Gamma$ (and $M$) point the bands No. 3 and 4 (No. 2
and 3) are degenerated and profiles: $\Gamma_{3}$ and $\Gamma_{4}$ ($M_{2}$
and $M_{3}$) have non-standard (for CLS) complementary form – i.e. their
combinations $\Gamma_{3}\pm i\Gamma_{4}$ ($M_{2}\pm iM_{3}$) gives NLS. To
obtain proper profiles of CLS, where the phase of procession flips around CLS
loop, we need to explore the vicinity of $\Gamma$ ($M$) point – see the grey
patches for the mode $\Gamma^{\leftarrow}_{4}$ ($M^{\leftarrow}_{2}$) with $+$
and $-$ signs.

The tight-binding model of Lieb-5 lattice predicts two flat bands with CLS:
the second (green) and fourth (magenta) band in the spectrum. The flat bands
in the tight-binding model are not intersected by Dirac cones but they are
degenerated at $\Gamma$ and $M$ point with the third band (red). These
features are reproduced in investigated magnonic Lieb-5 lattice (Fig. 2). The
dispersion relation for this system is presented in the Fig. 5(a). Also, we
have marked, with two rectangles (dark green and violet), the vicinities of
$\Gamma$ and $M$ points, where the flat bands (the fourth and second bands)
become degenerated with the third, dispersive band – Fig. 5(b). It is easy to
notice the essential frequency gaps ($\approx 33$ MHz and $\approx 84~{}$MHz
at $\Gamma$ and $M$ points, respectively), which qualitatively corresponds to
the prediction of the tight-binding model. It is worth noting that although
the low dispersion bands (the second and fourth band) are in general not
perfectly flat. Nevertheless, around the point $\Gamma$ and $M$ points the
bands are flattened and the $\Gamma-X$ and $X-M$ sections are very flat for
the fourth and second band, respectively.

The spin wave profiles of CLS at the high symmetry points: $\Gamma$ and $M$
are presented in Fig. 6. Exactly at $\Gamma$ and $M$ (the first and third
column), we can see the pairs of degenerated mods $\Gamma_{3}$, $\Gamma_{4}$
and $M_{2}$, $M_{3}$ which exhibit features of CLS predicated by the tight-
binding model (see the loops of sites on grey patches): (i) modes occupy only
the inclusions $B$ from majority sublattices, (ii) doublets of inclusions $B$
in the loops of CLS have opposite (the same) phases at $\Gamma$ ($M$) point.
The significant difference is that; once we switch one to another $B$-$B$
doublet, circulating the CLS loop the phase of precession charges by
$\pm\pi/2$ not by 0 or $\pi$. However, when we make combinations of
degenerated modes: $\Gamma_{3}\pm i\Gamma_{4}$ or $M_{2}\pm iM_{3}$, then we
obtain the NLS occupying the horizontal or vertical lines, where the
precession at exited $B$ inclusion will be in- or out-of-phase. The CLS modes
are clearly visible when we move slightly away from the high symmetry point
where the degeneracy occurs. In the proximity of $\Gamma$ and $M$ point, one
can see the CLS modes $\Gamma^{\leftarrow}_{4}$ and $M^{\leftarrow}_{2}$ for
which the phase of precession takes the relative values close to 0 or $\pi$.
The small discrepancies, are visible as a slight change in the colours
representing the phase, resulting from the fact that we are not exactly in
high symmetry points but shifted by 5% on the path $\Gamma-M$.

The extension of the presented analysis to magnonic Lieb - 7 lattice, where
the inclusions $A$ are liked by the chains composed of three inclusions $B$,
is presented in Supplementary Information A.

##  V  Conclusions

We proposed a possible realisation of the magnonic Lieb lattices where the
compact localized spin wave modes can be observed in flat bands. The presented
system qualitatively reproduces the spectral properties and the localization
features of the modes, predicted by the tight-binding model and observed for
photonic and electronic counterparts. The magnonic platform for the
experimental studies of Lieb lattices seems to be attractive due to the larger
flexibility in designing magnonic systems and the steering of its magnetic
configuration by external biases. The idea of the magnonic Lieb latices allows
considering many problems related to dynamics, localization and interactions
in flat band systems taking the advantage of the magnonic systems: presence
and possibility of tailoring of long-range interactions, intrinsic non-
linearity, etc.

##  VI  Acknowledgements

This work has received funding from National Science Centre Poland grants
UMO-2020/39/O/ST5/02110, UMO-2021/43/I/ST3/00550 and support from the Polish
National Agency for Academic Exchange grant BPN/PRE/2022/1/00014/U/00001.

## References

  * Davison and Steślicka [1996] S. Davison and M. Steślicka,  _Basic Theory of Surface States_ (Clarendon Press, 1996).
  * Abrahams [2010] E. Abrahams,  _50 Years of Anderson Localization_ (World Scientific, 2010).
  * Wannier [1962] G. H. Wannier, Dynamics of band electrons in electric and magnetic fields, Rev. Mod. Phys. 34, 645 (1962).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Sutherland [1986] B. Sutherland, Localization of electronic wave functions due to local topology, Phys. Rev. B 34, 5208 (1986).
  * Lieb [1989] E. H. Lieb, Two theorems on the Hubbard model, Phys. Rev. Lett. 62, 1201 (1989).
  * Ezawa [2020] M. Ezawa, Systematic construction of square-root topological insulators and superconductors, Phys. Rev. Res. 2, 033397 (2020).
  * Leykam  _et al._ [2018] D. Leykam, A. Andreanov, and S. Flach, Artificial flat band systems: from lattice models to experiments, ADV PHYS-X 3, 1473052 (2018).
  * Bhattacharya and Pal [2019] A. Bhattacharya and B. Pal, Flat bands and nontrivial topological properties in an extended Lieb lattice, Phys. Rev. B 100, 235145 (2019).
  * Zhang  _et al._ [2017] D. Zhang, Y. Zhang, H. Zhong, C. Li, Z. Zhang, Y. Zhang, and M. R. Belić, New edge-centered photonic square lattices with flat bands, Ann. Phys. 382, 160 (2017).
  * Flach  _et al._ [2014] S. Flach, D. Leykam, J. D. Bodyfelt, P. Matthies, and A. S. Desyatnikov, Detangling flat bands into fano lattices, EPL 105, 30001 (2014).
  * Aoki  _et al._ [1996] H. Aoki, M. Ando, and H. Matsumura, Hofstadter butterflies for flat bands, Phys. Rev. B 54, R17296 (1996).
  * Bergman  _et al._ [2008] D. L. Bergman, C. Wu, and L. Balents, Band touching from real-space topology in frustrated hopping models, Phys. Rev. B 78, 125104 (2008).
  * 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 and Yang [2019] J.-W. Rhim and B.-J. Yang, Classification of flat bands according to the band-crossing singularity of Bloch wave functions, Phys. Rev. B 99, 045107 (2019).
  * Rhim and Yang [2021] J.-W. Rhim and B.-J. Yang, Singular flat bands, ADV PHYS-X 6, 1901606 (2021).
  * Leykam  _et al._ [2017] D. Leykam, J. D. Bodyfelt, A. S. Desyatnikov, and S. Flach, Localization of weakly disordered flat band states, Eur Phys J B 90, 1 (2017).
  * Chalker  _et al._ [2010] J. T. Chalker, T. S. Pickles, and P. Shukla, Anderson localization in tight-binding models with flat bands, Phys. Rev. B 82, 104209 (2010).
  * Leykam  _et al._ [2013] D. Leykam, S. Flach, O. Bahat-Treidel, and A. S. Desyatnikov, Flat band states: Disorder and nonlinearity, Phys. Rev. B 88, 224203 (2013).
  * Häusler [2015] W. Häusler, Flat-band conductivity properties at long-range Coulomb interactions, Phys. Rev. B 91, 041102 (2015).
  * Kohmoto  _et al._ [1987] M. Kohmoto, B. Sutherland, and C. Tang, Critical wave functions and a cantor-set spectrum of a one-dimensional quasicrystal model, Phys. Rev. B 35, 1020 (1987).
  * Mieszczak  _et al._ [2022] S. Mieszczak, M. Krawczyk, and J. W. Kłos, Spin-wave localization on phasonic defects in a one-dimensional magnonic quasicrystal, Phys. Rev. B 106, 064430 (2022).
  * 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).
  * Drost  _et al._ [2017] R. Drost, T. Ojanen, A. Harju, and P. Liljeroth, Topological states in engineered atomic lattices, Nat. Phys. 13, 668 (2017).
  * Slot  _et al._ [2017] M. R. Slot, T. S. Gardenier, P. H. Jacobse, G. C. P. van Miert, S. N. Kempkes, S. J. M. Zevenhuizen, C. M. Smith, D. Vanmaekelbergh, and I. Swart, Experimental realization and characterization of an electronic Lieb lattice, Nat. Phys. 13, 672 (2017).
  * Shen  _et al._ [2010] R. Shen, L. B. Shao, B. Wang, and D. Y. Xing, Single Dirac cone with a flat band touching on line-centered-square optical lattices, Phys. Rev. B 81, 041410 (2010).
  * Taie  _et al._ [2015] S. Taie, H. Ozawa, T. Ichinose, T. Nishio, S. Nakajima, and Y. Takahashi, Coherent driving and freezing of bosonic matter wave in an optical Lieb lattice, Sci. Adv. 1, e1500854 (2015).
  * Swain and Karmakar [2020] N. Swain and M. Karmakar, Strain-induced superconductor-insulator transition on a Lieb lattice, Phys. Rev. Res. 2, 023136 (2020).
  * Xu  _et al._ [2021] F. Xu, L. Zhang, and L.-Y. Jiang, Temperature and doping dependent flat-band superconductivity on the Lieb-lattice, Chin. Phys. B 30, 067401 (2021).
  * Ma  _et al._ [2021] T.-X. Ma, Q.-S. Fan, C. Zhang, and Y.-S. Wang, Acoustic flatbands in phononic crystal defect lattices, J. Appl. Phys. 129, 145104 (2021).
  * 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).
  * Cracknell [1973] A. P. Cracknell, Van Hove singularities and zero-slope points in crystals, J. Phys. C: Solid State Phys. 6, 826 (1973).
  * Beugeling  _et al._ [2012] W. Beugeling, J. C. Everts, and C. Morais Smith, Topological phase transitions driven by next-nearest-neighbor hopping in two-dimensional lattices, Phys. Rev. B 86, 195129 (2012).
  * Jiang  _et al._ [2019a] W. Jiang, M. Kang, H. Huang, H. Xu, T. Low, and F. Liu, Topological band evolution between Lieb and kagome lattices, Phys. Rev. B 99, 125131 (2019a).
  * Beličev  _et al._ [2017] P. P. Beličev, G. Gligorić, A. Maluckov, M. Stepić, and M. Johansson, Localized gap modes in nonlinear dimerized Lieb lattices, Phys. Rev. A 96, 063838 (2017).
  * Ramachandran  _et al._ [2017] A. Ramachandran, A. Andreanov, and S. Flach, Chiral flat bands: Existence, engineering, and stability, Phys. Rev. B 96, 161104 (2017).
  * Jiang  _et al._ [2019b] W. Jiang, H. Huang, and F. Liu, A Lieb-like lattice in a covalent-organic framework and its Stoner ferromagnetism, Nat. Commun. 10, 2207 (2019b).
  * Poli  _et al._ [2017] C. Poli, H. Schomerus, M. Bellec, U. Kuhl, and F. Mortessagne, Partial chiral symmetry-breaking as a route to spectrally isolated topological defect states in two-dimensional artificial materials, 2D Mater. 4, 025008 (2017).
  * Khomeriki and Flach [2016] R. Khomeriki and S. Flach, Landau-Zener Bloch oscillations with perturbed flat bands, Phys. Rev. Lett. 116, 245301 (2016).
  * Chen  _et al._ [2014] L. Chen, T. Mazaheri, A. Seidel, and X. Tang, The impossibility of exactly flat non-trivial Chern bands in strictly local periodic tight binding models, J. Phys. A Math. Theor. 47, 152001 (2014).
  * Bergholtz and Liu [2013] E. J. Bergholtz and Z. Liu, Topological flat band models and fractional Chern insulators, Int J Mod Phys B 27, 1330017 (2013).
  * Vicencio Poblete [2021] R. A. Vicencio Poblete, Photonic flat band dynamics, ADV PHYS-X 6, 1878057 (2021).
  * Baba [2008] T. Baba, Slow light in photonic crystals, Nat. Photonics 2, 465 (2008).
  * Chumak  _et al._ [2022] A. V. Chumak _et al._ , Advances in magnetics roadmap on spin-wave computing, IEEE Trans. Magn. 58, 1 (2022).
  * Mahmoud  _et al._ [2020] A. Mahmoud, F. Ciubotaru, F. Vanderveken, A. V. Chumak, S. Hamdioui, C. Adelmann, and S. Cotofana, Introduction to spin wave computing, J. Appl. Phys. 128, 161101 (2020).
  * Tasaki [1998] H. Tasaki, From Nagaoka’s ferromagnetism to flat-band ferromagnetism and beyond: An introduction to ferromagnetism in the Hubbard model, Prog. Theor. Phys. 99, 489 (1998).
  * Cao  _et al._ [2015] X. Cao, K. Chen, and D. He, Magnon hall effect on the Lieb lattice, J. Phys. Condens. Matter 27, 166003 (2015).
  * Krawczyk and Grundler [2014] M. Krawczyk and D. Grundler, Review and prospects of magnonic crystals and devices with reprogrammable band structure, J. Phys. Condens. Matter 26, 123202 (2014).
  * Choudhury  _et al._ [2016] S. Choudhury, S. Saha, R. Mandal, S. Barman, Y. Otani, and A. Barman, Shape- and interface-induced control of spin dynamics of two-dimensional bicomponent magnonic crystals, ACS Appl. Mater. Interfaces 8, 18339 (2016).
  * Tacchi  _et al._ [2012] S. Tacchi, G. Duerr, J. W. Kłos, M. Madami, S. Neusser, G. Gubbiotti, G. Carlotti, M. Krawczyk, and D. Grundler, Forbidden band gaps in the spin-wave spectrum of a two-dimensional bicomponent magnonic crystal, Phys. Rev. Lett. 109, 137202 (2012).
  * Gubbiotti  _et al._ [2012] G. Gubbiotti, S. Tacchi, M. Madami, G. Carlotti, S. Jain, A. O. Adeyeye, and M. P. Kostylev, Collective spin waves in a bicomponent two-dimensional magnonic crystal, Appl. Phys. Lett. 100, 162407 (2012).
  * Böttcher  _et al._ [2022] T. Böttcher, M. Ruhwedel, K. O. Levchenko, Q. Wang, H. L. Chumak, M. A. Popov, I. V. Zavislyak, C. Dubs, O. Surzhenko, B. Hillebrands, A. V. Chumak, and P. Pirro, Fast long-wavelength exchange spin waves in partially compensated Ga:YIG, Appl. Phys. Lett. 120, 102401 (2022).
  * Gurevich and Melkov [1996] A. Gurevich and G. Melkov,  _Magnetization Oscillations and Waves_ (CRC Press., 1996).
  * Dechaumphai and Sucharitpwatskul [2019] P. Dechaumphai and S. Sucharitpwatskul,  _Finite Element Analysis with COMSOL_ (Alpha Science International, Limited, 2019).
  * Leykam  _et al._ [2012] D. Leykam, O. Bahat-Treidel, and A. S. Desyatnikov, Pseudospin and nonlinear conical diffraction in Lieb lattices, Phys. Rev. A 86, 031805 (2012).
  * Chen and Zhou [2016] R. Chen and B. Zhou, Finite size effects on the helical edge states on the Lieb lattice, Chin. Phys. B 25, 067204 (2016).
  * Dias and Gouveia [2015] R. G. Dias and J. D. Gouveia, Origami rules for the construction of localized eigenstates of the Hubbard model in decorated lattices, Sci. Rep. 5, 16852 (2015).
  * Liu  _et al._ [2020] J. Liu, X. Mao, J. Zhong, and R. A. Römer, Localization, phases, and transitions in three-dimensional extended Lieb lattices, Phys. Rev. B 102, 174207 (2020).
  * Mao  _et al._ [2020] X. Mao, J. Liu, J. Zhong, and R. A. Römer, Disorder effects in the two-dimensional Lieb lattice and its extensions, Physica E Low Dimens. Syst. Nanostruct. 124, 114340 (2020).
  * Graczyk and Krawczyk [2017] P. Graczyk and M. Krawczyk, Coupled-mode theory for the interaction between acoustic waves and spin waves in magnonic-phononic crystals: Propagating magnetoelastic waves, Phys. Rev. B 96, 024407 (2017).
  * Gallardo  _et al._ [2014] R. A. Gallardo, A. Banholzer, K. Wagner, M. Körner, K. Lenz, M. Farle, J. Lindner, J. Fassbender, and P. Landeros, Splitting of spin-wave modes in thin films with arrays of periodic perturbations: theory and experiment, New J. Phys. 16, 023015 (2014).

##  Supplementary Information

###  VI.1 Double extended Lieb lattice

We can generate further extensions of the magonic Lieb lattice by adding more
inclusions $B$, i.e. by introducing additional majority sublattices. We
considered here a doubly extended Lieb lattice (Lieb-7) to check to what
extent the magnonic system corresponds to the tight-binding model. The
mentioned lattice consists of seven nodes; six belong to majority sublattices
$B$ and one belongs to minority sublattice $A$ (Fig. 7). The magnetic
parameters were kept as for basic and Lieb-5 lattices, considered in the
manuscript. The geometrical parameters have changed only as a result of the
introduction of additional inclusions $B$. Therefore, the unit cell has
increased to the size 500x500 nm.

Figure 7: Doubly extended magnonic Lieb lattice: Lieb-7. Dimensions of the
ferromagnetic unit cell are equal to 500x500x59 nm. The unit cell contains
seven inclusions of 50 nm diameter. (a) The structure of extended Lieb
lattice, and (b) top view on Lieb-7 lattice unit cell where the node
(inclusion) from minority sublattice $A$ and two nodes (inclusions) from two
majority sublattices $B$ are marked.

In the case of a doubly extended Lieb lattice (Lieb-7), we expect (according
to the works [10, 59]) to obtain seven bands in the dispersion relation. The
tight-binding model predicts that the bands will be symmetric with respect to
the fourth band, exhibiting particle-hole symmetry. However, due to the
dipolar interaction, we did not expect such symmetry. Another feature which
one may deduce from the tight-binding model is that bands No. 2, 4 and 6
should be flat while bands No. 1, 3, 5 and 7 are considered dispersive.
Moreover, bands No. 3 and 5 suppose to form a Dirac cone intersecting flat
band No. 4 at the $\Gamma$ point.

We calculated the dispersion relation for magnonic Lieb-7 lattice (Fig. 8(a)),
which share many properties with those characteristic for the tight-binding
model [59]: (i) third and fifth bands form the Dirac cones which almost
intersect the flatter forth band at $\Gamma$ point; (ii) the third (and fifth)
band has a parabolic shape at $M$ point where it is degenerated with the
second (and six) band which is weakly dispersive. The mentioned regions of
dispersion are presented as 3D plots in Fig. 8(b). Also, we are going to
discuss shortly the profiles of spin wave eigenmodes (including CLS) in these
two regions of the dispersion relation, which are presented in Fig. 9.

Figure 8: Dispersion relation for the double extended magnonic Lieb lattice
(Lieb-7) containing seven inclusions in the unit cell: one inclusion $A$ from
minority sublattice and six inclusions $B$ from majority sublattices (see Fig.
2). (a) The dispersion relation is plotted along the high symmetry path
$\Gamma$-X-M-$\Gamma$ (see the inset). The first, third, fifth and seven bands
(dark blue, red, cyan and orange) are dispersive, while the second, fourth and
sixth bands (green, magenta and dark green bands) are the flatter bands,
supporting the magnonic CLS. Dirac cones occur at the $\Gamma$ point and
almost interact with the flatter fourth band, while at $M$ point, we observe
the degeneracy of the dispersive parabolic third (fifth) band with a flatter
second (six) band. (b) The zoomed vicinity of $\Gamma$ point (dark green
frame) and $M$ point (violet frame) regions are presented in 3D.  Figure 9:
The profiles of eigenmodes were obtained for magnonic Lieb-7. The modes are
presented for bands No. 3-5 at $\Gamma$ point and 5-7 at $M$ point. The modes
denoted as $\Gamma_{3}$ and $\Gamma_{4}$ are degenerated whereas the
$\Gamma_{5}$ is separated from them by extremely small gap $\approx~{}2$ MHz.
At $M$ point, we showed the profiles for bands No. 5, 6 and 7. The modes
$M_{5}$ and $M_{6}$ are degenerated and separated from $M_{7}$ by essential
gap – predicted by the tight-binding model.

Dirac cones appear at the $\Gamma$ point for bands No. 3 and 5. At this point,
as for the basic magnonic Lieb lattice (Fig. 3), there is a very narrow gap of
the width $\approx~{}2$ MHz. The profiles $\Gamma_{4}$ and $\Gamma_{5}$ (left
column in Fig. 9) represent the degenerated states originating from flat and
dispersive bands. Both of them do not occupy the inclusions $A$ and are more
focused on two inclusions $B$ arranged in horizontal ($\Gamma_{4}$) and
vertical lines ($\Gamma_{5}$) – see grey stripes. Therefore, their profiles
are similar to NLS, where the first and third inclusion $B$ in each three-
element chain, linking inclusions $A$, precesses out-of-phase and the second
(central) inclusion $B$ remains unoccupied.

At the $M$ point, the $M_{5}$ and $M_{6}$ bands are degenerated. For these
bands, the spin waves are localized in all inclusions $B$ and do not occupy
inclusions $A$ (see right column of Fig. 6) The first and third inclusion $B$
in each three-element chain, linking inclusions $A$, precess in-phase, whereas
the second (central) inclusion $B$ precesses out-of-phase with respect to the
first and third one. This pattern of occupation of inclusions and the phase
relations between them is similar to one observed for CLS (see grey patches
marking the loops of inclusions in the left column of Fig. 6), but has one
significant difference. The phase difference between successive three-element
chains of inclusions $B$, in the loop, is equal to $\pm\pi/2$. However, the
linear combination of the modes $M_{5}\pm iM_{6}$ produces, similarly to the
case of the Lieb-5 lattice, the NLS. To observe the proper profiles of CLS or
NLS, we need to shift slightly from the high symmetry points $\Gamma$ and $M$
to cancel the degeneracy.

###  VI.2  Realization of Lieb lattice
by shaping demagnetizing field

We have considered also an alternative realisation method for a magnonic Lieb
lattice in a ferromagnetic layer. This approach is based on shaping the
internal demagnetizing field. The structure under consideration is presented
in Fig. 10. It consists of a thin (28.5 nm) and infinite CoFeB layer on which
a Py antidot lattice (ADL), of 28.5 nm thickness, is deposited. The
cylindrical holes in ADL are arranged in shape of the basic Lieb lattice. The
size of the unit cell and diameter of holes remains the same as for the basic
Lieb lattice proposed in the main part of the manuscript (see Fig. 1). Due to
the absence of perpendicular magnetic anisotropy (PMA), we decided to apply a
much larger external magnetic field ($H_{0}=1500$ mT) to saturate the
ferromagnetic material in an out-of-plane direction.

Figure 10: Basic magnonic Lieb lattice where spin wave excitations in the
CoFeB layer are shaped by demagnetizing field from Py antidot lattice.
Dimensions of the ferromagnetic unit cell are equal to 250x250x59 nm and
contain 3 inclusions of 50 nm diameter. (a) structure of basic Lieb lattice,
(b) top view on basic Lieb lattice unit cell and differentiation to nodes of
sublattice $A$ and $B$.

We assumed the same gyromagentic ratio for both materials
$\gamma=187~{}{\rm\frac{GHz}{T}}$, the following values of material parameters
for CoFeB [60]: saturation magnetization - $M_{S}=1150~{}{\rm\frac{kA}{m}}$,
exchange stiffness constant - $A=15~{}{\rm\frac{pJ}{m}}$. For Py, we used
material parameters [61]: saturation magnetization -
$M_{S}=796~{}{\rm\frac{kA}{m}}$, exchange stiffiness constant -
$A=13~{}{\rm\frac{pJ}{m}}$.

The deposition of the ADL made of Py (material of lower $M_{S}$) above the
CoFeB layer (material of higher $M_{S}$) is critical for spin wave
localization in CoFeB below the exposed parts (holes) of the ADL. The
demagnetization field produced on CoFeB/Air interface creates wells partially
confining the spin waves. However, this pattern of internal demagnetizing
field becomes smoother with increasing distance from the ADL.

Figure 11: The dispersion relation obtained for basic Lieb lattice formed by
demagnetizing field of antidot lattice (see Fig. 10). (a) The dispersion
relation, (b) the 3D plot of dispersion relation in the region marked with the
green frame in (a). Results were obtained for $H_{0}=1500$ mT applied out-of-
plane.

The obtained dispersion relation is shown in Fig. 11. It is worth noting that
the lowest band is very dispersive, while the highest band is flattened more
than in the case of the structure presented in the main part of the manuscript
(see Fig. 3). The middle band, which suppose to support CLS, varies in extent
similar to the third band. For this structure, Dirac cones in the $M$ point
cannot be clearly unidentified.

###  VI.3 Lieb lattice formed by YIG inclusions
in non-magnetic matrix

The periodic arrangement of ferromagnetic cylinders surrounded by nonmagnetic
material (e.g. air) seems to be the simplest realization of the Lieb lattice.
To refer this structure to the bi-component system investigated in the main
part of the manuscript, we assumed the same material and geometrical
parameters for inclusions as for the structure presented in Fig. 1.

The advantage of this system is that the confinement of spin waves within the
areas of inclusions is ensured for arbitrarily high frequency. We are not
limited here by the FMR frequency of the matrix, as it was for bi-component
Lieb lattices (Figs. 1, 2). However, the coupling of magnetization dynamics
between the inclusions is here provided solely by the dynamical demagnetizing
field, i.e. the evanescent spin waves do not participate in the coupling.
Therefore, the interaction between inclusions is much smaller in general,
which leads to a significant narrowing of all magnonic bands (Fig. 11). The
widths of the second and third band can be even smaller than the gap
separating from the first bands – Fig. 11(b). Such strong modification of the
spectrum makes the applicability of the considered system for the studies of
magnonic CLS questionable.

Figure 12: Dispersion relations for basic Lieb lattice. (a) The results
obtained for YIG inclusions in Ga:YIG matrix (dashed lines) and YIG inclusions
without matrix (solid lines). (b) The zoomed dispersion relation obtained for
YIG inclusions without matrix, marked in (a) by the frame.

###  VI.4 Demagnetizing field
in YIG—Ga:YIG Lieb lattice

The difficulty in designing the magnonic system is not only due to the
adjustment of geometrical parameters of the system but also due to the shaping
of the internal magnetic field $\textbf{H}_{\rm eff}$.

The components of the effective magnetic field can be divided into long-range
and short-range. The realization of our model is inseparably linked to the
long-range dipole interactions through which the coupling between inclusions
is possible. This kind of interaction is sensitive to the geometry of the
ferromagnetic elements forming the magnonic system.

In Lieb lattice, the nodes of minority sublattice $A$ have four neighbours and
the nodes of majority sublattice $B$ have two. As a result, identical
inclusions (in terms of their shapes and material parameters) become
distinguishable, because of slightly different values of the internal
demagnetising field. This has consequences for the formation of a frequency
gap between Dirac cones at point $M$ in the dispersion relation obtained for
the basic Lieb lattice. In the literature, this phenomenon has been described
for the tight-binding model and is called node dimerisation of the lattice
[37].

In Fig. 13 we have shown the profile of the $z$-component of the demagnetising
field. For each inclusion through which the cut line passes, we have marked
the minimum value of the demagnetising field. The slightly lower value of
internal filed for inclusions $A$ is responsible for a tiny lowering of the
frequency for the mode $M_{1}$ (concentrated in inclusions $A$) respect the
degenerated modes $M_{2}$ and $M_{3}$ (confined in inclusions $B$).

Figure 13: Profile of static demagnetizing field plotted at cut through (a)
Lieb lattice unit cell. (b) The $z$-component of the demagnetizing field along
the cut line is shown in (a). In the plot, we have marked peaks for the areas
of inclusions $A$ and $B$. Please note the slightly different values of
demagnetizing in the centre of $A$ and $B$ inclusion due to different the
number neighboring of nodes: four for inclusion $A$, two for inclusion $B$.