Current and permanent address: ]Seagate Technology, 7801 Computer Ave.,
Bloomington, MN 55435

# Gænice: a general model for magnon band structure of artificial spin ices

Ghanem Alatteili  Victoria Martinez  Alison Roxburgh  Center for Magnetism and
Magnetic Nanostructures, University of Colorado Colorado Springs, Colorado
Springs, CO 80918, USA  Jack C. Gartside  Blackett Laboratory, Imperial
College London, Prince Consort Road, London SW7 2AZ, UK  Olle G. Heinonen  [
Materials Science Division, Argonne National Laboratory, Lemont, Illinois
60439, USA  Sebastian Gliga  Swiss Light Source, Paul Scherrer Institute, 5232
Villigen PSI, Switzerland  Ezio Iacocca<EMAIL_ADDRESS>Center for
Magnetism and Magnetic Nanostructures, University of Colorado Colorado
Springs, Colorado Springs, CO 80918, USA

###### Abstract

Arrays of artificial spin ices exhibit reconfigurable ferromagnetic resonance
frequencies that can be leveraged and designed for potential applications.
However, analytical and numerical studies of the frequency response of
artificial spin ices have remained somewhat limited due to the need of take
into account nonlocal dipole fields in theoretical calculations or by long
computation times in micromagnetic simulations. Here, we introduce Gænice, a
framework to compute magnon dispersion relations of arbitrary artificial spin
ice configurations. Gænice makes use of a tight-binding approach to compute
the magnon bands. It also provides the user complete control of the
interaction terms included, e.g., external field, anisotropy, exchange, and
dipole, making it useful also to compute ferromagnetic resonances for a
variety of structures, such as multilayers and ensembles of weakly or non-
interacting nanoparticles. Because it relies on a semi-analytical model,
Gænice is computationally inexpensive and efficient, making it an attractive
tool for the exploration of large parameter spaces.

††preprint: AIP/123-QED

##  I Introduction

Artificial spin ices (ASIs) are systems of structured nanomagnets arranged in
periodic patterns that are magneto-statically coupled. ASIs were originally
designed to mimic the behavior of natural spin ice materials [1] ,in order to
explore the fundamental principles of frustrated magnetism. Frustration arises
from competing magnetic interactions that cannot all be simultaneously
minimized [2], leading to highly degenerate states. ASIs can be also
considered as magnonic crystals [3, 4, 5] exhibiting reconfigurable magnonic
modes [6, 7, 8, 9, 5, 10, 11], nonlinear scattering [12], band structure [13,
14, 15], and hybrid modes [16, 17].

The arrangement of magnetic elements in a square lattice, known as square ice
[18],has been a test bed for the investigation of ASIs as magnonic crystals
because its relative simplicity allows for the understanding of the
fundamental physical phenomena. Analytically, square ices have proven
promising for reconfigurable magnonics because of the mode-dependent magnon
modes predicted [13, 14] as well as evidence of topological modes [19].
However, the study of similar effects in other geometries remains limited to
date. Experimentally, this is partly because of the large number of geometries
to explore [1] and the technical challenges to investigate wavevector- or
spatially-resolved magnons in a nanopatterned structure by, e.g. by Brillouin
light scattering [20]. From a numerical point of view, simulations using
micromagnetic modeling [21] are very time-consuming and often require large
memory allocations to investigate long-wavelength magnons that are easily
excited by microwave antennas in experiments. This means that more exotic
geometries, both in 2D [1] and 3D [22, 23, 24], have been slower to
materialize due to the lack of an efficient predictive tool for magnetization
dynamics.

Here, we present Gænice, a general formalism to compute the magnon dispersion
relation for arbitrary ASI geometries. The formalism is based on a Holstein-
Primakoff transformation [25] to obtain an eigenvalue problem that can be
solved numerically with little computational cost [13]. The main difference
between Gænice and other analytical methods is its generalization to arbitrary
nanomagnet orientations and magnetization states with an automatic
determination of the first Brillouin zone (FBZ). We expect that Gænice can
serve as a numerically efficient and computationally accurate tool to predict
magnonic functionality for ASIs and to direct more detailed studies of
promising geometries using micromagnetic simulations and experiments. In other
words, we envision Gænice as a tool to quickly explore the parameter space of
distinct ASI geometries and identify potentially interesting regimes that can
be then further explored with traditional computational and experimental
methods [26, 27, 28, 29, 30]

The remainder of the paper is organized as follows. In section II, we describe
the general formulation of the problem. The energy terms considered and their
implementation are detailed in section III. In section IV, we demonstrate the
functionality of Gænice by computing simple Kittel modes and ferromagnetic
resonance (FMR) modes in a variety of linear arrays of nanomagnets and the two
fundamental ASI configurations: square ice and Kagome ice.

##  II Generalized analytical model

We begin our description from the conservative Larmor torque equation

$\frac{\partial\mathbf{m}}{\partial
t}=-\gamma\mu_{0}\mathbf{m}\times\mathbf{H}_{\mathrm{eff}},$ (1)

where $\mathbf{m}$ is the normalized magnetization vector with
$|\mathbf{m}|=1$, $\gamma$ is the gyromagnetic ratio, and $\mu_{0}$ is the
vacuum permeability. The effective field $\mathbf{H}_{\mathrm{eff}}$ contains
physical terms and phenomena relevant to the magnetic material and interfaces
which are described within the context of our eigenvalue solver in section
III. Note that we neglect damping here given that we are interested in
resonant, propagating modes.

For small-amplitude excitations, such as magnons, the Larmor torque equation
can be rewritten as a Hamiltonian set of equations using a Holstein-Primakoff
transformation of the _complex_ small amplitude $a$ [25]

$a=\frac{m_{1}+im_{2}}{\sqrt{2(1+m_{3}^{2})}},$ (2)

where
$\mathbf{m}=m_{1}\hat{\mathbf{e}}_{1}+m_{2}\hat{\mathbf{e}}_{2}+m_{3}\hat{\mathbf{e}}_{3}$
is the magnetization vector expressed in a coordinate system where
$\hat{\mathbf{e}}_{3}$ defines the equilibrium orientation of the
magnetization vector and
$\hat{\mathbf{e}}_{1}\times\hat{\mathbf{e}}_{2}=\hat{\mathbf{e}}_{3}$. An
illustration of this basis is shown in Fig. 1(a).

From Eq. 2, we can relate the complex amplitudes to the magnetization vector
in the $\hat{\mathbf{e}}$ basis as

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle\sqrt{1-|a|^{2}}(a+a^{*})\approx\left(1-\frac{|a|^{2}}{2}\right)(a+a^{*})$
(3a) $\displaystyle m_{2}$ $\displaystyle=$ $\displaystyle
i\sqrt{1-|a|^{2}}(a-a^{*})\approx i\left(1-\frac{|a|^{2}}{2}\right)(a-a^{*})$
(3b) $\displaystyle m_{3}$ $\displaystyle=$ $\displaystyle(1-2|a|^{2})$ (3c)

Using the transformation of Eq. 2, we can approximately rewrite Eq. 1 as a
Hamiltonian system for the complex amplitude $a$

$\frac{\partial a}{\partial t}=-i\frac{\partial}{\partial
a^{*}}a\mathcal{H}a^{\dagger}$ (4)

and the Hamiltonian is defined over a magnetic volume as

$\mathcal{H}=-\mu_{0}M_{s}\int\mathbf{H}(\mathbf{m})\cdot{\mathbf{m}}dA,$ (5)

To describe the magnon band structure of an ensemble of nanomagnets in an ASI,
the auto-oscillator model can be generalized to an array of complex
amplitudes, as shown in Ref. 13. To account for the bending of the
magnetization at edges of the nanomagnets [31], we divide each nanomagnet into
3 _macrospins_. This is an important assumption in our model, making it valid
for magnetic elements with sizes of the order of hundreds of nanometers.
Therefore, given $N$ nanomagnets in the unit cell of the ASI, we define the
complex amplitude array $\underline{a}=[a_{1}~{}a_{2}~{}...~{}a_{3N}]$ and the
$2(3N)\times 2(3N)$ Hamiltonian matrix $\underline{\mathcal{H}}$ so that the
generalized Hamiltonian becomes

$\frac{d}{dt}\underline{a}=-i\frac{d}{d\underline{a}^{*}}\begin{bmatrix}\underline{a}&\underline{a}^{*}\end{bmatrix}\underline{\mathcal{H}}\begin{bmatrix}\underline{a}\\\
\underline{a}^{*}\end{bmatrix},$ (6)

The Hamiltonian matrix is further divided as

$\mathcal{H}=\begin{bmatrix}\mathcal{H}^{(1,2)}&\mathcal{H}^{(2,2)}\\\
\mathcal{H}^{(1,1)}&\mathcal{H}^{(2,1)}\end{bmatrix},$ (7)

where $\mathcal{H}^{(1,1)}=(\mathcal{H}^{(2,2)})^{*}$ and
$\mathcal{H}^{(1,2)}=(\mathcal{H}^{(2,1)})^{*}$ by symmetry of the Hamiltonian
equations. As further discussed below, this system describes bosonic
excitation (magnons) so that the eigenvalue problem can be solved using
Colpa’s grand dynamical matrix that ensures complex conjugate eigenvalues
[32].

Figure 1:  (a) The magnetization vector rotated by the polar and azimuthal
angles $(\theta_{m},\varphi_{m})$. The rotated frame
$(\hat{\mathbf{e}}_{1},\hat{\mathbf{e}}_{2},\hat{\mathbf{e}}_{3})$ is defined
relative to the equilibrium orientation of the magnetization vector
$\mathbf{m}$ (b) Nanomagnet orientation relative to the Cartesian coordinate
system. The unit vector $\hat{\mathbf{D}}$ indicates the direction of the
nanomagnet’s long axis and is defined by the angles $\theta_{d}$ and
$\varphi_{d}$.

###  II.1 Coordinate system

Our framework relies on a Cartesian coordinate system where the polar angle
$\theta=0$ and the azimuth angle $\varphi=0$ define the $z$-axis, as shown in
Fig. 1(a). Once this coordinate system is established, applying Eq. 1 requires
a rotation to the coordinate system defined by
$(\hat{\mathbf{e}}_{1},\hat{\mathbf{e}}_{2},\hat{\mathbf{e}}_{3})$, where the
direction of $\hat{\mathbf{e}}_{3}$ is parallel to the equilibrium orientation
of the magnetization vector at any given point in space.

This implies that coordinate transformations must be performed locally for
both the magnetization vector and the effective field for an arbitrary array
of magnetization vectors. We define the local rotation matrix

$R(\theta,\varphi)=\begin{bmatrix}\cos{\varphi}\cos\theta&\sin\varphi\cos\theta&-\sin\theta\\\
-\sin\varphi&\cos\varphi&0\\\
\cos\varphi\sin\theta&\sin\varphi\sin\theta&\cos\theta\end{bmatrix}.$ (8)

It is important to note that $R^{-1}(\theta,\varphi)=R(\theta,\varphi)$.

An arbitrary orientation of nanomagnets is considered, defined by a unit
vector $\hat{\mathbf{D}}$. The direction of $\hat{\mathbf{D}}$ is
parameterized by the polar and azimuth angles $\theta_{d}$ and $\varphi_{d}$
as shown in Fig. 1(b). When $\theta_{d}=0$ and $\varphi_{d}=0$, the nanomagnet
is aligned along the $z$-axis, and its thickness is aligned along the
$x$-axis. In this case, the angles $\theta_{d}$ and $\varphi_{d}$ represent
pitch and yaw, respectively.

###  II.2 Eigenvalue problem

The magnon dispersion relation $\omega(k)$ is obtained from Eq. 6 by invoking
Bloch’s theorem $a\rightarrow ae^{i\omega}$ and Colpa’s grand dynamical matrix
[32]

$\omega\underline{\Psi}\propto\begin{bmatrix}\mathcal{H}^{(1,2)}&-(\mathcal{H}^{(1,1)})^{*}\\\
\mathcal{H}^{(1,1)}&-(\mathcal{H}^{(1,2)})^{*}\end{bmatrix}\underline{\Psi},$
(9)

where $\underline{\Psi}$ is an array of eigenvectors. By introducing the
proportionality factor $\gamma/(2VM_{s})$, where $V$ is a volume and $M_{s}$
is the saturation magnetization, Eq. 9 becomes an equality. The volume and
saturation magnetization will be associated with a magnetization vector to
allow for maximal flexibility of the model. Therefore, we will use the
convention $\Omega=\gamma/(2VM_{s})\mathcal{H}$ to obtain

$\omega\underline{\Psi}=\begin{bmatrix}\underline{\Omega}^{(1,2)}&-(\underline{\Omega}^{(1,1)})^{*}\\\
\underline{\Omega}^{(1,1)}&-(\underline{\Omega}^{(1,2)})^{*}\end{bmatrix}\underline{\Psi},$
(10)

The eigenvalue problem of Eq. 10 can be solved numerically by standard
methods. This formulation has been implemented for square ice [13, 19] and it
has been referred to as semi-analytical because the Hamiltonian matrix
$\mathcal{H}$ is derived analytically and only the eigenvalue problem is
solved numerically. Here, we derive the Hamiltonian matrices for arbitrary ASI
configurations, so that the matrix is built in an automated way for any number
of nanomagnets in a unit cell.

##  III Effective field

The effective field is Gænice’s core, which gives rise to the Hamiltonian
matrix. In essence, the effective field is divided into two groups of physical
effects; local and non-local

$\mathbf{H}_{\mathrm{eff}}=\mathbf{H}_{\mathrm{l}}+\mathbf{H}_{\mathrm{nl}}.$
(11)

In its current implementation, Gænice includes a uniform external magnetic
field and anisotropy field as local contributions; and exchange interaction
and dipole-dipole interaction as non-local fields from a point of view that
the macrospins of the nanomagnets are coupled. In other words, these fields
lead to finite non-diagonal elements in the Hamiltonian block matrices. The
dipole-dipole contribution is fundamental to ASIs and is detailed below. In
the same manner, it is possible to extend Gænice to include other field
contributions, such as magnetocrystalline anisotropy, Dzyaloshinskii-Moriya
interaction, and RKKY exchange.

In the following subsections, we express these field contributions in the form
given in Eq. 10. All fields are defined in the Cartesian coordinate system and
rotated to the local magnetization coordinates described in section II.1.

###  III.1 External Field

A uniform external field, $\mathbf{H}_{0}$, leads to the energy

$E_{\mathrm{0}}=-V\mu_{0}M_{s}(R(\theta_{m},\varphi_{m})\cdot\mathbf{H}_{0})^{T}\cdot\mathbf{m}.$
(12)

The only quadratic term in $a$ in Eq. 12 is parallel to
$\hat{\mathbf{e}}_{3}$. Therefore, we can express the frequency contribution
due to an external field as

$\Omega_{0}=\gamma\mu_{0}|a|^{2}(R(\theta_{m},\varphi_{m})\cdot\mathbf{H}_{0})^{T}\cdot\hat{\mathbf{e}}_{3},$
(13)

with diagonal and off-diagonal blocks

$\displaystyle\underline{\Omega}_{0}^{(1,1)}$ $\displaystyle=$ $\displaystyle
0$ (14a) $\displaystyle\underline{\Omega}_{0}^{(1,2)}$ $\displaystyle=$
$\displaystyle\frac{\gamma\mu_{0}}{2}\left[(R(\theta_{m},\varphi_{m})\cdot\mathbf{H}_{0})^{T}\cdot\hat{\mathbf{e}}_{3}\right]\mathbf{I},$
(14b)

where $\mathbf{I}$ is the identity matrix.

###  III.2 Demagnetization Field

The demagnetization (demag) field is determined from the shape of the magnetic
element, which is an accurate approximation for soft magnets, such as
Permalloy. We consider a demag tensor $\underline{D}$ that we approximate with
diagonal demagnetizing factors $D_{1}<D_{2}<D_{3}$ such that

$\underline{D}=\begin{bmatrix}D_{3}&0&0\\\ 0&D_{2}&0\\\
0&0&D_{1}\end{bmatrix}.$ (15)

This approximation is consistent with the notion of macrospin elements, i.e.,
similar to the general ellipsoid [33].

This definition follows from the nanomagnet’s orientation in the Cartesian
coordinate system whereby the easy axis lies along the $z$-axis and the hard
axis along the $x$-axis. The demagnetizing factors can be found in a variety
of ways. Analytical expressions are available for oblate nanomagnets [33] and
for rectangular prisms [34]. Demag factors can also be obtained by fitting FMR
from simulated nanomagnets and empirical expressions can be found for a range
of aspect ratios [35] Gænice currently supports the oblate nanomagnet
analytical expressions given by.

$\displaystyle D_{1}$ $\displaystyle=$
$\displaystyle\frac{t\sqrt{1-e^{2}}(K-E)}{le^{2}},$ (16a) $\displaystyle
D_{2}$ $\displaystyle=$
$\displaystyle\frac{t\left(E-(1-e^{2})K\right)}{le^{2}\sqrt{1-e^{2}}},$ (16b)
$\displaystyle D_{3}$ $\displaystyle=$
$\displaystyle\frac{1-tE}{l\sqrt{1-e^{2}}},$ (16c)

where $K$ and $E$ are the complete elliptic integrals of the first and second
kind and $e=\sqrt{1-(w/l)^{2}}$. In the limit of a circular nanomagnet, one
must consider $D_{1}=D_{2}=0$ and $D_{3}=1$ to avoid a numerical singularity.

The anisotropy energy is expressed as

$E_{\mathrm{an}}=V\mu_{0}M_{s}^{2}\vec{m}\cdot\underline{D}\cdot\vec{m}^{T},$
(17)

We rotate the demagnetization tensor using Eq. 8 and the director vector
$\hat{\mathbf{D}}$ to align it with the magnetization direction in the
Cartesian reference frame

$\underline{C}=R(\theta_{m}-\theta_{d},\varphi_{m}-\varphi_{d})\cdot\underline{D}\cdot
R(\theta_{m}-\theta_{d},\varphi_{m}-\varphi_{d})^{-1},$ (18)

resulting in the nanomagnet-dependent Hamiltonian matrix

$\mathbf{\mathcal{H}}_{\mathrm{an}}=VM_{s}^{2}\vec{m}\cdot\underline{C}\cdot\vec{m}^{T}.$
(19)

We note that this matrix is a $3\times 3$ block that is defined for each
nanomagnet. Expressing Eq. 19 as a function of the complex amplitudes $a$,
ultimately results in the diagonal Hamiltonian block matrices

$\displaystyle\underline{\Omega}_{an}^{(1,1)}$ $\displaystyle=$
$\displaystyle\frac{\gamma\mu_{0}M_{s}}{2}\left[C_{11}-C_{22}+i(C_{12}+C_{21})\right]\mathbf{I},$
(20a) $\displaystyle\underline{\Omega}_{an}^{(1,2)}$ $\displaystyle=$
$\displaystyle\frac{\gamma\mu_{0}M_{s}}{2}\left[C_{11}-C_{22}-2C_{33}\right]\mathbf{I},$
(20b)

where the factors $C_{ij}$ are the coefficients of the matrix $\underline{C}$.

###  III.3 Exchange

We include exchange interaction as a minimal model for edge bending in the
magnetization of tightly packed nanomagnets [31, 13]. The nanomagnet is split
into three regions, and we use an effective exchange energy to parameterize
the exchange interaction. The nanomagnet splitting is shown in Fig. 2 for the
cases where the edge modes are (a) larger or (b) smaller than the stadium’s
semi-circular edges. It is assumed that the edge modes are symmetric in
volume. We refer to Appendix A for details. Here, we report the final form of
the block matrices used in the eigenvalue problem.

We consider that the nanomagnet is split in a bulk macrospin with volume
$V_{b}$ and two edge macrospins with identical volumes $V_{e}$, satisfying
$V=V_{b}+2V_{e}$. The volumes are uniquely determined by the parameter $\Delta
l$ defined as the length from the geometric center of the nanomagnet to the
center of the edge volume. The default value $\Delta l=(2l-w)/4$ is defined
when the edge volume is exactly contained at the semi-circular edges of
stadium-shaped nanomagnets. However, this parameter can be tuned.

The exchange energy is therefore defined as a pair-wise interaction between a
bulk macrospin $b$ and an edge macrospin $e+$ (upper) and $e-$ (lower)

$E_{\mathrm{ex}}^{(b,e)}=-\frac{J}{2}\mathbf{m}_{b}^{T}\cdot
R(\theta_{m_{e}}-\theta_{m_{b}},\varphi_{m_{e}}-\varphi_{m_{b}})\cdot\mathbf{m}_{e},$
(21)

where the exchange factor $J$ is given by

$J=\frac{2}{\Delta l^{2}}\left(V_{e}+\frac{V_{b}}{2}\right).$ (22)

Figure 2:  Symmetric splitting of stadium-shaped nanomagnets, where we discern
between a large edge volume (a) and a small edge volume (b), relative to the
semi-circular edges. The bulk and edge volumes, $V_{b}$ and $V_{e}$,
respectively, are uniquely determined by the parameter $\Delta l$.

Recognizing that the exchange interaction only occurs for neighboring
macrospins within a nanomagnet, we define the $3\times 3$ exchange energy
blocks per nanomagnet $N$

$\displaystyle\underline{E}_{ex,N}^{(1,1)}$ $\displaystyle=$
$\displaystyle-\frac{J}{2}\begin{bmatrix}0&E_{1}^{(b,e^{+})}/V_{e}&0\\\
E_{1}^{(e^{+},b)}/V_{b}&0&E_{1}^{(b,e^{-})}/V_{b}\\\
0&E_{1}^{(e^{-},b)}/V_{e}&0\end{bmatrix},$ (23a)
$\displaystyle\underline{E}_{ex,N}^{(1,2)}$ $\displaystyle=$
$\displaystyle-\frac{J}{2}\begin{bmatrix}-\Sigma
E^{+}&E_{2}^{(b,e^{+})}/V_{e}&0\\\ E_{2}^{(e^{+},b)}/V_{b}&-\Sigma
E^{+}-\Sigma E^{-}&E_{2}^{(b,e^{-})}/V_{b}\\\
0&E_{2}^{(e^{-},b)}/V_{b}&-\Sigma E^{-}\end{bmatrix},$ (23b)

where we define $\Sigma E^{e}=E_{3}^{b,e}/V_{e}+E_{3}^{e,b}/V_{b}$, and

$\displaystyle E_{1}^{b,e}$ $\displaystyle=$ $\displaystyle
R^{1,1}-R^{2,2}+i(R_{1,2}+R^{2,1}),$ (24a) $\displaystyle E_{2}^{b,e}$
$\displaystyle=$ $\displaystyle R^{1,1}+R^{2,2}-i(R_{1,2}-R^{2,1}),$ (24b)
$\displaystyle E_{3}^{b,e}$ $\displaystyle=$ $\displaystyle R^{3,3},$ (24c)

With these blocks, the frequency contribution due to exchange is written with
the block-diagonal matrices as

$\displaystyle\underline{\Omega}_{ex}^{(1,1)}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2M_{s}}\begin{bmatrix}\underline{E}_{ex,1}^{(1,1)}&&\\\
&\ddots&\\\ &&\underline{E}_{ex,N}^{(1,1)}\end{bmatrix}$ (25a)
$\displaystyle\underline{\Omega}_{ex}^{(1,2)}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2M_{s}}\begin{bmatrix}\underline{E}_{ex,1}^{(1,2)}&&\\\
&\ddots&\\\ &&\underline{E}_{ex,N}^{(1,2)}\end{bmatrix}$ (25b)

