# Skyrmion stacking in stray field-coupled ultrathin ferromagnetic multilayers

N J Dubicki, V V Slastikov, Anne Bernand-Mantel, C B Muratov

► To cite this version:

N J Dubicki, V V Slastikov, Anne Bernand-Mantel, C B Muratov. Skyrmion stacking in stray field-coupled ultrathin ferromagnetic multilayers. 2025. hal-05350918

HAL Id: hal-05350918

https://hal.science/hal-05350918v1

Preprint submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Skyrmion stacking in stray field-coupled ultrathin ferromagnetic multilayers

N. J. Dubicki*

V. V. Slastikov†

A. Bernand-Mantel‡

C. B. Muratov§¶

June 6, 2025

## Abstract

This paper explores the energy landscape of ferromagnetic multilayer heterostructures that feature magnetic skyrmions – tiny whirls of spins with non-trivial topology – in each magnetic layer. Such magnetic heterostructures have been recently pursued as possible hosts of room temperature stable magnetic skyrmions suitable for the next generation of low power information technologies and unconventional computing. The presence of stacked skyrmions in the adjacent layers gives rise to a strongly coupled nonlinear system, whereby the induced magnetic field plays a crucial stabilizing role. Starting with the micromagnetic modeling framework, we derive a general reduced energy functional for a fixed number of ultrathin ferromagnetic layers with perpendicular magnetocrystalline anisotropy. We next investigate this energy functional in the regime in which the energy is dominated by the intralayer exchange interaction and formally obtain a finite-dimensional description governed by the energy of a system of one skyrmion per layer as a function of the position, radius and the rotation angle of each of these skyrmions. For the latter, we prove that energy minimizers exist for all fixed skyrmion locations. We then focus on the simplest case of stray field-coupled ferromagnetic bilayers and completely characterize the energy minimizers. We show that the global energy minimizers exist and consist of two stray field-stabilized Néel skyrmions with antiparallel in-plane magnetization components. We also calculate the energy of two skyrmions of equal radius as a function of their separation distance.

## Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Micromagnetic model</td><td>2</td></tr><tr><td>3 Reduced model for ultrathin multilayers</td><td>5</td></tr><tr><td>4 Energy of N stacked skyrmions</td><td>9</td></tr><tr><td>5 The landscape of the energy function F<sub>N</sub></td><td>14</td></tr><tr><td>6 Application to the case of bilayers in the absence of DMI</td><td>18</td></tr><tr><td>A Appendix</td><td>25</td></tr><tr><td>References</td><td>29</td></tr></table>

<footnote>*Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, New Jersey 07102, USA.</footnote>

<footnote>†School of Mathematics, University of Bristol, Bristol BS8 1UG, United Kingdom.</footnote>

<footnote>‡CEMES, Université de Toulouse, CNRS, 29 Rue Jeanne Marvig, BP94347, 31055 Toulouse, France.</footnote>

<footnote>§Dipartimento di Matematica, Università di Pisa, Largo B. Pontecorvo, 5, 56127 Pisa, Italy.</footnote>

<footnote>¶Corresponding author: cyrill.muratov@unipi.it.</footnote>
# 1 Introduction

Magnetic heterostructures have recently garnered significant interest as promising platforms for hosting room-temperature stable magnetic skyrmions – tiny whirls of spins with non-trivial topology – making them ideal candidates for next-generation low-power information technologies and unconventional computing [19, 20, 39]. Magnetic skyrmions are particle-like topologically protected magnetic textures predicted to exist in ferromagnetic materials since the late 1980s [6–8]. They are typically stabilized by chiral Dzyaloshinskii-Moriya interaction (DMI) [15, 42] and have been observed in various magnetic materials, including chiral magnets, ferromagnetic ultrathin films and multilayers [24, 43, 47–49]. Magnetic skyrmions have been extensively investigated as potential information carriers for spintronic devices [18, 29, 45]. Such investigations are impossible without modeling and computational studies of these systems. Modeling of magnetic skyrmions often relies on two-dimensional models, a choice which may be justified in the case of bulk chiral materials [8], or for ultrathin magnetic films, when the film thickness is much smaller than the exchange length [4, 36].

In the last 10 years, there has been a growing body of skyrmion observations in multilayer heterostructures, whose intrinsic high tunability enables the observation of skyrmions at room temperature [41, 48]. These heterostructures may present a total thickness that can be larger than the exchange length, and, most importantly, the reduction of the interlayer exchange interaction due to the use of non-magnetic spacers facilitates the emergence of a thickness-dependent magnetization, leading to fully three-dimensional magnetization textures, as observed experimentally [14, 32, 40]. Such twisting of the magnetization in the direction perpendicular to the layers is a well known phenomena related to the influence of “magnetic charges” in ferromagnetic films with out-of-plane magnetocrystalline anisotropy. Indeed, the magnetic field induced by magnetic charges modifies the internal structure of domain walls, a phenomenon widely studied in the framework of bubble materials in the 1970’s [26, 46].

A similar phenomenon occurs in the case of multilayers with in-plane magnetocrystalline anisotropy in the absence of exchange coupling, i.e. in the case of purely stray field-coupled layers [22, 23]. For magnetic bilayers sufficiently thin to be in the Néel wall regime, the lower-energy state consists of two Néel walls with opposite rotations in both layers [23], partially cancelling the stray field. Similarly, in the case of magnetic bilayers with perpendicular magnetic anisotropy it has been observed experimentally that the Bloch wall, which is energetically favorable in a single layer, is replaced by two Néel walls with the opposite rotation senses [2]. This configuration can be further stabilized by the DMI, assuming a proper choice of asymmetric interfaces, as reported in the case of skyrmionic bubbles in bilayers with opposite DMI constants [25]. For larger number of repeats, the competition between the dipolar and the DMI interactions leads to twisted walls which present an asymmetry in the thickness direction with respect to the middle of the magnetic multilayer. This twist has been evidenced experimentally in the case of domain walls [17, 38] and skyrmionic bubbles [32]. The asymmetric twisted walls have also been the subject of detailed analytical modeling using domain wall and skyrmionic bubble *ansätze*, where the impact of twist on the spin-orbit torque dynamics of skyrmionic bubbles has been studied [33, 35].

In the present work, we focus on a study of compact magnetic skyrmions in ultrathin ferromagnetic multilayers, i.e., magnetic heterostructures whose total thickness lies below the Bloch wall width of the ferromagnetic material. We use the tools of asymptotic analysis to investigate multilayer systems with one skyrmion per layer in the case of purely stray field-coupled layers. We start with the micromagnetic energy functional for a three-dimensional multilayer system that includes the intralayer exchange, the magnetocrystalline anisotropy (both of bulk and interfacial origin), the Zeeman, the interfacial DMI and the full stray field interactions. We then derive a reduced energy functional in the case where each ferromagnetic layer is thin compared to the exchange length and the whole stack is thin compared to the characteristic length scale of variation of the magnetization in the film plane. We show that in addition to the usual local shape anisotropy term and the non-local dipolar interaction terms that are already present in the case of single ultrathin ferromagnetic layers [4, 13, 30, 44], the expansion of the stray field energy for multilayers leads to the appearance of a new local dipolar energy term corresponding to interlayer volume-surface interactions. This term cancels out in the case of identical magnetizations in the adjacent layers, but becomes equivalent to a stabilizing layer-dependent interfacial DMI when the in-plane magnetization components in two consecutive layers are opposite to one another.
The obtained reduced energy functional is investigated in the regime in which the energy is domi-
nated by the intralayer exchange interaction, which corresponds to the conformal limit studied in [5]
in the case of single ultrathin ferromagnetic layers. Following the arguments of [5], we formally ob-
tain a finite-dimensional description of the system of interacting compact magnetic skyrmions with
one skyrmion per layer as a function of the position, radius and the rotation angle of each of theses
skyrmions that is expected to be asymptotically exact in the considered limit. We then investigate
the energy landscape of the above system in the case of zero applied magnetic field and prove that
energy minimizers exist for all fixed skyrmion positions. The difficulty in obtaining such a result is
that a priori it is not clear whether the energy could not be reduced by some skyrmions shrinking to
zero radius. We prove that this phenomenon does not occur, if the skyrmion centers are fixed.

We then apply our results to the case of stray field-coupled ferromagnetic bilayers in the absence
of DMI, where we obtain a complete characterization of global energy minimizers in the intralayer
exchange-dominated regime. These minimizers consist of two concentric Néel skyrmions with anti-
parallel in-plane magnetization and a prescribed chirality. We also provide an expression for the
energy of two skyrmions of equal radius as a function of their separation distance and illustrate our
findings with a result of micromagnetic simulations. Notice that extending such a characterization to
the general case of multilayers presents a difficulty that the energy minimizing sequences could consist
of skyrmions in different layers moving far apart. For bilayers, we show that this cannot occur by an
explicit analysis of all the interaction terms, which, however, becomes intractable for higher numbers
of layers.

Our paper is organized as follows. In section 2, we specify the full three-dimensional micromag-
netic energy functional with all the relevant energy terms in the multilayer geometry and carry out its
non-dimensionalization. In section 3, we explicitly compute the energy of the magnetizations that are
independent of the thickness variable in each ferromagnetic layer. We then carry out an asymptotic
expansion of the energy as the ferromagnetic layer and interlayer thicknesses go to zero, with the num-
ber of layers fixed, to obtain a reduced two-dimensional variational model of ultrathin ferromagnetic
multilayers with perpendicular magnetic anisotropy. Then, in section 4 we introduce an ansatz in the
form of one truncated Belavin-Polyakov skyrmion in each ferromagnetic layer, which is asymptotically
valid in the intralayer exchange-dominated regime and asymptotically compute the finite-dimensional
energy function of such a skyrmion stack. In section 5, we then explore the basic properties of the
obtained energy function and prove that it admits a global energy minimizer for all fixed locations of
the skyrmion centers, see Theorem 1. Finally, in section 6 we completely characterize the global en-
ergy minimizer in the simplest non-trivial case of ultrathin stray field-coupled bilayers in the intralayer
exchange-dominated regime, see Theorem 2. We also corroborate our conclusions with the results of
micromagnetic simulations and characterize the interaction energy of two skyrmions of equal radius
separated by a prescribed distance.

**Acknowledgments.** N.J.D. and C.B.M. were partially supported by NSF via grant DMS-1908709. A. Bernand-Mantel was supported by France 2030 government investment plan managed by the French National Research Agency under grant reference PEPR SPIN [SPINTHEORY] ANR-22-EXSP- 0009 and grant NanoX ANR-17-EURE-0009 in the framework of the Programme des Investissements d'Avenir. C.B.M. is a member of INdAM-GNAMPA and acknowledges partial support by the MUR Excellence Department Project awarded to the Department of Mathematics, University of Pisa, CUP I57G22000700001, and by the PRIN 2022 PNRR Project P2022WJW9H. All analytical calculations in the paper involving special functions were carried out, using MATHEMATICA 14.2 software.

2 Micromagnetic model

Our starting point is the micromagnetic modeling framework (in the SI units), in which the observed
magnetization configurations are interpreted as local or global energy minimizers of a micromagnetic
energy functional evaluated on vector fields **M** : Ω → R³ of fixed length |**M**| = M_s [9, 27, 31]. Here
Ω ⊆ R³ is the spatial domain (an open set) occupied by a single centrosymmetric crystalline ferro-
magnetic material, **M** = **M**(**r**) is the magnetization vector at point **r** = (x, y, z) ∈ Ω, and M_s > 0
is the saturation magnetization (in A/m). The energy functional (in J) that contains the exchange, bulk magnetocrystalline anisotropy, Zeeman, stray field, interfacial magnetocrystalline anisotropy and interfacial DMI contributions, in that order, reads

$$ \mathcal{E}(\mathbf{M}) = \int_{\Omega} \left\{ \frac{A}{M_s^2} |\nabla \mathbf{M}|^2 + K_u \Phi_u \left( \frac{\mathbf{M}}{M_s} \right) - \mu_0 \mathbf{H}_a \cdot \mathbf{M} - \frac{\mu_0}{2} \mathbf{H}_d \cdot \mathbf{M} \right\} d^3 r \\ + \int_{\partial \Omega} \left\{ K_s \Phi_s \left( \frac{\mathbf{M}}{M_s}, \mathbf{r} \right) + \frac{D_s}{M_s^2} \left( M^{\|} \nabla_{\perp} \cdot \mathbf{M}^{\perp} - \mathbf{M}^{\perp} \cdot \nabla_{\perp} M^{\|} \right) \right\} d\mathcal{H}^2(\mathbf{r}). \quad (2.1) $$

Here $A > 0$ is the exchange stiffness (in J/m), $K_u \ge 0$ is the bulk magnetocrystalline anisotropy constant (in J/m³), $\Phi_u : \mathbb{S}^2 \to \mathbb{R}^+ = [0, \infty)$ specifies the dependence of the bulk anisotropy energy on the direction of $\mathbf{M}$, $\mu_0 = 4\pi \times 10^{-7}$ H/m is the vacuum permeability, $\mathbf{H}_a \in \mathbb{R}^3$ is the applied magnetic field in (A/m), and $\mathbf{H}_d : \mathbb{R}^3 \to \mathbb{R}^3$ is the demagnetizing field (in A/m) produced by $\mathbf{M}$ via the solution of the stationary Maxwell's equations

$$ \nabla \cdot (\mathbf{H}_d + \mathbf{M}) = 0, \quad \nabla \times \mathbf{H}_d = 0, \quad (2.2) $$

distributionally in $\mathbb{R}^3$, with $\mathbf{M}$ extended by zero in $\mathbb{R}^3 \setminus \bar{\Omega}$ [12]. The interfacial terms in the second line of (2.1) contain the interfacial DMI constant $D_s : \partial\Omega \to \mathbb{R}$ (in J/m), which may take different values depending on the adjacent non-magnetic material, the DMI energy term written in terms of $\mathbf{M} = (\mathbf{M}^{\perp}, M^{\|})$, where $M^{\|}$ is the component of $\mathbf{M}$ along the normal to $\partial\Omega$ and $\mathbf{M}^{\perp}$ is the tangential component to $\partial\Omega$, respectively, and $\nabla_{\perp}$ is the tangential gradient, $K_s : \partial\Omega \to \mathbb{R}^+$ is the interfacial magnetocrystalline anisotropy constant (in J/m²), which may also take different values depending on the adjacent non-magnetic material, and $\Phi_s : \mathbb{S}^2 \times \partial\Omega \to \mathbb{R}^+$ specifies the dependence of the interfacial anisotropy energy on the magnetization orientation relative to the crystalline axes in $\Omega$ and the normal to $\partial\Omega$. When $\Omega$ is unbounded, one often needs to subtract the contribution of a fixed reference configuration $\mathbf{M}^*$ from the integrand in (2.1) to make the resulting integrals convergent.

As usual, we introduce the normalized magnetization vector $\mathbf{m} : \Omega \to \mathbb{S}^2$, extended by zero to the rest of $\mathbb{R}^3$. In terms of $\mathbf{m}$ the energy becomes

$$ \mathcal{E}(M_s \mathbf{m}) = \int_{\Omega} (A |\nabla \mathbf{m}|^2 + K_u \Phi_u(\mathbf{m}) - 2K_d \mathbf{h} \cdot \mathbf{m} - K_d \mathbf{h}_d \cdot \mathbf{m}) d^3 r \\ + \int_{\partial \Omega} (K_s \Phi_s(\mathbf{m}, \mathbf{r}) + D_s(m^{\|} \nabla_{\perp} \cdot \mathbf{m}^{\perp} - \mathbf{m}^{\perp} \cdot \nabla m^{\|})) d\mathcal{H}^2(\mathbf{r}), \quad (2.3) $$

in which $K_d = \frac{1}{2}\mu_0 M_s^2$, $\mathbf{h} = \mathbf{H}_a/M_s$, $\mathbf{h}_d = \mathbf{H}_d/M_s$, $\mathbf{m} = (\mathbf{m}^{\perp}, m^{\|})$ on $\partial\Omega$ as before. Note that the demagnetizing field $\mathbf{h}_d$ satisfies [12]

$$ \mathbf{h}_d = -\nabla U, \quad \Delta U = \nabla \cdot \mathbf{m}, \quad (2.4) $$

distributionally in $\mathbb{R}^3$. In particular, the magnetostatic potential $U$ is given by

