# Thin cylindrical magnetic nanodots revisited: variational formulation,
accurate solution and phase diagram

Alexander Müller
Institute for Structural Mechanics
University of Stuttgart
70550 Stuttgart, Pfaffenwaldring 7
<EMAIL_ADDRESS>
& Manfred Bischoff
Institute for Structural Mechanics
University of Stuttgart
70550 Stuttgart, Pfaffenwaldring 7
<EMAIL_ADDRESS>& Marc-André Keip
Institute of Applied Mechanics
University of Stuttgart
70550 Stuttgart, Pfaffenwaldring 7
<EMAIL_ADDRESS>Corresponding author

###### Abstract

We investigate the variational formulation and corresponding minimizing
energies for the detection of energetically favorable magnetization states of
thin cylindrical magnetic nanodots. Opposed to frequently used heuristic
procedures found in the literature, we revisit the underlying governing
equations and construct a rigorous variational approach that takes both
exchange and demagnetization energy into account. Based on a combination of
Ritz’s method and a Fourier series expansion of the solution field, we are
able to pinpoint the precision of solutions, which are given by vortex modes
or single-domain states, down to an arbitrary degree of precision.
Furthermore, our model allows to derive an expression for the demagnetization
energy in closed form for the in-plane single-domain state, which we compare
to results from the literature. A key outcome of the present investigation is
an accurate phase diagram, which we obtain by comparing the vortex mode’s
energy minimizers with those of the single-domain states. This phase diagram
is validated with data of two- and three-dimensional models from literature.
By means of the phase diagram, we particularly find the critical radius at
which the vortex mode becomes unfavorable with machine precision. All relevant
data and codes related to the present contribution are available at [37].

_Keywords_ micromagnetics $\cdot$ variational formulation $\cdot$ Ritz method
$\cdot$ vortex mode $\cdot$ phase diagram

##  1 Introduction

Magnetic materials are usually characterized by the presence of magnetic
domains, that is, by regions with equally oriented magnetization. In thin
magnetic specimens, the formation of domains is often accompanied by the
creation of magnetic vortices [48, 25, 17, 45]. Such a vortex is characterized
by a swirl of magnetization that by default leads to an out-of-plane tilt of
magnetization in its core [7]. As a result of its circular shape, a favorable
specimen for the generation of a magnetic vortex is given by a cylindrical
magnetic nanodot. Due to the great stability of magnetic vortices in magnetic
nanodots, they are potential candidates for magnetic storage devices [14, 15].
Furthermore, magnetic nanodots are being explored for therapeutical and
diagnostical processes in biomedical engineering [33]. We refer to [47, 5, 24,
32, 10, 41, 20, 9, 3] for a variety of studies on the statics and dynamics of
magnetization in nanodots and to [21, 26] for a general overview of the topic.

In the present contribution, we are concerned with the analysis and the
prediction of magnetization states of thin cylindrical nanodots. In this
context, we are motivated by experiments [16, 44] and theoretical studies [47,
23]. On the theoretical side, we would like to highlight in particular the
contribution of Usov and Peschany [47], which has been serving as a foundation
for numerous approaches to study the stability as well as the associated
critical radius of the vortex configuration in the static setting. Here, the
critical radius is defined as the one at which the vortex mode becomes
unstable resulting in constant magnetization throughout the specimen, that is,
in a single-domain state. Theoretical upper and lower bounds of the critical
radius have been derived in the contributions [5, 35, 23], which are all based
on the aforementioned approach [47].111Deviating from [47], the approaches
documented in [4, 27] are based on an ansatz of the magnetization that is
given as a function of the thickness coordinate rather than the radial
coordinate. Associated approaches are related to a form of curling that is of
no interest in the present contribution.

In the present work we derive a variational principle based on the stray- and
exchange-energy density that provides the full set of micromagnetic governing
equations as Euler-Lagrange equations. We exploit this variational principle
to investigate energy minimizers with an appropriate Ritz ansatz for the
vortex mode. These minimizers are compared with the minimizers obtained from a
homogeneous in-plane magnetization. From the latter, we then derive the
magnetostatic energy in closed form for varying radii and thicknesses of a
nanodot. As a key outcome of the analysis, we eventually compute accurate
solutions for the critical radius of a stable vortex state.

Please note that in contrast to existing results found in literature, the
results presented herein are either exact or accurate up to machine precision.
In this way they perfectly suit as benchmarks for the validation of numerical
simulation schemes for magnetic nanodots and related problems of
magnetostatics.

##  2 Theory

###  2.1 Governing equations and constitutive law

In what follows, we consider a magnetic body in a configuration
${\mathcal{B}}\subset\mathbb{R}^{d}$ with boundary $\partial{\mathcal{B}}$,
where $\mathbb{R}^{d}$ is the $d$-dimensional Euclidean space. The body is
embedded in free space ${\Omega}\subset\mathbb{R}^{d}$ that is assumed to
include the space occupied by the magnetic body such that
${\mathcal{B}}\subset{\Omega}$. All space is spatially parameterized in the
coordinates ${\bm{x}}\in{\Omega}$. Then, in a purely magnetostatic setting,
the Maxwell equations read

$\operatorname{div}{\bm{B}}=0\quad\mbox{and}\quad\operatorname{curl}{\bm{H}}={\bm{0}}\quad\mbox{in}\
{\Omega},$ (1)

where ${\bm{B}}$ denotes the magnetic induction and ${\bm{H}}$ denotes the
magnetic field. The Maxwell equations are coupled through the micromagnetic
constitutive equation

${\bm{B}}=\mu_{0}({\bm{H}}+M_{s}{\bm{M}}),$ (2)

where $\mu_{0}$ and $M_{s}$ denote the magnetic permeability of vacuum and the
spontaneous magnetization, respectively. Furthermore, ${\bm{M}}$ is a unit-
vector field that encodes the orientation of magnetization, such that
${\bm{M}}\in{\mathcal{S}}^{2}$, where
${\mathcal{S}}^{d-1}=\\{{\bm{x}}\in\mathbb{R}^{d}\leavevmode\nobreak\
|\leavevmode\nobreak\ {\bm{x}}\cdot{\bm{x}}=1\\}$ is the $d$$-$$1$-dimensional
unit sphere embedded in $\mathbb{R}^{d}$. The corresponding static Landau-
Lifshitz equation will be derived as an Euler-Lagrange equation of the
variational formulation to be discussed next.

###  2.2 Variational formulation of the micromagnetic problem

####  2.2.1 Variational principle and Euler-Lagrange equations

As independent field variables of the micromagnetic problem we consider the
scalar magnetic potential $\Phi$ and the magnetization director ${\bm{M}}$ as