###  III.4 Dipole field

The dipole field is essential to compute the spin-wave band structure for
ASIs. We distinguish two contributions to the dipole field: a static
contribution originating from the equilibrium magnetization, and a dynamic
contribution originating from the long-range dynamics of macrospins. An
analogous way to phrase this, is that we consider a perturbation to the
dynamical matrix where the zeroth order term is the static stray field from
the magnetization and the first-order correction is the dipole-dipole
contribution.

####  III.4.1 Static contribution

To compute the static contribution of the dipole, we consider the stray field
from each nanomagnet in the ASI on a macrospin $i$. As an approximation, we
implemented the analytical expressions of the stray field from a rectangular
prism derived by R. Engel-Herbert and T. Hesjedal [36]. The resulting field
due to nanomagnet $n$ is $\mathbf{H}_{stray,N}^{(i)}$ and is computed as a
function of the center position of the nanomagnet $n$ and the position of the
_macrospin_ $i$. The analytical expressions derived in Ref. [36] are written
in Appendix B. Essentially, this computation provides a local field source for
macrospin $i$ so that it contributes to the Hamiltonian matrix as an external
magnetic field.

A subtle difference between a truly external field and the stray field is that
we need to scale the latter to the fractional volume of the macrospin it is
acting upon. In other words, we impose that the total energy on the target
nanomagnet due to nanomagnet $n$ is conserved