$$ U(\mathbf{r}) = - \int_{\mathbb{R}^3} \frac{\nabla \cdot \mathbf{m}(\mathbf{r}')}{4\pi |\mathbf{r} - \mathbf{r}'|} d^3 r', \quad \mathbf{r} \in \mathbb{R}^3. \quad (2.5) $$

The quantity $\rho_{\mathbf{m}} = -\nabla \cdot \mathbf{m}$ hence has the meaning of the magnetic charge density distributionally in $\mathbb{R}^3$, which includes the regular contribution of the volume charges from the absolutely continuous part of $\rho_{\mathbf{m}}$ in $\Omega$ and a singular contribution to $\rho_{\mathbf{m}}$ of the surface charges from the jump of $\mathbf{m}$ to zero across $\partial\Omega$.

We now specify the geometry of interest, which is that of a system of $N$ identical ferromagnetic layers of thickness $d$ separated by non-magnetic spacers of thickness $(a-1)d$ with $a > 1$, so that the total thickness of the magnetic layer plus the non-magnetic layer is $ad$:

$$ \Omega = \mathbb{R}^2 \times \bigcup_{n=1}^{N} [(n-1)ad, (n-1)ad + d], \quad (2.6) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematics of the geometry of a multilayer system. The heterostructure consists of N repeats of a sandwich of thickness ad in the form of a layer of one non-magnetic material (NM a), followed by a layer of a ferromagnet (FM) of thickness d, followed by a layer of another non-magnetic material (NM b), from the bottom to the top.</figcaption></figure>

see Fig. 1. The magnetic material is assumed to be uniaxial, with an easy axis perpendicular to the layers, so that $\Phi_u(\mathbf{m}) = |\mathbf{m}^{\perp}|^2$, where we now use the convention $\mathbf{m}(\mathbf{r}) = (\mathbf{m}^{\perp}(\mathbf{r}), m^{\parallel}(\mathbf{r}))$ with $\mathbf{m}^{\perp}(\mathbf{r}) \in \mathbb{R}^2$ and $m^{\parallel}(\mathbf{r}) \in \mathbb{R}$ being the in-plane and out-of-plane components of the magnetization, respectively, for all $\mathbf{r} \in \mathbb{R}^3$. Similarly, we assume that the interfacial anisotropy penalizes the tangential component of the magnetization. Hence, we take $\Phi_s(\mathbf{m}, \mathbf{r}) = |\mathbf{m}^{\perp}|^2$ as well. As is common in the spintronic multilayer materials [41, 47], the spacer is assumed to generally consist of two sublayers, so the upper surface of each magnetic layer is characterized by the interfacial DMI constant $D_s^+$ and the interfacial anisotropy constant $K_s^+$, while the bottom surface is characterized by the parameters $D_s^-$ and $K_s^-$, respectively. Notice that in the case of a single material spacer we simply set the parameters of the two surfaces equal to each other. Finally, we assume that the applied field is perpendicular to the film plane: $\mathbf{h} = \hat{h}\hat{z}$ for $\mathbf{h} \in \mathbb{R}$.

As can be easily seen from the above choices, for $K > K_d$, where

$$K = K_u + \frac{K_s^+ + K_s^-}{d} \qquad (2.7)$$

is the total magnetocrystalline anisotropy constant that takes into account the interfacial magnetocrystalline anisotropy, the uniform magnetization prefers to point out of the film plane. Indeed, if $\mathbf{m}$ is constant across $\Omega$, the solution of the Maxwell's equation in (2.2) is given by $\mathbf{h}_d = -\mathbf{m}$ in $\Omega$, resulting in the energy density $(K - K_d)|\mathbf{m}^{\perp}|^2$ minimized by $\mathbf{m} = \pm \hat{\mathbf{z}}$. Letting $\chi_n$ be the characteristic function of the n-th interval $I_n = ((n-1)ad, (n-1)ad+d) \subset \mathbb{R}$, we hence define the reference configuration in which the magnetization points downward in all the layers and is extended by zero to the whole of $\mathbb{R}^3$ as

$$\mathbf{m}^*(x, y, z) = -\sum_{n=1}^{N} \hat{\mathbf{z}}\chi_n(z), \qquad (2.8)$$

whose associated demagnetizing field is $\mathbf{h}_d^* = -\nabla U^* = -\mathbf{m}^*$, where $U^*$ is a bounded solution of (2.4) with $\mathbf{m} = \mathbf{m}^*$. Subtracting its contribution from the integrand in (2.3) and choosing the units of
energy and length to be given by Ad and the exchange length $\ell_{ex} = \sqrt{A/K_d}$, respectively, we can write the resulting renormalized energy functional (non-dimensionalized) as

$$ E(\mathbf{m}) = \frac{1}{\delta} \sum_{n=1}^{N} \int_{(n-1)a\delta}^{((n-1)a+1)\delta} \int_{\mathbb{R}^2} \left( |\nabla \mathbf{m}|^2 + Q_u |\mathbf{m}^{\perp}|^2 - 2h(m^{\parallel} + 1) + \mathbf{m} \cdot \nabla U - 1 \right) d^2r \, dz \\ + \sum_{n=1}^{N} \int_{\mathbb{R}^2} \left( Q_s^- |\mathbf{m}^{\perp}|^2 + 2\kappa^- \mathbf{m}^{\perp} \cdot \nabla m^{\parallel} \right) \Big|_{\{z=(n-1)a\delta\}} d^2r \\ + \sum_{n=1}^{N} \int_{\mathbb{R}^2} \left( Q_s^+ |\mathbf{m}^{\perp}|^2 - 2\kappa^+ \mathbf{m}^{\perp} \cdot \nabla m^{\parallel} \right) \Big|_{\{z=((n-1)a+1)\delta\}} d^2r, \qquad (2.9) $$

where we integrated the DMI by parts as in [4, 5] to make it more compact and mathematically well-behaved. Here

$$ \delta = \frac{d}{\ell_{ex}}, \quad Q_u = \frac{K_u}{K_d}, \quad Q_s^{\pm} = \frac{K_s^{\pm}}{dK_d}, \quad \kappa^{\pm} = \frac{D_s^{\pm}}{d\sqrt{AK_d}} \qquad (2.10) $$

are the dimensionless ferromagnetic layer thickness, the magnetocrystalline bulk and interfacial anisotropies' quality factors, and the dimensionless DMI strengths on the top and bottom surfaces of the film, respectively. Note that with this definition $E(\mathbf{m}^*) = 0$.

### 3 Reduced model for ultrathin multilayers

We now consider the situation in which each ferromagnetic layer is thin compared to the exchange length, i.e., the case $\delta \ll 1$, when the variations of the magnetization along the $z$-direction in each layer are highly penalized. In this case, following [21], it is appropriate to consider the magnetizations $\mathbf{m} = \tilde{\mathbf{m}}$ that do not vary across each layer and hence have the particular form

$$ \tilde{\mathbf{m}}(x, y, z) = \sum_{n=1}^{N} \mathbf{m}_n(x, y) \chi_n(z) = \sum_{n=1}^{N} \left( \mathbf{m}_n^{\perp}(x, y) + \hat{\mathbf{z}} m_n^{\parallel}(x, y) \right) \chi_n(z), \qquad (3.1) $$

where $\mathbf{m}_n : \mathbb{R}^2 \to \mathbb{S}^2$ for each $n = 1, \dots, N$, $\mathbf{m}_n^{\perp}$ is treated as an in-plane vector in $\mathbb{R}^3$, and now

$$ \chi_n(z) = \begin{cases} 1 & \text{if } z \in [(n-1)a\delta, ((n-1)a+1)\delta], \\ 0 & \text{otherwise.} \end{cases} \qquad (3.2) $$

Focusing on the renormalized stray field energy, which for configurations $\mathbf{m}_n(x, y)$ converging to $-\hat{\mathbf{z}}$ sufficiently fast as $x^2 + y^2 \to \infty$ can be rewritten, after integrations by parts, as

$$ E_d(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}^3} \left( |\nabla U|^2 - \sum_n \chi_n \right) d^3r. \qquad (3.3) $$

We define $\tilde{\mathbf{m}}_r = \tilde{\mathbf{m}} - \mathbf{m}^*$ and $U_r = U - U^*$, and observe that

$$ E_d(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}^3} \left( |\nabla U_r|^2 - 2 \sum_{n=1}^{N} \chi_n \partial_z U_r \right) d^3r = E'_d(\tilde{\mathbf{m}}) + E''_d(\tilde{\mathbf{m}}). \qquad (3.4) $$

Examining the second term, we note that

$$ \begin{aligned} \sum_{n=1}^{N} \int_{\mathbb{R}^3} \chi_n \partial_z U_r \, d^3r &= - \int_{\mathbb{R}^3} \nabla U^* \cdot \nabla U_r \, d^3r = \int_{\mathbb{R}^3} U^* \nabla \cdot \tilde{\mathbf{m}}_r \, d^3r \\ &= - \int_{\mathbb{R}^3} \tilde{\mathbf{m}}_r \cdot \nabla U^* \, d^3r = \sum_{n=1}^{N} \int_{\mathbb{R}^3} (m_n^{\parallel} + 1) \chi_n \, d^3r, \end{aligned} \qquad (3.5) $$
after integrations by parts and using (2.4). Therefore, using the fact that $|\mathbf{m}_n| = 1$ and hence $2(m_n^\parallel + 1) = |\mathbf{m}_n^\perp|^2 + |m_n^\parallel + 1|^2$ for all $n$, we obtain

$$E_d''(\tilde{\mathbf{m}}) = -\frac{2}{\delta} \sum_{n=1}^{N} \int_{\mathbb{R}^3} \chi_n \partial_z U_r \, d^3r = -\sum_{n=1}^{N} \int_{\mathbb{R}^2} (|\mathbf{m}_n^\perp|^2 + |m_n^\parallel + 1|^2) \, d^2r. \quad (3.6)$$

We now investigate

$$E_d'(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}^3} |\nabla U_r|^2 d^3r = \frac{1}{\delta} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \frac{\nabla \cdot \tilde{\mathbf{m}}_r(\mathbf{r}) \nabla \cdot \tilde{\mathbf{m}}_r(\mathbf{r}')}{4\pi |\mathbf{r} - \mathbf{r}'|} d^3r \, d^3r', \quad (3.7)$$

as can be seen from (2.5) and an integration by parts. Explicitly, this integral reads

$$E_d'(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}} \int_{\mathbb{R}} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{\nabla \cdot \tilde{\mathbf{m}}_r(\mathbf{r}, z) \nabla \cdot \tilde{\mathbf{m}}_r(\mathbf{r}', z')}{4\pi \sqrt{|\mathbf{r} - \mathbf{r}'|^2 + (z - z')^2}} \, d^2r \, d^2r' \, dz \, dz', \quad (3.8)$$

where from now on the variables of integration $\mathbf{r}, \mathbf{r}' \in \mathbb{R}^2$. Noting that

$$\partial_z \tilde{m}_r^\parallel(\mathbf{r}, z) = \sum_{n=1}^{n} (m_n^\parallel(\mathbf{r}) + 1)[\delta_{a(n-1)\delta}(z) - \delta_{(a(n-1)+1)\delta}(z)], \quad \forall \mathbf{r} \in \mathbb{R}^2, \quad (3.9)$$

where $\delta_\alpha(z)$ is the Dirac delta centered at $\alpha \in \mathbb{R}$, and integrating in $z$ and $z'$, we obtain

$$E_d'(\tilde{\mathbf{m}}) = \frac{1}{\delta} \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \left\{ \begin{aligned} & K_{vv}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \nabla \cdot \mathbf{m}_n^\perp(\mathbf{r}) \nabla \cdot \mathbf{m}_k^\perp(\mathbf{r}') \\ & + K_{vs}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \nabla \cdot \mathbf{m}_n^\perp(\mathbf{r}) (m_k^\parallel(\mathbf{r}') + 1) \\ & + K_{sv}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^\parallel(\mathbf{r}) + 1) \nabla \cdot \mathbf{m}_k^\perp(\mathbf{r}') \\ & + K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^\parallel(\mathbf{r}) + 1)(m_k^\parallel(\mathbf{r}') + 1) \end{aligned} \right\} d^2r \, d^2r', \quad (3.10)$$

where for $u \in \mathbb{R}$ we defined the volume-volume, volume-surface and surface-surface charge interaction kernels as

$$K_{vv}^{u}(r) = \frac{1}{4\pi} \int_{0}^{\delta} \int_{u\delta}^{(u+1)\delta} \frac{dz' dz}{\sqrt{r^2 + (z-z')^2}}, \quad (3.11)$$

$$K_{vs}^{u}(r) = \frac{1}{4\pi} \int_{0}^{\delta} \frac{dz}{\sqrt{r^2 + (z - u\delta)^2}} - \frac{1}{4\pi} \int_{0}^{\delta} \frac{dz}{\sqrt{r^2 + (z - (1 + u)\delta)^2}}, \quad (3.12)$$

$$K_{sv}^{u}(r) = \frac{1}{4\pi} \int_{0}^{\delta} \frac{dz'}{\sqrt{r^2 + (z' + u\delta)^2}} - \frac{1}{4\pi} \int_{0}^{\delta} \frac{dz'}{\sqrt{r^2 + (z' - (1 - u)\delta)^2}}, \quad (3.13)$$

$$K_{ss}^{u}(r) = \frac{1}{2\pi\sqrt{r^2 + u^2\delta^2}} - \frac{1}{4\pi\sqrt{r^2 + (u-1)^2\delta^2}} - \frac{1}{4\pi\sqrt{r^2 + (u+1)^2\delta^2}}, \quad (3.14)$$

We see clearly that $K_{vs}^u(r) = K_{sv}^{-u}(r)$, allowing these terms to be combined in the expression for $E_d'$. For simplicity, we express the energy as the sum of the interaction energies:

$$E_d'(\tilde{\mathbf{m}}) = \sum_{n=1}^{N} \sum_{k=1}^{N} \left\{ E_{vv}^{(nk)}(\tilde{\mathbf{m}}) + E_{vs}^{(nk)}(\tilde{\mathbf{m}}) + E_{ss}^{(nk)}(\tilde{\mathbf{m}}) \right\}, \quad (3.15)$$

where