$\Phi:\left\\{\begin{array}[]{l}{\Omega}\rightarrow\mathbb{R}\\\
{\bm{x}}\mapsto\Phi({\bm{x}})\end{array}\right.\quad\mbox{and}\quad{\bm{M}}:\left\\{\begin{array}[]{l}{\mathcal{B}}\rightarrow{\mathcal{S}}^{d-1}\\\
{\bm{x}}\mapsto{\bm{M}}({\bm{x}})\end{array}\right..$ (3)

Based on the magnetic potential, we define the magnetic field
${\bm{H}}=-\operatorname{grad}\Phi$, such that the Maxwell equation Eq. 12 is
identically satisfied.

In absence of externally applied fields and assuming isotropic behavior, the
micromagnetic energy reads

$\displaystyle{\Pi}(\Phi,{\bm{M}})=\int_{\mathcal{B}}\psi_{\text{ex}}(\operatorname{grad}{\bm{M}})\,{\mathrm{d}}V+\int_{\Omega}\psi_{\text{stray}}(\operatorname{grad}\Phi,{\bm{M}})\,{\mathrm{d}}V,$
(4)

where $\psi_{\text{ex}}$ is the exchange-energy density and
$\psi_{\text{stray}}$ is the stray-field (or demagnetization) energy density.
They are given by [36, 18, 29, 43, 19]

$\displaystyle\psi_{\text{ex}}(\operatorname{grad}{\bm{M}})$
$\displaystyle=\frac{A}{2}\|\operatorname{grad}{\bm{M}}\|^{2},$ (5)
$\displaystyle\psi_{\text{stray}}(\operatorname{grad}\Phi,{\bm{M}})$
$\displaystyle=\mu_{0}\widehat{M}_{s}{\bm{M}}\cdot\operatorname{grad}\Phi-\frac{\mu_{0}}{2}\|\operatorname{grad}\Phi\|^{2},$

where $A$ is the exchange-energy coefficient and

$\displaystyle\widehat{M}_{s}({\bm{x}})=\begin{cases}M_{s},&{\bm{x}}\in{\mathcal{B}},\\\
0,&\text{otherwise}.\end{cases}$ (6)

For a simpler notation, we introduce a dimensionless version of the total
energy. In that consequence, we also change the definition of the coordinates
$L{\bm{x}}\rightarrow{\bm{x}}$, where the length $L=\ell_{ex}/\sqrt{2}$ is
defined through the magnetostatic exchange length
$\ell_{ex}=\sqrt{2A/(\mu_{0}M_{s}^{2})}$ [1]. The dimensionless version of the
total energy is accomplished by using the relations

$\phi(L{\bm{x}})=\frac{\phi({\bm{x}})}{M_{s}L},\quad{\bm{M}}(L{\bm{x}})={\bm{m}}({\bm{x}}),\quad{\Pi}(\Phi,{\bm{M}})=\widehat{{\Pi}}(\phi,{\bm{m}})\frac{L^{2}}{A},\quad{\bm{H}}(L{\bm{x}})=\frac{{\bm{h}}({\bm{x}})}{M_{s}},\quad\widehat{M}_{s}=\widehat{m}_{s}M_{s},$
(7)

In what follows, we stick to the notation ’$\operatorname{curl}$’,
’$\operatorname{div}$’ and ’$\operatorname{grad}$’ for the differential
operators, although from now on they are related to the dimensionless
coordinates ${\bm{x}}$. The dimensionless total energy then reads

$\displaystyle\widehat{{\Pi}}(\phi,{\bm{m}})=\int_{\widehat{{\mathcal{B}}}}\frac{1}{2}\|\operatorname{grad}{\bm{m}}\|^{2}+\operatorname{grad}\phi\cdot(\widehat{m}_{s}{\bm{m}})\,{\mathrm{d}}v-\int_{\widehat{{\Omega}}}\frac{1}{2}\|\operatorname{grad}\phi\|^{2}\,{\mathrm{d}}v$
(8)

where $\widehat{{\Omega}}$ and $\widehat{{\mathcal{B}}}$ denote the
appropriately rescaled domains. The corresponding variational principle is
given as

$\displaystyle\boxed{\\{\phi^{\ast},{\bm{m}}^{\ast}\\}=\arg\left\\{\sup_{\phi\in
W^{k,p}(\widehat{{\Omega}},\mathbb{R})}\inf_{\vphantom{\textnormal{X}}{\bm{m}}\in
W^{l,q}(\widehat{{\Omega}},{\mathcal{S}}^{d-1})}{\Pi}(\phi,{\bm{m}})\right\\}}$
(9)

where $W^{k,p}$ denotes a Sobolev space. The corresponding Euler-Lagrange
equations read

$\displaystyle\updelta\phi:$
$\displaystyle\quad\operatorname{div}(-\operatorname{grad}\phi+{\bm{m}})\,$
$\displaystyle=0$ in $\displaystyle\widehat{{\mathcal{B}}}$
$\displaystyle\text{ in }\mathbb{R},$ (10) $\displaystyle\updelta\phi:$
$\displaystyle\quad\operatorname{div}(-\operatorname{grad}\phi)\,$
$\displaystyle=0$ in
$\displaystyle\widehat{{\Omega}}\setminus\widehat{{\mathcal{B}}}$
$\displaystyle\text{ in }\mathbb{R},$ (11) $\displaystyle\updelta\phi:$
$\displaystyle\quad(-\operatorname{grad}\phi+{\bm{m}})\cdot{\bm{n}}^{\mathcal{B}}\,$
$\displaystyle=0$ on $\displaystyle{\partial}\widehat{{\mathcal{B}}}$
$\displaystyle\text{ in }\mathbb{R},$ (12) $\displaystyle\updelta\phi:$
$\displaystyle\quad(-\operatorname{grad}\phi)\cdot{\bm{n}}^{\Omega}\,$
$\displaystyle=0$ on $\displaystyle{\partial}\widehat{{\Omega}}$
$\displaystyle\text{ in }\mathbb{R},$ (13) $\displaystyle\updelta{\bm{m}}:$
$\displaystyle\quad{\bm{m}}\times(\operatorname{div}(\operatorname{grad}{\bm{m}})+\operatorname{grad}\phi)\,$
$\displaystyle={\bm{0}}$ in $\displaystyle\widehat{{\mathcal{B}}}$
$\displaystyle\text{ in }\mathbb{R}^{3},$ (14)
$\displaystyle\updelta{\bm{m}}:$
$\displaystyle\quad{\bm{m}}\times(\operatorname{grad}{\bm{m}}\,{\bm{n}})\,$
$\displaystyle={\bm{0}}$ on $\displaystyle{\partial}\widehat{{\mathcal{B}}}$
$\displaystyle\text{ in }\mathbb{R}^{3},$ (15)

where ${\bm{n}}^{\mathcal{B}}$ and ${\bm{n}}^{\Omega}$ denote the unit outward
normal vectors on ${\partial}{\mathcal{B}}$ and ${\partial}{\Omega}$,
respectively. The first two equations correspond to Gauss’s law of
magnetostatics for magnetic and non-magnetic matter, respectively, see Eq. 1.
The third and fourth equation describe associated boundary conditions. The
last two equations are related to the conformation of the magnetization, which
needs to adhere the unit-sphere constraint $\|{\bm{m}}\|=1$. The latter is
usually accomplished by defining its variation as
$\updelta{\bm{m}}=\updelta{\bm{\omega}}\times{\bm{m}}$, where
$\updelta{\bm{\omega}}$ is a virtual axial vector of rotation, here given by
the virtual spin of the magnetization. This then results in the usual static
form of the Landau-Lifshitz equation
${\bm{m}}\times{\bm{h}}_{{\mathrm{e}}{\mathrm{f}}{\mathrm{f}}}={\bm{0}}$.
Therefore, the variation of the magnetization $\updelta{\bm{m}}$ is
constrained to lie in the tangent bundle $T{\mathcal{S}}^{d-1}$ of the unit
sphere ${\mathcal{S}}^{d-1}$ related to ${\bm{m}}$ [36]. Consequently, the
corresponding Euler-Lagrange equations need to be fulfilled only in the
tangent bundle. Therefore, and for the sake of brevity, we omit the usual
notation ${\bm{m}}\times\ldots$ of this constraint and simply demand that the
corresponding Euler-Lagrange equations

$\displaystyle\updelta{\bm{m}}:$ $\displaystyle\leavevmode\nobreak\
\operatorname{div}(\operatorname{grad}{\bm{m}})+\operatorname{grad}\phi\,$
$\displaystyle={\bm{0}}$ in $\displaystyle\widehat{{\mathcal{B}}}\text{ in
}T{\mathcal{S}}^{2},$ (16) $\displaystyle\updelta{\bm{m}}:$
$\displaystyle\leavevmode\nobreak\
\quad(\operatorname{grad}{\bm{m}})\,{\bm{n}}$ $\displaystyle={\bm{0}}$ on
$\displaystyle{\partial}\widehat{{\mathcal{B}}}\text{ in }T{\mathcal{S}}^{2},$
(17)

have to be fulfilled in the tangent bundle of ${\mathcal{S}}^{2}$.

Eqs. 12 and 13 can be recast using Eq. 6 at the common surface of
$\widehat{{\mathcal{B}}}$ and the free space
$\widehat{{\Omega}}\setminus\widehat{{\mathcal{B}}}$ as the jump condition

$\displaystyle{[\kern-1.49994pt[}-\operatorname{grad}\phi+\widehat{m}_{s}{\bm{m}}{]\kern-1.49994pt]}\cdot{\bm{n}}=0\
\text{ in }\ \mathbb{R},$ (18)

where ${\bm{n}}$ can be either ${\bm{n}}^{\widehat{{\mathcal{B}}}}$ or
${\bm{n}}^{\widehat{{\Omega}}\setminus\widehat{{\mathcal{B}}}}$ since
${\bm{n}}^{\widehat{{\mathcal{B}}}}=-{\bm{n}}^{\widehat{{\Omega}}\setminus\widehat{{\mathcal{B}}}}$.
Above, ${[\kern-1.49994pt[}(\cdot){]\kern-1.49994pt]}$ denotes the jump across
the interface separating the body from the surrounding free space.

From a given magnetization, the scalar potential can be directly calculated.
It is the solution of Eqs. 10, 11, 12 and 13. This solution for the scalar
potential $\phi$ is also denoted as the fundamental solution of Laplace’s
equation in $\mathbb{R}^{3}$, see [18], since it boils down to
$\Delta\phi=\widehat{m}_{s}\operatorname{div}{\bm{m}}$.222Strictly speaking
$\Delta\phi=\widehat{m}_{s}\operatorname{div}{\bm{m}}$ is a Poisson equation
since the right-hand side is non-zero. It is given as

$\displaystyle\phi({\bm{x}})=\frac{1}{4\pi}\left(-\int_{\widehat{{\mathcal{B}}}^{\prime}}\frac{\operatorname{div}^{\prime}{\bm{m}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}v^{\prime}+\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}\frac{{\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\right),$
(19)

where $\operatorname{div}^{\prime}$ denotes the divergence w.r.t.
${\bm{x}}^{\prime}$. For a source of Eq. 19, see [18],[6, Eq. 6.3.48] or [8,
Eq. 3.46].

The stray-field energy can then be calculated as

$\displaystyle\widehat{{\Pi}}_{\text{stray}}=\int_{\widehat{{\mathcal{B}}}}\widehat{\psi}_{\text{stray}}(\phi,{\bm{m}})\,{\mathrm{d}}v=\frac{1}{2}\left(\int_{\widehat{{\mathcal{B}}}}\rho_{v}\phi_{v}+2\rho_{v}\phi_{a}\,{\mathrm{d}}v+\int_{{\partial}\widehat{{\mathcal{B}}}}\sigma_{s}\phi_{a}\,{\mathrm{d}}a\right),$
(20)

see [39, Eq. 5]. Here, $\rho_{v}=-\operatorname{div}{\bm{m}}$ and
$\sigma_{s}={\bm{m}}\cdot{\bm{n}}$. These two quantities are denoted as
_volume charges_ and _surface charges_. Additionally, $\phi_{v}$ denotes the
volume contribution in Eq. 19 and $\phi_{a}$ denotes the surface contribution.

Thus, we have

$\displaystyle{\widehat{{\Pi}}}_{\text{stray}}=\frac{1}{8\pi}\Biggl{(}$
$\displaystyle\int_{\widehat{{\mathcal{B}}}}\int_{\widehat{{\mathcal{B}}}^{\prime}}\frac{\operatorname{div}{\bm{m}}({\bm{x}})\operatorname{div}^{\prime}{\bm{m}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}v^{\prime}\,{\mathrm{d}}v$
(21) $\displaystyle-$ $\displaystyle
2\int_{{\partial}\widehat{{\mathcal{B}}}}\int_{\widehat{{\mathcal{B}}}^{\prime}}\frac{\operatorname{div}^{\prime}{\bm{m}}({\bm{x}}^{\prime}){\bm{m}}({\bm{x}})\cdot{\bm{n}}({\bm{x}})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}v^{\prime}\,{\mathrm{d}}a+\int_{{\partial}\widehat{{\mathcal{B}}}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}\frac{{\bm{m}}({\bm{x}})\cdot{\bm{n}}({\bm{x}}){\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a\Biggr{)},$

where the double integral reflects the highly non-local nature of the stray-
field energy. Using the divergence theorem in the form $\int
g\operatorname{div}{\bm{f}}\,{\mathrm{d}}v=\int
g{\bm{f}}\cdot{\bm{n}}\,{\mathrm{d}}a-\int{\bm{f}}\cdot(\operatorname{grad}g)\,{\mathrm{d}}v$
yields the following representation of the stray energy

$\displaystyle{\widehat{{\Pi}}}_{\text{stray}}$
$\displaystyle=\frac{1}{8\pi}\Biggl{(}\int_{\widehat{{\mathcal{B}}}}\int_{\widehat{{\mathcal{B}}}^{\prime}}{\bm{m}}({\bm{x}})\cdot({\bm{x}}-{\bm{x}}^{\prime})\frac{\operatorname{div}^{\prime}{\bm{m}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}^{3}}\,{\mathrm{d}}v^{\prime}-\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}{\bm{m}}({\bm{x}})\cdot({\bm{x}}-{\bm{x}}^{\prime})\frac{{\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}^{3}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}v\Biggr{)},$
(22)

which can be compared to Eq. 19 and then translated to the well-known energy
expression

$\displaystyle{\widehat{{\Pi}}}_{\text{stray}}$
$\displaystyle=\frac{1}{2}\int_{\widehat{{\mathcal{B}}}}\operatorname{grad}\phi({\bm{m}})\cdot{\bm{m}}\,{\mathrm{d}}v=-\frac{1}{2}\int_{\widehat{{\mathcal{B}}}}{\bm{h}}({\bm{m}})\cdot{\bm{m}}\,{\mathrm{d}}v.$
(23)

The variation of the stray energy is

$\displaystyle\updelta{\widehat{{\Pi}}}_{\text{stray}}$
$\displaystyle=-\frac{1}{2}\int_{\widehat{{\mathcal{B}}}}\updelta{\bm{h}}({\bm{m}})\cdot{\bm{m}}+{\bm{h}}({\bm{m}})\cdot\updelta{\bm{m}}\,{\mathrm{d}}v$
(24)

and the first part of this variation reads

$\displaystyle\updelta{\bm{h}}({\bm{m}})$
$\displaystyle=\frac{1}{4\pi}\Biggl{(}\int_{\widehat{{\mathcal{B}}}^{\prime}}({\bm{x}}-{\bm{x}}^{\prime})\frac{\operatorname{div}^{\prime}\updelta{\bm{m}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}^{3}}\,{\mathrm{d}}v^{\prime}-\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}({\bm{x}}-{\bm{x}}^{\prime})\frac{\updelta{\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}^{3}}\,{\mathrm{d}}a^{\prime}\Biggr{)}.$
(25)

Exchanging the variables ${\bm{x}}$ and ${\bm{x}}^{\prime}$ introduces a minus
sign, so that we have
$\updelta{\bm{h}}({\bm{m}})={\bm{h}}({\bm{m}})\cdot\updelta{\bm{m}}$. This is
also a direct consequence of the _reciprocity theorem_ , see [9, Eq. 3-48].
The variation simplifies to

$\displaystyle\updelta{\widehat{{\Pi}}}_{\text{stray}}$
$\displaystyle=-\int_{\widehat{{\mathcal{B}}}}{\bm{h}}({\bm{m}})\cdot\updelta{\bm{m}}\,{\mathrm{d}}v.$
(26)

Thus, we have shown that the Euler-Lagrange equations are indeed the same as
in Eq. 14. We therefore have now

$\displaystyle\updelta{\bm{m}}:$
$\displaystyle\quad\operatorname{div}(\operatorname{grad}{\bm{m}})+\operatorname{grad}\phi({\bm{m}})\,$
$\displaystyle={\bm{0}}$ in $\displaystyle\widehat{{\mathcal{B}}}\text{ in
}T{\mathcal{S}}^{d-1},$ (27) $\displaystyle\updelta{\bm{m}}:$
$\displaystyle\quad(\operatorname{grad}{\bm{m}})\,{\bm{n}}$
$\displaystyle={\bm{0}}$ on
$\displaystyle{\partial}\widehat{{\mathcal{B}}}\text{ in
}T{\mathcal{S}}^{d-1}.$ (28)

In the following, these definitions are specified to the case of a nanodot.

####  2.2.2 Application to nanodots

The Euler-Lagrange equations derived above are valid for any shape of
${\mathcal{B}}$ and can usually not be solved analytically. As will be shown
next, it is however possible to derive accurate, semi-analytical solutions for
the magnetization field of cylindrical specimens such as magnetic nanodots.

##### Curling mode.

For certain specimen shapes, the energetic minimum is represented by a single,
rotationally symmetric vortex state. In the following, we present an ansatz
that is tailor-made to reproduce this phenomenon. It is encoded as

$\displaystyle{\bm{m}}=\begin{bmatrix}-\sin\theta\sin f(\rho)\\\
\cos\theta\sin f(\rho)\\\ \cos f(\rho)\end{bmatrix}.$ (29)

Here, $\rho\in[0,\widehat{R}]$ is the radial coordinate and
$\theta\in[0,2\pi]$ is the angular coordinate. The dimensionless radius of the
nanodot specimen is denoted by $\widehat{R}$. We will denote the dimensionless
thickness of the cylinder as $\widehat{H}$ and later also use
$\delta=\widehat{H}/2$. The vector components in Eq. 29 refer to the usual
Cartesian coordinates $x,y$ and $z$. Furthermore, $f(\rho)$ is an initially
unknown function for the rotationally symmetric ansatz that will be specified
in Section 3.1. The assumed slenderness of the nanodot also justifies the fact
that the magnetization given in Eq. 29 is constant across the thickness.

Introducing Eq. 29 into the energy Eq. 5, the exchange energy is obtained as

$\displaystyle\widehat{{\Pi}}_{\text{ex}}^{\text{CM}}(\operatorname{grad}{\bm{m}})=\int_{\widehat{{\Omega}}}\frac{1}{2}\|\operatorname{grad}{\bm{m}}\|^{2}\,{\mathrm{d}}V=2\pi\widehat{H}\int_{0}^{\widehat{R}}\frac{1}{2\rho}\left[{\left(\frac{{\partial}f(\rho)}{{\partial}\rho}\right)}^{2}\rho^{2}+\sin^{2}f(\rho)\right]\,{\mathrm{d}}\rho,$
(30)

where the superscript ’CM’ stands for ’curling mode’. The latter
representation of $\widehat{{\Pi}}_{\text{ex}}$ can also be found in, e.g.,
[9, 47]. Since $\widehat{{\Pi}}_{\text{ex}}$ is independent of the coordinates
$z$ and $\theta$, preintegration in these directions was carried out, reducing
the integral to only one dimension.

We now use Eq. 21 to calculate the stray-field energy. Due to the specific
form of the magnetization Eq. 29, we have $\operatorname{div}{\bm{m}}=0$, and
are left with the double integral over the surfaces. Since the direction of
magnetization according to Eq. 29 is tangential to the lateral surface, we
also have ${\bm{m}}\cdot{\bm{n}}=0$. Thus, only the contributions on the top
surface ${\mathcal{B}}^{\prime T}$ and the bottom surface
${\mathcal{B}}^{\prime B}$ at $z=\pm\delta=\pm\widehat{H}/2$ are nonzero. At
these surfaces we have ${\bm{m}}\cdot{\bm{n}}=\pm m_{z}=\pm\cos f(\rho)$ due
to ${\bm{n}}={[0,0,\pm 1]}^{T}$. From this it follows

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}({\bm{m}})$
$\displaystyle=\frac{1}{8\pi}\int_{{\partial}\widehat{{\mathcal{B}}}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}\frac{{\bm{m}}({\bm{x}})\cdot{\bm{n}}({\bm{x}}){\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a$
(31)
$\displaystyle=\frac{1}{8\pi}\left(\int_{{\partial}\widehat{{\mathcal{B}}}^{T}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime
T}}\frac{{\bm{m}}_{z}(\rho,\theta,\delta){\bm{m}}_{z}(\rho^{\prime},\theta^{\prime},\delta)}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a\right.$
$\displaystyle+\left.\int_{{\partial}\widehat{{\mathcal{B}}}^{B}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime
B}}\frac{(-{\bm{m}}_{z}(\rho,\theta,-\delta))(-{\bm{m}}_{z}(\rho^{\prime},\theta^{\prime},-\delta))}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a\right.$
$\displaystyle+2\left.\int_{{\partial}\widehat{{\mathcal{B}}}^{B}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime
T}}\frac{{\bm{m}}_{z}(\rho,\theta,\delta)(-{\bm{m}}_{z}(\rho^{\prime},\theta^{\prime},-\delta))}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a\right)$
$\displaystyle=\frac{1}{4\pi}\left(\int_{0}^{2\pi}\int_{0}^{\widehat{R}}\int_{0}^{2\pi}\int_{0}^{\widehat{R}}\frac{\cos
f(\rho)\cos
f(\rho^{\prime})}{d}\rho^{\prime}\,{\mathrm{d}}\rho^{\prime}\,{\mathrm{d}}\theta^{\prime}\rho\,{\mathrm{d}}\rho\,{\mathrm{d}}\theta\right.$
$\displaystyle+\left.\int_{0}^{2\pi}\int_{0}^{\widehat{R}}\int_{0}^{2\pi}\int_{0}^{\widehat{R}}\frac{-\cos
f(\rho)\cos
f(\rho^{\prime})}{\sqrt{d^{2}+4\delta^{2}}}\rho^{\prime}\,{\mathrm{d}}\rho^{\prime}\,{\mathrm{d}}\theta^{\prime}\rho\,{\mathrm{d}}\rho\,{\mathrm{d}}\theta\right),$

with
$d(\rho,\rho^{\prime},\theta,\theta^{\prime})^{2}=-2\rho\rho^{\prime}\cos(\theta-\theta^{\prime})+{(\rho^{\prime})}^{2}+\rho^{2}$.
This is derived by inserting the definitions
${\bm{x}}={[\rho\sin\theta,\rho\cos\theta,\pm\delta]}^{T}$ and
${\bm{x}}^{\prime}={[\rho^{\prime}\sin\theta^{\prime},\rho^{\prime}\cos\theta^{\prime},\pm\delta]}^{T}$.
Integration is accomplished with the usual Jacobian of polar coordinates, i.e.
$\int\int(\cdot)\,{\mathrm{d}}x\,{\mathrm{d}}y=\int\int(\cdot)\rho{\mathrm{d}}\theta{\mathrm{d}}\rho$.
Since the denominator can be made independent of $\theta$333By realizing
$\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{1}{\sqrt{-2\rho\rho^{\prime}\cos(\theta-\theta^{\prime})+{(\rho^{\prime})}^{2}+\rho^{2}}}\,{\mathrm{d}}\theta^{\prime}\,{\mathrm{d}}\theta=\int_{0}^{2\pi}\,{\mathrm{d}}\theta^{\prime}4\int_{0}^{\pi/2}\frac{1}{\sqrt{-4\rho\rho^{\prime}\sin^{2}(\theta)+({\rho^{\prime}}+\rho)^{2}}}\,{\mathrm{d}}\theta$.
, similar to [2, Eq. 19], one integration can be carried out directly and we
get

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}(f)$
$\displaystyle=2\left(\int_{0}^{\pi/2}\int_{0}^{\widehat{R}}\int_{0}^{\widehat{R}}\frac{\cos
f(\rho)\cos
f(\rho^{\prime})}{d}\rho^{\prime}\,{\mathrm{d}}\rho^{\prime}\rho\,{\mathrm{d}}\rho\,{\mathrm{d}}\theta-\frac{\cos
f(\rho)\cos
f(\rho^{\prime})}{\sqrt{d^{2}+4\delta^{2}}}\rho^{\prime}\,{\mathrm{d}}\rho^{\prime}\rho\,{\mathrm{d}}\rho\,{\mathrm{d}}\theta\right),$
(32)

with
$d(\rho,\rho^{\prime},\theta)^{2}=-4\rho\rho^{\prime}\sin^{2}(\theta)+{(\rho^{\prime}+\rho)}^{2}$.

Now, we can also integrate w.r.t. $\theta$. Both integrals in Eq. 32 have the
following form that can be integrated in terms of elliptic integrals,

$\displaystyle\int_{0}^{\pi/2}\frac{\,{\mathrm{d}}x}{\sqrt{a-b\sin^{2}(x)}}=\frac{1}{\sqrt{a}}\operatorname{K}\left(\frac{b}{a}\right),$
(33)

where $\operatorname{K}(m)$ is the _complete elliptic integral of the first
kind_ with the parameter $m=k^{2}$ and $k$ is the _elliptic modulus_. Thus, we
can simplify Eq. 32 as

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}(f)$
$\displaystyle=2\int_{0}^{\widehat{R}}\cos
f(\rho)\rho\int_{0}^{\widehat{R}}\cos
f(\rho^{\prime})\rho^{\prime}\left(\frac{1}{\widehat{a}}\operatorname{K}\left(\frac{b}{\widehat{a}^{2}}\right)-\frac{1}{\sqrt{\widehat{a}^{2}+4\delta^{2}}}\operatorname{K}\left(\frac{b}{\widehat{a}^{2}+4\delta^{2}}\right)\right){\mathrm{d}}\rho^{\prime}{\mathrm{d}}\rho,$
(34)

where we used $b=4\rho\rho^{\prime}$ and $\widehat{a}=\rho^{\prime}+\rho$.
Now, we further introduce the shortcut ${\mathcal{K}}$ for the symmetric
kernel

$\displaystyle{\mathcal{K}}(\rho,\rho^{\prime})=\frac{\rho\rho^{\prime}}{\widehat{a}}\operatorname{K}\left(\frac{b}{\widehat{a}^{2}}\right)-\frac{\rho\rho^{\prime}}{\sqrt{\widehat{a}^{2}+4\delta^{2}}}\operatorname{K}\left(\frac{b}{\widehat{a}^{2}+4\delta^{2}}\right).$
(35)

Finally, since $\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}(f)$ is symmetric
w.r.t. $\rho$ and $\rho^{\prime}$, we transform the integral from the
quadrilateral domain $[0,\widehat{R}]\times[0,\widehat{R}]$ to an integral
over the triangular domain $\rho\in[0,\widehat{R}],\rho^{\prime}\in[0,\rho]$
by multiplying with 2 and changing the interval of integration accordingly.
This yields

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}(f)$
$\displaystyle=4\int_{0}^{\widehat{R}}\cos f(\rho)\int_{0}^{\rho}\cos
f(\rho^{\prime}){\mathcal{K}}(\rho,\rho^{\prime}){\mathrm{d}}\rho^{\prime}{\mathrm{d}}\rho.$
(36)