$E=\mu_{0}M_{s}\frac{\sum_{i=1}^{3}V_{i}\mathbf{H}_{stray,n}^{(i)}\cdot\mathbf{m}_{i}}{V}.$
(26)

As a consequence, the contribution to the matrix becomes

$\displaystyle\underline{\Omega}_{stray}^{(1,1)}$ $\displaystyle=$
$\displaystyle 0,$ (27a) $\displaystyle\underline{\Omega}_{stray}^{(1,2)}$
$\displaystyle=$
$\displaystyle\frac{\gamma\mu_{0}}{2}\sum_{\tau_{1},\tau_{2}}^{\mathrm{ASI}}\left[\sum_{i=1}^{3}\frac{V_{i}}{V}R(\theta_{m},\varphi_{m})\cdot\mathbf{H}_{stray,n}^{(i)})^{T}\cdot\hat{\mathbf{e}}_{3}\right]\mathbf{I},$
(27b)

where $\tau_{1}$ and $\tau_{2}$ are integer numbers so that
$\tau_{1}=\tau_{2}=0$ is the unit cell. The maximum value of $\tau_{1}$ and
$\tau_{2}$ is capped so that the long-range dipole contributions converge with
sufficient numerical accuracy. Gænice permits to specify either the maximum
value for $\tau_{1}$ and $\tau_{2}$ or to expand the lattice until numerical
accuracy is achieved.

####  III.4.2 Dynamic contribution

The dipole field of macrospin $j$ acting on macrospin $i$ is calculated using
the following expression.

$\mathbf{H}_{\mathrm{d},ij}=\frac{V_{j}M_{s,j}}{4\pi}\left[\frac{3\mathbf{r}_{i,j}(\mathbf{r}_{i,j}\cdot\mathbf{m}_{j})}{|\mathbf{r}_{i,j}|^{5}}-\frac{\mathbf{m}_{j}}{|\mathbf{r}_{i,j}|^{3}}\right],$
(28)

where, $\mathbf{r}_{i,j}$ is the distance between the two macrospins $i$ and
$j$.

We adopt a tight-binding-like approach for periodic structures whereby the
dipole field is computed within and between unit cells. Therefore, the long-
range terms collapse into a single Hamiltonian matrix. Identifying each
macrospin’s spatial position within a nanomagnet is critical. Hence, we adopt
the convention that the bulk macrospin is located at the nanomagnet’s
geometric center, denoted by $\mathbf{X}_{b}$. The edge macrospins can be
computed by

$\mathbf{X}_{e^{\pm}}=\mathbf{X}_{b}\pm\Delta l\hat{\mathbf{D}}.$ (29)

We consider the periodic structure established by the translation vectors
$a_{1}$ and $a_{2}$. In general, the distance between macrospins $i$ and $j$
is given by

$\mathbf{r}_{ij}^{(\tau_{1},\tau_{2})}=\mathbf{X}_{i}-\mathbf{X}_{j}-(\tau_{1}\mathbf{a}_{1}+\tau_{2}\mathbf{a}_{2}).$
(30)

Therefore, the total nonlocal dipole field acting on a macrospin $i$ can be
written as

$\displaystyle\mathbf{H}_{\mathrm{d},ij}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\sum_{\tau_{1},\tau_{2}}^{\mathrm{ASI}}\sum_{j}^{\mathrm{U.C.}}V_{j}M_{s,j}\Big{[}\frac{3\mathbf{r}_{ij}^{(\tau_{1},\tau_{2})}(\mathbf{r}_{ij}^{(\tau_{1},\tau_{2})}\cdot\mathbf{m}_{j}^{(\tau_{1},\tau_{2})})}{|\mathbf{r}_{ij}^{(\tau_{1},\tau_{2})}|^{5}}$
(31) $\displaystyle-$
$\displaystyle\frac{\mathbf{m}_{j}}{|\mathbf{r}_{ij}^{(\tau_{1},\tau_{2})}|^{3}}\Big{]},$