$$E_{vv}^{(nk)}(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{vv}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \nabla \cdot \mathbf{m}_n^\perp(\mathbf{r}) \nabla \cdot \mathbf{m}_k^\perp(\mathbf{r}') \, d^2r \, d^2r', \quad (3.16)$$

$$E_{vs}^{(nk)}(\tilde{\mathbf{m}}) = -\frac{2}{\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{vs}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \mathbf{m}_n^\perp(\mathbf{r}) \cdot \nabla m_k^\parallel(\mathbf{r}') \, d^2r \, d^2r', \quad (3.17)$$

$$E_{ss}^{(nk)}(\tilde{\mathbf{m}}) = \frac{1}{\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^\parallel(\mathbf{r}) + 1)(m_k^\parallel(\mathbf{r}') + 1) \, d^2r \, d^2r', \quad (3.18)$$
and we integrated by parts in the second line.

We can explicitly evaluate the interaction kernels:

$$ K_{vv}^{u}(r) = \frac{1}{4\pi} \left[ 2\sqrt{r^2 + u^2\delta^2} - \sqrt{r^2 + (1+u)^2\delta^2} - \sqrt{r^2 + (1-u)^2\delta^2} + (1+u)\delta\sinh^{-1}\left(\frac{(1+u)\delta}{r}\right) + (1-u)\delta\sinh^{-1}\left(\frac{(1-u)\delta}{r}\right) - 2u\delta\sinh^{-1}\left(\frac{u\delta}{r}\right) \right], \quad (3.19) $$

$$ K_{vs}^{u}(r) = \frac{1}{4\pi} \left[ 2\sinh^{-1}\left(\frac{u\delta}{r}\right) - \sinh^{-1}\left(\frac{(u+1)\delta}{r}\right) - \sinh^{-1}\left(\frac{(u-1)\delta}{r}\right) \right]. \quad (3.20) $$

The obtained expressions for $K_{vv}^s$, $K_{vs}^u$ and $K_{ss}^u$ may be viewed as generalizations of those obtained in [21] for $u = 0$. Notice that by an explicit calculation we have

$$ \int_{\mathbb{R}^2} K_{ss}^{u}(|\mathbf{r}|) d^2r = \begin{cases} \delta(1-|u|), & |u| \le 1, \\ 0, & |u| > 1. \end{cases} \quad (3.21) $$

By a similar explicit calculation

$$ \int_{\mathbb{R}^2} K_{vs}^{u}(|\mathbf{r}|) d^2r = \frac{\delta^2}{2} \operatorname{sgn}(u), \quad |u| > 1. \quad (3.22) $$

Also notice that $K_{vs}^0(r) = 0$, which can also be seen immediately from the symmetry considerations. We can further simplify the surface-surface interaction term. First we observe that

$$ \begin{aligned} (m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}')) (m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}')) = & (m_n^{\parallel}(\mathbf{r}) + 1) (m_k^{\parallel}(\mathbf{r}) + 1) + (m_n^{\parallel}(\mathbf{r}') + 1) (m_k^{\parallel}(\mathbf{r}') + 1) \\ & - (m_n^{\parallel}(\mathbf{r}) + 1) (m_k^{\parallel}(\mathbf{r}') + 1) - (m_n^{\parallel}(\mathbf{r}') + 1) (m_k^{\parallel}(\mathbf{r}) + 1). \end{aligned} \quad (3.23) $$

Integrating with the kernel $K_{ss}^{a(k-n)}$, we find due to the invariance of $K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|)$ with respect to interchanging $\mathbf{r}$ and $\mathbf{r}'$ that

$$ \begin{aligned} & \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}')) (m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}')) d^2r \, d^2r' \\ & = 2 \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^{\parallel}(\mathbf{r}) + 1) (m_k^{\parallel}(\mathbf{r}) + 1) \, d^2r \, d^2r' \\ & - 2 \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^{\parallel}(\mathbf{r}) + 1) (m_k^{\parallel}(\mathbf{r}') + 1) \, d^2r \, d^2r'. \end{aligned} \quad (3.24) $$

Therefore, the surface-surface self-interaction energy for layer $n$ may be written as

$$ E_{ss}^{(nn)}(\tilde{\mathbf{m}}) = \int_{\mathbb{R}^2} |m_n^{\parallel} + 1|^2 d^2r - \frac{1}{2\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{0}(|\mathbf{r} - \mathbf{r}'|) (m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))^2 d^2r \, d^2r', \quad (3.25) $$

where we noted that by (3.21) we have $\int_{\mathbb{R}^2} K_{ss}^0(|\mathbf{r}|)d^2r = \delta$. At the same time, the surface-surface interaction energy for different layers $n \neq k$ is

$$ E_{ss}^{(nk)}(\tilde{\mathbf{m}}) = -\frac{1}{2\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}')) (m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}')) \, d^2r \, d^2r', \quad (3.26) $$

where again we used (3.21) and noted that $\int_{\mathbb{R}^2} K_{ss}^u(|\mathbf{r}|)d^2r = 0$ for all $u > 1$, recalling that $a > 1$ and $|k - n| \ge 1$. Then, combining all three interactions in (3.15) with (3.6) we get an exact expression for
the energy of $\tilde{\mathbf{m}}$:

$$ E_d(\tilde{\mathbf{m}}) = -\sum_{n=1}^{N} \int_{\mathbb{R}^2} |\mathbf{m}_n^{\perp}|^2 d^2r + \frac{1}{\delta} \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{vv}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \nabla \cdot \mathbf{m}_n^{\perp}(\mathbf{r}) \nabla \cdot \mathbf{m}_k^{\perp}(\mathbf{r}') d^2r \, d^2r' \\ -\frac{2}{\delta} \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{vs}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) \mathbf{m}_n^{\perp}(\mathbf{r}) \cdot \nabla \mathbf{m}_k^{\parallel}(\mathbf{r}') d^2r \, d^2r' \\ -\frac{1}{2\delta} \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} K_{ss}^{a(k-n)}(|\mathbf{r} - \mathbf{r}'|) (\mathbf{m}_n^{\parallel}(\mathbf{r}) - \mathbf{m}_n^{\parallel}(\mathbf{r}')) (\mathbf{m}_k^{\parallel}(\mathbf{r}) - \mathbf{m}_k^{\parallel}(\mathbf{r}')) d^2r \, d^2r'. \quad (3.27) $$

We now simplify this cumbersome expression in such a way that it is valid asymptotically as $\delta \to 0$ with $\tilde{\mathbf{m}}$ and all the other parameters fixed. Since the layer displacement parameter, $u$, is understood to be fixed with respect to $\delta \to 0$, one expands to find the following asymptotic behaviors for $\delta \ll 1$:

$$ K_{vv}^{u}(r) = \frac{\delta^2}{4\pi r} + o(\delta^2), \quad (3.28) $$

$$ K_{vs}^{u}(r) = \frac{u\delta^3}{4\pi r^3} + o(\delta^3), \quad (3.29) $$

$$ K_{ss}^{u}(r) = \frac{\delta^2}{4\pi r^3} + o(\delta^3). \quad (3.30) $$

Therefore, for $\mathbf{m}_n$ sufficiently smooth the expression in (3.25) may be rendered asymptotically for $\delta \ll 1$ as

$$ E_{ss}^{(nn)}(\tilde{\mathbf{m}}) = \int_{\mathbb{R}^2} |m_n^{\parallel} + 1|^2 d^2r - \delta \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{(m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))^2}{8\pi |\mathbf{r} - \mathbf{r}'|^3} d^2r \, d^2r' + o(\delta), \quad (3.31) $$

where we noted that the strong singularity of the kernel is partially cancelled by $(m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))^2 = O(|\mathbf{r} - \mathbf{r}'|^2)$ for $m_n^{\parallel} \in C^1(\mathbb{R}^2)$, making the integral in the right-hand side of (3.31) convergent. Similarly, we have

$$ E_{ss}^{(nk)}(\tilde{\mathbf{m}}) \approx -\delta \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{(m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))(m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}'))}{8\pi |\mathbf{r} - \mathbf{r}'|^3} d^2r \, d^2r' + o(\delta) \quad (3.32) $$

for all $n \neq k$. Meanwhile, using (3.28) the volume-volume interactions can be asymptotically expressed as

$$ E_{vv}^{(nk)}(\tilde{\mathbf{m}}) = \delta \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{\nabla \cdot \mathbf{m}_n^{\perp}(\mathbf{r}) \nabla \cdot \mathbf{m}_k^{\perp}(\mathbf{r}')}{4\pi |\mathbf{r} - \mathbf{r}'|} d^2r \, d^2r' + o(\delta). \quad (3.33) $$

As can be seen from (3.29), we have $\delta^{-2}K_{vs}^{u}(r) \to 0$ as $\delta \to 0$ for all $r > 0$. Care, however, is needed in passing to the limit in the integral, as the strong singularity of the kernel in (3.29) precludes passing the pointwise limit to the value of $E_{vs}^{(nk)}(\tilde{\mathbf{m}})$ at $O(\delta)$. In fact, from (3.22) and the above observation it is clear that for all $|u| > 1$ we have

$$ \frac{2}{\delta^2} K_{vs}^{u}(|\mathbf{r}|) \to \operatorname{sgn}(u)\delta(\mathbf{r}) \quad \text{as } \delta \to 0, \quad (3.34) $$

in the sense of distributions, where $\delta(\mathbf{r})$ is a Dirac delta centered at the origin in $\mathbb{R}^2$. As a consequence, by an argument similar to the one in the case of the surface-surface interactions the volume-surface interaction energies admit the following asymptotic expansion:

$$ E_{vs}^{(nk)}(\tilde{\mathbf{m}}) = -\operatorname{sgn}(k-n)\delta \int_{\mathbb{R}^2} \mathbf{m}_n^{\perp} \cdot \nabla m_k^{\parallel} d^2r + o(\delta). \quad (3.35) $$
Thus, combining all the terms as in (3.27), we obtain that to within $o(\delta)$ accuracy for $\delta \ll 1$ and with all other parameters fixed, the total energy is asymptotically $E(\tilde{\mathbf{m}}) \simeq E_N(\{\mathbf{m}_n\})$, where

$$ 
\begin{aligned}
E_N(\{\mathbf{m}_n\}) = & \sum_{n=1}^{N} \int_{\mathbb{R}^2} (|\nabla \mathbf{m}_n|^2 + (Q-1)|\mathbf{m}_n^{\perp}|^2 - 2h(m_n^{\parallel} + 1) - 2\kappa \mathbf{m}_n^{\perp} \cdot \nabla m_n^{\parallel}) d^2r \\
& - \delta \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \int_{\mathbb{R}^2} (\mathbf{m}_n^{\perp} \cdot \nabla m_k^{\parallel} - \mathbf{m}_k^{\perp} \cdot \nabla m_n^{\parallel}) d^2r \\
& + \delta \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \left( \frac{\nabla \cdot \mathbf{m}_n^{\perp}(\mathbf{r}) \nabla \cdot \mathbf{m}_k^{\perp}(\mathbf{r}')}{4\pi|\mathbf{r}-\mathbf{r}'|} - \frac{(m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))(m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}'))}{8\pi|\mathbf{r}-\mathbf{r}'|^3} \right) d^2r \, d^2r',
\end{aligned} 
\quad (3.36) $$

where we defined the material quality factor and the dimensionless DMI strength

$$ Q = \frac{K}{K_d}, \quad \kappa = \kappa^+ - \kappa^-, \quad (3.37) $$

respectively. We note that, as can be seen from the above derivation, the asymptotic formula in (3.36) is rigorously valid in the limit $\delta \to 0$ at least for every $\mathbf{m}_n : \mathbb{R}^2 \to \mathbb{S}^2$ such that $\mathbf{m}_n + \hat{\mathbf{z}} \in C_c^{\infty}(\mathbb{R}^2; \mathbb{R}^3)$ for every $n = 1, \dots, N$, and by an approximation argument also extends to the natural class of functions $\mathbf{m}_n + \hat{\mathbf{z}} \in H^1(\mathbb{R}^2; \mathbb{R}^3)$.

The obtained energy functional in (3.36) generalizes the one obtained in [3,4] for a single layer to the case of multilayers. Notice that in addition to the non-local dipolar interactions present in the case of a single-layer, a peculiar new local term appears in the second line of (3.36) that corresponds to the interaction of the out-of-plane magnetic moments with the approximately vertical magnetic field formed by the volume charges immediately above and below the layers (see also [35]). This term vanishes, as expected, in the case when all layers have identical magnetizations, in which case the energy is equivalent to that of the magnetization of a single layer of thickness $N\delta$. However, when the in-plane magnetizations in two layers are opposite to one another: $\mathbf{m}_n^{\perp} = -\mathbf{m}_k^{\perp}$, this term works as an effective interfacial DMI term, favoring Néel rotation in these layers.

The reduced energy is applicable in the situation in which the characteristic scale of variation of $\mathbf{m}_n(x, y)$ exceeds the thickness $aN\delta$ of the entire stack. Note that in this regime the energy is independent of $a$ and, hence, does not see the presence of the non-magnetic spacers. This property is known to be violated once the magnetization configurations acquire the lateral size comparable to the stack thickness [32], in which case the full kernels $K_{vv}^u$, $K_{vs}^u$ and $K_{ss}^u$ need to be utilized, resulting in a much more cumbersome model that only permits a numerical treatment [10, 33, 35].

# 4 Energy of N stacked skyrmions

We have demonstrated that for ultrathin multilayers the micromagnetic energy functional of a magnetization configuration $\mathbf{m}: \bar{\Omega} \to \mathbb{S}^2$ obeys $E(\mathbf{m}) \simeq E_N(\{\mathbf{m}_n\})$, where $E_N(\{\mathbf{m}_n\})$ is given by (3.36). This energy may be further simplified by taking advantage of the scaling properties of different terms in the energy. For $Q > 1$, introducing a rescaling and the new parameters:

$$ \mathbf{r} \to \frac{\mathbf{r}}{\sqrt{Q-1}}, \quad \bar{\delta} = \frac{\delta}{\sqrt{Q-1}}, \quad \bar{h} = \frac{h}{Q-1}, \quad \bar{\kappa} = \frac{\kappa}{\sqrt{Q-1}}, \quad (4.1) $$

we see that $E_N(\{\mathbf{m}_n(\cdot/\sqrt{Q-1})\}) = \bar{E}_N(\{\mathbf{m}_n\})$, where

$$ 
\begin{aligned}
\bar{E}_N(\{\mathbf{m}_n\}) = & \sum_{n=1}^{N} \int_{\mathbb{R}^2} (|\nabla \mathbf{m}_n|^2 + |\mathbf{m}_n^{\perp}|^2 - 2\bar{h}(m_n^{\parallel} + 1) - 2\bar{\kappa} \mathbf{m}_n^{\perp} \cdot \nabla m_n^{\parallel}) d^2r \\
& - \bar{\delta} \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \int_{\mathbb{R}^2} (\mathbf{m}_n^{\perp} \cdot \nabla m_k^{\parallel} - \mathbf{m}_k^{\perp} \cdot \nabla m_n^{\parallel}) d^2r \\
& + \bar{\delta} \sum_{n=1}^{N} \sum_{k=1}^{N} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \left( \frac{\nabla \cdot \mathbf{m}_n^{\perp}(\mathbf{r}) \nabla \cdot \mathbf{m}_k^{\perp}(\mathbf{r}')}{4\pi|\mathbf{r}-\mathbf{r}'|} - \frac{(m_n^{\parallel}(\mathbf{r}) - m_n^{\parallel}(\mathbf{r}'))(m_k^{\parallel}(\mathbf{r}) - m_k^{\parallel}(\mathbf{r}'))}{8\pi|\mathbf{r}-\mathbf{r}'|^3} \right) d^2r \, d^2r'.
\end{aligned} 
\quad (4.2) $$
Hence in the following we focus our attention on the study of the energy $\bar{E}_N$ for the magnetization configurations consisting of a single skyrmion in each ferromagnetic layer. The latter may be specified by prescribing the topological degree +1 to the magnetization in each layer:

$$ \mathcal{N}(\mathbf{m}_n) = \frac{1}{4\pi} \int_{\mathbb{R}^2} \mathbf{m}_n \cdot (\partial_x \mathbf{m}_n \times \partial_y \mathbf{m}_n) \, d^2r = 1 \quad \forall n = 1, \dots, N, \qquad (4.3) $$

which together with some additional technical assumptions should ensure existence of local minimizers of the energy [4, 5].

For $N = 1$, it is known that for $\bar{h} = 0$ the energy $\bar{E}_N$ admits minimizers in an appropriate function class for all $0 < \bar{\delta} < \bar{\delta}_0$, with $\bar{\delta}_0 > 0$ universal [4]. Moreover, the energy-minimizing profiles admit a complete asymptotic characterization in the conformal limit $\bar{\kappa}, \bar{\delta} \to 0$ [5]. After a suitable translation, dilation and rotation, these profiles approach the canonical Belavin-Polyakov (BP) profile

$$ \mathbf{m}_{\infty}(\mathbf{r}) = \left( -\frac{2\mathbf{r}}{1 + |\mathbf{r}|^2}, \frac{1 - |\mathbf{r}|^2}{1 + |\mathbf{r}|^2} \right). \qquad (4.4) $$

This is due to the fact that for $\bar{\delta} \ll 1$ the energy of a skyrmion in a single layer gets close to the value of the Dirichlet energy of harmonic maps with degree +1, whose minimizers are well-known [1, 5, 16]. Furthermore, the rigidity estimates for almost harmonic maps of degree +1 and a fine analysis of the tail of the skyrmion profiles yield closeness of the profiles to those of the form $\mathbf{m}(\mathbf{r}) = R\mathbf{m}_{\infty}((\mathbf{r}-\mathbf{r}_0)/\rho)$, where $R \in SO(3)$ is a rotation around the z-axis [5]. As the anisotropy energy evaluated on $\mathbf{m}_{\infty}$ diverges, in order to determine the skyrmion radius $\rho$ and the rotation angle $\theta$ as functions of $\bar{\kappa}$ and $\bar{\delta}$ one needs to consider a truncated BP profile, whose energy yields the asymptotic dependence of the skyrmion characteristics on $\bar{\kappa}, \bar{\delta} \to 0$ [5].

To that end, we define

$$ f(r) = \frac{2r}{1 + r^2}, \quad r \ge 0, \qquad (4.5) $$

and its truncated version

$$ f_L(r) = \begin{cases} f(r), & \text{if } r \le \sqrt{L}, \\ \frac{f(\sqrt{L})}{K_1(1/\sqrt{L})} K_1(r/L), & \text{if } r > \sqrt{L}, \end{cases} \qquad (4.6) $$

where $K_1(x)$ is the modified Bessel function of the second kind. This choice of the truncation is motivated by the asymptotic decay of the skyrmion solution at infinity to the leading order in $\bar{\kappa}, \bar{\delta} \ll 1$ determined by the exchange and anisotropy terms [28]. For $L > 1$ we then define

$$ \mathbf{m}_{\rho, \theta, L, \mathbf{r}_0}(\mathbf{r}) = \left( -f_L\left(\frac{|\mathbf{r} - \mathbf{r}_0|}{\rho}\right) \frac{R_{\theta}(\mathbf{r} - \mathbf{r}_0)}{|\mathbf{r} - \mathbf{r}_0|}, \operatorname{sgn}(\rho - |\mathbf{r} - \mathbf{r}_0|)\sqrt{1 - f_L^2\left(\frac{|\mathbf{r} - \mathbf{r}_0|}{\rho}\right)} \right), \qquad (4.7) $$

where $R_{\theta} \in SO(2)$ is a counter-clockwise rotation by angle $\theta \in [-\pi, \pi)$. Using the above expression as an ansatz, the task of determining the values of $\rho$ and $\theta$ for a skyrmion amounts to minimizing the energy of $\mathbf{m}_{\rho, \theta, L, \mathbf{r}_0}$ in $\rho, \theta$ and $L$ to the leading order in $\bar{\kappa}, \bar{\delta} \ll 1$.

It is clear that the above considerations should remain valid for several layers, each containing a single skyrmion. Therefore, we now consider an asymptotic expansion of the energy $\bar{E}_N$ for the profiles of the form (see Fig. 2 for an illustration):

$$ \mathbf{m}_n = \mathbf{m}_{\rho_n, \theta_n, L_n, \mathbf{r}_n}, \quad n = 1, \dots, N, \qquad (4.8) $$

and study the energy landscape of $\bar{E}_N(\{\mathbf{m}_{\rho_n, \theta_n, L_n, \mathbf{r}_n}\})$ in terms of its dependence on $\{\rho_n, \theta_n, L_n, \mathbf{r}_n\}$ for $\bar{\delta}, |\bar{\kappa}| \ll 1$. Notice that despite having reduced our problem to a finite-dimensional one, we still need to study a strongly nonlinear, fully coupled system whose analysis is a significant challenge. In particular, it is not a priori clear whether the minimum of $\bar{E}_N(\{\mathbf{m}_{\rho_n, \theta_n, L_n, \mathbf{r}_n}\})$ is attained, as it may be energetically favorable for a skyrmion in one of the layers to collapse, which would correspond to $\rho_n \to 0$, violating (4.3) in the limit.
<figure><img src="image_3.png" /><figcaption>Figure 2: An example of a skyrmion configuration from (4.8) with three distinct radii and centers in a ferromagnetic trilayer. The skyrmion in the bottom layer is of Néel type ($\theta_1 = 0$), the skyrmion in the top layer is of Bloch type ($\theta_3 = \pi/2$), and the skyrmion in the middle layer is of mixed type ($\theta_2 = \pi/4$).</figcaption></figure>

For simplicity of notation, we can write

$$ \bar{E}_N(\{\mathbf{m}_{\rho_n, \theta_n, L_n, \mathbf{r}_n}\}) = \sum_{n=1}^{N} (\bar{E}_{\text{ex}}(\mathbf{m}_n) + \bar{E}_{\text{an}}(\mathbf{m}_n) + \bar{E}_{\text{Z}}(\mathbf{m}_n) + \bar{E}_{\text{DMI}}(\mathbf{m}_n)) \\ + \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \bar{E}_{\text{vs}}(\mathbf{m}_n, \mathbf{m}_k) + \sum_{n=1}^{N} \sum_{k=1}^{N} (\bar{E}_{\text{vv}}(\mathbf{m}_n, \mathbf{m}_k) + \bar{E}_{\text{ss}}(\mathbf{m}_n, \mathbf{m}_k)), \quad (4.9) $$

where each term in the sum corresponds to the respective one in (4.2). As was shown in [5], for all $L_n \ge L_0$, with some $L_0 > 1$ universal, one can carry out an expansion of each term in (4.9) that becomes asymptotically exact as $L_0 \to \infty$. Utilizing the results from [5, Lemma A.6] and [3], we obtain that for $L_0 \gg 1$ we have

$$ \bar{E}_{\text{ex}}(\mathbf{m}_{n, \rho_n, \theta_n, L_n}) - 8\pi \approx \frac{4\pi}{L_n^2}, \quad (4.10) $$

$$ \bar{E}_{\text{an}}(\mathbf{m}_{n, \rho_n, \theta_n, L_n}) \approx 4\pi\rho_n^2 \ln\left(\frac{4L_n^2}{e^{2(1+\gamma)}}\right), \quad (4.11) $$

$$ \bar{E}_{\text{Z}}(\mathbf{m}_{n, \rho_n, \theta_n, L_n}) \approx -4\pi\bar{h}\rho_n^2 \ln\left(\frac{4L_n^2}{e^{1+2\gamma}}\right), \quad (4.12) $$

$$ \bar{E}_{\text{DMI}}(\mathbf{m}_{n, \rho_n, \theta_n, L_n}) \approx -8\pi\bar{\kappa}\rho_n \cos \theta_n. \quad (4.13) $$

We next introduce the Fourier transform of $\mathbf{m}_{r,n} = \mathbf{m}_n + \hat{\mathbf{z}} \in H^1(\mathbb{R}^2; \mathbb{R}^3)$:

$$ \hat{\mathbf{m}}_{r,n}(\mathbf{q}) = \int_{\mathbb{R}^2} e^{-i\mathbf{q} \cdot \mathbf{r}} \mathbf{m}_{r,n}(\mathbf{r}) \, d^2\mathbf{r} \quad (4.14) $$
and write the stray field-mediated terms as [37]

$$ \bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) = \frac{\bar{\delta}}{2} \int_{\mathbb{R}^2} \frac{(\mathbf{q} \cdot \overline{\hat{\mathbf{m}}}_{r,n}^{\perp})(\mathbf{q} \cdot \overline{\hat{\mathbf{m}}}_{r,k}^{\perp})}{|\mathbf{q}|} \frac{d^2 q}{(2\pi)^2}, \quad (4.15) $$

$$ \bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) = -\frac{\bar{\delta}}{2} \int_{\mathbb{R}^2} |\mathbf{q}| (\overline{\hat{\mathbf{m}}}_{r,n}^{\parallel}) \overline{\hat{\mathbf{m}}}_{r,k}^{\parallel} \frac{d^2 q}{(2\pi)^2}, \quad (4.16) $$

$$ \bar{E}_{vs}(\mathbf{m}_n, \mathbf{m}_k) = -i\bar{\delta} \int_{\mathbb{R}^2} \mathbf{q} \cdot \left( \overline{\hat{\mathbf{m}}}_{r,n}^{\perp} \overline{\hat{\mathbf{m}}}_{r,k}^{\parallel} - \overline{\hat{\mathbf{m}}}_{r,k}^{\perp} \overline{\hat{\mathbf{m}}}_{r,n}^{\parallel} \right) \frac{d^2 q}{(2\pi)^2}, \quad (4.17) $$

where the overline denotes complex conjugate.

We now assume that the magnetization $\mathbf{m}_n$ takes the form (4.8). Using the fact that

$$ R_{\theta_n} \frac{\mathbf{r} - \mathbf{r}_n}{|\mathbf{r} - \mathbf{r}_n|} = \frac{\mathbf{r} - \mathbf{r}_n}{|\mathbf{r} - \mathbf{r}_n|} \cos \theta_n + \frac{(\mathbf{r} - \mathbf{r}_n)^{\perp}}{|\mathbf{r} - \mathbf{r}_n|} \sin \theta_n, \quad (4.18) $$

where $(\mathbf{r}-\mathbf{r}_n)^{\perp}$ denotes the 90° counter-clockwise rotation of $\mathbf{r}-\mathbf{r}_n$, we can split the in-plane component of $\mathbf{m}_n$ into two terms and note that the second term has zero divergence and will not factor into the volume charge energy. Using the results from [5, Lemmas A.5 and A.6], we obtain asymptotically for $L_0 \gg 1$:

$$ \bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) \approx 2\bar{\delta}\rho_n^2 \rho_k^2 \cos \theta_n \cos \theta_k \int_{\mathbb{R}^2} e^{i\mathbf{q} \cdot (\mathbf{r}_n - \mathbf{r}_k)} |\mathbf{q}| K_1(\rho_n |\mathbf{q}|) K_1(\rho_k |\mathbf{q}|) d^2 q, \quad (4.19) $$

$$ \bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) \approx -2\bar{\delta}\rho_n^2\rho_k^2 \int_{\mathbb{R}^2} e^{i\mathbf{q} \cdot (\mathbf{r}_n - \mathbf{r}_k)} |\mathbf{q}| K_0(\rho_n |\mathbf{q}|) K_0(\rho_k |\mathbf{q}|) d^2 q, \quad (4.20) $$

$$ \begin{aligned} \bar{E}_{vs}(\mathbf{m}_n, \mathbf{m}_k) \approx & -4\bar{\delta}\rho_n^2\rho_k^2 \cos \theta_n \int_{\mathbb{R}^2} e^{i\mathbf{q} \cdot (\mathbf{r}_n - \mathbf{r}_k)} |\mathbf{q}| K_1(\rho_n |\mathbf{q}|) K_0(\rho_k |\mathbf{q}|) d^2 q \\ & + 4\bar{\delta}\rho_n^2\rho_k^2 \cos \theta_k \int_{\mathbb{R}^2} e^{i\mathbf{q} \cdot (\mathbf{r}_k - \mathbf{r}_n)} |\mathbf{q}| K_0(\rho_n |\mathbf{q}|) K_1(\rho_k |\mathbf{q}|) d^2 q. \end{aligned} \quad (4.21) $$

Integrating in polar coordinates with $q = |\mathbf{q}|$ and $\varphi$ being the polar angle between the vectors $\mathbf{r}_n - \mathbf{r}_k$ and $\mathbf{q}$, we can further reduce the above integrals with the help of the well-known formula

$$ \frac{1}{2\pi} \int_{0}^{2\pi} e^{iq|\mathbf{r}_n - \mathbf{r}_k| \cos \varphi} d\varphi = J_0(q|\mathbf{r}_n - \mathbf{r}_k|), \quad (4.22) $$

where $J_0(x)$ is the Bessel function of the first kind. Introducing the new variables

$$ \alpha = \sqrt{\frac{\rho_k}{\rho_n}}, \quad \beta = \sqrt{\rho_n \rho_k}, \quad \lambda = \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}}, \quad (4.23) $$