The total energy, only in terms of $f$, is

$\displaystyle\widehat{{\Pi}}^{\text{CM}}(f)=4\int_{0}^{\widehat{R}}\cos
f(\rho)\int_{0}^{\rho}\cos
f(\rho^{\prime}){\mathcal{K}}(\rho,\rho^{\prime})\,{\mathrm{d}}\rho^{\prime}+2\pi\widehat{H}\frac{1}{2\rho}\left[{\left(\frac{{\partial}f(\rho)}{{\partial}\rho}\right)}^{2}\rho^{2}+\sin^{2}f(\rho)\right]\,{\mathrm{d}}\rho.$
(37)

This gives rise to the minimization principle

$\displaystyle\boxed{f^{\ast}=\arg\left\\{\inf_{f\in
W^{k,p}({\mathcal{P}},\mathbb{R})}\widehat{{\Pi}}^{\text{CM}}(f)\right\\}}$
(38)

with
$\widehat{{\Pi}}^{\text{CM}}(f)=\widehat{{\Pi}}^{\text{CM}}_{\text{ex}}(f)+\widehat{{\Pi}}^{\text{CM}}_{\text{stray}}(f)$
and ${\mathcal{P}}=[0,\widehat{R}]$.

The first variation of the individual energy contributions reads

$\displaystyle{\delta}\widehat{{\Pi}}_{\text{stray}}^{\text{CM}}(f)$
$\displaystyle=-8\int_{0}^{\widehat{R}}\sin f(\rho)\int_{0}^{\rho}\cos
f(\rho^{\prime}){\mathcal{K}}(\rho,\rho^{\prime}){\mathrm{d}}\rho^{\prime}{\delta}f(\rho){\mathrm{d}}\rho$
(39)