Because the distances in Eq. (30) are computed in the natural Cartesian
coordinates, we need to rotate into the basis of each macrospin to compute the
products as a function of coupled complex amplitudes. For this we define two
distances in the rotated reference frame

$\displaystyle\mathbf{\rho}_{ij}$ $\displaystyle=$ $\displaystyle
R(\theta_{j},\phi_{j})\cdot\mathbf{r_{ij}},$ (32a)
$\displaystyle\mathbf{\alpha}_{ij}$ $\displaystyle=$ $\displaystyle
R(\theta_{j}-\theta_{i},\phi_{j}-\phi_{i})\cdot\mathbf{r_{ij}}.$ (32b)

Therefore, the net field on macrospin $i$, expressed in the basis of $i$, is

$\displaystyle\mathbf{H}_{\mathrm{d},ij}^{(i)}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\sum_{\tau_{1},\tau_{2}}^{\mathrm{ASI}}\sum_{j}^{\mathrm{U.C.}}V_{j}M_{s,j}\Big{[}\frac{\mathbf{\alpha}_{ij}(\mathbf{\rho}_{ij}\cdot\mathbf{m_{j}})}{r_{ij}^{5}}$
(33)
$\displaystyle-\frac{R(\theta_{j}-\theta_{i},\phi_{j}-\phi_{i})\cdot\mathbf{m_{j}}}{r_{ij}^{3}}\Big{]}$

The last step to collapse the sums into a single Hamiltonian matrix is to
incorporate a tight-binding approach. We invoke Bloch’s theorem but we make
the assumption that the phase between the complex amplitudes is solely given
by the translation vectors between unit cells. This is the main approximation
in our model and ensures that the resulting band structure is periodic within
the FBZ. If macrospin-to-macrospin phases were to be included, then length
scales smaller than the FBZ would be resolved, which is outside the model’s
scope. Therefore, we apply Bloch’s theorem as

$\mathbf{m}_{j}^{(\tau_{1},\tau_{2})}=\mathbf{m}_{j}e^{\Phi}=\mathbf{m}_{j}e^{-i\mathbf{a}_{ij}^{(\tau_{1},\tau_{2})}\cdot\mathbf{k}}\,$
(34)

and compute the dipole energy,

$\mathbf{\mathcal{H}^{(i)}}=\mu_{0}M_{s,i}\mathbf{m}_{i}^{T}\cdot\mathbf{H}_{\mathrm{d},ij}^{(i)}.$
(35)

Rewriting the energy as a function of the complex amplitudes $a$ and rescaling
to units of frequency, we obtain the block Hamiltonian matrices

$\displaystyle\underline{\Omega}_{d}^{1,1}$ $\displaystyle=$
$\displaystyle\sum_{\tau_{1},\tau_{2}}^{\mathrm{ASI}}e^{\Phi}\begin{bmatrix}0&C_{12}&0\\\
C_{21}&0&C_{23}\\\ 0&C_{32}&0\end{bmatrix}$ (36a)
$\displaystyle\underline{\Omega}_{d}^{1,2}$ $\displaystyle=$
$\displaystyle\sum_{\tau_{1},\tau_{2}}^{\mathrm{ASI}}e^{\Phi}\begin{bmatrix}G_{11}&D_{12}&0\\\
D_{21}&G_{22}&C_{23}\\\ 0&D_{32}&G_{33}\end{bmatrix},$ (36b)

where,

$\displaystyle C_{ij}$ $\displaystyle=$
$\displaystyle\frac{3}{|\mathbf{r}_{ij}|^{5}}[\alpha_{1}\rho_{1}+i(\alpha_{1}\rho_{2}+\alpha_{2}\rho_{1})-\alpha_{2}\rho_{2}]$
(37a)
$\displaystyle-\frac{1}{|\mathbf{r}_{ij}|^{3}}[R^{(1,1)}-R^{(2,2)}+i(R^{(1,2)}+R^{(2,1)})],$
$\displaystyle D_{ij}$ $\displaystyle=$
$\displaystyle\frac{3}{|\mathbf{r}_{ij}|^{5}}[\alpha_{1}\rho_{1}+i(\alpha_{1}\rho_{2}-\alpha_{2}\rho_{1})+\alpha_{2}\rho_{2}]$
(37b)
$\displaystyle-\frac{1}{|\mathbf{r}_{ij}|^{3}}[R^{(1,1)}+R^{(2,2)}-i(R^{(1,2)}-R^{(2,1)})],$
$\displaystyle G_{ij}$ $\displaystyle=$
$\displaystyle-\frac{3}{|\mathbf{r}_{ij}|^{5}}[\alpha_{3}\rho_{3}]-\frac{2}{|\mathbf{r}_{ij}|^{5}}R^{(3,3)},$
(37c)

and we have used a slightly shorthand notation where
$\rho_{i,j}=(\rho_{1},\rho_{2},\rho_{3})$,
$\alpha_{i,j}=(\alpha_{1},\alpha_{2},\alpha_{3})$, $R$ represent the
components of the $3\times 3$ matrix
$R(\theta_{j}-\theta_{i},\phi_{j}-\phi_{i})$, and the sums over the ASI modify
the phase $\Phi$.

##  IV Validation

Gænice is implemented in MATLAB and it can be obtained from
http://doi.org/10.17605/OSF.IO/YUNHD as well as a script that reproduces the
results presented below.

Figure 3:  Comparison of numerical computation of ferromagnetic resonance
(FMR) and Kittel’s equation to validate the geometry and implementation of the
rotation matrix. (a) Field-dependent FMR of a perpendicularly magnetized easy-
plane ferromagnet. The blue circles represent calculations for a magnet in the
$x-y$ plane and magnetized along the $z$ direction. The calculations agree
with Kittel’s equation (38). The gold asterisks are obtained when the magnet
and the magnetization are oriented along $\theta_{d}=\theta=33$ deg and
$\varphi_{d}=\varphi=233$ deg. The field is perpendicular to this orientation,
and we recover the same field-dependent frequency. (b) Angle dependence FMR of
an in-plane magnetic saturated at $H=1,000$ kA/m. Both numerical computations
(blue circles) and Kittel’s equation (39) agree. (c) Validation of the
anisotropy field implementation. Field-dependent FMR for magnets of different
sizes (colored symbols) and the corresponding Kittel’s equation (41).

###  IV.1 Local fields

The validity of the implementation of local fields can be verified my means of
the Kittel equation. For the purposes of the model presented here, it is
imperative to verify the field magnitude and angle dependent ferromagnetic
resonance as well as its independence from the coordinate system.

We first model a circular thin-film which can be considered as a single
macrospin due to the fact that only the hard axis contributes to the demag
tensor, i.e., $D_{3}=1$ and $D_{1}=D_{2}=0$. We use a saturation magnetization
of $M_{s}=800$ kA/m. Kittel’s equation as a function of the magnetic field
amplitude $H$ is thus:

$\omega=\gamma\mu_{0}M_{s}\left(\frac{H}{M_{s}}-1\right),$ (38)

valid for $H>M_{s}$. To model this scenario, we use an applied field oriented
along the $z$-axis with magnitude in the range $800$ kA/m $<H<1,200$ kA/m. The
magnetization is parallel to the applied field so that $\theta_{m}=0$ and
$\varphi_{m}=0$. The magnetic film must be rotated so that the hard axis is
also oriented along the $z$-axis. In other words, $\theta_{d}=\pi/2$ and
$\varphi_{d}=0$. The numerical results are shown in Fig. 3(a) by blue circles.
The solution of the Kittel equation is shown at the top by a solid black line.
The agreement is within numerical error ($<4\times 10^{-15}$).

Validation of the rotation matrix is achieved by computing the same field
dependence when the magnetization, external field, and nanomagnet are rotated
by arbitrary polar and azimuth angles. For example, selecting a rotation
$\theta=33$ deg and $\varphi=233$ deg, we recover the correct solution within
numerical error ($<9\times 10^{-15}$), shown in Fig. 3(a) by gold asterisks.

We now set the external field magnitude to $H=1,000$ kA/m and vary its angle,
$\theta_{0}$. The frequency as a function of angle is obtained from the Kittel
equation expressed as

$\omega=\gamma\mu_{0}\sqrt{H_{i}\left(H_{i}+M_{s}\cos^{2}{(\theta_{0})}\right)},$
(39)

where $H_{i}$ is the internal magnetic field magnitude obtained by solving the
magnetostatic equations

$\displaystyle(H_{i}+M_{s})\cos{\theta_{i}}$ $\displaystyle=$ $\displaystyle
H\cos{(\theta_{0})},$ (40a) $\displaystyle H_{i}\sin{(\theta_{i}})$
$\displaystyle=$ $\displaystyle H\sin{(\theta_{0})}.$ (40b)

The magnetization vector is oriented along the internal magnetic field angle
for a saturating field, $\theta_{m}=\theta_{i}$ and $\varphi_{m}=0$. We define
$\theta_{0}=0$ as the out-of-plane component so that $\theta_{d}=\pi/2$ and
$\varphi_{d}=0$ for all cases. The results shown in Fig. 3(b) further validate
the implementation of the external field and demag fields.