and the rescaled variable of integration $\xi = q\sqrt{\rho_n \rho_k}$, we can write the energies above as

$$ \bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) \approx 4\pi\bar{\delta}\beta \cos \theta_n \cos \theta_k \int_0^\infty J_0(\lambda\xi) K_1(\alpha\xi) K_1(\xi/\alpha) \xi^2 d\xi, \quad (4.24) $$

$$ \bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) \approx -4\pi\bar{\delta}\beta \int_0^\infty J_0(\lambda\xi) K_0(\alpha\xi) K_0(\xi/\alpha) \xi^2 d\xi, \quad (4.25) $$

$$ \begin{aligned} \bar{E}_{vs}(\mathbf{m}_n, \mathbf{m}_k) \approx & -8\pi\bar{\delta}\beta \cos \theta_n \int_0^\infty J_0(\lambda\xi) K_0(\alpha\xi) K_1(\xi/\alpha) \xi^2 d\xi \\ & + 8\pi\bar{\delta}\beta \cos \theta_k \int_0^\infty J_0(\lambda\xi) K_1(\alpha\xi) K_0(\xi/\alpha) \xi^2 d\xi. \end{aligned} \quad (4.26) $$

Now define the functions

$$ F_{vv}(\alpha, \lambda) = \frac{32}{3\pi^2} \int_0^\infty \xi^2 J_0(\lambda\xi) K_1(\alpha\xi) K_1(\xi/\alpha) d\xi, \quad (4.27) $$

$$ F_{ss}(\alpha, \lambda) = \frac{32}{\pi^2} \int_0^\infty \xi^2 J_0(\lambda\xi) K_0(\alpha\xi) K_0(\xi/\alpha) d\xi, \quad (4.28) $$

$$ F_{vs}(\alpha, \lambda) = 2 \int_0^\infty \xi^2 J_0(\lambda\xi) K_0(\alpha\xi) K_1(\xi/\alpha) d\xi, \quad (4.29) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: The plots of $F_{vv}(\alpha, 0)$ (a), $F_{ss}(\alpha, 0)$ (b), and $F_{vs}(\alpha, 0)$ (c).</figcaption></figure>

normalized so that $F_{vv}(1, 0) = F_{ss}(1, 0) = F_{vs}(1, 0) = 1$. Notice that all these three functions are uniformly bounded. Indeed, since $K_{0,1}(t) > 0$ for all $t > 0$ and $|J_0(t)| \le 1$, we have $|F_{vv}(\alpha, \lambda)| \le F_{vv}(\alpha, 0)$, $|F_{ss}(\alpha, \lambda)| \le F_{ss}(\alpha, 0)$ and $|F_{vs}(\alpha, \lambda)| \le F_{vs}(\alpha, 0)$, with the equality achieved only for $\lambda = 0$ (see Lemmas 1 and 2 in the Appendix). At the same time, for $0 < \alpha < 1$ there holds

$$F_{vv}(\alpha, 0) = \frac{16\alpha ((\alpha^4 + 1) E(1 - \alpha^4) - 2\alpha^4 K(1 - \alpha^4))}{3\pi (\alpha^4 - 1)^2}, \quad (4.30)$$

$$F_{ss}(\alpha, 0) = \frac{16\alpha^3 ((\alpha^4 + 1) K(1 - \alpha^4) - 2E(1 - \alpha^4))}{\pi (\alpha^4 - 1)^2}, \quad (4.31)$$

where $K(m)$ and $E(m)$ are the complete elliptic integrals of the first and second kind, respectively.<fnref n="1" /> For $\alpha > 1$ one can infer the values of the above functions via identities $F_{vv}(\alpha, 0) = F_{vv}(\alpha^{-1}, 0)$ and $F_{ss}(\alpha, 0) = F_{ss}(\alpha^{-1}, 0)$. We also can explicitly compute for any $\alpha > 0$

$$F_{vs}(\alpha, 0) = \frac{2\alpha^3 (\alpha^4 - 4 \ln \alpha - 1)}{(\alpha^4 - 1)^2}, \quad (4.32)$$

extending the last expression by continuity to $F_{vs}(\alpha, 0) = 1$ at $\alpha = 1$. The graphs of the above functions are presented in Fig. 3.

With these definitions, we have the following representation for the energies:

$$\bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) \approx \frac{3\pi^3}{8} \bar{\delta}\beta \cos \theta_n \cos \theta_k F_{vv}(\alpha, \lambda), \quad (4.33)$$

$$\bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) \approx -\frac{\pi^3}{8} \bar{\delta}\beta F_{ss}(\alpha, \lambda), \quad (4.34)$$

$$\bar{E}_{vs}(\mathbf{m}_n, \mathbf{m}_k) \approx -4\pi\bar{\delta}\beta \cos \theta_n F_{vs}(\alpha, \lambda) + 4\pi\bar{\delta}\beta \cos \theta_k F_{vs}(\alpha^{-1}, \lambda). \quad (4.35)$$

<footnote n="1">¹We use the convention $K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$ and $E(m) = \int_0^{\pi/2} \sqrt{1-m\sin^2\theta} d\theta$.</footnote>
Coming back to the original set of parameters {$\rho_n, L_n, \theta_n, \mathbf{r}_n$}, it is clear that

$$ \bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) \approx \frac{3\pi^3}{8} \bar{\delta} \sqrt{\rho_n \rho_k} \cos \theta_n \cos \theta_k F_{vv} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right), \quad (4.36) $$

$$ \bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) \approx -\frac{\pi^3}{8} \bar{\delta} \sqrt{\rho_n \rho_k} F_{ss} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right), \quad (4.37) $$