and

$\displaystyle\updelta\widehat{{\Pi}}_{\text{ex}}^{\text{CM}}(f)$
$\displaystyle=2\pi\widehat{H}\int_{0}^{\widehat{R}}\left(\frac{\sin
f(\rho)\cos
f(\rho)}{\rho}-\frac{{\partial}f(\rho)}{{\partial}\rho}-\frac{{\partial}^{2}f(\rho)}{{\partial}\rho^{2}}\rho\right){\delta}f(\rho){\mathrm{d}}\rho.$
(40)

The corresponding Euler-Lagrange equation is a non-autonomous, nonlinear
second order _integro-differential_ equation that reads

$\displaystyle\updelta$ $\displaystyle f:\sin
f(\rho)\left(-4\int_{0}^{\widehat{R}}\cos
f(\rho^{\prime}){\mathcal{K}}(\rho,\rho^{\prime})\,{\mathrm{d}}\rho^{\prime}+2\pi\widehat{H}\frac{\cos
f(\rho)}{\rho}\right)-2\pi\widehat{H}\left(\frac{{\partial}f(\rho)}{{\partial}\rho}+\frac{{\partial}^{2}f(\rho)}{{\partial}\rho^{2}}\rho\right)=0\quad$
$\displaystyle\forall\rho\in{\mathcal{P}},$ (41) $\displaystyle\updelta$
$\displaystyle f:\frac{{\partial}f(\rho)}{{\partial}\rho}=0\quad$
$\displaystyle\text{at }\rho=\widehat{R}.$