Finally, we vary the size of the magnetic element so that all three
demagnetizing factors are computed. We consider three different oblate
spheres: “Large” ($10,000~{}\mathrm{nm}\times 1,000~{}\mathrm{nm}\times 5$
nm), “Medium” ($1,000~{}\mathrm{nm}\times 100~{}\mathrm{nm}\times 5$ nm), and
“Small” ($100~{}\mathrm{nm}\times 10~{}\mathrm{nm}\times 5$ nm). The field is
once again considered to be oriented along the $z$-axis and its magnitude is
varied between $800$ kA/m $<H<1,200$ kA/m. The frequency dependence as a
function of field is given by Kittel’s equation

$\omega=\gamma\mu_{0}M_{s}\sqrt{\left(\frac{H}{M_{s}}+D_{1}-D_{3}\right)\left(\frac{H}{M_{s}}+D_{2}-D_{3}\right)}.$
(41)

The results shown in Fig. 3(c) validate the demagnetization field and its
nanomagnet-dependent implementation because all three nanomagnets are
concurrently simulated. This also shows that Gænice can be used as a tool to
quickly compute FMR for an ensemble of uncoupled nanomagnets.

###  IV.2 Nonlocal field: exchange

The dynamic contribution of the exchange energy introduces the splitting of
the resonant frequencies within a single nanomagnet. As a test case, we set a
stadium-shaped nanomagnet with dimensions $l=280$ nm, $w=100$ nm, and $t=10$
nm. The nanomagnet is oriented along the $x$ axis, $\theta_{d}=\pi/2$ and
$\varphi_{d}=0$.

Figure 4:  Frequencies computed by adding the exchange Hamiltonian. (a) A
single macrospin is azimuthally rotated, resulting in a sizeable variation of
the frequency in one band. The minima occur at $\varphi_{1}=0$ deg and
$\varphi_{1}=90$ deg, consistent with a dynamic coupling mediated only by the
magnetization’s $z$ component. (b) Frequency variation as a function of
$\Delta l$, showing divergence as either the bulk or edge volumes tend to
zero. The frequencies in the vicinity of the default value $\Delta l=(2l-w)/4$
(the transition between the white and gray areas) are approximately constant.
Each band is displayed in different colors for clarity. Figure 5:  Frequency
dependence on the nanomagnet’s size. The aspect ratio is maintained in the
calculations so that the length is a representative metric. The solid lines
are frequencies obtained from computations on a single nanomagnet. The circles
are obtained from a single computation including five non-interacting
nanomagnets. Each band is displayed in different colors for clarity.

We first explore the effect of the magnetization’s relative angles. For this,
we set the magnetization parallel to the nanomagnet orientation, and we vary
the azimuth angle of the magnetization at one extremum, $\varphi_{1}$. The
computed frequencies are shown in Fig. 4, where different colors and dashed
curves were used for each branch for clarity. One frequency branch exhibits a
sinusoidal variation, consistent with one magnetization being rotated and
modifying the exchange energy. The maximum occurs at $45$ deg, implying that
the maximum exchange contribution occurs when the adjacent magnetization
vectors dynamically couple in both $x$ and $y$. Indeed, when $\varphi_{1}=0$
deg or $90$ deg, the $z$ components of the magnetization vectors are coupled,
leading to identical energy contributions to the eigenmodes. Note that this is
different than the static exchange energy computed in Eq. (43).

We now explore the influence of the bulk and edge volume ratios. From Eq.
(23), the exchange energy diverges as either the bulk or edge volume tends to
zero. This is expected because of the underlying assumption that the
nanomagnet is separated in three macrospins. In other words, such a divergence
has no physical origin. The frequencies computed as a function of $\Delta l$
when all the magnetization vectors are aligned with the nanomagnet are shown
in Fig. 4(b). Clearly, the frequencies diverge when the bulk and edge volumes
tend to zero towards the left and right extrema of the figure, respectively.
The frequencies are relatively constant close to the default distance $\Delta
l=(2l-w)/4=110$ nm.

As the size of the nanomagnet increases, the effect of the exchange
interaction in the frequencies must necessarily decrease insofar as the
nanomagnet is split into three macrospins. We compute this test scenario by
locking $\Delta l$ to its default value and varying the nanomagnet’s size. We
maintain the aspect ratio of the nanomagnet so that the length is
representative of the nanomagnet’s volume. The results are shown by solid and
dashed curves in Fig. 5, where the x axis is shown in natural logarithmic
scale and the colors represent different branches for clarity. As expected,
the frequency splitting diminishes as the nanomagnet’s size increases.

A final test for the exchange interaction, is to verify that its
implementation is independent of the number of nanomagnets. For this, we
specify five nanomagnets with lengths $100$, $200$, $300$, $400$, and $1000$
nm and dimensions consistent with the aspect ratio of the test case considered
in this section. The resulting eigenvalue problem requires solving for a
matrix of dimension $30\times 30$. The frequencies are shown by circles in
Fig. 5, color-coded according to the branches of the single-nanomagnet
calculations. We note that the frequencies are not automatically sorted for
each nanomagnet: only the computation of the eigenvectors can return such type
of sorting which is not currently computed in our implementation, as discussed
in the conclusions. In this case, the frequencies were manually sorted. The
results are in agreement with the calculations done for each nanomagnets,
validating that the exchange interaction is nonlocal but intrinsic to each
nanomagnet, i.e, there is no coupling between nanomagnets.

Figure 6:  Frequencies as a function of distance between two identical
nanomagnets interacting via (a) nonlocal dipole field and (b) both nonlocal
dipole field and exchange interaction within each nanomagnet. The magnets are
parallel to one another, and we distinguish the relative magnetization being
parallel (solid black curves) and anti-parallel (dashed red curves). In all
cases, the nonlocal dipole field becomes negligible at large distances and the
bands become degenerate, as expected for non-interacting identical
nanomagnets.  Figure 7:  Geometrical variations between two interacting
nanomagnets. The nanomagnet with director $\hat{D}_{1}$ is located at a
distance $d=300$ nm along the $\hat{y}$ direction. The director’s orientation
is varied by the polar and azimuth angles $\theta_{1}$ and $\varphi_{1}$,
respectively.  Figure 8:  Computed frequencies for a pair of interacting
nanomagnets when one of the nanomagnets is rotated about (a) the polar angle
$\theta_{1}$ and (b) the azimuth $\varphi_{1}$. The angles are shown in the
schematic Fig. 7.

###  IV.3 Nonlocal field: dipole

We now investigate the interaction between two identical nanomagnets of
dimensions $l=280$ nm, $w=100$ nm, and $t=10$ nm, as used in the previous
section. We focus here on collective excitation, so that $|\mathbf{k}|=0$ and
the phase contributions in Eqs. (36) simplify to 1.

The first test ensures that the nonlocal dipole field’s strength depends on
the distance between the nanomagnets. For this, we consider a varying distance
$d$ along the $y$ axis ranging from $100$ nm to $1,000$ nm. The computed
frequencies considering only nonlocal dipole fields are shown in Fig. 6(a).
The relative magnetization orientation is parallel for the solid black curves
and antiparallel for the dashed red curves. In both cases, the modes are
degenerate at long distances. This is a clear indication that nonlocal dipole
field does not affect the internal modes of non-interacting (or weakly
interacting) nanomagnets. Modes are visibly split under a distance of $\approx
400$ nm. Red and blue-shifts are observed for the parallel and antiparallel
cases, respectively, in agreement with the static dipole energy for each.
Including exchange energy,shown in Fig. 6(b), naturally leads to larger split
bands because of the additional energy. As expected in all cases, the modes
converge towards degenerate values at large distances indicating a negligible
interaction mediated by the nonlocal dipole field.

We next explore the frequency dependence on the relative orientation between
the two nanomagnets. For this, we consider a nanomagnet located at the global
origin of the Cartesian coordinate and a second nanomagnet located at a
distance of $d=300$ nm in the $y$ direction with a varying unit vector
$\hat{\mathbf{D}_{1}}$. We consider both polar and azimuth rotations
parametrized by the angles $\theta_{1}$ and $\varphi_{1}$, respectively, as
shown Fig. 7. Note that these angles are measured relative to the orientation
of the fixed nanomagnet.

The computed frequencies are shown in Fig. 8. The frequency variation as a
function of the polar angle $\theta_{1}$ is shown in (a). In these
computations, we disabled the exchange interaction to focus on the symmetry of
the static dipole field. There is a modest change in the frequency that is
maximal at $\theta_{1}=180$ deg. The symmetry is also consistent with the fact
that $90$ deg and $270$ deg are degenerate. The frequency variation as a
function of the azimuth angle $\varphi_{1}$ is shown in (b). There are again
clear symmetries consistent with the rotation of the nanomagnet despite the
increased number of modes originating from the non-collinear magnetization
orientations. Notably, at $90$ deg and $270$ deg, the rotated nanomagnet is
perpendicular to the fixed nanomagnet and the spacing is just $110$ nm. Strong
variations are observed close to these conditions. As expected, the computed
frequencies are periodic for both $\theta_{1}$ and $\varphi_{1}$.

This concludes the verification of the static dipole field, which follows the
qualitative expectations of decay with distance and symmetries due to
different types of rotations.

###  IV.4 Band structure

####  IV.4.1 Nanomagnet chain

A one-dimensional chain of nanomagnets is modeled by a single nanomagnet with
dimensions $l=280$ nm, $w=100$ nm, and $t=10$ nm subject to a translation
vector $\mathbf{a}_{1}$ oriented at an azimuth $\varphi_{a}$ and lattice
constant $|\mathbf{d}|=300$ nm. Magnons with wave vectors $\mathbf{k}$
oriented at an azimuth $\varphi_{k}$ are computed, as shown schematically in
Fig. 9.