$$ \begin{aligned} \bar{E}_{vs}(\mathbf{m}_n, \mathbf{m}_k) \approx & -4\pi\bar{\delta}\sqrt{\rho_n\rho_k}\cos\theta_n F_{vs} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right) \\ & + 4\pi\bar{\delta}\sqrt{\rho_n\rho_k}\cos\theta_k F_{vs} \left( \sqrt{\frac{\rho_n}{\rho_k}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right), \end{aligned} \quad (4.38) $$

and from the definitions of $F_{vv}$ and $F_{ss}$ we observe that

$$ \bar{E}_{vv}(\mathbf{m}_n, \mathbf{m}_k) \approx \frac{3\pi^3}{8} \bar{\delta} \rho_n \cos^2 \theta_n, \quad (4.39) $$

$$ \bar{E}_{ss}(\mathbf{m}_n, \mathbf{m}_k) \approx -\frac{\pi^3}{8} \bar{\delta} \rho_n. \quad (4.40) $$

Using the above expressions, we may asymptotically express the full energy, to leading order in $L_0 \gg 1$ as $\bar{E}_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\}) - 8\pi N \approx F_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\})$, where

$$ \begin{aligned} F_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\}) = & \sum_{n=1}^{N} \left[ \frac{4\pi}{L_n^2} + 4\pi\rho_n^2 \ln\left(\frac{4L_n^2}{e^{2(1+\gamma)}}\right) - 4\pi\bar{h}\rho_n^2 \ln\left(\frac{4L_n^2}{e^{1+2\gamma}}\right) - 8\pi\bar{\kappa}\rho_n \cos\theta_n + \bar{\delta}\frac{\pi^3}{8}\rho_n(3\cos^2\theta_n - 1) \right] \\ & + \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \bar{\delta}\frac{\pi^3}{4}\sqrt{\rho_n\rho_k} \left[ 3\cos\theta_n \cos\theta_k F_{vv} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right) - F_{ss} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right) \right] \\ & - \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} 4\pi\bar{\delta}\sqrt{\rho_n\rho_k} \left[ \cos\theta_n F_{vs} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right) - \cos\theta_k F_{vs} \left( \sqrt{\frac{\rho_n}{\rho_k}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n\rho_k}} \right) \right]. \end{aligned} \quad (4.41) $$

This is the reduced energy function that determines the energy landscape of an N-skyrmion stack in terms of the skyrmion parameters.

# 5 The landscape of the energy function $F_N$

We now investigate the energy landscape governed by the function $F_N$ in the regime of its applicability to stacked magnetic skyrmions in ferromagnetic multilayers. To simplify the discussion, from now on we assume that there is no applied external magnetic field, $\bar{h} = 0$, and omit the Zeeman contribution to the energy from now on. The analysis below can be easily extendable to the case $\bar{h} < 0$, when the field is applied opposite to the magnetization direction in the skyrmion core. At the same time, for fields applied along the magnetization in the core, $\bar{h} > 0$, the energy landscape becomes more complex due to skyrmion bursting [3], a new phenomenon whose study goes beyond the present paper.

As we are interested in the minimization of $F_N$, we introduce several auxiliary functions which are obtained by partial minimization. With a slight abuse of notation, we still utilize the symbol $F_N$ to denote those functions, which now depend on fewer variables. For example, we introduce

$$ F_N(\{\rho_n, \theta_n, \mathbf{r}_n\}) = \min_{L_n>0} F_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\}), \quad (5.1) $$

which is obtained by minimizing the energy function in all $L_n$. An explicit calculation shows that $F_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\})$ is minimized for $L_n = \rho_n^{-1}$, which is also the unique critical point of $F_N(\{\rho_n, \theta_n, L_n, \mathbf{r}_n\})$
in $L_n$, resulting in

$$
\begin{align}
F_N(\{\rho_n, \theta_n, \mathbf{r}_n\})
&= \sum_{n=1}^{N} \left[ -4\pi \rho_n^2 \ln \left( \frac{e^{1+2\gamma}}{4} \rho_n^2 \right) - 8\pi \bar{\kappa} \rho_n \cos \theta_n + \bar{\delta} \frac{\pi^3}{8} \rho_n (3 \cos^2 \theta_n - 1) \right] \nonumber \\
&\quad + \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \bar{\delta} \frac{\pi^3}{4} \sqrt{\rho_n \rho_k} \left[ 3 \cos \theta_n \cos \theta_k F_{vv} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) - F_{ss} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) \right] \tag{5.2} \\
&\quad - \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} 4\pi \bar{\delta} \sqrt{\rho_n \rho_k} \left[ \cos \theta_n F_{vs} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) - \cos \theta_k F_{vs} \left( \sqrt{\frac{\rho_n}{\rho_k}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) \right]. \nonumber
\end{align}
$$

Note, however, that this value of $L_n$ is admissible only when $L_n > L_0$, where $L_0 \gg 1$ is the parameter that measures the applicability of the energy function $F_N$ to the system of $N$ stacked skyrmions (see section 4). This means that for consistency we need to restrict the admissible values of $\rho_n$ to $\rho_n < L_0^{-1}$ for a fixed $L_0 > 0$ sufficiently large. We thus define, for $\{\mathbf{r}_n\}_{n=1}^N \subset \mathbb{R}^2$, an admissible class

$$ \mathcal{A}_N(\{\mathbf{r}_n\}) = \{\rho_n \in (0, L_0^{-1}), \theta_n \in [-\pi, \pi), \mathbf{r}_n\}_{n=1}^N. \quad (5.3) $$

Without loss of generality, from now on we may assume that $L_0 > \bar{L}_0$, where

$$ \bar{L}_0 = \frac{1}{2} e^{2+\gamma}, \quad (5.4) $$

which ensures that the first term in the first line of the right-hand side of (5.2) is strictly convex in $\rho_n$.

We next investigate the energy $F_N(\{\mathbf{r}_n\})$ obtained by fixing the positions of the skyrmions and minimizing with respect to angles and radii. We note that as the admissible set of the radii $0 < \rho_n < L_0^{-1}$ is not closed, the existence of solutions for this minimization problem is not a priori clear, since some skyrmions may prefer to collapse, $\rho_n \to 0$, or burst, $\rho_n \to L_0^{-1}$, in the course of minimization. As a first step, we prove that the minimization of $F_N(\{\rho_n, \theta_n, \mathbf{r}_n\})$ in $\theta_n$ and $\rho_n$ is well defined, i.e., that minimizers of this problem exist. The obtained minimal energy can be used to understand the interaction between magnetic skyrmions at different locations. It also allows to estimate from below the energy of the saddle points that separate the basins of attraction of different skyrmion spatial arrangements. We will subsequently illustrate the explicit solution of this problem in the simplest case of stray field-coupled ferromagnetic bilayers with no DMI.

**Theorem 1.** Let $N \in \mathbb{N}$ and $L_0 > \bar{L}_0$. There exists $\bar{\delta}_0 > 0$ such that for any fixed $\{\mathbf{r}_n\} \subset \mathbb{R}^2$ and all $\bar{\delta}, |\bar{\kappa}| < \bar{\delta}_0$ there exists a minimizer of the problem

$$ F_N(\{\mathbf{r}_n\}) = \min_{\mathcal{A}_N(\{\mathbf{r}_n\})} F_N(\{\rho_n, \theta_n, \mathbf{r}_n\}). \quad (5.5) $$

*Proof. Step 1.* We first minimize in angles $\theta_n$, keeping $\{\rho_n, \mathbf{r}_n\}$ fixed, with $1 \le n \le N$. It is clear that if we fix $\mathbf{r}_n \in \mathbb{R}^2$ and $\rho_n \in (0, L_0^{-1})$, the function $F_N(\{\rho_n, \theta_n, \mathbf{r}_n\})$ is continuous and periodic in $\theta_n$ with period $2\pi$ and, therefore, there exists a global minimizer $\{\theta_n^*\} \subset (-\pi, \pi)^N$. We define the resulting minimal energy as $F_N(\{\rho_n, \mathbf{r}_n\}) = \min_{\theta_n} F_N(\{\rho_n, \theta_n, \mathbf{r}_n\})$, where

$$
\begin{align}
F_N(\{\rho_n, \mathbf{r}_n\}) ={}& \sum_{n=1}^{N} \left[ -4\pi \rho_n^2 \ln \left( \frac{e^{1+2\gamma}}{4} \rho_n^2 \right) - 8\pi \bar{\kappa} \rho_n \cos \theta_n^* + \bar{\delta} \frac{\pi^3}{8} \rho_n (3 \cos^2 \theta_n^* - 1) \right] \nonumber \\
& + \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} \bar{\delta} \frac{\pi^3}{4} \sqrt{\rho_n \rho_k} \left[ 3 \cos \theta_n^* \cos \theta_k^* F_{vv} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) - F_{ss} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) \right] \nonumber \\
& - \sum_{n=1}^{N-1} \sum_{k=n+1}^{N} 4\pi \bar{\delta} \sqrt{\rho_n \rho_k} \left[ \cos \theta_n^* F_{vs} \left( \sqrt{\frac{\rho_k}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) - \cos \theta_k^* F_{vs} \left( \sqrt{\frac{\rho_n}{\rho_k}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n \rho_k}} \right) \right]. \tag{5.6}
\end{align}
$$
**Step 2.** Now we show that $F_N(\{\rho_n, \mathbf{r}_n\})$ is bounded from below. From Lemmas 1 and 2 in the Appendix, we know that $F_{ss}(\alpha, \lambda)$, $F_{vv}(\alpha, \lambda)$, and $F_{vs}(\alpha, \lambda)$ are bounded functions. Then since $2\sqrt{\rho_n\rho_k} \le \rho_n + \rho_k$, it follows that for some $C > 0$ universal we have

$$F_N(\{\rho_n, \mathbf{r}_n\}) \ge \sum_{n=1}^{N} \left[ -4\pi \rho_n^2 \ln \left( \frac{e^{1+2\gamma}}{4} \rho_n^2 \right) - C(|\bar{\kappa}| + |\bar{\delta}|) \rho_n \right] \ge -C_1 \quad (5.7)$$

for some $C_1 > 0$ and all $0 < \rho_n < L_0^{-1}$.

Therefore, for $1 \le n \le N$ there are minimizing sequences $(\rho_{n,l})_l$ such that as $l \to \infty$

$$F_N(\{\rho_{n,l}, \mathbf{r}_n\}) \to \inf_{\rho_n} F_N(\{\rho_n, \mathbf{r}_n\}) > -\infty. \quad (5.8)$$

Since $\rho_{n,l} \in (0, L_0^{-1})$, up to extraction of a subsequence (not relabeled) they converge:

$$\lim_{l \to \infty} \rho_{n,l} = \rho_n^* \in [0, L_0^{-1}]. \quad (5.9)$$

**Step 3.** We now show that in a minimizing sequence $\rho_{n,l}$ do not converge to 0. Assume this is not true and, without loss of generality, there is a minimizing sequence with $\rho_{N,l} \to 0$, while for $1 \le n \le N-1$ we have $\rho_{n,l} \to \rho_n^* \ge 0$ as $l \to \infty$. Let us assume $\rho_n^* > 0$ for all $1 \le n \le N-1$, as the other case is simpler and will follow in a similar fashion. We also assume that $\mathbf{r}_N \ne \mathbf{r}_n$ for all $1 \le n \le N-1$, as the other case is also simpler. We now observe that the infimum of the energy is achieved on $\{\rho_n^*\}_{n=1}^{N-1} \cup \{0\}$:

$$\begin{aligned} \inf_{\rho_n} F_N(\{\rho_n, \mathbf{r}_n\}) &= \sum_{n=1}^{N-1} \left[ -4\pi(\rho_n^*)^2 \ln \left( \frac{e^{1+2\gamma}}{4} (\rho_n^*)^2 \right) - 8\pi\bar{\kappa}\rho_n^* \cos\theta_n^* + \bar{\delta}\frac{\pi^3}{8}\rho_n^*(3\cos^2\theta_n^* - 1) \right] \\ &\quad + \sum_{n=1}^{N-2} \sum_{k=n+1}^{N-1} \bar{\delta}\frac{\pi^3}{4}\sqrt{\rho_n^*\rho_k^*} \left[ 3\cos\theta_n^* \cos\theta_k^* F_{vv} \left( \sqrt{\frac{\rho_k^*}{\rho_n^*}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n^*\rho_k^*}} \right) - F_{ss} \left( \sqrt{\frac{\rho_k^*}{\rho_n^*}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n^*\rho_k^*}} \right) \right] \\ &\quad - \sum_{n=1}^{N-2} \sum_{k=n+1}^{N-1} 4\pi\bar{\delta}\sqrt{\rho_n^*\rho_k^*} \left[ \cos\theta_n^* F_{vs} \left( \sqrt{\frac{\rho_k^*}{\rho_n^*}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n^*\rho_k^*}} \right) - \cos\theta_k^* F_{vs} \left( \sqrt{\frac{\rho_n^*}{\rho_k^*}}, \frac{|\mathbf{r}_n - \mathbf{r}_k|}{\sqrt{\rho_n^*\rho_k^*}} \right) \right] \\ &= F_{N-1}(\{\rho_n^*, \mathbf{r}_n\}). \end{aligned} \quad (5.10)$$

Therefore, we only need to show that there exists $\rho_N \ne 0$ such that with $\rho_n = \rho_n^*$ for $1 \le n \le N-1$ we have

$$F_N(\{\rho_n, \mathbf{r}_n\}) < F_{N-1}(\{\rho_n^*, \mathbf{r}_n\}). \quad (5.11)$$

Since $F_N(\{\rho_n, \mathbf{r}_n\})$ is obtained by minimizing in $\theta_n$, from (5.10) we know that by taking $\theta_n = \theta_n^*$ for $1 \le n \le N-1$ and $\theta_N = \frac{\pi}{2}$ we have

$$F_N(\{\rho_n, \mathbf{r}_n\}) \le F_N(\{\rho_n, \theta_n, \mathbf{r}_n\}) = F_{N-1}(\{\rho_n^*, \mathbf{r}_n\}) - G(\{\rho_n, \mathbf{r}_n\}), \quad (5.12)$$

where

$$G(\{\rho_n, \mathbf{r}_n\}) = 4\pi\rho_N^2 \ln \left( \frac{e^{1+2\gamma}}{4} \rho_N^2 \right) + \bar{\delta}\frac{\pi^3}{8}\rho_N + \bar{\delta}\frac{\pi^3}{4}\sum_{n=1}^{N-1} \sqrt{\rho_n\rho_N} F_{ss} \left( \sqrt{\frac{\rho_N}{\rho_n}}, \frac{|\mathbf{r}_N - \mathbf{r}_n|}{\sqrt{\rho_N\rho_n}} \right) \quad (5.13)$$

$$+ 4\pi\bar{\delta} \sum_{n=1}^{N-1} \sqrt{\rho_n \rho_N} \cos \theta_n^* F_{vs} \left( \sqrt{\frac{\rho_N}{\rho_n}}, \frac{|\mathbf{r}_N - \mathbf{r}_n|}{\sqrt{\rho_n \rho_N}} \right). \quad (5.14)$$

We now need to show that $G(\{\rho_n, \mathbf{r}_n\}) > 0$ as then we will have a contradiction with minimality. We will take $\rho_N > 0$ small and investigate $F_{ss}(\alpha, \lambda) = F_{ss}(\alpha^{-1}, \lambda)$ defined in (4.28), where $\alpha = \sqrt{\frac{\rho_N}{\rho_n}}$, $\lambda = \frac{|\mathbf{r}_N - \mathbf{r}_n|}{\sqrt{\rho_N \rho_n}}$. We can use the change of variables $t = \xi\lambda$ to obtain

$$F_{ss}(\alpha, \lambda) = \frac{32}{\pi^2 \lambda^3} \int_0^\infty t^2 J_0(t) K_0(\alpha t / \lambda) K_0(t / (\alpha \lambda)) dt. \quad (5.15)$$
Next we observe that for $\rho_N \ll 1$ the value of $\beta = \frac{\alpha}{\lambda} = \frac{\rho_N}{|\mathbf{r}_N - \mathbf{r}_n|} > 0$ is small, while $\frac{1}{\alpha\lambda} = \frac{\rho_n}{|\mathbf{r}_N - \mathbf{r}_n|}$ is fixed. We have estimates on the Bessel function:

$$0 < K_0(t) < h(t) = \max\{-\ln t + 1, 1\}, \quad \forall t > 0. \qquad (5.16)$$

It is now clear that with the help of $|J_0(t)| \le 1$ we can estimate