Here, the interval of the inner integral is again $\widehat{R}$ instead of
$\rho$, since Eq. 41 is not symmetric in $\rho$ and $\rho^{\prime}$ anymore.
This means that the variation of the demagnetization contains the factor $4$
instead of $8$. Parts of this ordinary differential equation can also be found
in the literature, specifically in [9, p. 97, eq. 6–19]. In the latter
reference, the stray energy is missing, but anisotropic effects and an
external applied field are taken into account. We are not aware of a closed-
form solution of the derived Eq. 41.

##### In-plane single-domain state.

The in-plane single-domain state of the nanodot is encoded in a homogeneous
magnetization throughout the nanodot. Note that this is not a realistic
assumption for cylindrical nanodots, since a perfectly homogeneous
magnetization can only be observed in ellipsoidal bodies [40]. Nevertheless,
we assume that for a thin specimen it is a reasonable approximation, see [30,
49] and [6, Ch. 6.1.2] for reference. Due to isotropy, we could assume any
orientation of magnetization in the $x$-$y$-plane. In what follows, we choose
${\bm{m}}={\bm{e}}_{x}$. Due to homogeneity, we further have
$\operatorname{grad}{\bm{m}}={\bm{0}}$, so that the exchange energy vanishes
identically and only the stray energy Eq. 21 has to be evaluated. Since
${\bm{m}}$ is tangential to the top and bottom surface, we only have to
evaluate Eq. 21 at the lateral surface of the cylinder. Again, the volume
terms vanish due to $\operatorname{div}{\bm{m}}=0$. Hence, with
${\bm{m}}\cdot{\bm{n}}=\cos\theta$ the stray energy at the lateral surface is

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{IP}}$
$\displaystyle=\frac{1}{8\pi}\int_{{\partial}\widehat{{\mathcal{B}}}}\int_{{\partial}\widehat{{\mathcal{B}}}^{\prime}}\frac{{\bm{m}}({\bm{x}})\cdot{\bm{n}}({\bm{x}}){\bm{m}}({\bm{x}}^{\prime})\cdot{\bm{n}}({\bm{x}}^{\prime})}{\|{\bm{x}}-{\bm{x}}^{\prime}\|_{2}}\,{\mathrm{d}}a^{\prime}\,{\mathrm{d}}a$
(42)
$\displaystyle=\frac{1}{8\pi}\int_{-\delta}^{\delta}\int_{0}^{2\pi}\int_{-\delta}^{\delta}\int_{0}^{2\pi}\frac{\cos(\theta)\cos(\theta^{\prime})}{\sqrt{{(z-z^{\prime})}^{2}-2\widehat{R}^{2}\cos(\theta-\theta^{\prime})+2\widehat{R}^{2}}}\widehat{R}^{2}\,{\mathrm{d}}z^{\prime}\,{\mathrm{d}}\theta^{\prime}\,{\mathrm{d}}\theta\,{\mathrm{d}}z,$

where ’IP’ stands for ’in plane’. Furthermore,
${\bm{x}}=[\widehat{R}\sin\theta,\widehat{R}\cos\theta,z]^{T}$,
${\bm{x}}^{\prime}=[\widehat{R}\sin\theta^{\prime},\widehat{R}\cos\theta^{\prime},z^{\prime}]^{T}$,
and the infinitesimal area element is
$\|\frac{\partial{\bm{n}}}{\partial\theta}\times\frac{\partial{\bm{n}}}{\partial
z}\|_{2}\,{\mathrm{d}}\theta\,{\mathrm{d}}z=\widehat{R}\,{\mathrm{d}}\theta\,{\mathrm{d}}z$.
Additionally, for the transformation of the denominator, we used
$\sin(\theta^{\prime})\sin(\theta)+\cos(\theta^{\prime})\cos(\theta)=\cos(\theta^{\prime}-\theta)$.

The integrals w.r.t. the thickness coordinates $z$ and $z^{\prime}$ can be
evaluated analytically. Thus, we get

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{IP}}$
$\displaystyle=-\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}\cos\theta\cos\theta^{\prime}(\delta\ln(-2\delta+q)-\delta\ln(2\delta+q)+q-\sqrt{-4\delta^{2}+q^{2}})\,{\mathrm{d}}\theta^{\prime}\,{\mathrm{d}}\theta,$
(43)

with $q^{2}=4\delta^{2}+2\widehat{R}^{2}(1-\cos(\theta-\theta^{\prime}))$. The
single integration w.r.t. $z^{\prime}$ was also shown in [11]. However, in
contrast to our case, the authors of [11] derived $\phi({\bm{x}})$ by
integrating the equivalent of Eq. 19 w.r.t. $z^{\prime}$ and
$\theta^{\prime}$.

Since we are interested in the energy, two further integrations remain.

We start with the integration w.r.t. $\theta^{\prime}$. First, we perform a
substitution of variables by using $x=\pi-(\theta^{\prime}-\theta)$ and get

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{IP}}$
$\displaystyle=-\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\cos\theta\int_{-\pi+\theta}^{\pi+\theta}\cos(\theta-x)\Bigg{(}\underbrace{\ln(2\delta+\sqrt{2\widehat{R}^{2}\cos
x+2\widehat{R}^{2}+4\delta^{2}})\delta}_{(1)}$ (44)
$\displaystyle\underbrace{-\ln(-2\delta+\sqrt{2\widehat{R}^{2}\cos
x+2\widehat{R}^{2}+4\delta^{2}})\delta}_{(2)}\underbrace{-\sqrt{2\widehat{R}^{2}\cos(x)+2\widehat{R}^{2}+4\delta^{2}}}_{(3)}\underbrace{+\sqrt{2}\widehat{R}\sqrt{\cos
x+1}}_{(4)}\Bigg{)}\,{\mathrm{d}}x\,{\mathrm{d}}\theta,$

which yields a better separation of variables between $x$ and $\theta$. Next,
we evaluate this integral part by part, starting with the fourth underbraced
term. Here, direct integration succeeds with the computer algebra system Maple
[34],

$\displaystyle-\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\cos\theta\int_{-\pi+\theta}^{\pi+\theta}\cos(\theta-x)\sqrt{2}\widehat{R}\sqrt{\cos
x+1}\,{\mathrm{d}}x\,{\mathrm{d}}\theta=-\frac{2}{3\pi}\widehat{R}^{3}\int_{0}^{2\pi}\cos^{2}\theta\,{\mathrm{d}}\theta=-\frac{2}{3}\widehat{R}^{3}.$
(45)

The third term can also be calculated with Maple as

$\displaystyle\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\cos\theta\int_{-\pi+\theta}^{\pi+\theta}\cos(\theta-x)\sqrt{2\widehat{R}^{2}\cos(x)+2\widehat{R}^{2}+4\delta^{2}}\,{\mathrm{d}}x\,{\mathrm{d}}\theta$
(46)
$\displaystyle=\frac{2}{3\pi}\left((R^{2}+2\delta^{2})\operatorname{E}(-\beta^{2})-2\operatorname{K}(-\beta^{2})(R^{2}+\delta^{2})\right)\cos^{2}\theta\delta,$

where we introduced the factor of slenderness $\beta=\widehat{R}/\delta$.
Additionally, we used

$\displaystyle\operatorname{E}(m)=\int_{0}^{\pi/2}\sqrt{1-m\sin^{2}(x)}\,{\mathrm{d}}x,$
(47)

where $\operatorname{E}$ is called _complete elliptic integral of the second
kind_. For the evaluation of the second term in Eq. 44, we need to apply
integration by parts, where we use $\cos(\theta-x)$ as one factor, which can
be integrated in a straight-forward manner. This yields

$\displaystyle\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\cos\theta\int_{-\pi+\theta}^{\pi+\theta}\cos(\theta-x)\ln(-2\delta+d(x))\delta\,{\mathrm{d}}x\,{\mathrm{d}}\theta=-\int_{0}^{2\pi}\int_{-\pi+\theta}^{\pi+\theta}\frac{\cos(\theta)\widehat{R}^{4}\sin(\theta-x)\sin
x\delta}{4\pi
d(x)\left(-2\delta+d(x)\right)}\,{\mathrm{d}}x\,{\mathrm{d}}\theta,$ (48)

with $d(x)=\sqrt{2\widehat{R}^{2}\cos x+2\widehat{R}^{2}+4\delta^{2}}$. Here,
the boundary terms vanish, since $\sin(-(-\pi+\theta)+\theta)=0$ and
$\sin(-(\pi+\theta)+\theta)=0$.

Now, by expanding the numerator with the angle difference identity
$\sin(\theta-x)=\sin\theta\cos x-\cos\theta\sin x$, we arrive at two terms
that have separated variables