The magnon dispersion is computed for cases where we set the wavevector
parallel to the $x$-axis and rotating the translation vector. In other words,
$\mathbf{a}_{1}=d\hat{\mathbf{a}}_{1}=d[\cos{(\varphi_{a})}\hat{x}+\sin{(\varphi_{a})}\hat{y}]$.
The resulting dispersion relations as a function of
$k=2\pi/d(\mathbf{k}\cdot\hat{\mathbf{a}}_{1})$ within the first Brillouin
zone for selected azimuths are shown in Fig. 10(a). Three bands are observed
because a single nanomagnet composes the unit cell of the chain. The band
structures show a pronounced periodic behavior in the FBZ that is symmetric
relative to $\varphi_{a}=90$ deg, stemming from the product
$\mathbf{k}\cdot\mathbf{a}_{1}$. This symmetry validates the implementation of
the sums performed in Eqs. (36). It is also shown in the mid panel of Fig.
10(a) that the band structure perpendicular to the chain orientation at
$\varphi_{a}=90$ deg, is flat. This is because the phase in Eqs. (36) is
exactly zero when the translation vector and wavevector are perpendicular.

Figure 9:  A chain of nanomagnets is modeled by a single nanomagnet (green)
upon which the nonlocal dipole field from an infinite chain of nanomagnets
(blue) acts. The chain can be defined along an arbitrary in-plane direction by
setting the translation vector $\mathbf{a}_{1}$ oriented at an azimuth
$\varphi_{1}$. The magnon dispersion can be computed for arbitrary in-plane
wavevectors $\mathbf{k}$ given the azimuth $\varphi_{k}$.  Figure 10:
Dispersion relation for a 1D nanomagnet chain upon (a) setting the wavevector
along $\hat{x}$ and rotating $\mathbf{a}_{1}$ and (b) setting $\mathbf{a}_{1}$
along the $\hat{x}$ direction and rotating the wavevector. The bands exhibit
the most changes when the wavevector and translation vector are parallel and
are flat when these are orthogonal. This is in agreement with our tight-
binding definition of the phase in Eqs. (36). Figure 11:  (a) Gænice
representation of the 2D array of nanomagnet chains, depicting the translation
vectors. (b) Resulting band structure where the FBZ is shown under the band
structure. The FBZ is directly computed by Gænice from the translation vectors
and the high-symmetry points are also identified. The band structure is
obtained by performing a Delaunay triangulation over the FBZ and evaluating
the resulting wavevectors. The color scale represents the frequency and it is
also shown in the vertical axis. (b) Irreducible path in the FBZ exhibiting
the asymmetry of this geometry as well as the periodicity achieved by the
tight-binding method. Figure 12:  Gænice representation of the square ice
geometry for the (a) vortex and (b) remanent states. The respective band
structure for each case in shown in (c) and (d) while the irreducible path in
the FBZ are shown in (e) and (f).

It is also possible to consider another case of $\varphi_{a}\neq\varphi_{k}$.
We set the translation vector along the $x$ axis, i.e.,
$\mathbf{a}_{1}=d\hat{x}$ and we vary $\varphi_{k}$, such that $\mathbf{k}$=
$|\mathbf{k}|[\cos{(\varphi_{k})}\hat{x}+\sin{(\varphi_{k})}\hat{y}]$. The
dispersion relations for $k=|\mathbf{k}|$ up to the FBZ are shown in Fig.
10(b). As for Fig. 10(a), the expected symmetries are respected, e.g. the band
structure perpendicular to the chain orientation is flat; see $\varphi_{k}=90$
and $270$ deg, for the same reasons outlined above.. Note that here we extend
the rotation of $\varphi_{k}$ to a full cycle.

####  IV.4.2 Multiple interacting chains

We now calculate the band structure of interacting nanomagnet chains. Each
nanomagnet is oriented at $\phi_{d}=0$ with respect to the x-axis and the
array is generated from the single nanomagnet unit cell due to translation
vectors $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ with a lattice constant of
$d=300$ nm. A visualization of this configuration is produced by Gænice to
ensure the correct geometry definition, shown in Fig. 11(a). The magnon band
structure is computed by an automatic determination of the FBZ and its
subsequent Delaunay triangulation to produce an array of wavevectors
$\mathbf{k}=k_{x}\hat{x}+k_{y}\hat{y}$. This feature allows to optimally map
the FBZ and produce band surfaces, as shown in Fig. 11(b).

As for the 1D nanomagnet chain, the FBZ also shows three bands because the
unit cell consists of a single nanomagnet. By examining the band structure
depicted in Fig. 11(b), the frequencies calculated along $k_{x}\hat{x}$
exhibit pronounced variations while it is predominantly flat along
$k_{y}\hat{y}$. This is consistent with our tight-binding approach whereby the
dynamic dipole coupling depends on the gap distance between nanomagnets.

Figure 13:  Gænice representation of the Kagome ice geometry for (a) identical
nanomagnets and (b) anisotropy modified nanomagnets. The respective band
structure for each case in shown in (c) and (d) while the irreducible path in
the FBZ exhibiting are shown in (e) and (f).

The irreducible path in the FBZ can be also directly computed in Gænice. We
observe that the band structure is different when the path is taken towards
the $X$ and $X^{\prime}$ points. This is because the array is asymmetric, such
that the dipole field is different along $\hat{x}$ and $\hat{y}$ directions.

####  IV.4.3 Square ice

We now use Gænice to compute the magnon band structure for square ASI, where
four nanomagnets are placed around a vertex, each at an angle of 90 degrees to
one another and equidistant from the vertex. We maintain the previously used
nanomagnet dimensions $l=280$ nm, $w=100$ nm, and $t=10$ nm but now set a
center-to-center distance of $d=430$ nm. For simplicity, we investigate the
band structure for states where the magnetization is in a homogeneous (onion)
state. Edge bending in the magnetization state leads to $S$ and $C$ states
that are known to modify the band structure [13].

We investigate both the vortex (type-I) and remanent (type-II) configuration.
The vortex state has four nanomagnets in the unit cell and is defined by the
translation vectors $\mathbf{a}_{1}=2\hat{x}$ and
$\mathbf{a}_{2}=\hat{x}+\hat{y}$. The remanent state has two nanomagnets in
the unit cell and is defined by $\mathbf{a}_{1}=\hat{x}$ and
$\mathbf{a}_{2}=\hat{y}$. These configurations are shown in Fig 12(a) and (b),
respectively.

The band structures in the FBZ are shown in Fig. 12(c) and 12(d), for the
vortex and remanent states, respectively. By dividing the nanomagnet into
three macrospins, the vortex state has twelve bands. The band structure
exhibits little dispersion, which may be expected by the fact that the stray
fields are largely compensated in a type-I configuration. The two high-
frequency bands are bulk modes while the low-frequency bands are edge modes,
as shown previously [13]. This band separation is clearly seen in the
irreducible path shown in Fig. 12(e), exhibiting a band-gap of about $\approx
5$ GHz.

In the remanent state, there are six bands. In this case, the bands are very
close together, with a visible dip at the $\Gamma$ point. It is also evident
that the band structure is skewed, which is a consequence of the likewise
skewed static dipole field in this configuration.

The results in this section are in agreement with previous calculations [13]
demonstrating the reconfigurability of the magnon band structure for square
ices. However, the improved dipole field implementation in Gænice showcases
more subtleties in the band structure as well as asymmetries that could in
principle indicate directional magnon propagation, as recently surmised in a
combined experimental and micromagnetic study [12]. We also emphasize that we
have only explored here the onion state, but it is well-known that the
magnetization tilts at the edges of the nanomagnets due to stray fields. This
adds an additional degree of freedom for tuning the band structure.

####  IV.4.4 Kagome ice

We now explore the band structure for Kagome ASI. The Kagome unit cell
comprises three nanomagnets with lattice constant $d=800$ nm which we define
as twice the radius of the circle in which the hexagonal structure is
embedded. Considering the center of the triad of nanomagnets as the origin, we
define the translation vectors $\mathbf{a}_{1}=\hat{x}$ and
$\mathbf{a}_{2}=(1/2)\hat{x}+(1/\sqrt{3}+1/4)\hat{y}$.

We consider two cases: a “regular” Kagome ice where the nanomagnets have
identical dimensions $l=280$ nm, $w=100$ nm, and $t=10$ nm; and an anisotropy
modified Kagome ice inspired by the work by T. Dion et al. [9], where we use
three different widths $w=100$ nm, $w=180$ nm, and $w=60$ nm for the
nanomagnets in the unit cell. The geometries are shown in Fig. 13(a) and Fig.
13(b). In both cases, the array is in a degenerate ground state where the unit
cell triad has a 2-in/1-out vertex.

In the “regular” Kagome ice we find a modest band structure with all nine
bands contributing to the band structure, shown in Fig. 13(c). However, the
anisotropy modified Kagome ice exhibits only four bands, as shown in Fig.
13(d) with other five softened to exactly zero. In the context of our
framework, a zero-frequency band entails a real, evanescent solution. The
bands in the irreducible paths in Figs. 13(e) and (f) further confirm that the
bands are relatively flat in all cases. An important distinction is the
anisotropy modified Kagome ice exhibits band-gaps which is consistent with the
different FMR for each nanomagnet, i.e., different demagnetization factors.
While this is certainly not an in-depth investigation of the frequency
response of anisotropy modified Kagome ices, it showcases the functionality of
Gænice to compute the band structure of relatively complex geometries with
ease.

##  V Conclusions