$$\begin{aligned} \left| \int_0^\infty t^2 J_0(t) K_0(\alpha t/\lambda) K_0(t/(\alpha\lambda)) dt \right| &\le \int_0^\infty t^2 K_0(t/(\alpha\lambda)) h(\beta t) dt \\ &= \ln \beta^{-1} \int_0^{\beta^{-1}} t^2 K_0(t/(\alpha\lambda)) dt + \int_0^{\beta^{-1}} t^2 K_0(t/(\alpha\lambda)) |\ln t| dt \qquad (5.17) \\ &\quad + \int_0^\infty t^2 K_0(t/(\alpha\lambda)) dt \le C |\ln \rho_N|, \end{aligned}$$

for some $C = C(\rho_n, |\mathbf{r}_N - \mathbf{r}_n|) > 0$ and all $\rho_N > 0$ small enough. Since $\lambda^{-3} = (\frac{\sqrt{\rho_N \rho_n}}{|\mathbf{r}_N - \mathbf{r}_n|})^3$, we obtain

$$\bar{\delta} \frac{\pi^3}{4} \sum_{n=1}^{N-1} \sqrt{\rho_N \rho_n} F_{ss} \left( \sqrt{\frac{\rho_n}{\rho_N}}, \frac{|\mathbf{r}_N - \mathbf{r}_n|}{\sqrt{\rho_N \rho_n}} \right) \ge -C \bar{\delta} \rho_N^2 |\ln \rho_N| \qquad (5.18)$$

for some $C = C(\rho_n, |\mathbf{r}_N - \mathbf{r}_n|) > 0$ and all $\rho_N > 0$ small enough.

We now use the same approach to estimate $F_{vs}$ defined in (4.29):

$$F_{vs}(\alpha, \lambda) = 2 \int_{0}^{\infty} k^2 J_0(\lambda k) K_0(\alpha k) K_1(k/\alpha) dk = \frac{2}{\lambda^3} \int_{0}^{\infty} t^2 J_0(t) K_0(\alpha t/\lambda) K_1(t/(\alpha\lambda)) dt. \qquad (5.19)$$

Using the fact that

$$|F_{vs}(\alpha, \lambda)| \le \frac{2}{\lambda^3} \int_0^\infty t^2 K_1(t/(\alpha\lambda)) h(\beta t) dt \le \frac{C |\ln \rho_N|}{\lambda^3}, \qquad (5.20)$$

where $C = C(\rho_n, |\mathbf{r}_N - \mathbf{r}_n|) > 0$, for all $\rho_N > 0$ sufficiently small, we deduce

$$4\pi\bar{\delta} \sum_{n=1}^{N-1} \sqrt{\rho_n \rho_N} \cos \theta_n^* F_{vs} \left( \sqrt{\frac{\rho_N}{\rho_n}}, \frac{|\mathbf{r}_n - \mathbf{r}_N|}{\sqrt{\rho_n \rho_N}} \right) \ge -C\bar{\delta}\rho_N^2 |\ln \rho_N|. \qquad (5.21)$$

It follows that

$$G(\{\rho_n, \mathbf{r}_n\}) \ge 4\pi\rho_N^2 \ln \left( \frac{e^{1+2\gamma}}{4} \rho_N^2 \right) + \bar{\delta} \frac{\pi^3}{8} \rho_N - C\bar{\delta}\rho_N^2 |\ln \rho_N| > 0 \qquad (5.22)$$

for any $\bar{\delta} > 0$ and all $\rho_N > 0$ small enough.

**Step 4.** It remains to show that $\rho_n^* < L_0^{-1}$ with $L_0 > \bar{L}_0$ for all $1 \le n \le N$. Assume this is not true and, without loss of generality, that $\rho_N^* = L_0^{-1}$. In this case it is not difficult to show that for all $\bar{\delta} < \bar{\delta}_0$ with some $\bar{\delta}_0 > 0$ small enough depending only on $N$ and $L_0$ we have

$$\inf_{\rho_n} F(\{\rho_n, \mathbf{r}_n\}) > \lim_{l \to 0} F(\{\rho_{n,l}, \mathbf{r}_n\}), \qquad (5.23)$$

where $\rho_{n,l} = \rho_n^*$ for all $1 \le n \le N - 1$ and $\rho_{N,l} \to 0$ as $l \to \infty$. Therefore, $\rho_N = L_0^{-1}$ cannot be a minimizer and $\rho_n^* \in (0, L_0^{-1})$ for all $1 \le n \le N$. Finally, as the function $F_N(\{\rho_n, \mathbf{r}_n\})$ is continuous for $\rho_n \in (0, L_0^{-1})$, the minimum is attained in $A_N(\{\mathbf{r}_n\})$ for all $\bar{\delta} < \bar{\delta}_0$. $\square$

Having established existence of minimizers of $F_N(\{\rho_n, \theta_n, \mathbf{r}_n\})$ over $A_N(\{\mathbf{r}_n\})$, it would next be interesting to understand the nature of the minimizers of $F_N(\{\mathbf{r}_n\})$ with respect to the positions $\mathbf{r}_n$ of the skyrmions in each layer. This, however, is in general a daunting task, as the interaction energy describes a fully coupled, strongly nonlinear system, in which the dependence on positions arises, in addition to the direct interaction term, due to the implicit and a priori unknown dependence of the skyrmion radii $\rho_n$ and angles $\theta_n$ on the skyrmion positions $\{\mathbf{r}_n\}$. In particular, it is not a priori clear whether minimizers of $F_N(\{\mathbf{r}_n\})$ over the positions exist, as the interaction between different skyrmions could be repulsive, leading to failure of compactness of the minimizing sequences and, as a consequence, to failure of existence of minimizers.
## 6 Application to the case of bilayers in the absence of DMI

Instead of treating the problem in its full generality, in the remainder of this paper we focus on the simplest particular case of stray field-coupled ferromagnetic bilayers in the absence of DMI. Here we can obtain a complete characterization of the energy minimizers. We show that minimizers indeed exist, signifying an attractive interaction of the skyrmions in the adjacent layers. We find that the energy minimizers consist of pairs of concentric identical skyrmions of Néel type, except that their in-plane magnetizations are anti-parallel. Furthermore, these skyrmion pairs are chiral, as the sense of the magnetization rotation in each of the layers is fixed by the direction of the magnetization at infinity. This is in contrast with the case of ferromagnetic monolayers, which are known to support stray field-stabilized skyrmions of Bloch type that can have two chiralities.

We now define the energy by setting $N = 2$ and $\bar{\kappa} = 0$ in (5.2):

$$ F_2(\{\rho_n, \theta_n, \mathbf{r}_n\}) = \sum_{n=1}^{2} \left[ -4\pi\rho_n^2 \ln\left(\frac{e^{1+2\gamma}}{4}\rho_n^2\right) + \bar{\delta}\frac{\pi^3}{8}\rho_n(3\cos^2\theta_n - 1) \right. \\ \left. + \bar{\delta}\frac{\pi^3}{4}\sqrt{\rho_1\rho_2} \left[ 3\cos\theta_1\cos\theta_2 F_{vv}\left(\sqrt{\frac{\rho_2}{\rho_1}}, \frac{|\mathbf{r}_2 - \mathbf{r}_1|}{\sqrt{\rho_1\rho_2}}\right) - F_{ss}\left(\sqrt{\frac{\rho_2}{\rho_1}}, \frac{|\mathbf{r}_2 - \mathbf{r}_1|}{\sqrt{\rho_1\rho_2}}\right) \right] \right. \\ \left. - 4\pi\bar{\delta}\sqrt{\rho_1\rho_2} \left[ \cos\theta_1 F_{vs}\left(\sqrt{\frac{\rho_2}{\rho_1}}, \frac{|\mathbf{r}_2 - \mathbf{r}_1|}{\sqrt{\rho_1\rho_2}}\right) - \cos\theta_2 F_{vs}\left(\sqrt{\frac{\rho_1}{\rho_2}}, \frac{|\mathbf{r}_2 - \mathbf{r}_1|}{\sqrt{\rho_1\rho_2}}\right) \right]. \quad (6.1) $$

and would like to investigate the energy minimizing configurations of stacked skyrmion pairs in a bilayer. We have the following result that gives a complete characterization of the minimizers of this problem.

**Theorem 2.** Let $L_0 > \bar{L}_0$. Then there exists $\bar{\delta}_0 > 0$ such that for every $\bar{\delta} < \bar{\delta}_0$ the minimizers of the energy $F_2(\{\rho_n, \theta_n, \mathbf{r}_n\})$ among $\rho_{1,2} \in (0, L_0^{-1})$, $\theta_{1,2} \in [-\pi, \pi)$ and $\mathbf{r}_{1,2} \in \mathbb{R}^2$ exist and satisfy

i) $\mathbf{r}_1 = \mathbf{r}_2$;

ii) $\theta_1 = 0, \theta_2 = -\pi$;

iii) $\rho_1 = \rho_2 = \rho$, where

$$ \rho = \frac{(16 + \pi^2)\bar{\delta}}{-64W_{-1}\left(-\frac{16+\pi^2}{128}e^{1+\gamma}\bar{\delta}\right)}, \quad (6.2) $$

and $W_{-1}(t)$ is the Lambert W function.

*Proof.* We first observe that by Lemmas 1 and 2 in the Appendix the absolute values of the functions $F_{vv}(\alpha, \lambda)$, $F_{ss}(\alpha, \lambda)$ and $F_{vs}(\alpha, \lambda)$ are maximized at $\lambda = 0$ for fixed $\alpha > 0$. Moreover, at $\lambda = 0$ all these functions are positive. It follows, that

$$ F_2(\{\rho_n, \theta_n, \mathbf{r}_n\}) \geq \sum_{n=1}^{2} \left[ -4\pi\rho_n^2 \ln\left(\frac{e^{1+2\gamma}}{4}\rho_n^2\right) + \bar{\delta}\frac{\pi^3}{8}\rho_n(3\cos^2\theta_n - 1) \right. \\ \left. - \bar{\delta}\frac{\pi^3}{4}\sqrt{\rho_1\rho_2} \left[ 3|\cos\theta_2| |\cos\theta_1| F_{vv}\left(\sqrt{\frac{\rho_2}{\rho_1}}, 0\right) + F_{ss}\left(\sqrt{\frac{\rho_2}{\rho_1}}, 0\right) \right] \right. \\ \left. - 4\pi\bar{\delta}\sqrt{\rho_1\rho_2} \left[ |\cos\theta_1| F_{vs}\left(\sqrt{\frac{\rho_2}{\rho_1}}, 0\right) + |\cos\theta_2| F_{vs}\left(\sqrt{\frac{\rho_1}{\rho_2}}, 0\right) \right], \quad (6.3) $$

and equality is achieved if and only if $\mathbf{r}_1 = \mathbf{r}_2$, $\cos\theta_1 \ge 0$ and $\cos\theta_2 \le 0$.

We now observe by completing the square that the function

$$ h_0(\theta_1, \theta_2) = \frac{3\pi^3}{8}\bar{\delta}\left(\sum_{n=1}^{2} \rho_n \cos^2 \theta_n - 2\sqrt{\rho_1\rho_2} F_{vv}\left(\sqrt{\frac{\rho_2}{\rho_1}}, 0\right) |\cos \theta_1| |\cos \theta_2|\right) \quad (6.4) $$

is non-negative due to the fact that $F_{vv}(\alpha, 0) \le 1$ with equality achieved only at $\alpha = 1$ (see Lemma 1 in the Appendix). Moreover, for $\rho_n > 0$ the function $h_0(\theta_1, \theta_2) = 0$ if and only if one of the two alternatives below holds:
<figure><img src="image_5.png" /><figcaption>Figure 4: Plot of $F(\rho)$ when $\rho \in (0, \rho_0)$ with $\rho_0 = \bar{L}_0^{-1}$ for $\bar{\delta} = 0.25$.</figcaption></figure>

i) $\rho_1 \neq \rho_2$ and $|\theta_{1,2}| = \frac{\pi}{2}$;

ii) $\rho_1 = \rho_2$ and $|\cos \theta_1| = |\cos \theta_2|$.

We also observe that

$$ -|\cos \theta_1| F_{vs} \left( \sqrt{\frac{\rho_2}{\rho_1}}, 0 \right) - |\cos \theta_2| F_{vs} \left( \sqrt{\frac{\rho_1}{\rho_2}}, 0 \right) \geq -F_{vs} \left( \sqrt{\frac{\rho_2}{\rho_1}}, 0 \right) - F_{vs} \left( \sqrt{\frac{\rho_1}{\rho_2}}, 0 \right), \quad (6.5) $$

with equality achieved at $\theta_{1,2} = 0$ or $\theta_{1,2} = \pi$. Minimizing the right-hand side of the above expression in $\rho_1, \rho_2$ we obtain $\rho_1 = \rho_2$ and hence the following lower bound achievable if $\rho_1 = \rho_2$:

$$ -F_{vs} \left( \sqrt{\frac{\rho_2}{\rho_1}}, 0 \right) - F_{vs} \left( \sqrt{\frac{\rho_1}{\rho_2}}, 0 \right) \geq -2, \quad (6.6) $$

where we used Lemma 2 in the Appendix. Recalling that $0 < F_{ss}(\alpha, 0) \leq 1$ with equality achievable if and only if $\alpha = 1$ (see Lemma 1 of the Appendix) and combining our findings above, we obtain

$$ \begin{aligned} F_2(\{\rho_n, \theta_n, \mathbf{r}_n\}) &\geq \sum_{n=1}^{2} \left[ -4\pi\rho_n^2 \ln \left( \frac{e^{1+2\gamma}}{4}\rho_n^2 \right) - \bar{\delta}\frac{\pi^3}{8}\rho_n \right] - \bar{\delta}\frac{\pi^3}{4}\sqrt{\rho_1\rho_2} - 8\pi\bar{\delta}\sqrt{\rho_1\rho_2} \\ &\geq \sum_{n=1}^{2} \left[ -4\pi\rho_n^2 \ln \left( \frac{e^{1+2\gamma}}{4}\rho_n^2 \right) - \bar{\delta}\frac{\pi^3}{4}\rho_n - 4\pi\bar{\delta}\rho_n \right] \geq 2F(\rho), \end{aligned} \quad (6.7) $$

with equality achieved if and only if $\mathbf{r}_1 = \mathbf{r}_2$, $\theta_1 = 0$ and $\theta_2 = \pi$, and $\rho_1 = \rho_2 = \rho$, where

$$ F(\rho) = -4\pi\rho^2 \ln \left( \frac{e^{1+2\gamma}}{4}\rho^2 \right) - \bar{\delta}\frac{\pi^3}{4}\rho - 4\pi\bar{\delta}\rho. \quad (6.8) $$

The graph of $F(\rho)$ is illustrated in Fig. 4.

It is now clear that

$$ \inf_{\rho_n} \min_{\theta_n, \mathbf{r}_n} F_2(\{\rho_n, \theta_n, \mathbf{r}_n\}) = 2 \inf_{\rho} F(\rho). \quad (6.9) $$

We can minimize $F(\rho)$ in the admissible interval of $\rho \in (0, L_0^{-1})$. It is easy to check that for $L_0 > \bar{L}_0$ and $\bar{\delta}_0 > 0$ sufficiently small depending only on $L_0$ the minimum of $F(\rho)$ is attained at the point $\rho \in (0, L_0^{-1})$ satisfying

$$ 0 = F'(\rho) = -8\pi\rho \left[ \ln \left( \rho^2 \frac{e^{1+2\gamma}}{4} \right) + 1 \right] - \frac{\bar{\delta}\pi^3}{4} - 4\pi\bar{\delta}. \quad (6.10) $$

After a little algebra, one can see that the solution of this equation in the interval $(0, \rho_0)$, where $\rho_0 = \bar{L}_0^{-1}$, is given by (6.2). This concludes the proof. $\square$
To summarize, we have demonstrated that the global energy minimizers of $F_2(\{\rho_n, \theta_n, \mathbf{r}_n\})$ are characterized by a certain symmetry that makes the in-plane components of the magnetization in (4.8) anti-parallel. It is instructive to see what this assumption would lead to on the level of the original energy in (4.2), before introducing the truncated BP ansatz. Setting $(\mathbf{m}_1^{\perp}, m_1^{\|}) = (-\mathbf{m}_2^{\perp}, m_2^{\|}) = \mathbf{m}$, we see that the energy $\bar{E}_2$ of the configuration $\{\mathbf{m}_1, \mathbf{m}_2\}$ becomes $\bar{E}_2(\{\mathbf{m}_1, \mathbf{m}_2\}) = 2\bar{E}^{\pm}(\mathbf{m})$, where