$\displaystyle-\int_{0}^{2\pi}\int_{-\pi+\theta}^{\pi+\theta}\frac{\cos\theta\widehat{R}^{4}\sin(\theta-x)\sin
x\delta}{4\pi
d(x)(-2\delta+d(x))}\,{\mathrm{d}}x\,{\mathrm{d}}\theta=-\frac{\widehat{R}^{4}\delta}{4\pi}\int_{0}^{2\pi}\int_{-\pi+\theta}^{\pi+\theta}\frac{\cos\theta\sin\theta\cos
x\sin
x-\cos^{2}\theta\sin^{2}x}{d(x)\left(-2\delta+d(x)\right)}\,{\mathrm{d}}x\,{\mathrm{d}}\theta.$
(49)

The first term in the numerator vanishes, since the integral of this part is
only a function $\cos(x)$ which yields the same value for $x=-\pi+\theta$ and
$x=\pi+\theta$. The second part however is non-zero and yields, again
evaluated with Maple,

$\displaystyle\frac{\widehat{R}^{4}\delta}{4\pi}\int_{0}^{2\pi}\int_{-\pi+\theta}^{\pi+\theta}\frac{\cos^{2}\theta\sin^{2}x}{d(x)\left(-2\delta+d(x)\right)}\,{\mathrm{d}}x\,{\mathrm{d}}\theta=\frac{\delta}{4\pi}(4(\widehat{R}^{2}+\delta^{2})\operatorname{K}(-\beta^{2})-4\operatorname{E}(-\beta^{2})\delta^{2}+\pi\widehat{R}^{2})\int_{0}^{2\pi}\cos^{2}\theta\,{\mathrm{d}}\theta$
(50)
$\displaystyle=\delta\left((\widehat{R}^{2}+\delta^{2})\operatorname{K}(-\beta^{2})-\operatorname{E}(-\beta^{2})\delta^{2}+\frac{\pi\widehat{R}^{2}}{4}\right).$

The first term can be evaluated similarly to the second one, i.e., via
integration by parts and by realizing the vanishing boundary terms. Then,
expanding the integral yields two terms. After integration, one of these terms
is a function of $\cos x$, so that it vanishes in consideration of the
integration limits. The remaining term is again evaluated with the help of
computer algebra [34], which yields

$\displaystyle-\frac{\widehat{R}^{2}}{4\pi}\int_{0}^{2\pi}\cos\theta\int_{-\pi+\theta}^{\pi+\theta}\cos(\theta-x)\ln(2\delta+\sqrt{2\widehat{R}^{2}\cos
x+2\widehat{R}^{2}+4\delta^{2}})\delta\,{\mathrm{d}}x\,{\mathrm{d}}\theta$
(51)
$\displaystyle=\frac{\delta}{4\pi}(4(\widehat{R}^{2}+\delta^{2})\operatorname{K}(-\beta^{2})-4\operatorname{E}(-\beta^{2})\delta^{2}-\pi\widehat{R}^{2})\int_{0}^{2\pi}\cos^{2}\theta\,{\mathrm{d}}\theta$
$\displaystyle=\delta\left((\widehat{R}^{2}+\delta^{2})\operatorname{K}(-\beta^{2})-\operatorname{E}(-\beta^{2})\delta^{2}-\frac{\pi\widehat{R}^{2}}{4}\right).$

In total, evaluation of the stray energy, Eq. 42, provides

$\displaystyle\widehat{{\Pi}}_{\text{stray}}^{\text{IP}}$
$\displaystyle=-\frac{2}{3}((\beta^{3}-\beta)\operatorname{E}(-\beta^{2})-(\beta^{3}+\beta)\operatorname{K}(-\beta^{2})+1)R^{3},$
(52)

which, to our best knowledge, cannot be found in this particular form in the
literature. We would however like to mention the demagnetization energy
derived in [46], which has a different representation, but which gives
numerically equivalent results. In contrast to
LABEL:{eq:totalEnergyStrayNanoDot_final}, the representation in [46] involves
hypergeometric functions but can be used for any homogeneous orientation of
magnetization. Furthermore, there are certain similarities to the
demagnetization factors calculated in [28] and to the scalar potential
calculated in [11].

##### Out-of-plane single-domain state.

The energy of the out-of-plane single-domain state can be directly derived
from the curling mode Eq. 29 by setting $f(\rho)=0$. Then, we end up with the
energy

$\displaystyle\widehat{{\Pi}}^{\text{OOP}}(f)=4\int_{0}^{\widehat{R}}\int_{0}^{\rho}{\mathcal{K}}(\rho,\rho^{\prime})\,{\mathrm{d}}\rho^{\prime}\,{\mathrm{d}}\rho,$
(53)

where ’OOP’ stands for ’out of plane’. The first part of the kernel, Eq. 35,
can be evaluated analytically,

$\displaystyle\widehat{{\Pi}}^{\text{OOP}}_{\text{stray}}=\frac{4R^{3}}{3}-\int_{0}^{\widehat{R}}\int_{0}^{\rho}\frac{4\rho\rho^{\prime}}{\sqrt{\widehat{a}^{2}+4\delta^{2}}}\operatorname{K}\left(\frac{b}{\widehat{a}^{2}+4\delta^{2}}\right)\,{\mathrm{d}}\rho^{\prime}\,{\mathrm{d}}\rho,$
(54)

and the second part is left for numerical integration. In [46], this term is
derived analytically using hypergeometric functions,

$\displaystyle\widehat{{\Pi}}^{\text{OOP}}_{\text{stray}}=\frac{R^{3}}{6}\left[-4+3\pi\sqrt{\tau^{2}+1}\,_{2}{\mathrm{F}}_{1}\left(-\frac{1}{2},\frac{3}{2};2;\frac{1}{\tau^{2}+1}\right)\right],$
(55)

where $\,{}_{2}{\mathrm{F}}_{1}$ is the hypergeometric function and
$\tau=1/\beta$.

##  3 Results and discussion

###  3.1 Profiles of the out-of-plane magnetization

We now use Eq. 52 and Eq. 37 to determine at what point the curling mode
becomes unstable. For this, we have to take three different modes into
account. These three modes are shown in Fig. 1.

Figure 1: The three minimizing configurations depending on the thickness-to-
radius ratio. On the left, the magnetization is homogeneously distributed
throughout the specimen and points in the direction of the cylinder’s axis. In
the middle, the curling mode is shown with the magnetization pointing out of
plane in the center of the vortex as discussed in Section 2.2.2. On the right,
the homogeneous, in-plane magnetization as discussed in Section 2.2.2 is
shown.

The left and the center one are entirely captured by the curling mode Eq. 29.
For the right one, with in-plane constant magnetization, the magnetostatic
energy was derived in Section 2.2.2. The (sudden) switch between these three
states is of particular interest and arises as an instability phenomenon.

To study the minimizers of the curling mode, we use Ritz’s method. This
distinguishes our approach from classical models in the literature, from which
many are based on the approach suggested by [47] or slight modifications
thereof. For completeness, we write down the ansatz for the out-of-plane
magnetization given in [47]

$\displaystyle
m_{z}(\rho)=\begin{cases}\sqrt{1-\frac{2a\rho}{a^{2}+\rho^{2}}},&0\leq\rho\leq
a,\\\ 0,&a\leq\rho\leq\widehat{R},\end{cases}$ (56)

where the size of the domain with non-zero out-of-plane magnetization $a$ is
derived by minimizing the energy w.r.t. $a$. We note that the above function
for the magnetization is non-differentiable at $\rho=a$. In other words, it
has a non-physical kink precisely at this position.

Since we want to arrive at the exact minimizers of Eq. 37, at least
asymptotically, we use a Fourier series expansion given by

$\displaystyle\widehat{f}(\rho)=\sum_{i=1}^{n}a_{i}\sin\left(\rho\frac{(2i-1)\pi}{\widehat{R}/2}\right)$
(57)

and assume that $\widehat{f}(\rho)\in V\subset
W^{k,p}({\mathcal{P}},\mathbb{R})$. The series expansion includes only terms
that satisfy the boundary conditions $f(0)=0$ and
$\frac{{\partial}f(\rho)}{{\partial}\rho}|_{\rho=\widehat{R}}=0$ for any
$a_{i}$.The first condition follows from circular symmetry, requiring the
magnetization to point in $z$-direction at $\rho=0$. The second condition
follows from the Euler-Lagrange equation at the boundary, see Eq. 15.

For the numerical solution of the nonlinear problem of energy minimization,
the trust region method [13] is used. The gradient and the Hessian of the
energy w.r.t. $a_{i}$ are calculated using automatic differentiation in our
C++-implementation using [31].

The obtained profiles of the out-of-plane magnetization $m_{z}=\cos f(\rho)$
are shown in Fig. 2 for several nanodot radii and a fixed thickness
$\widehat{H}=3\sqrt{2}$. For all computations, we increase the number of
coefficients $a_{i}$ until the change in the minimizing energy is below
$1e-6$. For reference, $\widehat{{\Pi}}^{\text{CM}}$ varies from $33.54$ to
$153.91$ for geometries $\\{\widehat{R};\widehat{H}\\}$ given by
$\\{3.28;5.95\\}$ and $\\{14.14;14.14\\}$, respectively, where the smaller
energy is associated with the separated red dot in the phase diagram shown in
Fig. 3 (please refer to the zoomed region for a detailed view). For the raw
data, see [37]. The energy decrease of $1e-6$ is usually captured with $i$
between 3 and 25. The results are in qualitative agreement with the results
obtained from micromagnetic simulations documented in [23, Fig. 2] and [42,
Fig. 3]. In particular, the overshooting to the negative regime $M_{z}<0$ in
the periphery of the vortex is captured. Furthermore, the decay
$M_{z}(\widehat{R})\rightarrow 0$ for $\widehat{R}\rightarrow\infty$ is
foreseeable.