We have presented Gænice, a computational tool to compute the dispersion
relation of arbitrary artificial spin ice geometries. The theoretical
framework of Gænice relies on the excitation of small-amplitude perturbations
and produces the dispersion relation by computing both static and dynamic
dipole contributions to the Hamiltonian matrices. Our framework also relies on
a tight-binding approach to ensure the periodicity of solutions within the
FBZ, which composes the main simplification of the model.

Gænice can be also used for FMR computations of relatively complicated
geometries. For example, Gænice has been recently applied for square ices
based on trilayers and exhibited remarkable agreement with experiments and
micromagnetic simulations of field-dependent FMR [20]. Because both the
exchange and dipole interactions can be toggled, Gænice can be easily be used
to study the FMR of ensembles of interacting or non-interacting nanoparticles
and extended to 3D structures.

There are three main limitations to Gænice in its current form. First, the
computations are accurate for nanoparticles and nanomagnets because of the
assumption of three macrospins. Larger nanomagnets possess higher degrees of
freedom that will reduce the relative energy contributions. Therefore, Gænice
is likely to overestimate the frequency split when nanomagnets are brought
very close together. A possible solution to this issue is to further split the
nanomagnets into more macrospins, with the caveat that the number of
macrospins should be kept to a minimum to maintain a computational advantage
over micromagnetic simulations. Another way to solve this issue is to compute
the energy of the system to actively modify the magnetization’s edge bending
due to stray fields, as recently shown in Ref. [37]. Second, the wavefunctions
are not currently computed. As is well known, the linearization of the
Hamiltonian leads to wavefunctions with large errors. This limitation will be
resolved in future work. Third, a 3D band structure is not currently
supported. However, the basic framework is written and a generalization in 3D
will compose a simple expansion of the dipole phases in the tight-binding
approximation, a method that is well-known in solid-state physics.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
under Grant No. 2205796. AR and EI acknowledge support from the UCCS Committee
on Research and Creative Works (CRCW). VM and EI acknowledge the Department of
Physics and Energy Science at UCCS for the use of their facilities and
equipment. This work was supported by the Royal Academy of Engineering
Research Fellowships, awarded to JCG. JCG was supported by EPSRC grant
EP/X015661/1. Work by OGH at Argonne was supported by the US Department of
Energy, Basic Energy Sciences Division of Materials Sciences and Engineering.

## Data Availability Statement

The data that support the findings of this study are openly available in Open
Science Framework (OSF) at http://doi.org/10.17605/OSF.IO/YUNHD.

##  Appendix A Derivation of exchange energy

To compute the exchange energy, we use a simple quasi-one-dimensional spin
chain model to estimate the energy along the chain and relate it to the
nanomagnet’s regions and their volume. Consider a chain of length $l$ where
the magnetization vector is linearly rotated, so that

$\mathbf{m}=\cos{(k_{0}x)}\hat{\mathbf{x}}+\sin{(k_{0}x)}\hat{\mathbf{y}}.$
(42)

It can be shown that the exchange energy is given by
$E_{\mathrm{ex}}=AVk_{0}^{2}+E_{0}$, where $A$ is the exchange constant in
units of pJ/m, $V$ is the volume of the quasi-1D chain, and $E_{0}$ is a
constant of integration.

We consider now a nanomagnet of length $l$, width $w$, and thickness $d$,
split in three unequal pieces with boundaries at $l_{1}$ and $l_{2}$ so that
their volumes are $V_{1}=wtl_{1}$, $V_{2}=wt(l_{2}-l_{1})=wt\Delta l_{1,2}$
and $V_{3}=wt(l_{3}-l_{2})=wt\Delta l_{2,3}$. This scenario is schematically
shown in Fig. 14. The total exchange energy is

$E_{\mathrm{ex}}=J_{1}\left(\mathbf{m}_{1}\cdot\mathbf{m}_{2}\right)+J_{2}\left(\mathbf{m}_{2}\cdot\mathbf{m}_{3}\right),$
(43)

where the magnetization vectors are taken in the geometric center of each
piece. This leads to $\mathbf{m}_{1}\cdot\mathbf{m}_{2}=\cos{(k_{0}\Delta
l_{1,2})}$ and $\mathbf{m}_{2}\cdot\mathbf{m}_{3}=\cos{(k_{0}\Delta
l_{2,3})}$. Expanding the cosine to first order in Eq. (43) and equating to
the continuum solution, we obtain

$E_{\mathrm{ex}}=\left(J_{1}+J_{2}\right)+\left(J_{1}\frac{\Delta
l^{2}_{1,2}}{2}+J_{2}\frac{\Delta
l^{2}_{2,3}}{2}\right)k_{0}^{2}=E_{0}+AVk_{0}^{2},$ (44)

Given that the exchange constant is uniform in the nanomagnet, we can set
$J_{1}=C_{1,2}A$ and $J_{2}=C_{2,3}A$. From geometry, it can be shown that

$\displaystyle C_{1,2}$ $\displaystyle=$ $\displaystyle\frac{2}{\Delta
l^{2}_{1,2}}\left(V_{1}+\frac{V_{2}}{2}\right),$ (45a) $\displaystyle C_{2,3}$
$\displaystyle=$ $\displaystyle\frac{2}{\Delta
l^{2}_{2,3}}\left(V_{3}+\frac{V_{2}}{2}\right).$ (45b)

In the case of stadium-shaped nanomagnets, one can consider a symmetric
splitting so that $V_{1}=V_{3}=V_{e}$ and $V_{2}=V_{b}$, leading to the
expression shown in Eq. (22).

The edge volume $V_{e}$ and the bulk volume $V_{b}$ can be computed as a
function of $\Delta l$. We have two cases.

####  A.0.1 Case $l/4<\Delta l<(2l-w)/4$

This corresponds to the situation where an edge mode occupies more than the
half-circle in the stadium’s edge at the expense of the bulk mode. The edge
and bulk volumes are

$\displaystyle V_{e}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[wtl-w^{2}t\left(1-\frac{\pi}{4}\right)-V_{b}\right],$
(46a) $\displaystyle V_{b}$ $\displaystyle=$ $\displaystyle(4\Delta l-l)wt.$
(46b)

####  A.0.2 Case $(2l-w)/4\leq\Delta l<l/2$

This corresponds to the situation where an edge mode is confined to the half-
circle in the stadium’s edge. Computing the cone angle

$\theta=2\arccos{\left(1-\frac{2l-4\Delta l}{w}\right)},$ (47)

the edge and bulk volumes are

$\displaystyle V_{e}$ $\displaystyle=$
$\displaystyle\frac{\theta-\sin{(\theta)}}{8}w^{2}t,$ (48a) $\displaystyle
V_{b}$ $\displaystyle=$ $\displaystyle
wtl-w^{2}t\left(1-\frac{\pi}{4}\right)-2V_{e}.$ (48b)

The edge and bulk volumes as a function of $\Delta l$ are shown in Fig. 15.
The limiting case $\Delta l=(2l-w)/4$ is considered to be the default.

Figure 14:  Toy model for a quasi-1D spin chain splitted into unequal pieces
to estimate the exchange energy.  Figure 15:  Fractional ratio between edge
and volume modes. The default is considered at the edge of the gray area in
the limiting case $\Delta l=(2l-w)/4$.  Figure 16:  Stray fields computed with
Eqs. (49). We show the computed field under several rotations within Gænice
and in the position of neighboring macrospins. The rotation angles are
displayed for each figure.

##  Appendix B Derivation of exchange energy

In Ref. [36], the authors considered a rectangular prisms with its geometric
center at the origin of the Cartesian reference frame and sides
$2x_{b}>2y_{b}>2z_{b}$. The resulting expressions for the stray field are:

$\displaystyle H_{x}(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{M_{s}}{4\pi}\sum_{k,l,m=1}^{2}(-1)^{k+l+m}\mathrm{ln}\left\\{z+(-1)^{m}z_{b}+\sqrt{L(k,l,m)}\right\\},$
(49a) $\displaystyle H_{y}(x,y,z)$ $\displaystyle=$
$\displaystyle-\frac{M_{s}}{4\pi}\sum_{k,l,m=1}^{2}(-1)^{k+l+m}\frac{\left[y+(-1)^{l}y_{b}\right]\left[x+(-1)^{k}x_{b}\right]}{|y+(-1)^{l}y_{b}||x+(-1)^{k}x_{b}|}\times\mathrm{arctan}\left\\{\frac{|x+(-1)^{k}x_{b}|\left[z+(-1)^{m}z_{b}\right]}{|y+(-1)^{l}y_{b}|L(k,l,m)}\right\\},$
(49b) $\displaystyle H_{z}(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{M_{s}}{4\pi}\sum_{k,l,m=1}^{2}(-1)^{k+l+m}\mathrm{ln}\left\\{x+(-1)^{k}x_{b}+\sqrt{L(k,l,m)}\right\\},$
(49c)

where

$L(k,l,m)=\left[x+(-1)^{k}x_{b}\right]^{2}+\left[y+(-1)^{l}y_{b}\right]^{2}+\left[z+(-1)^{m}z_{b}\right]^{2}$
(50)

We note that the assumed orientation of the rectangular prism in Ref. [36] is
different than that assumed in Gænice. For this reason, we rotate the
expressions of Eq. (49) such that the easy axis of the rectangular prism
aligns with the $z$ axis. In Fig. 16 we show the calculated stray field from
rectangular prisms with different director vectors. In all cases, it is seen
by inspection that the stray field is computed correctly.

## References

  * Skjærvø  _et al._ [2019] S. H. Skjærvø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, “Advances in artificial spin ice,” Nature Reviews Physics (2019), 10.1038/s42254-019-0118-3.
  * Heyderman and Stamps [2013] L. J. Heyderman and R. L. Stamps, “Artificial ferroic systems: Novel functionality from structure, interactions and dynamics,” Journal of Physics: Condensed Matter 25, 363201 (2013).
  * Gliga  _et al._ [2013] S. Gliga, A. Kákay, R. Hertel, and O. G. Heinonen, “Spectral analysis of topological defects in an artificial spin-ice lattice,” Phys. Rev. Lett. 110, 117205 (2013).
  * Gliga, Iacocca, and Heinonen [2020] S. Gliga, E. Iacocca, and O. G. Heinonen, “Dynamics of reconfigurable artificial spin ice: Toward magnonic functional materials,” APL Materials 8, 040911 (2020), https://doi.org/10.1063/1.5142705 .
  * Lendinez and Jungfleisch [2019] S. Lendinez and M. B. Jungfleisch, “Magnetization dynamics in artificial spin ice,” J. Phys.: Condens. Matter 32, 013001 (2019).
  * Gartside  _et al._ [2018] J. C. Gartside, D. M. Arroo, D. M. Burn, V. L. Bemmer, A. Moskalenko, L. F. Cohen, and W. R. Branford, “Realization of ground state in artificial kagome spin ice via topological defect-driven magnetic writing,” Nature Nanotechnology 13, 53?58 (2018).
  * Mamica  _et al._ [2018] S. Mamica, X. Zhou, A. Adeyeye, M. Krawczyk, and G. Gubbiotti, “Spin-wave dynamics in artificial anti-spin-ice systems: Experimental and theoretical investigations,” Phys. Rev. B 98, 054405 (2018).
  * Arroo, Gartside, and Branford [2019] D. M. Arroo, J. C. Gartside, and W. R. Branford, “Sculpting the spin-wave response of artificial spin ice via microstate selection,” Phys. Rev. B 100, 214425 (2019).
  * Dion  _et al._ [2019a] T. Dion, D. M. Arroo, K. Yamanoi, T. Kimura, J. C. Gartside, L. F. Cohen, H. Kurebayashi, and W. R. Branford, “Tunable magnetization dynamics in artificial spin ice via shape anisotropy modification,” Phys. Rev. B 100, 054433 (2019a).
  * Iacocca, Gliga, and Heinonen [2020] E. Iacocca, S. Gliga, and O. G. Heinonen, “Tailoring spin-wave channels in a reconfigurable artificial spin ice,” Phys. Rev. Applied 13, 044047 (2020).
  * Micaletti and Montoncello [2023] P. Micaletti and F. Montoncello, “Dynamic footprints of the specific artificial spin ice microstate on its spin waves,” Magnetochemistry 9 (2023), 10.3390/magnetochemistry9060158.
  * Lendinez  _et al._ [2021] S. Lendinez, M. T. Kaffash, O. G. Heinonen, S. Gliga, E. Iacocca, and M. B. Jungfleisch, “Nonlinear multi-magnon scattering in artificial spin ice,” Nature Communications 14, 3419 (2021).
  * Iacocca  _et al._ [2016] E. Iacocca, S. Gliga, R. L. Stamps, and O. Heinonen, “Reconfigurable wave band structure of an artificial square ice,” Phys. Rev. B 93, 134420 (2016).
  * Lasnier and Wysin [2020] T. D. Lasnier and G. M. Wysin, “Magnetic oscillation modes in square-lattice artificial spin ice,” Phys. Rev. B 101, 224428 (2020).
  * Montoncello  _et al._ [2023] F. Montoncello, M. T. Kaffash, H. Carfagno, M. F. Doty, G. Gubbiotti, and M. B. Jungfleisch, “A Brillouin light scattering study of the spin-wave magnetic field dependence in a magnetic hybrid system made of an artificial spin-ice structure and a film underlayer,” Journal of Applied Physics 133, 083901 (2023).
  * Graczyk  _et al._ [2018] P. Graczyk, M. Krawczyk, S. Dhuey, W.-G. Yang, H. Schmidt, and G. Gubbiotti, “Magnonic band gap and mode hybridization in continuous permalloy films induced by vertical dynamic coupling with an array of permalloy ellipses,” Phys. Rev. B 98, 174420 (2018).
  * Negrello  _et al._ [2022] R. Negrello, F. Montoncello, M. T. Kaffash, M. B. Jungfleisch, and G. Gubbiotti, “Dynamic coupling and spin-wave dispersions in a magnetic hybrid system made of an artificial spin-ice structure and an extended nife underlayer,” APL Materials 10, 091115 (2022).
  * Wang  _et al._ [2006] R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, “Artificial spin ice in a geometrically frustrated lattice of nanoscale ferromagnetic islands,” Nature 439, 303–306 (2006).
  * Iacocca and Heinonen [2017] E. Iacocca and O. Heinonen, “Topologically nontrivial magnon bands in artificial square spin ices with dzyaloshinskii-moriya interaction,” Phys. Rev. Applied 8, 034015 (2017).
  * Dion  _et al._ [2019b] T. Dion, K. D. Stenning, A. Vanstone, H. H. Holder, R. Sultana, G. Alatteili, V. Martinez, M. T. Kaffash, T. Kimura, H. Kurebayashi, W. R. Branford, E. Iacocca, B. M. Jungfleisch, and J. C. Gartside, “Ultrastrong magnon-magnon coupling and chiral symmetry breaking in a 3d magnonic metamaterial,” arXiv:2306.16159 (2019b).
  * Abert [2019] C. Abert, “Micromagnetics and spintronics: models and numerical methods,” The European Physical Journal B 92, 120 (2019).
  * May  _et al._ [2019] A. May, M. Hunt, A. Van Den Berg, A. Hejazi, and S. Ladak, “Realisation of a frustrated 3d magnetic nanowire lattice,” Communication Physics 2, 13 (2019).
  * May  _et al._ [2021] A. May, M. Saccone, A. van den Berg, J. Askey, M. Hunt, and S. Ladak, “Magnetic charge propagation upon a 3d artificial spin ice,” Nature Communications 12, 3217 (2021).
  * Sahoo  _et al._ [2021] S. Sahoo, A. May, A. van Den Berg, A. K. Mondal, S. Ladak, and A. Barman, “Observation of coherent spin waves in a three-dimensional artificial spin ice structure,” Nano Letters 21, 4629–4635 (2021).
  * Slavin and Tiberkevich [2009] A. Slavin and V. Tiberkevich, “Nonlinear auto-oscillator theory of microwave generation by spin-polarized current,” Magnetics, IEEE Transactions on 45, 1875 –1918 (2009).
  * Grimsditch  _et al._ [2004] M. Grimsditch, L. Giovannini, F. Montoncello, F. Nizzoli, G. K. Leaf, and H. G. Kaper, “Magnetic normal modes in ferromagnetic nanoparticles: A dynamical matrix approach,” Phys. Rev. B 70, 054409 (2004).
  * Neusser  _et al._ [2011] S. Neusser, G. Duerr, S. Tacchi, M. Madami, M. L. Sokolovskyy, G. Gubbiotti, M. Krawczyk, and D. Grundler, “Magnonic minibands in antidot lattices with large spin-wave propagation velocities,” Phys. Rev. B 84, 094454 (2011).
  * Rychły  _et al._ [2015] J. Rychły, J. W. Kłos, M. Mruczkiewicz, and M. Krawczyk, “Spin waves in one-dimensional bicomponent magnonic quasicrystals,” Phys. Rev. B 92, 054414 (2015).
  * Gubbiotti  _et al._ [2018] G. Gubbiotti, X. Zhou, Z. Haghshenasfard, M. G. Cottam, and A. O. Adeyeye, “Reprogrammable magnonic band structure of layered permalloy/cu/permalloy nanowires,” Phys. Rev. B 97, 134428 (2018).
  * Lisiecki  _et al._ [2019] F. Lisiecki, J. Rychły, P. Kuświk, H. Głowiński, J. W. Kłos, F. Groß, N. Träger, I. Bykova, M. Weigand, M. Zelent, E. J. Goering, G. Schütz, M. Krawczyk, F. Stobiecki, J. Dubowik, and J. Gräfe, “Magnons in a quasicrystal: Propagation, extinction, and localization of spin waves in fibonacci structures,” Phys. Rev. Applied 11, 054061 (2019).
  * Gliga  _et al._ [2015] S. Gliga, A. Kákay, L. J. Heyderman, R. Hertel, and O. G. Heinonen, “Broken vertex symmetry and finite zero-point entropy in the artificial square ice ground state,” Phys. Rev. B 92, 060413 (2015).
  * Colpa [1978] J. H. P. Colpa, “Diagonalization of the quadratic boson hamiltonian,” Physica A: Statistical Mechanics and its Applications 93, 327 – 353 (1978).
  * Osborn [1945] J. A. Osborn, “Demagnetizing factors of the general ellipsoid,” Phys. Rev. 67, 351–357 (1945).
  * Aharoni [1998] A. Aharoni, “Demagnetizing factors for rectangular ferromagnetic prisms,” Journal of Applied Physics 83, 3432–3434 (1998).
  * Martinez and Iacocca [2023] V. Martinez and E. Iacocca, “A numerical method to determine demagnetization factors of stadium-shaped nanoislands,” IEEE Magnetics Letters xx, x–x (2023).
  * Engel-Herbert and Hesjedal [2005] R. Engel-Herbert and T. Hesjedal, “Calculation of the magnetic stray field of a uniaxial magnetic domain,” Journal of Applied Physics 97, 074504 (2005).
  * Saccone  _et al._ [2023] M. Saccone, J. Carter-Gartside, K. Stenning, W. R. Branford, and F. Caravelli, “From vertices to vortices in magnetic nanoislands,” Physics of Fluids 35, 017101 (2023).