$$ \bar{E}^{\pm}(\mathbf{m}) = \int_{\mathbb{R}^2} (|\nabla \mathbf{m}|^2 + |\mathbf{m}^{\perp}|^2 - \bar{\delta}\mathbf{m}^{\perp} \cdot \nabla m^{\|}) d^2r - 2\bar{\delta} \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{(m^{\|}(\mathbf{r}) - m^{\|}(\mathbf{r}'))^2}{8\pi|\mathbf{r} - \mathbf{r}'|^3} d^2r d^2r'. \quad (6.11) $$

One can observe that for these configurations the volume-surface interactions act as an effective interfacial DMI term, favoring the Néel rotation of the magnetization with a particular rotation sense. At the same time, the volume-surface and surface-surface interactions act constructively to stabilize the skyrmion pair, while the penalizing volume-volume charge interaction is absent. The reason for the stabilizing action of the volume-surface interaction may be seen from Fig. 5. For skyrmions with anti-parallel in-plane magnetizations in the two layers and counter-clockwise rotation in the bottom layer the volume-surface interaction energy is lower than that of all other possible skyrmion configurations with anti-parallel in-plane magnetizations and the same out-of-plane component. Notice that the stray field energy due to the volume-volume interactions may be forced to be zero by choosing the Bloch rotation in both layers instead (not necessarily anti-parallel). However, in that case the DMI-like term due to the volume-surface interaction does not contribute to the energy of the skyrmion pair, making the Bloch rotation as compared to the Néel rotation with anti-parallel in-plane components of the magnetization less favorable. Finally, notice that the configuration in Fig. 5(f) that corresponds to the lowest energy is reminiscent of a flux closure structure in bulk ferromagnets.

To verify the predictions of our asymptotic analysis, we carried out a numerical study of skyrmion profiles in stray field-coupled ferromagnetic bilayers, using MuMax3 software [34]. For the material parameters, we chose those corresponding to a ferrimagnetic material such as GdCo with the material parameters $A = 20$ pJ/m, $M_s = 10^5$ A/m, $K_u = 6700$ J/m$^3$ [11]. We consider two 5 nm-thick layers of such a material with a non-magnetic separator of negligible thickness. These parameters give an exchange length $\ell_{ex} \simeq 56.4$ nm and a small dimensionless layer thickness $\delta \simeq 0.089$, justifying the use of the reduced model appropriate for ultrathin films. The quality factor associated with the uniaxial magnetocrystalline anisotropy is $Q \simeq 1.066$, giving the parameter $\bar{\delta} \simeq 0.344$ characterizing the strength of the stray field interaction that is also within the validity range of the asymptotic theory in $\bar{\delta} \ll 1$. For this value of $\bar{\delta}$, the formula in (6.2) predicts $\rho \simeq 0.0938$, which corresponds to the dimensional skyrmion radius of 20.5 nm. The resulting profiles obtained, using the minimize function of MuMax3 on a 2048 x 2048 x 2 grid with the in-plane discretization steps $\Delta x = \Delta y = 0.5$ nm, the out-of-plane step equal to $d = 5$ nm and the number of repeats in $(X, Y, Z)$ set to $(5, 5, 0)$ to approximate the periodic boundary conditions, are presented in Fig. 6. In a very good agreement with the theoretical prediction, the obtained profiles consist of a pair of concentric Néel skyrmions that are very close to the Belavin-Polyakov profiles with radius 20 nm obtained by fitting the numerical profiles to the Belavin-Polyakov profiles.

We conclude our study by attempting to calculate the function $F_2(\{\mathbf{r}_1, \mathbf{r}_2\})$ from Theorem 1 that determines the interaction energy of two skyrmions in the adjacent layers of a bilayer. By the translational and rotational symmetries of the problem this function depends only on $r = |\mathbf{r}_1 - \mathbf{r}_2|$ and is obtained by minimizing the function $F_2(\{\rho_n, \theta_n, \mathbf{r}_n\})$ in four parameters $\rho_{1,2}$ and $\theta_{1,2}$. Nevertheless, the problem is still intractable analytically, since no manageable closed form analytical expressions for the functions $F_{vv}(\alpha, \lambda)$, $F_{ss}(\alpha, \lambda)$ and $F_{vs}(\alpha, \lambda)$ are available for $\lambda > 0$ and general $\alpha > 0$. It would seem plausible, and is supported by the numerical evaluation of the energy at a few points of the parameter space that the minimizers from Theorem 1 in the case of bilayers, $N = 2$, would exhibit a symmetry for the skyrmion radius with respect to the layer position, which is also certainly true for the global energy minimizers by Theorem 2. Thus, it is reasonable to conjecture that in a minimizer from Theorem 1 at fixed $r > 0$ we have $\rho_1 = \rho_2 = \rho$, and thus $F_2(\{\mathbf{r}_1, \mathbf{r}_2\}) = F_2^{\text{sym}}(r)$, where $F_2^{\text{sym}}(r)$
<figure><img src="image_6.png" /><figcaption>Figure 5: An illustration of the stray field interactions between the volume and surface charges in a bilayer for a skyrmion with anti-parallel in-plane magnetization components (i.e., with $m_1^\parallel = m_2^\parallel$ and $m_1^\perp = -m_2^\perp$): (a,c,e) clockwise rotation in the bottom layer; (b,d,f) anti-clockwise rotation in the bottom layer. The volume charge density $\rho_m^{\text{vol}} = -\nabla_\perp \cdot \mathbf{m}^\perp$ in one layer is indicated in blue (negative) and red (positive), and its associated magnetic field lines are shown by the lines with arrows going from red to blue regions. Only the out-of-plane component $m_n^\parallel$ of the magnetization in the other layer that contributes to the volume-surface interaction is shown in (a-d), while the corresponding full magnetization profiles are shown in (e,f). In (a,c), the out-of-plane magnetic moments of one layer point agains the field lines from the other layer, resulting in a higher energy. In (b,d), the out-of-plane magnetic moments of one layer point along the field lines from the other layer, resulting in a lower energy. The total surface-surface interaction energies are identical in both cases, and the total volume-volume interaction energy is asymptotically zero.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Stray field-stabilized Néel skyrmion in an ultrathin ferromagnetic bilayer. (a) The out-of-plane components $m_{1,2}^{\parallel}$ of the magnetization in the bottom and top layers, respectively, along the horizontal line through the skyrmion center. (b) The in-plane components $\hat{\mathbf{x}} \cdot \mathbf{m}_{1,2}^{\perp}$ in the bottom and top layers, respectively, along the horizontal line through the skyrmion center. (c) and (d) The top view of the magnetization in the bottom and top layers, respectively. Results of the MUMax3 simulations on the 2048 × 2048 × 2 grid with the in-plane discretization steps $\Delta x = \Delta y = 0.5$ nm, the out-of-plane discretization step $\Delta z$ equal to the single layer thickness $d = 5$ nm, and periodic boundary conditions in the plane. The material parameters are $A = 20$ pJ/m, $M_s = 10^5$ A/m, $K_u = 6700$ J/m$^3$, corresponding to a ferrimagnetic material, with a non-magnetic spacer of negligible thickness. In all the panels, only a 400 nm × 400 nm region around the skyrmion center is shown. In (a) and (b), the Néel Belavin-Polyakov profiles with radius 20 nm are also shown by thin red lines.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Plot of $-\ln[-F_2^{\text{sym}}(\rho, r)]$ for $\bar{\delta} = 0.25$. The choice of the function plotted helps to visualize the part of the parameter space in which the energy is negative.</figcaption></figure>

can now be calculated in closed form. Indeed, we have

$$ F_{vv}(1, \lambda) = \frac{32}{3\pi^2\lambda^2(\lambda^2+4)^{3/2}} \left[ \begin{aligned} & \left(\lambda^2\sqrt{\lambda^2+4} + 4\sqrt{\lambda^2+4} + 8\right) K^2\left(\frac{1}{2} - \frac{\sqrt{\lambda^2+4}}{4}\right) \\ & -4\left(\lambda^2\sqrt{\lambda^2+4} + 2\sqrt{\lambda^2+4} + 4\right) K\left(\frac{1}{2} - \frac{\sqrt{\lambda^2+4}}{4}\right) E\left(\frac{1}{2} - \frac{\sqrt{\lambda^2+4}}{4}\right) \\ & +4(\lambda^2+2)\sqrt{\lambda^2+4} E^2\left(\frac{1}{2} - \frac{\sqrt{\lambda^2+4}}{4}\right) \end{aligned} \right], \quad (6.12) $$

where, again, $K(m)$ and $E(m)$ are the complete elliptic integrals of the first and second kind, respectively. Similarly

$$ F_{ss}(1, \lambda) = -\frac{32}{\pi^2 \lambda^2 (\lambda^2 + 4)^{3/2}} \left[ \begin{aligned} & \left(\left(\sqrt{\lambda^2 + 4} + 4\right) \lambda^2 + 4\left(\sqrt{\lambda^2 + 4} + 2\right)\right) K^2\left(\frac{1}{2} - \frac{\sqrt{\lambda^2 + 4}}{4}\right) \\ & -8\left(\lambda^2 + \sqrt{\lambda^2 + 4} + 2\right) K\left(\frac{1}{2} - \frac{\sqrt{\lambda^2 + 4}}{4}\right) E\left(\frac{1}{2} - \frac{\sqrt{\lambda^2 + 4}}{4}\right) \\ & +8\sqrt{\lambda^2 + 4} E^2\left(\frac{1}{2} - \frac{\sqrt{\lambda^2 + 4}}{4}\right) \end{aligned} \right], \quad (6.13) $$

and

$$ F_{vs}(1, \lambda) = 2 \left( \frac{1}{\lambda^2 + 4} + \frac{4 \sinh^{-1}\left(\frac{\lambda}{2}\right)}{\lambda (\lambda^2 + 4)^{3/2}} \right). \quad (6.14) $$

With these expressions for $F_{vv}(\alpha, \lambda)$, $F_{ss}(\alpha, \lambda)$ and $F_{vs}(\alpha, \lambda)$, we can proceed to minimize the energy in (6.1) with respect to the angles $\theta_{1,2}$ to obtain a closed form expression for $F_2^{\text{sym}}(\rho, r) = \min_{\theta_1, \theta_2} F_2(\rho_1, \theta_1, \rho_2, \theta_2, r)$ with $\rho_1 = \rho_2 = \rho$. The plot of this function for a particular value of $\bar{\delta}$ is presented in Fig. 7.

The function $F_2^{\text{sym}}(\rho, r)$ can finally be minimized numerically for a given value of $r > 0$. The results for a particular choice $\bar{\delta} = 0.25$ are presented in Fig. 8, where the minimal energy, the optimal
<figure><img src="image_9.png" /><figcaption>Figure 8: The result of the global numerical minimization of $F_2^{\text{sym}}(\rho, r)$ in $\rho$ for $\bar{\delta} = 0.25$: (a) the minimum energy $F_2^{\text{sym}}(r)$; (b) zooming in on the minimum energy $F_2^{\text{sym}}(r)$ to show the repulsion at large separation distances; (c) optimal value of $\rho$ as a function of $r$; (d) optimal value of $\cos\theta_1$ as a function of $r$.</figcaption></figure>

value of $\rho$ and the optimal value of $\cos\theta_1 = -\cos\theta_2 > 0$ (the latter is due to the fact that from $0 < F_{vv}(1, \lambda) < 1$ for $\lambda > 0$ follows that the energy is a strictly convex function of $\cos\theta_1$ and $\cos\theta_2$) are plotted. We observe that the minimum of the energy is indeed attained at $r = 0$, as it should be by Theorem 2. As the value of $r$ increases, the equilibrium radius $\rho$ also slightly increases, while the optimal angles remain those of the minimizer: $\theta_1 = 0$ and $\theta_2 = -\pi$, in Theorem 2. We note that as can be seen from Fig. 7, the energy landscape in the $(\rho, r)$ plane consists of a wide valley for not too big values of $\rho \sim r$ and a narrow gorge providing an escape path to $r \to \infty$ and corresponding to an almost constant value of $\rho$ and increasing values of $r$. At $r \gg 1$ the latter corresponds to a pair of non-interacting skyrmions in each layer.

When the value of $r$ is increased from zero, at first the minimum is found in the wide valley, but as the parameters push the energy towards the steep wall at a certain critical value of $r$ (close to 0.054 for $\bar{\delta} = 0.25$) the minimizer jumps into the gorge and continues to follow along it. Notice that for this value of $\bar{\delta}$ the angles switch abruptly from those corresponding to the Néel skyrmions in a global energy minimizer to $|\theta_{1,2}| \simeq \pi/2$ corresponding to a pair of Bloch skyrmions. Also notice that the radius of the skyrmions in the gorge is an order of magnitude smaller than that of the global minimizer, and the absolute value of its energy is two orders of magnitude lower, respectively. In sum, the skyrmions remain of Néel type and exhibit a strongly attractive interaction within a certain range of separation distances. Beyond that range they abruptly change their nature and behave as a pair of weakly interacting Bloch skyrmions, and at large separations these skyrmions exhibit a weak repulsive interaction dominated by the dipolar forces. The strongly nonlinear attractive interaction of a pair of skyrmions at sufficiently short distances thus may provide an important stabilization mechanism that is highly desirable for spintronic applications.
# A Appendix

In this section, we establish a few basic properties of the functions $F_{vv}(\alpha, \lambda)$, $F_{ss}(\alpha, \lambda)$ and $F_{vs}(\alpha, \lambda)$ necessary in proving our theorems. An impatient reader may quickly convince oneself about these facts by simply plotting those functions evaluated numerically from their respective integral definitions.

**Lemma 1.** For $\alpha \in (0, 1)$, let $F_{vv}(\alpha, 0)$ and $F_{ss}(\alpha, 0)$ be defined by (4.30) and (4.31), respectively. Then $0 < F_{vv}(\alpha, 0) < 1$ and $0 < F_{ss}(\alpha, 0) < 1$. Furthermore, for all $\alpha > 0$ and $\lambda \ge 0$

$$|F_{vv}(\alpha, \lambda)| \le F_{vv}(\alpha, 0), \quad |F_{ss}(\alpha, \lambda)| \le F_{ss}(\alpha, 0), \qquad (A.1)$$

where $F_{vv}(\alpha, \lambda)$ and $F_{ss}(\alpha, \lambda)$ are defined in (4.27) and (4.28), and $F_{vv}(\alpha, \lambda)$ and $F_{ss}(\alpha, \lambda)$ are uniquely maximized by $(\alpha, \lambda) = (1, 0)$, with $\max F_{vv} = \max F_{ss} = 1$.

*Proof.* The statements about $F_{vv}(\alpha, \lambda)$ and $F_{ss}(\alpha, \lambda)$ follow immediately from the first part of the lemma, since from the definitions of $F_{vv}(\alpha, \lambda)$ and $F_{ss}(\alpha, \lambda)$, and from the fact that $|J_0(\lambda\xi)| < 1$ for all $\lambda > 0$ and $\xi > 0$ we have $|F_{vv}(\alpha, \lambda)| < F_{vv}(\alpha, 0)$ and $|F_{ss}(\alpha, \lambda)| < F_{ss}(\alpha, 0)$ for all $\alpha > 0$ and $\lambda > 0$, together with the facts that $F_{vv}(\alpha, 0) = F_{vv}(\alpha^{-1}, 0)$ and $F_{ss}(\alpha, 0) = F_{ss}(\alpha^{-1}, 0)$.

We now prove that $F_{ss}(\alpha, 0)$ attains its maximum only at $\alpha = 1$. Since $F_{ss}(\alpha, 0) = F_{ss}(\alpha^{-1}, 0)$ and $F_{ss}(1, 0) = 1$, we only need to investigate the interval of $\alpha \in (0, 1)$. Recalling the definition of $F_{ss}(\alpha, 0)$ and using the change of variables $t = 1 - \alpha^4$ with $t \in (0, 1)$, we note that

$$F_{ss}((1-t)^{\frac{1}{4}}, 0) = \frac{16(1-t)^{\frac{3}{4}}}{\pi t^2} ((2-t)K(t) - 2E(t)). \qquad (A.2)$$

We now need to show that $f(t) = F_{ss}((1-t)^{\frac{1}{4}}, 0)$ has a unique maximum at $t = 0$. We can differentiate to obtain

$$f'(t) = -\frac{4}{\pi t^3 (1-t)^{\frac{1}{4}}} ((3t^2 - 16t + 16)K(t) + 8(t-2)E(t)). \qquad (A.3)$$

If we show that $f_1(t) = (3t^2 - 16t + 16)K(t) + 8(t-2)E(t) > 0$ on $(0, 1)$ then it follows that $f'(t) < 0$ on $(0, 1)$, and hence the maximum of $F_{ss}(\alpha, 0)$ is achieved at $\alpha = 1$.

We note that $f_1(0) = 0$ and

$$f'_{1}(t) = -\frac{3}{2(1-t)}((3t^2 - 11t + 8)K(t) + (7t - 8)E(t)). \qquad (A.4)$$

Hence we need to show that $f_2(t) = (3t^2 - 11t + 8)K(t) + (7t - 8)E(t) < 0$. Again, we note that $f_2(0) = 0$ and

$$f'_{2}(t) = \frac{9}{2}(t-2)K(t) + 9E(t). \qquad (A.5)$$

Hence we need to show that $f_3(t) = \frac{9}{2}(t-2)K(t) + 9E(t) < 0$. We note that $f_3(0) = 0$ and

$$f'_{3}(t) = -\frac{9}{4(1-t)}(K(t)(t-1) + E(t)). \qquad (A.6)$$

Finally, we need to show that $f_4(t) = K(t)(t-1) + E(t) > 0$. We note that $f_4(0) = 0$ and

$$f'_{4}(t) = \frac{1}{2}K(t) > 0. \qquad (A.7)$$

The result is proved.

We now prove that $F_{vv}(\alpha, 0)$ attains its maximum only at $\alpha = 1$. Similarly, we consider only the interval of $\alpha \in (0, 1)$. Recalling the definition of $F_{vv}(\alpha, 0)$ and using the change of variables $t = 1 - \alpha^4$ with $t \in (0, 1)$, we obtain

$$F_{vv}((1-t)^{\frac{1}{4}}, 0) = \frac{16(1-t)^{\frac{1}{4}}}{3\pi t^2} ((2-t)E(t) - 2(1-t)K(t)). \qquad (A.8)$$
We now need to show that $g(t) = F_{vv}((1-t)^{\frac{1}{4}}, 0)$ has maximum at $t = 0$. We can differentiate to obtain

$$g'(t) = \frac{4}{3\pi t^3(1-t)^{\frac{3}{4}}} (8(t^2 - 3t + 2)K(t) - (t^2 - 16t + 16)E(t)). \quad (A.9)$$

If we show that $g_1(t) = 8(t^2 - 3t + 2)K(t) - (t^2 - 16t + 16)E(t) < 0$ on $(0, 1)$ then it follows that $g'(t) < 0$ on $(0, 1)$, and, hence, the maximum of $F_{vv}(\alpha, 0)$ is achieved at $\alpha = 1$.

Arguing as in the case of $f(t)$, we note that $g_1(0) = 0$ and

$$g'_{1}(t) = -\frac{5}{2}((8 - 5t)K(t) + (t - 8)E(t)). \quad (A.10)$$

We now need to show that $g_2(t) = -((8 - 5t)K(t) + (t - 8)E(t)) < 0$. Again, we note that $g_2(0) = 0$ and

$$g'_{2}(t) = -\frac{3}{2(t-1)}(2(1-t)K(t) + (t-2)E(t)). \quad (A.11)$$

Hence, we now need to show that $g_3(t) = -(2(1-t)K(t) + (t-2)E(t)) > 0$. We note that $g_3(0) = 0$ and

$$g'_{3}(t) = \frac{3}{2}(K(t) - E(t)) > 0. \quad (A.12)$$

The result is proved. □

**Lemma 2.** For $\alpha > 0$, let $F_{vs}(\alpha, 0)$ be defined by (4.32). Then $0 < F_{vs}(\alpha, 0) \le 2$, and for all $\lambda \ge 0$

$$|F_{vs}(\alpha, \lambda)| \le F_{vs}(\alpha, 0), \quad (A.13)$$

where the function $F_{vs}(\alpha, \lambda)$ is defined in (4.29). Furthermore,

$$|F_{vs}(\alpha, \lambda)| + |F_{vs}(\alpha^{-1}, \lambda)| \le 2, \quad (A.14)$$

with equality achieved only at $(\alpha, \lambda) = (1, 0)$.

*Proof.* As in the case with $F_{vv}(\alpha, \lambda)$ and $F_{ss}(\alpha, \lambda)$, the statements about $F_{vs}(\alpha, \lambda)$ follow, once we demonstrate the desired properties of $F_{vs}(\alpha, 0)$, since by the same argument $|F_{vs}(\alpha, \lambda)| \le F_{vs}(\alpha, 0)$. From the definition of the latter, it is clear that the function $F_{vs}(\alpha, 0)$ is positive and continuous for all $\alpha > 0$, including at $\alpha = 1$, since $\lim_{\alpha \to 1} F_{vs}(\alpha, 0) = 1$. Therefore, since $F_{vs}(\alpha, 0) \to 0$ as $\alpha \to 0$ or $\alpha \to \infty$, the function $F_{vs}(\alpha, 0)$ is uniformly bounded.

To show the inequality in (A.14), we observe that

$$\begin{aligned} F_{vs}(\alpha, 0) + F_{vs}(\alpha^{-1}, 0) &= \frac{2\alpha(\alpha^2 - 1)(\alpha^4 + 4\alpha^2 \ln \alpha - 1)}{(\alpha^4 - 1)^2} \\ &= \sech \omega(1 + \omega \sech \omega \operatorname{csch} \omega), \end{aligned} \quad (A.15)$$

where $\omega = \ln \alpha$, which is manifestly maximized at $\omega = 0$ among all $\omega \in \mathbb{R}$. Finally, the upper bound on $F_{vs}(\alpha, 0)$ follows directly from (A.14). □

## References

[1] A. A. Belavin and A. M. Polyakov. Metastable states of two-dimensional isotropic ferromagnets. *JETP Lett.*, 22:245–247, 1975.

[2] A. Bellec, S. Rohart, M. Labrune, J. Miltat, and A. Thiaville. Domain wall structure in magnetic bilayers with perpendicular anisotropy. *Europhys. Lett.*, 91:17009, 2010.
[3] A. Bernand-Mantel, A. Fondet, S. Barnova, T. M. Simon, and C. B. Muratov. Theory of magnetic field-stabilized compact skyrmions in thin film ferromagnets. *Phys. Rev. B*, 108:L161405, 2023.

[4] A. Bernand-Mantel, C. B. Muratov, and T. M. Simon. Unraveling the role of dipolar versus Dzyaloshinskii-Moriya interactions in stabilizing compact magnetic skyrmions. *Phys. Rev. B*, 101:045416, 2020.

[5] A. Bernand-Mantel, C. B. Muratov, and T. M. Simon. A quantitative description of skyrmions in ultrathin ferromagnetic films and stability of degree ±1 harmonic maps from $\mathbb{R}^2$ to $\mathbb{S}^2$. *Arch. Rat. Mech. Anal.*, 239:219-299, 2021.

[6] A. Bogdanov and A. Hubert. Thermodynamically stable magnetic vortex states in magnetic crystals. *J. Magn. Magn. Mater.*, 138:255-269, 1994.

[7] A. N. Bogdanov, M. V. Kudinov, and D. A. Yablonskii. Theory of magnetic vortices in easy-axis ferromagnets. *Sov. Phys. - Solid State*, 31:1707-1710, 1989.

[8] A. N. Bogdanov and D. A. Yablonskii. Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets. *Sov. Phys. - JETP*, 68:101-103, 1989.

[9] W. F. Brown. Micromagnetics. Interscience Tracts of Physics and Astronomy 18. Interscience Publishers (Wiley & Sons), 1963.

[10] F. Büttner, I. Lemesh, and G. S. D. Beach. Theory of isolated magnetic skyrmions: From fundamentals to room temperature applications. *Scientific Reports*, 8:4464, 2018.

[11] L. Caretta, L. Mann, F. Büttner, K. Ueda, B. Pfau, C. M. Günther, P. Hessing, A. Churikova, C. Klose, M. Schneider, D. Engel, C. Marcus, D. Bono, K. Bagschik, S. Eisebitt, and G. S. D. Beach. Fast current-driven domain walls and small skyrmions in a compensated ferrimagnet. *Nature Nanotechnol.*, 13:1154-1160, 2018.

[12] G. Di Fratta, C. B. Muratov, F. N. Rybakov, and V. V. Slastikov. Variational principles of micromagnetics revisited. *SIAM J. Math. Anal.*, 52:3580-3599, 2020.

[13] G. Di Fratta, C. B. Muratov, and V. V. Slastikov. Reduced energies for thin ferromagnetic films with perpendicular anisotropy. *Math. Models Methods Appl. Sci.*, 34:1861-1904, 2024.

[14] Y. Dovzhenko, F. Casola, S. Schlotter, T. X. Zhou, F. Büttner, R. L. Walsworth, G. S. D. Beach, and A. Yacoby. Magnetostatic twists in room-temperature skyrmions explored by nitrogen-vacancy center spin texture reconstruction. *Nature Commun.*, 9:2712, 2018.

[15] I. Dzyaloshinskii. A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics. *J. Phys. Chem. Solids*, 4:241-255, 1958.

[16] J. Eells and L. Lemaire. Another report on harmonic maps. *Bull. Lond. Math. Soc.*, 20:385-524, 1988.

[17] K. Fallon, S. McVitie, W. Legrand, F. Ajejas, D. Maccariello, S. Collin, V. Cros, and N. Reyren. Quantitative imaging of hybrid chiral spin textures in magnetic multilayer systems by Lorentz microscopy. *Phys. Rev. B*, 100:214431, 2019.

[18] A. Fert, V. Cros, and J. Sampaio. Skyrmions on the track. *Nature Nanotechnol.*, 8:152-156, 2013.

[19] A. Fert, N. Reyren, and V. Cros. Magnetic skyrmions: advances in physics and potential applications. *Nat. Rev. Mater.*, 2:17031, 2017.

[20] G. Finocchio, M. Di Ventra, K. Y. Camsari, K. Everschor-Sitte, P. Khalili Amiri, and Z. Zeng. The promise of spintronics for unconventional computing. *J. Magn. Magn. Mater.*, 521:167506, 2021.
[21] C. J. Garcia-Cervera. *Magnetic Domains and Magnetic Domain Walls*. PhD thesis, New York University, 1999.

[22] C. J. Garcia-Cervera. Structure of the Bloch wall in multilayers. *Proc. Roy. Soc. A*, 461:1911–1926, 2005.

[23] Carlos J. Garcia-Cervera. Néel walls in low anisotropy symmetric double layers. *SIAM J. Appl. Math.*, 65:1726–1747, 2005.

[24] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blugel. Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions. *Nature Phys.*, 7:713–718, 2011.

[25] A. Hrabec, J. Sampaio, M. Belmeguenai, I. Gross, R. Weil, S. M. Chérif, A. Stashkevich, V. Jacques, A. Thiaville, and S. Rohart. Current-induced skyrmion generation and dynamics in symmetric bilayers. *Nature Commun.*, 8:15765, 2017.

[26] A. Hubert. Statics and dynamics of domain walls in bubble materials. *J. Appl. Phys.*, 46:2276–2287, 1975.

[27] A. Hubert and R. Schäfer. *Magnetic Domains*. Springer, Berlin, 1998.

[28] B. Ivanov, V. Stephanovich, and A. Zhmudskii. Magnetic vortices: The microscopic analogs of magnetic bubbles. *J. Magn. Magn. Mater.*, 88:116–120, 1990.

[29] W. Jiang, G. Chen, K. Liu, J. Zang, S. G.E. te Velthuis, and A. Hoffmann. Skyrmions in magnetic multilayers. *Physics Reports*, 704:1–49, 2017.

[30] H. Knüpfer, C. B. Muratov, and F. Nolte. Magnetic domains in thin ferromagnetic films with strong perpendicular anisotropy. *Arch. Rat. Mech. Anal.*, 232:727–761, 2019.

[31] L. D. Landau and E. M. Lifshitz. *Course of Theoretical Physics*, volume 8. Pergamon Press, London, 1984.

[32] W. Legrand, J.-Y. Chauleau, D. Maccariello, N. Reyren, S. Collin, K. Bouzehouane, N. Jaouen, V. Cros, and A. Fert. Hybrid chiral domain walls and skyrmions in magnetic multilayers. *Sci. Adv.*, 4:eaat0415, 2018.

[33] W. Legrand, N. Ronceray, N. Reyren, D. MacCariello, V. Cros, and A. Fert. Modeling the shape of axisymmetric skyrmions in magnetic multilayers. *Phys. Rev. Appl.*, 10:064042, 2018.

[34] J. Leliaert, M. Dvornik, J. Mulkers, J. De Clercq, M. V. Milošević, and B. Van Waeyenberge. Fast micromagnetic simulations on GPU—recent advances made with Mumax3. *J. Phys. D: Appl. Phys.*, 51:123002, 2018.

[35] I. Lemesh and G. S. D. Beach. Twisted domain walls and skyrmions in perpendicularly magnetized multilayers. *Phys. Rev. B*, 98:104402, 2018.

[36] A. O. Leonov, T. L. Monchesky, N. Romming, A. Kubetzka, A. N. Bogdanov, and R. Wiesendanger. The properties of isolated chiral skyrmions in thin magnetic films. *New J. Phys.*, 18:065003, 2016.

[37] E. H. Lieb and M. Loss. *Analysis*. American Mathematical Society, Providence, RI, 2010.

[38] J. Lucassen, M. J. Meijer, O. Kurnosikov, H. J. M. Swagten, B. Koopmans, R. Lavrijsen, F. Kloodt-Twesten, R. Frömter, and R. A. Duine. Tuning magnetic chirality by dipolar interactions. *Phys. Rev. Lett.*, 123:157201, 2019.

[39] S. Luo and L. You. Skyrmion devices for memory and logic applications. *APL Materials*, 9:050901, 2021.
[40] S. A. Montoya, S. Couture, J. J. Chess, J. C. T. Lee, N. Kent, D. Henze, S. K. Sinha, M.-Y. Im, S. D. Kevan, P. Fischer, B. J. McMorran, V. Lomakin, S. Roy, and E. E. Fullerton. Tailoring magnetic energies to form dipole skyrmions and skyrmion lattices. *Phys. Rev. B*, 95:024415, 2017.

[41] C. Moreau-Luchaire, C. Moutafis, N. Reyren, J. Sampaio, C. A. F. Vaz, N. Van Horne, K. Bouzehouane, K. Garcia, C. Deranlot, P. Warnicke, P. Wohlhüter, J.-M. George, M. Weigand, J. Raabe, V. Cros, and A. Fert. Additive interfacial chiral interaction in multilayers for stabilization of small individual skyrmions at room temperature. *Nature Nanotechnol.*, 11:444–448, 2016.

[42] T. Moriya. Anisotropic superexchange interaction and weak ferromagnetism. *Phys. Rev.*, 120:91–98, 1960.

[43] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni. Skyrmion lattice in a chiral magnet. *Science*, 323:915–919, 2009.

[44] C. B. Muratov and V. V. Slastikov. Domain structure of ultrathin ferromagnetic elements in the presence of Dzyaloshinskii-Moriya interaction. *Proc. R. Soc. A*, 473:20160666, 2017.

[45] N. Nagaosa and Y. Tokura. Topological properties and dynamics of magnetic skyrmions. *Nature Nanotechnol.*, 8:899–911, 2013.

[46] E. Schlömann. Domain walls in bubble films. I. General theory of static properties. *J. Appl. Phys.*, 44:1837–1849, 1973.

[47] A. Soumyanarayanan, M. Raju, A. L. Gonzalez Oyarce, A. K. C. Tan, M.-Y. Im, A. P. Petrovic, P. Ho, K. H. Khoo, M. Tran, C. K. Gan, F. Ernult, and C. Panagopoulos. Tunable room temperature magnetic skyrmions in Ir/Fe/Co/Pt multilayers. *Nature Mat.*, 16:898–904, 2017.

[48] S. Woo, K. Litzius, B. Kruger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand, P. Agrawal, I. Lemesh, M.-A. Mawass, P. Fischer, M. Kläui, and G. S. D. Beach. Observation of room-temperature magnetic skyrmions and their current-driven dynamics in ultrathin metallic ferromagnets. *Nature Mat.*, 15:501–506, 2016.

[49] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura. Real-space observation of a two-dimensional skyrmion crystal. *Nature*, 465:901–904, 2010.