Figure 2: Profiles of the $z$-component of the magnetization ${\bm{m}}$ for
different nanodot radii and a fixed height $\widehat{H}=3\sqrt{2}$ at the
minimizing energy configuration. $\widehat{H}=3\sqrt{2}$ corresponds to three
times the exchange length $\ell_{\text{ex}}$. The plots were obtained using
Eq. 29. The horizontal axis is scaled such that the unit is in terms of the
exchange length $\ell_{\text{ex}}$.

###  3.2 Phase diagram

To show the favourable minimizers in dependence of the nanodot’s thickness and
radius, we sampled the minimizer of Eq. 37 with our ansatz Eq. 57 and compared
it with Eq. 52. The results of this procedure were then assembled in a phase
diagram taking into account the three different configurations given in Fig.
1, see Fig. 3. Here, the colors indicate the respective mode that provides the
lowest energy minimum. The similarity of our results with [10, Fig. 2] w.r.t.
the border between the vortex configuration and the in-plane single-domain
state is obvious. However, [10] takes the full three-dimensional problem into
account.

Motivated by [44], we also show the sampled data in a recast format in the
same figure. Comparing the latter representation with [44, Fig. 6] indicates
that in particular the border between the out-of-plane and the in-plane
configuration as well as the border between the vortex and the in-plane
configuration are nicely captured.

We note that the border between the out-of-plane configuration and the in-
plane configuration can be determined exactly by comparing the corresponding
energies Eq. 55 and Eq. 52 such that

$\displaystyle\frac{\widehat{{\Pi}}^{\text{IP}}_{\text{stray}}}{\widehat{{\Pi}}^{\text{OOP}}_{\text{stray}}}\stackrel{{\scriptstyle!}}{{=}}1.$
(58)

This leads to the equation

$\displaystyle
3\pi\sqrt{\tau^{2}+1}\,_{2}{\mathrm{F}}_{1}\left(-\frac{1}{2},\frac{3}{2};2;\frac{1}{\tau^{2}+1}\right)-2\pi\tau=4,$
(59)

which can be solved numerically for $\tau=\widehat{H}/(2\widehat{R})=1/\beta$,
providing $\tau=0.9064\,7615\,0066\,0646$. This result was also derived in [3]
by comparing the corresponding demagnetization factors of the two
configurations. The value for $\tau$ can be also found explicitly or
implicitly in, e.g., [46, 44, 22]. Note that the value for $\tau$ can also be
nicely seen in Fig. 3, where it separates the in-plane and the out-of-plane
single-domain states.

For completeness, we mention that the border between the out-of-plane single-
domain state and the vortex mode is quite different to the one reported by
Ross et al. [44]. This can be explained by the fact that in [44] a three-
dimensional micromagnetic simulation is carried out, which captures the out-
of-plane flower state. In addition, the experimental data documented in [44]
has been obtained from arrays of almost cylindrical nanodots that are densely
packed. Since we take into account individual nanodots with perfectly
cylindrical shape, a quantitative comparison seems implausible. As a final
note, we would like to mention the apparent similarities with the simulation
results shown in [12, Fig. 9a].

Figure 3: Results obatined by minimizing Eq. 37 with the ansatz Eq. 57 and by
comparing the minimizing energy to the energy of the in-plane state Eq. 52. We
sampled w.r.t. the radius and the height by using a grid of $100\times 100$
data points. Additionally, we sub-sampled the at a height of $4.2$ and a
radius of $2.2$ to capture the character of the boundary separating the
individual states with higher precision. Each of the samples is colored
depending on its minimizing configuration type given in Fig. 1. Here, green
refers to a homogeneous in-plane magnetization, red refers to a curling mode,
and yellow refers to a homogeneous, out-of-plane magnetization.

##  4 Summary and outlook

We investigated the variational formulation and corresponding minimizing
energies for typical magnetization states of thin cylindrical magnetic
nanodots. For that, we considered both the exchange and the demagnetization
energy and used Ritz’s method in conjunction with a Fourier series expansion
to calculate the energy minimizers. By comparing the minimizers of the vortex
mode to the energies of in- and out-of-plane single-domain states, we derived
a phase diagram of the nanodot and compared it to data obtained from two- and
three-dimensional models from the literature. Our results allowed us to
determine the critical radius at which the vortex mode becomes unfavorable
with arbitrary precision. Current research is devoted to establishing a
sophisticated, three-dimensional, micromagnetic, numerical simulation
framework that will be benchmarked against the results that have been obtained
in the present study.

## Data availability statement

The data and the scripts for creating the data are available via [37]. The
software Ikarus is available at [38].

## Acknowledgements

We gratefully acknowledge the support for this work by Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC 2075 – 390740016.

## References

  * Abo et al. [2013] Abo, G.S., Hong, Y.K., Park, J., Lee, J., Lee, W., Choi, B.C., 2013.  Definition of Magnetic Exchange Length.  IEEE Transactions on Magnetics 49, 4937–4939.  doi:10.1109/TMAG.2013.2258028.
  * Aharoni [1983] Aharoni, A., 1983.  Magnetostatic energy of a ferromagnetic cylinder.  Journal of Applied Physics 54, 488–492.  doi:10.1063/1.332100.
  * Aharoni [1989] Aharoni, A., 1989.  Single-domain ferromagnetic cylinder.  IEEE Transactions on Magnetics 25, 3470–3472.  doi:10.1109/20.42338.
  * Aharoni [1990a] Aharoni, A., 1990a.  Magnetostatics of curling in a finite cylinder.  Journal of Applied Physics 68, 255–258.  doi:10.1063/1.347125.
  * Aharoni [1990b] Aharoni, A., 1990b.  Upper bound to a single‐domain behavior of a ferromagnetic cylinder.  Journal of Applied Physics 68, 2892–2900.  URL: https://doi.org/10.1063/1.346422, doi:10.1063/1.346422, arXiv:https://doi.org/10.1063/1.346422.
  * Aharoni [1996] Aharoni, A., 1996.  Introduction to the Theory of Ferromagnetism.  Oxford : Clarendon Press ; New York : Oxford University Press.
  * Anirban [2021] Anirban, A., 2021.  Stable magnetic vortices.  Nature Reviews Physics 3, 4--4.  doi:10.1038/s42254-020-00270-6.
  * Bertotti [1998] Bertotti, G., 1998.  Maxwell’s Equations in Magnetic Media, in: Hysteresis in Magnetism. Elsevier, pp. 73--102.  doi:10.1016/B978-012093270-2/50052-0.
  * Brown [1963] Brown, W.F., 1963.  Micromagnetics.  18, interscience publishers.
  * Buda et al. [2002] Buda, L.D., Prejbeanu, I.L., Ebels, U., Ounadjela, K., 2002\.  Micromagnetic simulations of magnetisation in circular cobalt dots.  Computational Materials Science 24, 181--185.  doi:10.1016/S0927-0256(02)00184-2.
  * Caciagli et al. [2018] Caciagli, A., Baars, R.J., Philipse, A.P., Kuipers, B.W., 2018\.  Exact expression for the magnetic field of a finite cylinder with arbitrary uniform magnetization.  Journal of Magnetism and Magnetic Materials 456, 423--432.  doi:10.1016/j.jmmm.2018.02.003.
  * Chung et al. [2010] Chung, S.H., McMichael, R.D., Pierce, D.T., Unguris, J., 2010\.  Phase diagram of magnetic nanodisks measured by scanning electron microscopy with polarization analysis.  Phys. Rev. B 81, 024410\.  URL: https://link.aps.org/doi/10.1103/PhysRevB.81.024410, doi:10.1103/PhysRevB.81.024410.
  * Conn et al. [2000] Conn, A.R., Gould, N.I.M., Toint, P.L., 2000.  Trust Region Methods.  Society for Industrial and Applied Mathematics.  URL: https://epubs.siam.org/doi/abs/10.1137/1.9780898719857, doi:10.1137/1.9780898719857, arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9780898719857.
  * Cowburn [2002] Cowburn, R.P., 2002.  Magnetic nanodots for device applications.  Journal of Magnetism and Magnetic Materials 242--245, 505--511.  doi:10.1016/S0304-8853(01)01086-1.
  * Cowburn [2007] Cowburn, R.P., 2007.  Change of direction.  Nature Materials 6, 255--256.  doi:10.1038/nmat1877.
  * Cowburn et al. [1999] Cowburn, R.P., Koltsov, D.K., Adeyeye, A.O., Welland, M.E., Tricker, D.M., 1999.  Single-Domain Circular Nanomagnets.  Physical Review Letters 83, 1042--1045.  doi:10.1103/PhysRevLett.83.1042.
  * Dao et al. [2001] Dao, N., Whittenburg, S.L., Cowburn, R.P., 2001.  Micromagnetics simulation of deep-submicron supermalloy disks.  Journal of Applied Physics 90, 5235--5237.  doi:10.1063/1.1412838.
  * Di Fratta et al. [2020] Di Fratta, G., Muratov, C.B., Rybakov, F.N., Slastikov, V.V., 2020\.  Variational principles of micromagnetics revisited.  SIAM Journal on Mathematical Analysis 52, 3580--3599.  URL: https://doi.org/10.1137/19M1261365, doi:10.1137/19M1261365, arXiv:https://doi.org/10.1137/19M1261365.
  * Dorn and Wulfinghoff [2022] Dorn, C., Wulfinghoff, S., 2022\.  Computing magnetic noise with micro-magneto-mechanical simulations.  IEEE Transactions on Magnetics , 1--1doi:10.1109/TMAG.2022.3212764.
  * Gouva et al. [1989] Gouva, M.E., Wysin, G.M., Bishop, A.R., Mertens, F.G., 1989\.  Vortices in the classical two-dimensional anisotropic Heisenberg model.  Physical Review B 39, 11840--11849.  doi:10.1103/PhysRevB.39.11840.
  * Guslienko [2008] Guslienko, K.Y., 2008.  Magnetic Vortex State Stability, Reversal and Dynamics in Restricted Geometries.  Journal of Nanoscience and Nanotechnology 8, 2745--2760.  doi:10.1166/jnn.2008.18305.
  * Guslienko et al. [2000] Guslienko, K.Y., Choe, S.B., Shin, S.C., 2000.  Reorientational magnetic transition in high-density arrays of single-domain dots.  Applied Physics Letters 76, 3609--3611.  URL: https://doi.org/10.1063/1.126722, doi:10.1063/1.126722, arXiv:https://doi.org/10.1063/1.126722.
  * Guslienko and Novosad [2004] Guslienko, K.Y., Novosad, V., 2004\.  Vortex state stability in soft magnetic cylindrical nanodots.  Journal of Applied Physics 96, 4451--4455.  doi:10.1063/1.1793327.
  * Guslienko et al. [2001] Guslienko, K.Y., Novosad, V., Otani, Y., Shima, H., Fukamichi, K., 2001.  Magnetization reversal due to vortex nucleation, displacement, and annihilation in submicron ferromagnetic dot arrays.  Physical Review B 65, 024414\.  doi:10.1103/PhysRevB.65.024414.
  * Hehn et al. [1996] Hehn, M., Ounadjela, K., Bucher, J.P., Rousseaux, F., Decanini, D., Bartenlian, B., Chappert, C., 1996.  Nanoscale magnetic domains in mesoscopic magnets.  Science 272, 1782--1785.  URL: https://www.science.org/doi/abs/10.1126/science.272.5269.1782, doi:10.1126/science.272.5269.1782, arXiv:https://www.science.org/doi/pdf/10.1126/science.272.5269.1782.
  * Hubert and Schäfer [1998] Hubert, A., Schäfer, R., 1998\.  Magnetic Domains: The Analysis of Magnetic Microstructures.  Springer, Berlin ; New York.  URL: https://doi.org/10.1007/978-3-540-85054-0, doi:10.1007/978-3-540-85054-0.
  * Ishii and Sato [1989] Ishii, Y., Sato, M., 1989.  Magnetization curling in a finite cylinder.  Journal of Applied Physics 65, 3146--3150.  doi:10.1063/1.342712.
  * Joseph [1966] Joseph, R.I., 1966.  Ballistic demagnetizing factor in uniformly magnetized cylinders.  Journal of Applied Physics 37, 4639--4643.  URL: https://doi.org/10.1063/1.1708110, doi:10.1063/1.1708110, arXiv:https://doi.org/10.1063/1.1708110.
  * Keip and Sridhar [2019] Keip, M.A., Sridhar, A., 2019\.  A variationally consistent phase-field approach for micro-magnetic domain evolution at finite deformations.  Journal of the Mechanics and Physics of Solids 125, 805--824.  doi:10.1016/j.jmps.2018.11.012.
  * Kobayashi and Ishikawa [1992] Kobayashi, M., Ishikawa, Y., 1992\.  Surface magnetic charge distributions and demagnetizing factors of circular cylinders.  IEEE Transactions on Magnetics 28, 1810--1814.  doi:10.1109/20.141290.
  * Leal et al. [2018] Leal, A.M.M., et al., 2018.  Autodiff, a modern, fast and expressive C++ library for automatic differentiation.  https://autodiff.github.io.  URL: https://autodiff.github.io.
  * Lee et al. [2008] Lee, K.S., Kim, S.K., Yu, Y.S., Choi, Y.S., Guslienko, K.Y., Jung, H., Fischer, P., 2008\.  Universal Criterion and Phase Diagram for Switching a Magnetic Vortex Core in Soft Magnetic Nanodots.  Physical Review Letters 101, 267206\.  doi:10.1103/PhysRevLett.101.267206.
  * Manzin et al. [2021] Manzin, A., Ferrero, R., Vicentini, M., 2021.  From micromagnetic to in silico modeling of magnetic nanodisks for hyperthermia applications.  Advanced Theory and Simulations 4, 2100013.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/adts.202100013, doi:10.1002/adts.202100013, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adts.202100013.
  * Maplesoft, a division of Waterloo Maple Inc. [2018] Maplesoft, a division of Waterloo Maple Inc., 2018\.  Maple.  URL: https://hadoop.apache.org.
  * Metlov and Guslienko [2002] Metlov, K.L., Guslienko, K.Y., 2002\.  Stability of magnetic vortex in soft magnetic nano-sized circular cylinder.  Journal of Magnetism and Magnetic Materials 242--245, 1015--1017.  doi:10.1016/S0304-8853(01)01360-9.
  * Miehe and Ethiraj [2012] Miehe, C., Ethiraj, G., 2012\.  A geometrically consistent incremental variational formulation for phase field models in micromagnetics.  Computer Methods in Applied Mechanics and Engineering 245-246, 331--347.  URL: https://www.sciencedirect.com/science/article/pii/S0045782512000977, doi:https://doi.org/10.1016/j.cma.2012.03.021.
  * Müller [2023] Müller, A., 2023.  Thin cylindrical magnetic nanodots revisited: Scripts and data.  URL: https://doi.org/10.18419/darus-3103, doi:10.18419/darus-3103.
  * Müller and Vinod Kumar Mitruka [2023] Müller, A., Vinod Kumar Mitruka, T.K.M., 2023\.  Ikarus v0.3.  URL: https://doi.org/10.18419/darus-3303, doi:10.18419/darus-3303.
  * Nonaka et al. [1985] Nonaka, K., Hirono, S., Hatakeyama, I., 1985.  Magnetostatic energy of magnetic thin-film edge having volume and surface charges.  Journal of Applied Physics 58, 1610--1614.  doi:10.1063/1.336049.
  * Osborn [1945] Osborn, J.A., 1945.  Demagnetizing factors of the general ellipsoid.  Phys. Rev. 67, 351--357.  URL: https://link.aps.org/doi/10.1103/PhysRev.67.351, doi:10.1103/PhysRev.67.351.
  * Pigeau [2012] Pigeau, B., 2012.  Magnetic vortex dynamics nanostructures.  Theses. Université Paris Sud - Paris XI.  URL: https://theses.hal.science/tel-00779597.
  * Raabe et al. [2000] Raabe, J., Pulwey, R., Sattler, R., Schweinböck, T., Zweck, J., Weiss, D., 2000\.  Magnetization pattern of ferromagnetic nanodisks.  Journal of Applied Physics 88, 4437\.  doi:10.1063/1.1289216.
  * Reichel et al. [2022] Reichel, M., Xu, B.X., Schröder, J., 2022.  A comparative study of finite element schemes for micromagnetic mechanically coupled simulations.  Journal of Applied Physics 132, 183903\.  URL: https://doi.org/10.1063/5.0105613, doi:10.1063/5.0105613, arXiv:https://doi.org/10.1063/5.0105613.
  * Ross et al. [2002] Ross, C.A., Hwang, M., Shima, M., Cheng, J.Y., Farhoud, M., Savas, T.A., Smith, H.I., Schwarzacher, W., Ross, F.M., Redjdal, M., Humphrey, F.B., 2002.  Micromagnetic behavior of electrodeposited cylinder arrays.  Physical Review B 65, 144417\.  doi:10.1103/PhysRevB.65.144417.
  * Schneider et al. [2000] Schneider, M., Hoffmann, H., Zweck, J., 2000.  Lorentz microscopy of circular ferromagnetic permalloy nanodisks.  Applied Physics Letters 77, 2909--2911.  doi:10.1063/1.1320465.
  * Tandon et al. [2004] Tandon, S., Beleggia, M., Zhu, Y., De Graef, M., 2004\.  On the computation of the demagnetization tensor for uniformly magnetized particles of arbitrary shape. Part I: Analytical approach.  Journal of Magnetism and Magnetic Materials 271, 9--26.  doi:10.1016/j.jmmm.2003.09.011.
  * Usov and Peschany [1993] Usov, N.A., Peschany, S.E., 1993\.  Magnetization curling in a fine cylindrical particle.  Journal of Magnetism and Magnetic Materials 118, L290--L294.  doi:10.1016/0304-8853(93)90428-5.
  * Wachowiak et al. [2002] Wachowiak, A., Wiebe, J., Bode, M., Pietzsch, O., Morgenstern, M., Wiesendanger, R., 2002\.  Direct Observation of Internal Spin Structure of Magnetic Vortex Cores.  Science (New York, N.Y.) 298, 577--80.  doi:10.1126/science.1075302.
  * Wysin [2015] Wysin, G.M., 2015.  Magnetic Excitations and Geometric Confinement Theory and Simulations.  IOP Publishing.  doi:10.1088/978-0-7503-1074-1.