# Monte Carlo Studies on Geometrically Confined Skyrmions in Nanodots:
Stability and Morphology under Radial Stresses

G. Diguet<EMAIL_ADDRESS>Micro System Integration Center,
Tohoku University, Sendai, Japan  B. Ducharne  INSA Lyon, Universite de Lyon,
Villeurbanne Cedex, France  ELyTMaX, CNRS-Universite de Lyon-Tohoku
University, Sendai, Japan  S. El Hog  Universit$\acute{\rm e}$ de Monastir
(LMCN), Monastir, Tunisie  F. Kato  National Institute of Technology (KOSEN),
Ibaraki College, Hitachinaka, Japan.  H. Koibuchi<EMAIL_ADDRESS><EMAIL_ADDRESS>National Institute of Technology (KOSEN), Ibaraki College,
Hitachinaka, Japan.  T. Uchimoto  Institute of Fluid Science (IFS), Tohoku
University, Sendai, Japan  ELyTMaX, CNRS-Universite de Lyon-Tohoku University,
Sendai, Japan  H. T. Diep<EMAIL_ADDRESS>CY Cergy Paris University, Cergy-
Pontoise, France

###### Abstract

We numerically study the stability and morphology of geometrically confined
skyrmions in nanodots using Finsler geometry (FG) modeling technique. The FG
model dynamically implements anisotropies in ferromagnetic interaction,
Dzyaloshinskii-Moriya interaction, and magneto-elastic coupling in response to
mechanical stresses. Without the stresses, there exists a geometrically
confined effect originating from the surface effect of small nanodots, in
which skyrmions are stabilized under a low external magnetic field. This
surface effect is enhanced by radial stresses, which significantly reduce the
surface DMI compared to the bulk DMI. The radial stresses also alter the
interactions to be anisotropic. Owing to these position- and direction-
dependent interactions, incomplete skyrmions emerge at the center of the
nanodots under the tensile stress. In addition to the incomplete skyrmions,
target skyrmions are observed under the compressive stress. Our numerical
results indicate that the strain-enhanced surface effect and the strain-
induced interaction anisotropies suitably explain the skyrmion stability in
nanodots with zero magnetic field.

##  I Introduction

A topologically stable spin configuration called skyrmion emerges in chiral
magnetic materials under competition between ferromagnetic interaction (FMI)
and Dzyaloshinskii-Moriya interaction (DMI), and many studies have been
conducted Romming-etal-Science2013 ; Fert-etal-NatReview2017 ; Zhang-JPCM2020
; Tokura-Kanazawa-ACS2020 ; Gobel-etal-PhysRep2021 ; Wang-etal-JMMM2022 ; Li-
etal-Wiley2022  since its experimental discovery Uchida-etal-SCI2006 ;
Pfleiderer-etal-Science2009 ; Yu-etal-Nature2010 . It is well known that
skyrmions are stabilized by an external magnetic field in Zeeman energy and
mechanical strains from magnetoelastic coupling (MEC) and magnetic anisotropy
Bogdanov-PRL2001 ; Butenko-etal-PRB2010 . The confinement of skyrmions, called
geometric confinement (GC), has also been proposed as a stabilization
technique and emphasized to be caused by surface or edge boundary conditions
for spin configurations in nanodot Rohart-Thiaville-PRB2013 , and several
experimental and numerical studies on GC and strain effects have been
conducted HDu-etal-NatCom2015 ; Matsumoto-etal-NanoLett2018 ; CJin-etal-
NatCom2017 ; ZHou-etal-AcsNano2019 ; PHo-etal-PRAp2019 ; Seki-etal-PRB2017 ;
YWang-etal-NatCom2020 ; Li-etal-PRAP2020 .

As a GC effect inside nanowires, thermodynamically stable unusual skyrmions
were theoretically investigated Leonov-etal-EPJConf2014 . One is called
”target skyrmion” surrounded by multiple circular stripes. The target
skyrmions can be identified with ”helix” numerically investigated in the
presence of dipolar coupling Du-etal-PRBL2013  and were experimentally
observed in nanodisks of FeGe Zheng-etal-PRL2017  and Permalloy Kent-etal-
APL2019 . The other is of a vortex-like structure called ”incomplete
skyrmions” with the topological charge or skyrmion number $N_{\rm sk}$ less
than 1. The incomplete skyrmions were numerically studied for their responses
to the variations of external magnetic field and nanodot diameter in the
ground states Beg-etal-SRep2015 ; Beg-etal-PRB2017 .

Y. Wang et al. demonstrated the application of GC of skyrmions to electric
switching and confirmed a controllability of skyrmion on nanodots of diameter
350 nm without a magnetic field, in combination with GC and strain effects
YWang-etal-NatCom2020 . Their results show that the strain effects on magnetic
anisotropy and DMI are crucial for the morphological changes. The skyrmion
stabilization was numerically confirmed on nanodots without an external
magnetic field under a specific boundary condition Li-etal-PRAP2020 , and the
effects of magnetic anisotropy were studied numerically Mehmood-etal-JMMM2022
. A GC effect for stabilizing skyrmion vortex was also observed in tetrahedral
nanoparticles Niitsu-etal-Natmatt2022 . However, the origins of GC and the
effects of strain are unclear and need to be studied further.

In this study, we apply Finsler geometry (FG) modeling to clarify the strain
effect on morphological changes in spin configurations in nanodots using the
Monte Carlo (MC) simulation technique Diguet-etal-ICFD2023 . Dynamical
anisotropies in the FMI, DMI, and MEC emerging under radial stresses are
numerically evaluated to explain the skyrmion stability. For the
stabilization, FG modeling implements anisotropic effects of strains on
interactions by internal degrees of freedom, such as the strain direction,
which dynamically modifies the interactions to be direction-dependent and
position-dependent. For this simple geometric treatment of interactions, FG
modeling allows us to perform theoretical studies on skyrmions based on the
same input-output relation as experimental studies, leading to a better
understanding of the surface effect of skyrmions in nanodots. This is the main
advantage of FG modeling over standard modeling, in which interactions are
manually controlled by fixing the coefficients to be anisotropic.

This study is motivated by the experimental study in Ref. [20]. However, the
targeted material is simplified for modeling: for example a spontaneous in-
plane magnetic anisotropy assumed in YWang-etal-NatCom2020  is not necessary
in our model. Furthermore, our model is based on a single layer that may not
be as thin as experimental samples, and the experimental technique for the
application of radial strains to nanodots is not the same as in our model.
Nevertheless, we believe that these differences do not change the main effects
of the radial strains observed in YWang-etal-NatCom2020  and the present
numerical study.

##  II 3D cylindrical lattices and Hamiltonian

###  II.1 3D cylindrical lattices and radial stresses

Figure 1:  (a) A cylindrical lattice of radius $d$ and height $h$, (b)
illustrations of tensile ($f\\!>\\!0$) and compressive ($f\\!<\\!0$) stresses
radially applied to cylindrical lattice, (c) strain direction $\vec{\tau}_{1}$
at vertex 1 and its component $|\vec{\tau}_{1}\cdot\vec{e}_{12}|$, where
$\|\vec{e}_{12}\|\\!=\\!1$, along a local coordinate axis $x^{1}$ of a
tetrahedron with vertices 1, 2, 3 and 4. The variables ${\bf s},\vec{\tau}$
are defined at vertices and $\zeta$ on bonds, and these are introduced in the
following subsection. Free boundary condition is assumed for all surfaces.

Cylindrical lattices constructed by the Voronoi tessellation technique
Friedberg-Ren-NPB1984  are used to simulate skyrmions in nanodots YWang-etal-
NatCom2020  (Fig. 1(a)). The total number of vertices is $N\\!=\\!1533$, 2083,
5430, 8465, and 11932, and the ratio $R\\!=\\!d/h$ of diameter $d$ and a fixed
height $h$ is approximately given by $R\\!=\\!1.01$, $1.20$, $2.00$, $2.51$,
and $R\\!=\\!3.00$ (see Appendix A for further details of the lattice
structure). Mechanical stresses, $f\\!>\\!0$ and $f\\!<\\!0$ applied along the
radial direction (Fig. 1(b)), cause a variation of strain variable
$\vec{\tau}$ (Fig. 1(c)) leaving the lattice structure non-deformed, and a
magnetic field $B$ is applied along the $z$ direction. The direction of
magnetic field and that of stresses are the same as those in YWang-etal-
NatCom2020 . A free boundary condition is assumed for all surfaces, including
the lower disk, which contacts with a polymer substrate for strain transfer in
the reported experiment in YWang-etal-NatCom2020 ; hence, the boundary
condition in the present study is not exactly the same as that of the
experiment.

###  II.2 Hamiltonian and Monte Carlo

Figure 2:  (a) Illustration of bond $a$ or $ij$ and the connected bonds. The
total number of the connected bonds at vertex $i$ ($j$) is $q_{i}\\!-\\!1$
($q_{j}\\!-\\!1$), where $q_{i}$ ($q_{j}$) is the coordination number at
vertex $i$ ($j$). Illustrations of (b) $\sum_{ij(\Delta)}$; the sum over bonds
$ij(\Delta)$ of tetrahedron $\Delta$, and (c) $\sum_{\Delta(ij)}$; the sum
over tetrahedrons $\Delta(ij)$ sharing bond $ij$. $\sum_{ij(\Delta)}1\\!=\\!6$
for all $\Delta$, while the number $n_{ij}\\!=\\!\sum_{\Delta(ij)}1$ depends
on bond $ij$, and we have the identity $6N_{\rm
tet}\\!=\\!\sum_{\Delta}\sum_{ij(\Delta)}1\\!=\\!\sum_{ij}\sum_{\Delta(ij)}1\\!=\\!\sum_{ij}n_{ij}$,
where $N_{\rm tet}$ is the total number of tetrahedrons.

The discrete Hamiltonian is

$\displaystyle H({\bf s},\vec{\tau},\zeta)=\lambda H_{\rm FM}+DH_{\rm
DM}-H_{B}-\alpha H_{\rm ME}-H_{f}-\delta H_{\zeta},$ (1)

where ${\bf s}(\in S^{2}:{\rm unit\;sphere})$ and $\vec{\tau}(\in S^{2}/2:{\rm
half\;sphere})$ denote the spin and strain variables, respectively, defined at
every lattice vertex, and $\zeta(\in\\{1,-1\\})$ on the bonds (Fig. 1(c));
$\zeta$ is a variable to implement a dynamical anisotropy in MEC. The symbols
$\lambda,D,\alpha,\delta$ on the right-hand side are the interaction
coefficients; $\lambda$, $D$ and $\alpha$ are strengths of FMI, DMI and MEC,
respectively, and $\delta$ denotes a strength of $\zeta$ interaction for
Finsler length of MEC (Appendix B).

The terms on the right hand side are given as follows:

$\displaystyle\begin{split}&H_{\rm
FM}=\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})\left(1-{\bf
s}_{i}\cdot{\bf s}_{j}\right),\quad H_{\rm
DM}=\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j},\\\ &H_{B}=\sum_{i}{\bf
s}_{i}\cdot\vec{B},\;\vec{B}=(0,0,B),\quad H_{\rm
ME}=f\sum_{\Delta}\sum_{ij(\Delta)}\Omega_{ij}(\vec{\tau},\zeta)\left({\bf
s}_{i}\cdot{\bf s}_{j}\right)^{2},\\\ &H_{f}={\rm
sgn}(f)\sum_{i}\left(\vec{\tau}_{i}\cdot\vec{f}\right)^{2},\quad\vec{f}={f}\vec{e}^{\;r},\quad{\rm
sgn}(f)=\left\\{\begin{array}[]{@{\,}ll}1\quad(f>0)\\\
-1\quad(f<0)\end{array}\right.,\\\
&H_{\zeta}=\sum_{(ab)}\zeta_{a}\zeta_{b}\end{split}$ (2)

The first term $H_{\rm FM}$ is an FMI energy that is modified to have a
$\vec{\tau}$-dependent interaction coefficient $\Gamma_{ij}(\vec{\tau})$,
which depends on a small number $v_{0}$ playing a role in the strength of
anisotropy (see Appendix B for FG modeling details). The same
$\Gamma_{ij}(\vec{\tau})$ is included in the second term $H_{\rm DM}$ for DMI
energy. The third term $H_{B}$ is the Zeeman energy, and the fourth term
$H_{\rm ME}$ is the energy for the ME coupling quadratic with respect to ${\bf
s}$ with a coefficient $\Omega_{ij}(\vec{\tau},\zeta)$, which is dynamically
variable depending on $\zeta$ and, in this sense, different from
$\Gamma_{ij}(\vec{\tau})$, which varies depending only on $\vec{\tau}$, in
$H_{\rm FM}$ and $H_{\rm DM}$ (Appendix B). The fifth term $H_{f}$ denotes the
energy required for the response of strain $\vec{\tau}$ to an external stress
$\vec{f}$ along the radial direction $\vec{e}^{\;r}$. The final term
$H_{\zeta}$, which is not included in the model in Ref. Diguet-etal-ICFD2023 ,
describes the dynamical Finsler length for an anisotropic MEC in $H_{\rm ME}$.
$H_{\zeta}$ is defined as the sum of connected bonds $a$ and $b$ denoted by
$(ab)$ (Fig. 2(a)). We note that $H_{\zeta}$ is equivalent to the Potts model
Hamiltonian, $H_{\rm P}\\!=\\!\sum_{ab}\delta_{\zeta_{a},\zeta_{b}}$, where
$\delta_{\zeta_{a},\zeta_{b}}\\!=\\!1\;(\zeta_{a}\\!=\\!\zeta_{b})$ and
$\delta_{\zeta_{a},\zeta_{b}}\\!=\\!0\;(\zeta_{a}\\!\not=\\!\zeta_{b})$.
$H_{\rm FM}$, $H_{\rm DM}$ and $H_{\rm ME}$ are commonly defined by
$\sum_{\Delta}\sum_{ij(\Delta)}$, the sum over tetrahedrons $\Delta$, where
$\sum_{ij(\Delta)}$ denotes the sum over bonds $ij(\Delta)$ of tetrahedron
$\Delta$ (Fig. 2(b)).

In the case of a zero stress $f\\!=\\!0$, the spin configurations are
determined only by the first three terms. For this reason, no spontaneous
magnetic anisotropy is assumed in the model; therefore, the material studied
in the present paper is not exactly the same as that experimentally studied in
Ref. YWang-etal-NatCom2020  in addition to a difference in the application
technique of radial stress on nanodots. Moreover, the materials studied in
Ref. YWang-etal-NatCom2020  are layered ones and a Neel-type skyrmion is
assumed in their simulations; in contrast, the materials in the present paper
are single-layer and not always thin, in which a Bloch-type skyrmion is
assumed.

The partition function to be simulated by MC technique is given by

$\displaystyle Z=\sum_{{\bf s}}\sum_{\vec{\tau}}\sum_{\zeta}\exp\left[-H({\bf
s},\vec{\tau},\zeta)/T\right],$ (3)

where $T$ is the temperature and $\sum_{s,\tau,\zeta}$ denotes the sum of all
possible configurations of ${\bf s},\vec{\tau}$ and $\zeta$. The partition
cannot be numerically obtained, and instead, the sum over configurations,
defined by $\sum_{{\bf
s}}\sum_{\vec{\tau}}\\!=\\!(\int\Pi_{i=1}^{N}d\vec{s}_{i})(\int\Pi_{i=1}^{N}d\vec{\tau}_{i})$
the multiple integrals on the unit sphere and
$\sum_{\zeta}\\!=\\!\Pi_{i=1}^{N}\sum_{\zeta_{i}\in\\{-1,1\\}}$, is simulated
by Metropolis importance MC sampling Metropolis-JCP-1953 ; Landau-PRB1976 ,
where $\vec{\tau}_{i}$ is identified with $-\vec{\tau}_{i}$ in the
Hamiltonians. The mean value of physical quantity $P_{q}$ is defined by the
ensemble average $\langle
P_{q}\rangle\\!=\\!\sum_{s,\tau,\zeta}P_{q}\exp(-H/T)/Z$, which is calculated
using Metropolis MC technique Metropolis-JCP-1953 ; Landau-PRB1976  with the
sample average $\langle P_{q}\rangle\\!=\\!\sum_{i}P_{q}(i)/\sum_{i}1$, where
$P_{q}(i)$ is $i$-th sample evaluated by the lattice average. A new variable
${\bf s}^{\prime}$ is randomly generated, independent of the present ${\bf s}$
and updated with probability ${\rm Min}[1,\exp-(H_{\rm new}\\!-\\!H_{\rm
old})]$. The same procedure is applied to the variable $\vec{\tau}$. When
updating $\zeta(\in\\!\\{-1,1\\})$, we fix $\zeta^{\prime}$ to the opposite to
the present $\zeta$ and accept $\zeta^{\prime}$ with the same procedure. $N$
consecutive updates for each of the three variables are examined, and a
uniform random number is used to select $N$ from $N_{B}$ $\zeta$-variables
defined on the bonds, where $N_{B}$ is the total number of bonds (Appendix A).

The total number of MC sweeps (MCS) for the calculation of physical quantities
is $5.25\\!\times\\!10^{7}$ for the lattices of $N\\!\leq\\!8465$ and
$(1.75\sim 3.5)\\!\times\\!10^{8}$ for the $N\\!=\\!11932$ lattice. The
samples are calculated for every 500 MCS after the thermalization iterations
of $7.5\\!\times\\!10^{6}$ $\sim$ $12.5\\!\times\\!10^{6}$ MCS for the
lattices of $N\\!\leq\\!8465$ and $5\\!\times\\!10^{7}$ MCS for the
$N\\!=\\!11932$ lattice starting with random initial configurations.

##  III Results

###  III.1 Size effect without strains

Figure 3:  Phase diagram of $B$ vs. $d$ for confined skyrmion ($\bullet$) and
non-confined or stripe ($\circ$) observed on four different lattices, where
$d$ is the diameter of lattice in the unit of the edge length ($a\\!=\\!1$) of
a regular triangle (Appendix A). Spins of $s^{z}\\!<\\!0$ are plotted with
small cones in the snapshots. Small cylinders or lines, shown only on the
surface for clear visualization, denote $\vec{\tau}$, which is isotropic and
does not influence the spins.

First, we demonstrate a size effect that enables skyrmions to emerge on
smaller nanodots with a small external magnetic field. The assumed parameters
are shown in Fig. 3. The temperature is fixed to $T\\!=\\!0.5$ in all
simulations, because it is low enough to see skyrmions.
$v_{0}\\!=\\!(0.1,0.1)$ denotes that $v_{0}\\!=\\!0.1$ for both the FMI and
DMI in Eq. (15). $v_{0}$ plays a role in the strength of anisotropy in the
cases that strains $\vec{\tau}$ are uniformly aligned; however, it does not
play any role in anisotropy of interactions in the case of $f\\!=\\!0$, where
$\vec{\tau}$ is isotropic. Plots in Fig. 3 represent configurations of
confined skyrmion ($\bullet$) and non-confined or stripe ($\circ$) observed on
four different lattices. The color code represents $s^{z}$ and the spins of
$s^{z}\\!<\\!0$ are plotted with cones, and small cylinders or lines represent
the direction of $\vec{\tau}$.

All the plotted data are obtained by visually evaluating whether the skyrmions
are confined using snapshots. The ”non-confined” means (i) stripe
configurations or (ii) skyrmions are located at the boundary, and the
”confined” means (iii) skyrmions are inside the boundary. However, these
criteria are not strict, and the phase diagram shown in Fig. 3 is a rough
estimate. The dashed line in red (blue) shows the upper (lower) limit of $B$
for the non-confined (confined) phase. The region between these two dashed
lines is similar to the range of two-phase coexistence; however, it does not
always imply a first-order phase transition because no clear jump is observed
in the physical quantities, including the topological charge between the
confined skyrmion and non-confined skyrmion phases. Interestingly, the
decrease of $B$ for the confined region with decreasing $d$ is qualitatively
consistent with experimental data reported in Ref. YWang-etal-NatCom2020 .

Notably, the confinement mechanism in the model originates from a surface
effect, in which the absolute value of the surface DMI coefficient is smaller
than that of the bulk DMI Diguet-etal-JMMM2023 . For this reason, skyrmions
hardly appear on the periphery region. A static part of this difference in the
DMI strength for $f\\!=\\!0$ is intuitively understood from the discrete form
of $H_{\rm DM}$ in Eq. (2) (Appendix B), and a dynamical enhancement of the
difference for $f\\!>\\!0$ will be shown below.

###  III.2 Strain effect

Figure 4:  The snapshots obtained for (a) $f\\!=\\!-4.8$ (compression), (b)
$f\\!=\\!0$ , and (c) $f\\!=\\!2.28$ (tension) on a small lattice of size
$N\\!=\\!2803$ for clear visualization. We find the target skyrmions (target
Sk), stripes, and incomplete skyrmion (iSk) phases for negative, zero, and
positive stresses. These three phases correspond to vortex, stripe and
skyrmion in the thin-layered case of the reported experimental data YWang-
etal-NatCom2020 . Small cylinders for $\vec{\tau}$ are shown only on the
surface for clear visualization. The parameters $T,\lambda,D,\alpha,\delta$
and $v_{0}$ shown in the figure are used in all simulations below.

Next, we show the effects of the strain caused by radial tension ($f>0$) and
compression ($f<0$) under zero external magnetic field ($B\\!=\\!0$). To
visualize the configurations of ${\bf s}$ and $\vec{\tau}$, a small cylinder
of size $N\\!=\\!2803$ is used in this subsection. The assumed parameters are
shown in Fig. 4; $v_{0}\\!=\\!(0.1,0.1,0.2)$ denotes that $v_{0}\\!=\\!0.1$
for both FMI and DMI in Eq. (15) and $v_{0}\\!=\\!0.2$ for MEC in Eq. (18),
which are fixed in this study.

Snapshots are shown in Figs. 4(a)–(c), where top and side views are plotted in
(a) $f\\!=\\!-4.8$, (b) $f\\!=\\!0$, and (c) $f\\!=\\!2.28$. Small cylinders
representing $\vec{\tau}$ can be viewed as (a) a spiral along the $z$
direction, (b) isotropic, and (c) radially aligned. These observations are
consistent with expectations because $\vec{\tau}$ represents the elongation
direction owing to the radial stress $f$.

The configuration shown in Fig. 4(a) is a skyrmion enclosed by circular
stripes corresponding to the target skyrmion denoted by ”target Sk” or simply
by ”target” Leonov-etal-EPJConf2014 , which is not exactly the same as
”vortex” observed in Ref. YWang-etal-NatCom2020  for the thin-layered case.
The configuration in Fig. 4(c) is also a skyrmion, which is called ”incomplete
skyrmion” in Leonov-etal-EPJConf2014  denoted by iSK in the figure. The
direction of ${\bf s}$ at the periphery is not completely opposite to that at
the center Leonov-etal-EPJConf2014 ; Beg-etal-SRep2015 ; Beg-etal-PRB2017 .
However, in their experiment on the thick-layered case, stripes were observed,
which are consistent with our numerical data. We checked that these
morphological changes under radial stress can also be observed in Neel-type
skyrmions. Thus, we consider that these morphological changes under radial
stresses correspond to vortices, stripes, and skyrmions in the thin-layered
case of reported experimental data YWang-etal-NatCom2020 , although the vortex
and skyrmions in YWang-etal-NatCom2020  are not exactly corresponding to the
target and incomplete skyrmions in this paper, respectively.

It should be noted that the spin direction at the center of the target
skyrmion shown in Fig. 4(a) is spontaneously determined and opposite to that
shown in Fig. 3. This spontaneous direction occurs because $B\\!=\\!0$, in
contrast with the case in Ref. YWang-etal-NatCom2020 , where a small magnetic
field is applied for atomic microscopy measurements, even with zero external
magnetic field. For the same reason, the polarity at the center of the
incomplete skyrmion in Fig. 4(c) is also spontaneously determined. Note also
that the radial stress, for example, compression, is well defined only in
3-dimensional models because of the induced deformation of the nanodots along
the $z$ direction indicated by the $\vec{\tau}$ direction.

Figure 5:  The snapshots of target skyrmions on the lattices of size (a)
$N\\!=\\!1533$, (b) $N\\!=\\!5430$ and (c) $N\\!=\\!11932$, and (d) incomplete
skrymion on the $N\\!=\\!1533$ lattice, (e) an isolated skrymion and (f)
incomplete skrymion on the $N\\!=\\!11932$ lattice. The polarity variation
between the center and periphery is less than $\pi$ in (d), (f), whereas it is
more than $\pi$ in (e).

To observe the lattice size dependence, snapshots of the target, incomplete,
and isolated skyrmions are shown in Figs. 5(a)–(f) for lattices of size
$N\\!=\\!1533$, $N\\!=\\!5430$ and $N\\!=\\!11932$. Note that the total number
of circular stripes of $N\\!=\\!5430$ and $N\\!=\\!11932$ lattices is larger
than that of the $N\\!=\\!1533$ lattice because the lattice diameter $d$
approximately increases $1.97$ and $2.97$ times, respectively (Table 1 in
Appendix A). Such a change in the number of morphology can also be observed in
the stripe phase. Therefore, the number of circular stripes is determined by
the scale $\lambda/D$ of the FMI and DMI coefficients in Eq. (1).

However, the total number (=1) of incomplete skyrmions remains unchanged and
is independent of $d$ on the lattices of $N\\!=\\!1533$ in Fig. 5(d) and
$N\\!=\\!5430$ (not shown in Fig. 5), in sharp contrast to the case of nonzero
magnetic field with $f\\!=\\!0$ in Fig. 3; the incomplete skyrmion looks like
a vortex Shinjo-etal-science2000 ; Okuno-etal-JMMM2002 . Nevertheless, when
the diameter $d$ increases further, an isolated skyrmion, denoted by
”isolated”, emerges in the positive region of $f$ (Fig. 5(e)), in which the
polarity between the center and periphery changes more than $\pi$ in contrast
to the case of incomplete skyrmion. When $f$ is further increased to
$f\\!=\\!2.56$, incomplete skyrmions appear even on the $N\\!=\\!11932$
lattice (Fig. 5(f)).

The emergence of incomplete skyrmions is expected to be owing to the strong
MEC effect for a relatively large $f(>0)$, where a strain configuration
($\Leftrightarrow\\!\vec{\tau}\\!\propto\\!\vec{e}^{\;r}\\!=\\!\vec{r}/r$)
shares the center with the skyrmion for the rotational symmetry independent of
$N$. This rotationally symmetric configuration of $\vec{\tau}$ enhances the
surface effect of GC. The $\vec{\tau}$ direction is almost perpendicular to
the bond directions, and the FMI and DMI are expected to be very small on the
side surface owing to the FG modeling prescription. Therefore, there is no
chance for incomplete skyrmions to appear independently of the primary
incomplete skyrmions on the nanodots, which emerge as a ground state without
strain under zero magnetic field Beg-etal-SRep2015 ; Beg-etal-PRB2017 . This
problem will be further discussed below.

###  III.3 Order parameters and energies

Figure 6:  (a) Nonpolar order parameters $Q_{s}^{\mu}(\mu\\!=\\!r,\theta,z)$
and $Q_{s}^{\rm in}(=\\!Q_{s}^{r}\\!+\\!Q_{s}^{\theta})$ vs. $f$, where the
indices represent the directions $\mu$ and spin $s$, and ”in” denotes ”in-
plane”. (b)–(i) Snapshots for $-5.6\leq f\leq 2.64$ obtained on the
$N\\!=\\!5430$ lattice under $B\\!=\\!0$, and (j) FMI, DMI, MEC, and stress-
strain energies in Eq. (2) vs. $f$. We find three different phases; target
skyrmion (target Sk) and stripe, and incomplete skyrmion (iSk), and the stripe
phase can be divided into two phases; rounded stripe (”r-stripe”) and stripe.
We find target skyrmions in (f) and (g), elongated incomplete skyrmion in (h)
and rotationally symmetric incomplete skyrmion in (i), where ${s}^{z}$ at the
periphery is $s_{z}\\!>\\!-1$. Note that $-\alpha H_{\rm ME}$ and $-H_{f}$ are
included in the total Hamiltonian of Eq. (1).

In Fig. 6(a), we plot nonpolar order parameters of spins defined by

$\displaystyle\begin{split}Q_{s}^{\mu}=&\frac{3}{2}\left(\langle({\bf
s}\cdot\vec{e}_{i}^{\;\mu})^{2}\rangle-\frac{1}{3}\right),\\\ &(-0.5\leq
Q_{s}^{\mu}\leq 1),\quad(\mu=r,\theta,z),\\\ Q_{s}^{\rm
in}=&Q_{s}^{r}+Q_{s}^{\theta}\end{split}$ (4)

where $\vec{e}_{i}^{\;r}$ and $\vec{e}_{i}^{\;\theta}$ denote unit vectors
along the radial and tangential directions at vertex $i$, respectively
(Appendix C). Order parameters of this type are always used for liquid crystal
molecules, which are nonpolar Andrienko-JML2018 . The reason why nonpolar
order parameters are calculated for ${\bf s}$ is that $Q_{s}^{\mu}$ are useful
to find the phase boundaries as shown below.

The parameters $Q_{s}^{r}$ ($\circ$), $Q_{s}^{\theta}$ ($\triangle$) and
$Q_{s}^{\rm in}$ ($\bullet$) discontinuously change at $f\\!\simeq\\!1$
indicating that the incomplete skyrmion phase is discontinuously separated
from the stripe phase. The jump of $Q_{s}^{\rm in}$ at the phase boundary is
owing to those in $Q_{s}^{r,\theta}$. The jump in $Q_{s}^{\theta}$ implies an
anisotropic alignment of ${\bf s}_{i}$ in the $\pm{\vec{e}}^{\;\theta}$
directions, and this in-plane magnetic anisotropy as a response to $f(>\\!0)$
is a non-trivial effect of FMI, DMI, and MEC implemented in $H_{\rm
FM,DM,ME}$. We find an intermediate phase denoted by ”r-stripe” (rounded
stripe) between the target skyrmion and stripe phases has continuous but cusp-
like changes in the variation of $Q_{s}^{r}$ ($\circ$) and $Q_{s}^{\theta}$
($\triangle$) vs. $f$. The incomplete skyrmion and stripe phases are also
characterized by $Q_{s}^{\rm in}$ ($\bullet$), although the difference in
$Q_{s}^{\rm in}$ between the two phases is not very large. Note that behaviors
of $Q_{s}^{\mu}$ depend on the lattice size $N$; the transition point between
the incomplete skyrmion and stripe phases depends on $N$. Here, we only
describe the results obtained on the $N\\!=\\!5430$ lattice, which is not
sufficiently large. However, this lattice size is reasonable for studying the
surface effect.

Here, we discuss the differences in $Q_{s}^{r}$ between the target skyrmion,
r-stripe, stripe, and incomplete skyrmion phases. Note that
$Q_{s}^{r}\\!\to\\!-0.5$ if all ${\bf s}$ satisfy ${\bf
s}\\!\perp\\!\vec{e}^{\;r}$, $Q_{s}^{r}\\!\to\\!0$ for random ${\bf s}$ with
respect to the direction $\vec{e}^{\;r}$, and $Q_{s}^{r}\\!\to\\!1$ for ${\bf
s}\\!\parallel\\!\vec{e}^{\;r}$. These properties of $Q_{s}^{\mu}$ describe
directional- or anisotropic-property of magnetization independent of
$\pm$-direction. In the target skyrmion phase (Figs. 6(f),(g)), $Q_{s}^{r}$
remains monotonically decreasing as $f$ increases but almost constant;
$Q_{s}^{r}\\!\simeq\\!-0.2$, indicating that ${\bf s}\\!\perp\\!\vec{e}^{\;r}$
for many ${\bf s}$. In the r-stripe phase (Fig. 6(b)), $Q_{s}^{r}$ starts to
increase implying that the condition ${\bf s}\\!\perp\\!\vec{e}^{\;r}$ is not
always satisfied as $f$ increases. This increase in $Q_{s}^{r}$ terminates and
decreases again as $f$ moves to the right in the stripe phase (Figs. 6(c),(d)
and (e)), and it discontinuously decreases at $f\\!\simeq\\!1$ illuminating
the phase boundaries.

In the incomplete skyrmion phase (Figs. 6(h),(i)), we find
$Q_{s}^{r}\\!\to\\!-0.5$ for large $f(>\\!0)$ corresponding to ${\bf
s}\\!\perp\\!\vec{e}^{\;r}$. This ${\bf s}\\!\perp\\!\vec{e}^{\;r}$ implies
that ${\bf s}$ aligns along the $z$-direction or $\theta$-direction and is
reasonably understood from the incomplete skyrmion configuration. The
alignment in the $\theta$-direction can be confirmed in $Q_{s}^{\theta}$ in
Fig. 6(a), and the alignment along the $z$-direction will be shown below.
Thus, we find a strain-induced magnetic anisotropy along the
$\theta$-direction and $z$-direction in the incomplete skyrmion phase for
increasing tensile stress.

The sign of $s^{z}$ at the periphery of the skyrmions that emerge under
tensile stress is not always completely opposite to that at the center, as
shown in Figs. 6(h),(i). Therefore, these skyrmions resemble vortex
configurations, as mentioned above Shinjo-etal-science2000 ; Okuno-etal-
JMMM2002 . However, the opposite sign of $s^{z}$ at the periphery is apparent
when $f$ is small; that is, $f\\!=\\!1.68$ in Fig. 6(i). Therefore, the
configurations shown in Figs. 6(h),(i) are considered to be incomplete
skyrmions.

FMI, DMI and MEC energies per tetrahedron plotted in Fig. 6(j), where $H_{\rm
FM}$, $H_{\rm DM}$, $H_{\rm ME}$ and $H_{f}$ are given in Eq. (2) and $N_{\rm
tet}$ is the total number of tetrahedrons. We find that $H_{\rm FM}/N_{\rm
tet}$ ($H_{\rm DM}/N_{\rm tet}$) in the incomplete skyrmion phase is smaller
(negatively larger) than that in the other phases, implying that the
incomplete skyrmion phase is energetically stable for the FMI and DMI,
although $H_{\rm DM}/N_{\rm tet}$ discontinuously increases when $f$ increases
at $f\\!\simeq\\!1$. $H_{\rm ME}/N_{\rm tet}$ and $H_{f}/N$ increase almost
linearly with increasing $f$ for $f>0$ as expected from the definitions of
$H_{\rm ME}$ and $H_{f}$ in Eq. (2). These results also imply that $H_{\rm
ME}$ and $H_{f}$ are small or negatively large in the incomplete skyrmion
phase because $-\alpha H_{\rm ME}$ and $-H_{f}$ are included in the total
Hamiltonian of Eq. (1). The shape of $H_{f}/N_{\rm tet}$ is almost constant
zero for $f\\!<\\!0$ implying that
$(\vec{\tau}_{i}\cdot\vec{f})^{2}\\!\to\\!0$ for negatively large $f$. This
result is consistent with the observation that
$\vec{\tau}_{i}\perp\vec{e}^{\;r}$ in the target skyrmion phase (Figs. 4(a)
and 5(a)-(c)).

To see effects of strains on magnetization, we plot the absolute of polar
order parameters

$\displaystyle M^{z}=\frac{1}{N}|\sum_{i}s_{i}^{z}|,\quad
M^{xy}=\frac{1}{N}\sqrt{(\sum_{i}s_{i}^{x})^{2}+(\sum_{i}s_{i}^{y})^{2}}$ (5)

in Fig. 7(a), where the dashed vertical lines represent the phase boundaries
plotted in Fig. 6(a). In the incomplete skyrmion phase, we find that $M^{z}$
increases from $M^{z}\\!\simeq\\!0$ to $M^{z}\\!>\\!0$ as $f$ increases. The
small increment of $M^{z}$ with increasing $f$ is considered a signal of ”the
strain-induced alignment of ${\bf s}$ along the $\pm z$ direction”. This
result is consistent with ”the strain-mediated increase of in-plane magnetic
anisotropy” observed under increasing compression in Ref. YWang-etal-
NatCom2020 , because we can read the plotted $M^{z}$ such that $M^{z}$
decreases with decreasing $f$ in the incomplete skyrmion phase. The result
$M^{xy}\\!\simeq\\!0$ does not imply that the tensile stress $f$ decreases the
in-plane magnetization but is simply owing to the cancellation in $\pm
s_{i}^{x}$ and $\pm s_{i}^{y}$. This is the reason why the non-polar order
parameters $Q_{s}^{\mu}$ plotted in Fig. 6(a) are useful to observe the
$\pm$-direction independent alignment. The plotted data $M^{z}\\!\to\\!0$ in
the stripe and target skyrmion phases are also reasonable because ${\bf s}$
that align along the $\pm z$ direction cancel each other.

Figure 7:  (a) Polar order parameter $M^{z}$ and $M^{xy}$ vs. $f$, and (b)
order parameters $Q_{\zeta}$ and $Q_{\tau}^{\,\mu}(\mu=r,\theta,z)$ vs. $f$.

An order parameter $Q_{\zeta}$ for $\zeta$ defined by

$\displaystyle Q_{\zeta}=(1/N_{B})\sum_{a}\zeta_{a},$ (6)

is plotted in Fig. 7(b). The value $Q_{\zeta}\\!=\\!1$ ($Q_{\zeta}\\!=\\!-1$)
corresponds to the cosine (sine) type $v_{ij}$ in Eq. (18). An abrupt change
in $Q_{\zeta}$ at $f\\!=\\!0$ does not always indicate an internal phase
transition caused by external mechanical stress $f$, because no abrupt change
is observed in energies including $H_{\zeta}$. However, this change in
$Q_{\zeta}$ causes a discontinuity in the effective coupling constant $\Omega$
as shown below. Notably, we verified that no target skyrmion configuration
would emerge if
$Q_{\zeta}\\!=\\!1(\Leftrightarrow\\!v_{ij}\\!=\\!|\cos\theta_{ij}|\\!+\\!v_{0})$
for $f\\!<\\!0$ by manually fixing $\zeta$ to be $\zeta\\!=\\!1$. Therefore,
this transition in $Q_{\zeta}$ (Fig 7(b)) is expected to be crucial for the
morphological change on nanodots under radial stress. From the direction-
dependent nonpolar order parameters
$Q_{\tau}^{\mu}\\!=\\!\frac{3}{2}\left(\langle(\vec{\tau}\cdot\vec{e}_{i}^{\;\mu})^{2}\rangle-\frac{1}{3}\right),\;(\mu=r,\theta,z)$
in Fig. 7(b), we find that the strain $\vec{\tau}$ is at random
($\Leftrightarrow\\!Q_{\tau}^{\mu}\to\\!0,(\mu\\!=\\!r,\theta,z)$) for
$f\\!\to\\!0$ and aligns such that $\tau^{r}\\!\to\\!0$
($\Leftrightarrow\\!Q_{\tau}^{r}\to\\!-0.5$) [($\tau^{\theta,z}\\!\to\\!0$
($\Leftrightarrow\\!Q_{\tau}^{\theta,z}\to\\!-0.5$)] in the target
[incomplete] skyrmion phase, as expected.

###  III.4 Direction-dependent interaction coefficients with effective
Hamiltonian and topological charge

Figure 8:  Position-dependent energies. Bulk and side energies per bond for
(a) FMI, (b) DMI, and (c) MEC. The coefficients
$(\lambda,D,\alpha)\\!=\\!(1,0.7,0.3)$ are not included in these energies. We
find a surface effect that ${H}_{\rm FM}({\rm side})/N_{B}\\!\to\\!0$,
${H}_{\rm DM}({\rm side})/N_{B}\\!\to\\!0$ and ${H}_{\rm ME}({\rm
side})/N_{B}\\!\to\\!0$ in the incomplete skyrmion (iSk) phase. We also find a
competing nature of FMI and DMI in the iSk phase, such that ${H}_{\rm FM}({\rm
bulk})/N_{B}$ decreases, whereas ${H}_{\rm DM}({\rm bulk})/N_{B}$ increases
with increasing $f$. Note that $-\alpha H_{\rm ME}$ is included in the total
Hamiltonian of Eq. (1).

First, we plot the FMI and DMI Hamiltonians on the side surface and inside or
bulk to show the surface effect. From the identity
$\sum_{\Delta}\sum_{ij(\Delta)}\\!=\\!\sum_{ij}\sum_{\Delta(ij)}$ for the
summation convention, energies per tetrahedron can also be written as energies
per bond such that

$\displaystyle\begin{split}&H_{\rm
FM}/N_{B}=\frac{1}{N_{B}}\dot{\sum_{ij}}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\left(1-{\bf
s}_{i}\cdot{\bf s}_{j}\right),\\\ &H_{\rm
DM}/N_{B}=\frac{1}{N_{B}}\dot{\sum_{ij}}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j},\\\ &H_{\rm
ME}/N_{B}=\frac{f}{N_{B}}\dot{\sum_{ij}}\sum_{\Delta(ij)}\Omega_{ij}(\vec{\tau},\zeta)\left({\bf
s}_{i}\cdot{\bf s}_{j}\right)^{2},\end{split}$ (7)

where $\dot{\sum}_{ij}$ denotes that the sum over bond $ij$ is limited to a
part of the lattice, and $N_{B}(=\\!\dot{\sum}_{ij}1)$ is the total number of
bonds included in that part. The expression of energy per bond is more
convenient for observing the difference in energies between the surface and
bulk, because it is apparent whether a bond $ij$ is on the surface or not.
Note that, when $H_{\rm DM}/N_{B}$ is limited to the side surface, $N_{B}$ is
the total number of bonds on the side surface. The circular edges shared by
the side surface and disks are one-dimensional objects; however, the edges are
included on the side surface, for simplicity.

The data obtained on the $N\\!=\\!5430$ lattice are plotted in Figs. 8(a)–(c).
We find a large difference between ${H}_{\rm DM}({\rm side})/N_{B}$ and
${H}_{\rm DM}({\rm bulk})/N_{B}$. One reason for this difference is that (i)
there exists a sum over tetrahedrons $\dot{\sum}_{\Delta(ij)}$ inside the sum
over bonds $\dot{\sum}_{ij}$ such that $\dot{\sum}_{ij}\sum_{\Delta(ij)}$ in
$H_{\rm DM}$, and this ${\sum}_{\Delta(ij)}$ depends on whether bond $ij$ is
on the surface or inside; the large difference is not obtained if we use the
sum over bonds $H_{\rm DM}\\!=\\!\dot{\sum}_{ij}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j}$. The other reason is a dynamical one that (ii) the
intensive part of energy $\Gamma_{ij}(\vec{\tau})$ included in $H_{\rm DM}$
depends on position and direction. This second part (ii) is a main topic in
this paper.

Figure 9:  (a) Position-dependent $\Gamma({\rm bulk,side,disk})$ in Eq. (9);
$\Gamma$ ($\bullet$) is calculated using the whole bonds. (b) Position-
dependent $\Omega({\rm bulk,side,disk})$ for ME. (c) Direction-dependent
coefficients $\Gamma^{\mu}({\rm bulk})$, and (d) the corresponding bulk
Hamiltonian $h_{\rm DM}^{\mu}$ in Eq. (10) and (e) $h_{\rm DM}^{\mu}({\rm
side})$ on the side surface. (f) The bulk coefficient of the magneto-elastic
coupling $\Omega^{\mu}$, and (g) the corresponding bulk Hamiltonian $h_{\rm
ME}^{\mu}({\rm bulk})$, (h) FMI bulk Hamiltonian $h_{\rm FM}^{\mu}({\rm
bulk})$ and side Hamiltonian $h_{\rm FM}^{\mu}({\rm side})$.

To see this second reason (ii), we plot the effective DMI (and FMI) coupling
constant defined by

$\displaystyle\Gamma=\frac{1}{N_{B}}\dot{\sum}_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})$
(8)

and $\Omega$ defined by the same expression in Figs. 9(a),(b). In this lattice
average, the direction dependence of $\Gamma_{ij}(\vec{\tau})$ is integrated
out or summed over. To extract the direction dependence of
$\Gamma_{ij}(\vec{\tau})$, we define the direction-dependent coefficients

$\displaystyle\begin{split}&\Gamma^{\mu}=\frac{1}{N_{B}}\dot{\sum}_{ij}\sum_{\Delta(ij)}|\Gamma_{ij}^{\mu}(\vec{\tau})|,\;(\mu=r,\theta,z)\end{split}$
(9)

which are plotted in Fig. 9(c), where $\Gamma_{ij}^{\mu}(\vec{\tau})$ is given
by Eq. (28). The corresponding direction-dependent energies $h_{\rm DM}^{\mu}$
plotted in Fig. 9(d) are defined by

$\displaystyle\begin{split}&h_{\rm DM}^{\mu}=H_{\rm DM}^{\mu}/N_{B},\\\
&H_{\rm
DM}^{\mu}=\frac{\dot{\sum}_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;\mu})^{2}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf
s}_{j}}{\frac{1}{N_{B}}\dot{\sum}_{ij}\sum_{\Delta(ij)}|\Gamma_{ij}^{\mu}(\vec{\tau})|},\;(\mu=r,\theta,z)\end{split}$
(10)

where $\sum_{\Delta}\sum_{ij(\Delta)}$ is replaced by
$\dot{\sum}_{ij}\sum_{\Delta(ij)}$ in $H_{\rm DM}^{\mu}$ of Eq. (30) and
$N_{B}$ is given by $N_{B}\\!=\\!\dot{\sum}_{ij}1$ as in Eq. (7).

The reason for $\Gamma^{\mu}\\!\simeq\\!0.5$ at $f\\!\to\\!0$ in Fig. 9(c) is
that the normalization factor $\hat{\Gamma}$ of $\Gamma_{ij}$ in Eq. (22) is
defined such that $\Gamma_{ij}\\!\to\\!1$ for isotropic configuration of
$\vec{\tau}$ and that the lattices are constructed isotropic such that the
mean value of $|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|$ is
$|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|\\!\simeq\\!0.5$. The large difference
between $\Gamma^{\rm r}({\rm bulk})$ and $\Gamma^{\theta,z}({\rm bulk})$ in
the incomplete skyrmion phase shown in Fig. 9(c) indicates a strong helical
order along the $\vec{e}^{\;r}$ direction than along the
$\vec{e}^{\;\theta,z}$ directions. Discontinuities of $\Omega$ and
$\Omega({\rm bulk})$ at $f\\!\to\\!0$ observed in Fig. 9(b) are owing to the
difference between $v_{ij}(\zeta_{ij}\\!=\\!1)$ and
$v_{ij}(\zeta_{ij}\\!=\\!-1)$ in Eq. (18) corresponding to the discontinuous
change in $Q_{\zeta}$ plotted in Fig. 7(b).

The surface DMI ($=\\!\Gamma(\rm side)$ and $\Gamma^{\theta,z}({\rm side})$)
plotted in the lower part of Figs. 9(a) and 9(c) are remarkably smaller than
the bulk DMI ($=\\!\Gamma({\rm bulk})$ and $=\\!\Gamma^{\theta,z}({\rm
bulk})$) indicated by ($\updownarrow$) in the incomplete skyrmion phase; the
$r$ component is not defined on the side surface because
$\vec{e}_{ij}\perp\vec{e}^{\;r}$ is expected on the side surface. Note also
that the direction dependence or anisotropy and the position dependence of
$\Gamma^{\mu}$ are apparent only for $f\\!>\\!0$, where
$\vec{\tau}\\!\propto\\!\vec{e}^{\;r}$ (Figs. 4(c), 5(f)). This dynamical
enhancement of the surface effects is considered to be an origin of the GC
effect on the skyrmions in nanodots and compatible with the assumption of zero
surface DMI in a confinement model of Ref. Diguet-etal-JMMM2023 . By contrast,
the direction and position dependence of $\Gamma^{\mu}$ disappears for
negative large $f(<0)$, where $\vec{\tau}$ is spiral along the $z$ direction
(Figs. 4(a), 5(a)–(c)).

No apparent difference is observed between the surface DMI ($=\\!\Gamma^{\rm
r}({\rm disk})$) on the upper/lower disk and the bulk DMI ($=\\!\Gamma^{\rm
r}({\rm bulk})$). The other component $\Gamma^{\rm\theta}({\rm disk})$ that is
not plotted is comparable to the bulk component $\Gamma^{\rm\theta}({\rm
bulk})$. Note that $\Gamma^{\mu},(\mu=r,\theta,z)$ are considered the $\mu$
components of a DMI vector Fert-etal-NatReview2017 ; Zhang-JPCM2020 ; Tokura-
Kanazawa-ACS2020 ; Gobel-etal-PhysRep2021 ; Wang-etal-JMMM2022  (Appendix D).
It is also noteworthy that $\Gamma^{\mu}({\rm bulk}),(\mu\\!=\\!r,\theta,z)$
increase with increasing $f$ in the whole range of $f$ and this behavior is
qualitatively consistent with the reported result that $D_{\rm ave}$ decreases
with increasing strain in Ref. YWang-etal-NatCom2020 , where $-D_{\rm ave}$
corresponds to $\Gamma$ plotted in Fig. 9(a) with the symbol ($\bullet$).

The direction dependence of $h_{\rm DM}^{\mu}({\rm bulk})$ (Fig. 9(d)) is
apparent in the target and incomplete skyrmion phases. The direction-dependent
energies $h_{\rm DM}^{\mu}({\rm side})$ plotted in Fig. 9(e) are comparable
with the bulk ones $h_{\rm DM}^{\mu}({\rm bulk})$ in 9(d). This result implies
that the reason for a large difference between $h_{\rm DM}({\rm side})$ and
$h_{\rm DM}({\rm bulk})$ in Fig. 8(b) is mainly due to the difference between
$\Gamma^{\mu}({\rm side})$ and $\Gamma^{\mu}({\rm bulk})$, because the DMI
energies plotted in Fig. 8(b) can also be expressed by using the direction-
dependent quantities in Eq. (32) such that

$\displaystyle H_{\rm DM}=\Gamma^{r}H_{\rm DM}^{r}+\Gamma^{\theta}H_{\rm
DM}^{\theta}+\Gamma^{z}H_{\rm DM}^{z}.$ (11)

Thus, we should emphasize that the reason for $|H_{\rm DM}({\rm
side})/N_{B}|\\!\ll\\!|H_{\rm DM}({\rm bulk})/N_{B}|$ in Fig. 8(b) is that the
intensive part of energy satisfies

$\Gamma^{\theta,z}({\rm side})\ll\Gamma^{\theta,z}({\rm bulk})$

as shown in Fig. 9(c). Note that a small discontinuity in $h_{\rm
DM}^{\theta}({\rm bulk})$ at the phase boundary between the skyrmion and
stripe, denoted by ”gap” in Fig. 9(d), indicates a weak first-order
transition.

It is also interesting to observe that the variations of $\Omega^{\mu}$ in
Fig. 9(f) are almost the same as those of $\Gamma^{\mu}$; the definition of
$\Omega_{ij}$ is considerably different from that of $\Gamma_{ij}$ (Appendix
B). $\Omega^{\theta,z}({\rm side})$ are very small compared with
$\Omega^{\theta,z}({\rm bulk})$ in the incomplete skyrmion phase (Fig. 9(f)).
The corresponding $h_{\rm ME}^{\mu}$ plotted in Fig. 9(g) are direction-
dependent only in the target and incomplete skyrmion phases similar to $h_{\rm
DM}^{\mu}$. The symbol $\Omega_{\rm eff}^{z}$ denotes an effective coupling
constant

$\displaystyle\Omega_{\rm eff}^{z}=\alpha f\Omega^{z},$ (12)

where the constant $\alpha$ is included. This $\Omega_{\rm eff}^{z}$
corresponds to the magnetic anisotropy $K_{\rm ave}$ in Ref. YWang-etal-
NatCom2020 , however, the $z$ component of $H_{\rm ME}$ does not always play a
role in the magnetic anisotropy in our model. Nevertheless, we find
qualitative consistency between the observation that $\Omega_{\rm eff}^{z}$
increases with increasing $f$ and the behavior of $K_{\rm ave}$ versus strain
in Ref. YWang-etal-NatCom2020 . Note that $\Omega_{\rm eff}^{z}\\!\to\\!0$ for
$f\\!\to\\!0$, and therefore no spontaneous magnetic anisotropy is implemented
in the model, as mentioned in Section II.2; all anisotropies are dynamically
generated and can be evaluated.

We should note that a strong anisotropy in $\Omega^{\mu}({\rm
bulk}),(\mu\\!=\\!r,\theta,z)$, which is considered a strain-induced
interaction anisotropy, appears only in the region of $f\\!>\\!0$ in the
proposed model as confirmed from Fig. 9(f). In the region of $f\\!>\\!0$, we
find $\Omega^{r}({\rm bulk})>\Omega^{\theta,z}({\rm bulk})$, which indicates
that the tensile stress makes the spin correlation $({\bf s}_{i}\cdot{\bf
s}_{j})^{2}$ along the $r$-axis to be stronger than those along the other
axes. This strong and anisotropic correlation implemented in $\Omega^{r}({\rm
bulk})$ of MEC, in addition to the enhanced effective FMI, competes with the
enhanced DMI along the $r$-axis. Consequently, the spin correlation along the
$r$ direction becomes stronger than that without MEC. This strain-enhanced
competition modified by ME along the $r$-axis increases the skyrmion size and
is considered the reason for the appearance of incomplete skyrmions. Note that
the inequality $\Omega^{r}({\rm bulk})>\Omega^{\theta,z}({\rm bulk})$ does not
imply anisotropy of ${\bf s}_{i}$ along the $r$ direction because
$Q_{s}^{r}\\!\to\\!-0.5$ in Fig. 6(a). In contrast, strain-mediated magnetic
anisotropies along the $\theta,z$ directions are expected in the incomplete
skyrmion phase owing to the observation that $\Gamma^{\theta,z}({\rm bulk})$
and $\Omega^{\theta,z}({\rm bulk})$ in the incomplete skyrmion phase are
larger than those in the stripe and target skyrmion phases (Figs. 9(c),(f)).
This expectation of strain-mediated magnetic anisotropy is consistent with the
increase in both $M^{z}$ (Fig. 7(a)) and $Q_{s}^{\theta}$ (Fig. 6(a)) with
increasing $f$ in the incomplete skyrmion phase.

We also find that $h_{\rm DM}^{\mu}({\rm bulk})$ increase as shown in Fig.
9(d), while $h_{\rm FM}^{\mu}({\rm bulk})$ decrease as shown in Fig. 9(h), as
$f$ increases in the incomplete skyrmion phase. This behavior represents the
competing nature of the interactions necessary for skyrmion emergence and is
consistent with the observation in the full Hamiltonians $H_{\rm FM}({\rm
bulk})/N_{B}$ and $H_{\rm DM}({\rm bulk})/N_{B}$ in Figs. 8(a),(b).

Figure 10:  (a) A triangle $ijk$ and the spins ${\bf s}_{i}$, ${\bf s}_{j}$
and ${\bf s}_{k}$ for the calculation of the topological charge $N_{\rm sk}$
on the upper disk. (b) $|N_{\rm sk}|$ vs. $f$ and (c) absolute topological
charge $N_{\rm sk}^{\rm ab}$ vs. $f$. The dashed vertical lines are
discontinuous phase boundaries between incomplete skyrmion and stripe phases
($N\\!=\\!1533$, $N\\!=\\!5430$), and a continuous phase boundary between
incomplete skyrmion and isolated skyrymion phases ($N\\!=\\!11932$). The
dashed vertical lines represent the phase boundaries depending on $N$.

The topological charge $N_{\rm sk}\\!=\\!\frac{1}{4\pi}\int d^{2}x\,{\bf
s}\\!\cdot\\!\frac{\partial{\bf s}}{\partial
x^{1}}\\!\times\\!\frac{\partial{\bf s}}{\partial x^{2}}$ is calculated using
the discrete expression

$\displaystyle N_{\rm sk}=\frac{1}{4\pi}\sum_{\it\Delta(ijk)}{\bf
s}_{i}\cdot\left({\bf s}_{j}-{\bf s}_{i}\right)\times\left({\bf s}_{k}-{\bf
s}_{i}\right)$ (13)

on the upper disk of the lattice, where $\sum_{\it\Delta(ijk)}$ denotes the
sum of triangles with vertices $i$, $j$ and $k$ (Fig. 10(a)). We observe a
topological phase transition between the skyrmion and stripe phases (Fig.
10(b)) at least for $N\\!\leq\\!5430$, although $|N_{\rm sk}|$ is smaller than
$|N_{\rm sk}|\\!=\\!1$ in the skyrmion phase owing to the nature of incomplete
skyrmions that the polarity change between the center and the periphery is
less than $\pi$. Another reason for $|N_{\rm sk}|\\!<\\!1$ is the numerical
error originating from discretization. The order of the transition is
considered first-order if we employ the definition of the first-order
transition such that there exists a physical quantity that is discontinuously
changing. Moreover, $|N_{\rm sk}|$ in the skyrmion phase for the
$N\\!\geq\\!5430$ lattices is larger than $|N_{\rm sk}|(=\\!0.5)$ of the
vortex, implying that the skyrmion configurations under tensile stress are
different from those of vortices, as discussed in Section III.3. Note that an
isolated skyrmion phase appears on the $N\\!=\\!11932$ lattice between the
incomplete skyrmion and stripe phases in the region of $f$ indicated by
($\leftrightarrow$) in Fig. 10(a). However, surface effects are expected to be
weak due to the increasing bulk volume on large nanodots, and for this reason,
a discontinuous transition is observed only on smaller nanodots.

The absolute of topological charge defined by

$\displaystyle N_{\rm sk}^{\rm ab}=\frac{1}{4\pi}\sum_{\it\Delta(ijk)}|{\bf
s}_{i}\cdot\left({\bf s}_{j}-{\bf s}_{i}\right)\times\left({\bf s}_{k}-{\bf
s}_{i}\right)|$ (14)

is convenient to see the total topological excitation (Fig. 10(c)). This
quantity is called absolute skyrmion number in Ref. Beg-etal-SRep2015 ; Beg-
etal-PRB2017 . We find that $N_{\rm sk}^{\rm ab}$ increases with increasing
$N$ and smoothly varies with respect to $f$.

##  IV Concluding remarks

In this paper, we numerically study the skyrmion stability in nanodots with
and without a magnetic field using the Finsler geometry (FG) modeling
technique and find that a single incomplete and/or isolated skyrmion is
stabilized at the center of the nanodots under radial tension. Metropolis
Monte Carlo technique is used for the simulations on 3D tetrahedral lattices
of cylindrical shape. Our Hamiltonian is defined with the variables of
electric spins and strain direction, and an Ising-like variable to dynamically
treat the anisotropy of the magneto-elastic coupling (MEC). In such FG
modeling, position- and direction-dependent coupling constants emerge
dynamically in response to the the tensile stress.

As a GC effect without stress, we confirm that the minimal external magnetic
field for skyrmions decreases when the nanodot diameter decreases. This
observation is consistent with previously reported experimental data by Y.
Wang et al. YWang-etal-NatCom2020  for materials with a spontaneous in-plane
magnetic anisotropy. When radial stresses are applied under zero magnetic
field, we confirm three different phases: incomplete (or isolated) skyrmions,
stripes, and target skyrmions, by varying the stress from positive (tension)
to negative (compression). In the incomplete skyrmion phase, we find that the
effective DMI coefficient increases with increasing tensile stress, and there
is a large difference in the effective coupling constants, including FMI and
MEC, between the surface and bulk in the nanodots. This modification of the
coupling constants is considered a strain-enhanced surface effect. We also
find a strain-induced $z$-direction alignment of spin configurations in the
polar-order parameter. These observations are consistent with the experimental
results reported in Ref. YWang-etal-NatCom2020  that the DMI coefficient and
spontaneous in-plane magnetic anisotropy are modified by strains. Moreover, we
observe that the incomplete skyrmions close to the stripe phase change to
isolated skyrmions when the disk diameter increases.

Thus, we find that a stable incomplete (or isolated) skyrmion emerges at the
center of the nanodot owing to the dynamical enhancement of the surface effect
by radially applied tensile stresses. In summary, the origins of skyrmion
stability in nanodots without a magnetic field are (i) dynamically enhanced or
strain-enhanced surface effects for position-dependent interactions and (ii)
dynamically generated or strain-induced interaction anisotropies of the FMI,
DMI, and MEC by radial stresses.

##  V Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

###### Acknowledgements.

This work is supported in part by Collaborative Research Project J23Ly07 of
the Institute of Fluid Science (IFS), Tohoku University. Numerical simulations
were performed on the supercomputer system AFI-NITY at the Advanced Fluid
Information Research Center, IFS, Tohoku University.

##  Appendix A Lattice construction

Figure 11:  Cylindrical lattices of size (a) $N\\!=\\!2803$, (b) $N\\!=\\!5430$ and (c) $N\\!=\\!11932$ composed of tetrahedrons. The diameter $d$ and the height $h$ are written in the plots; $h$ is fixed to $h\\!=\\!12b\\!=\\!6\sqrt{3}a$, where $a$ is the edge length of regular triangle on the side cylinder and $b$ is the height. The diameter $d\\!=\\!a/\sin(\pi/n_{e})$ is that of circumscribed circle of regular $n_{e}$-gon with $n_{e}$ given in Table 1, because the upper and lower disks are identified with the regular $n_{e}$-gon. (d) Distribution $h(\ell)$ of bond length $\ell$, where both $h(\ell)$ and $\ell$ are normalized. The peak at $\ell\\!\simeq\\!0.6$ corresponds to the bond of the cylindrical side surface. The height shares the symbol $h$ with the distribution $h(\ell)$, however, no confusion is expected.  Table 1: Lattices are characterized by the numbers $N,N_{B},N_{T},N_{\rm tet}$, which are the total number of vertices, bonds, triangles, and tetrahedrons. The Euler number $\chi(=\\!N\\!-\\!N_{B}\\!+\\!N_{T}\\!-\\!N_{\rm tet})$ satisfies $\chi\\!=\\!1$, implying that the lattices are filled with tetrahedrons. $n_{e}$ is the total number of edges of the upper/lower disk, which is a regular $n_{e}$-gon. The diameter $d$ and ratio $R(=\\!d/h)$ characterize the lattice shape. The height $h$ of the cylinders is $h\\!=\\!12b\\!=\\!6\sqrt{3}a$ (Fig. 11).  $N$ |  | $N_{B}$ |  | $N_{T}$ |  | $N_{\rm tet}$ |  | $\chi$ |  | $n_{e}$ |  | $d$ |  | $R(=d/h)$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1533 |  | 9970 |  | 16315 |  | 7877 |  | 1 |  | 33 |  | $10.5a$ |  | $1.02$
2083 |  | 13686 |  | 22501 |  | 10897 |  | 1 |  | 39 |  | $12.4a$ |  | $1.20$
5430 |  | 36720 |  | 61093 |  | 29802 |  | 1 |  | 65 |  | $20.7a$ |  | $1.99$
8465 |  | 57882 |  | 96700 |  | 47282 |  | 1 |  | 82 |  | $26.1a$ |  | $2.51$
11932 |  | 82097 |  | 137492 |  | 67326 |  | 1 |  | 98 |  | $31.2a$ |  | $3.00$

The height $h$ of the cylinder is fixed to $h\\!=\\!12b$ in the unit of
regular triangle height $b(=(\sqrt{3}/2)a)$ with the edge length $a$ (Fig.
11). The upper and lower disks are identified with a regular $n_{e}$-gon,
where $n_{e}$ is the total number of vertices or edges of the disk, and the
diameter $d$ is given by the circumscribed circle diameter
$d\\!=\\!a/\sin(\pi/n_{e})$. Therefore, the cylinders are also characterized
by the ratio $R\\!=\\!d/h$. We plot three cylinders of $R$ approximately given
by $R\\!=\\!1.20$, $R\\!=\\!1.99$ and $R\\!=\\!3.00$ in Figs. 11(a)–(c). The
cylindrical surface is composed of a regular triangle. The total number of
vertices on each of the top and bottom disks, including the edge, is given by
$\pi(d/2)^{2}/(\sqrt{3}/2)$, where $\sqrt{3}/2$ is the area of two regular
triangles in the case of $a\\!=\\!1$, and the vertices are randomly
distributed as uniformly as possible. The distribution of vertices on the
upper disk is the same as that on the bottom disk. The vertices inside the
surface are also randomly distributed as uniformly as possible, and the total
number is calculated using the volume of vertices inside the surface, in which
the space is assumed to be filled with regular tetrahedrons. Vertex volume is
calculated as $(\sqrt{2}/48)\\!\times\\!22.8$, where $(\sqrt{2}/48$ is a
quarter of the regular tetrahedron volume and $22.8$ is an approximate number
of the regular tetrahedrons emanating from the inner vertices calculated by
$4\pi/\Omega$ with the solid angle $\Omega$ of a vertex of the tetrahedron.
The peak of the distributions $h(\ell)$ of the bond length $\ell$ plotted in
11(d) corresponds to the bonds on the cylindrical surface. Table 1 lists the
numbers characterizing the lattices.

The mean value of the total number of tetrahedrons
$n_{ij}(=\\!\sum_{\Delta(ij)}1)$ sharing bond $ij$ are shown in Table 2. We
find that
$\frac{\overline{n}(\rm{side})}{\overline{n}(\rm{bulk})}\\!\simeq\\!0.51$.
This ratio is one of the origins of the surface effect on GC that the surface
DMI on the side surface is smaller than the bulk DMI. The surface effect will
be described in Appendix B. To show that the direction $\vec{e}_{ij}$ of bond
$ij$ is isotropic, the lattice averages
$(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}$,
$(\vec{e}_{ij}\cdot\vec{e}^{\;\theta})^{2}$ and
$(\vec{e}_{ij}\cdot\vec{e}^{\;z})^{2}$ are listed. We find from
$(\vec{e}_{ij}\cdot\vec{e}^{\;\mu})^{2}\\!\simeq\\!1/3,(\mu\\!=\\!r,\theta,z)$
that the bond direction $\vec{e}_{ij}$ is isotropically distributed.
Interaction anisotropy comes only from a dynamical reason in contrast to the
surface effect.

Table 2: The mean value of the total number of tetrahedrons sharing bond $ij$ are calculated using $\overline{n}=\\!\frac{\sum_{ij}\sum_{\Delta(ij)}1}{\sum_{ij}1}$ on upper and lower disks, side surface and inside, denoted by $\overline{n}(\rm{disk})$, $\overline{n}(\rm{side})$ and $\overline{n}(\rm{bulk})$, respectively. $\overline{n}$ is the mean value for all bonds $ij$. $(\vec{e}_{ij}\cdot\vec{e}^{\;\mu})^{2},(\mu\\!=\\!r,\theta,z)$ are the lattice averages of squared $\mu$ components of bond direction $\vec{e}_{ij}$, where $\vec{e}^{\;\mu}(\mu\\!=\\!r,\theta,z)$ denote the $(r,\theta,z)$ direction at the center of bond $ij$ (see Appendix C).  $N$ |  | $\overline{n}(\rm{disk})$ |  | $\overline{n}(\rm{side})$ |  | $\overline{n}(\rm{bulk})$ |  | $\overline{n}$ |  | $(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}$ |  | $(\vec{e}_{ij}\cdot\vec{e}^{\;\theta})^{2}$ |  | $(\vec{e}_{ij}\cdot\vec{e}^{\;z})^{2}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1533 |  | 2.87 |  | 2.64 |  | 5.15 |  | 4.74 |  | 0.334 |  | 0.337 |  | 0.326
2083 |  | 2.85 |  | 2.65 |  | 5.16 |  | 4.78 |  | 0.340 |  | 0.336 |  | 0.330
5430 |  | 2.80 |  | 2.64 |  | 5.17 |  | 4.87 |  | 0.333 |  | 0.330 |  | 0.338
8465 |  | 2.80 |  | 2.65 |  | 5.17 |  | 4.90 |  | 0.329 |  | 0.329 |  | 0.340
11932 |  | 2.80 |  | 2.65 |  | 5.17 |  | 4.92 |  | 0.330 |  | 0.329 |  | 0.342

##  Appendix B Discretization of direction-dependent interaction coefficient

Figure 12:  (a) Trajectory of unit Finsler length $v_{ij}$ at $i$ is not
always spherical or symmetric in Euclidean space, (b) the trajectory for
$v_{ij}$ assumed in this paper, where $v_{0}\\!=\\!0$ is assumed for
simplicity. In (b),
$v_{ij}\\!=\\!\cos\theta_{ij}\\!>\\!v_{ij}\\!=\\!\sin\theta_{ij}$ for the
given strain direction $\vec{\tau}_{i}(\in S^{2}/2)$. Finsler length $v_{ij}$
does not correspond to the real length in materials but it is used to define
direction-dependent interaction anisotropy between spins ${\bf s}$ and strains
$\vec{\tau}$.

Finsler geometry framework provides interaction anisotropy between two
particles via a direction-dependent Finsler length $v_{ij}$ (Fig. 12(a)). The
Finsler length $v_{ij}$ for $H_{\rm FM},H_{\rm DM}$ is fixed to be cosine
type, and $v_{ij}$ for $H_{\rm ME}$ is defined to be dynamically convertible
to cosine type or sine type depending on $\zeta_{ij}$ ($\in\\{-1,1\\}$ on bond
$ij$):

$\displaystyle v_{ij}=|\cos\theta_{ij}|+v_{0},\quad(v_{0}=0.1;\;H_{\rm
FM},H_{\rm DM})$ (15) $\displaystyle
v_{ij}=\left\\{\begin{array}[]{@{\,}ll}|\cos\theta_{ij}|+v_{0},\;(\zeta_{ij}=1)\\\
\sin\theta_{ij}+v_{0},\;(\zeta_{ij}=-1)\end{array}\right.,(v_{0}=0.2;\;H_{\rm
ME}),$ (18)

where $\cos\theta_{ij}$ and $\sin\theta_{ij}$ (Fig. 12(b)) are

$\displaystyle\cos\theta_{ij}=\vec{\tau}_{i}\cdot\vec{e}_{ij},\quad\sin\theta_{ij}=\sqrt{1-|\vec{\tau}_{i}\cdot\vec{e}_{ij}|^{2}}\;(\geq
0),$ (19)

and $v_{0}$ is a small cutoff; the smaller $v_{0}$, the larger interaction
anisotropy.

Using these $v_{ij}$, the Finsler metric is defined by

$\displaystyle g_{ab}=\begin{pmatrix}v_{12}^{-2}&0&0\\\ 0&v_{13}^{-2}&0\\\
0&0&v_{14}^{-2}\end{pmatrix},\quad\sqrt{g}=\sqrt{\det
g_{ab}}=v_{12}^{-1}v_{13}^{-1}v_{14}^{-1},\quad g^{ab}=(g_{ab})^{-1}.$ (20)

These expressions are used in the continuous Hamiltonians

$\displaystyle\begin{split}&H_{{\rm
FM}}=\frac{1}{2}\int\sqrt{g}d^{3}xg^{ab}\frac{\partial{\bf s}}{\partial
x^{a}}\cdot\frac{\partial{\bf s}}{\partial x^{b}},\\\ &H_{{\rm
DM}}=\int\sqrt{g}d^{3}xg^{ab}\frac{\partial{\vec{r}}}{\partial x^{a}}\cdot{\bf
s}\times\frac{\partial{\bf s}}{\partial x^{b}},\\\ &H_{{\rm
ME}}=\frac{1}{4}\int\sqrt{g}d^{3}x\left(g^{ab}\frac{\partial{\bf s}}{\partial
x^{a}}\cdot\frac{\partial{\bf s}}{\partial x^{b}}\right)^{2},\end{split}$ (21)

where the factor $f$ is removed from $H_{\rm ME}$ for simplicity. From these
continuous Hamiltonians, we obtain the discrete ones in Eq. (2) by replacing
the differentials with differences such that $\partial_{1}{\bf s}\\!\to\\!{\bf
s}_{2}\\!-\\!{\bf s}_{1}$, $\partial_{2}{\bf s}\\!\to\\!{\bf s}_{3}\\!-\\!{\bf
s}_{1}$ and $\partial_{3}{\bf s}\\!\to\\!{\bf s}_{4}\\!-\\!{\bf s}_{1}$ for a
local coordinate system with the origin at the vertex $1$, and by replacing
$\int\sqrt{g}d^{3}x$ with $\sum_{\it\Delta}\sum_{ij(\Delta)}$; the sum over
tetrahedrons $\sum_{\it\Delta}$ and the sum over bonds $\sum_{ij(\Delta)}$ of
the tetrahedron $\Delta$, and by summing over all symmetric expressions
obtained by the cyclic replacements
$1\\!\to\\!2,2\\!\to\\!3,3\\!\to\\!4,4\\!\to\\!1$ (Fig. 1(c)). For the
discretization of $H_{\rm ME}$, only quadratic terms such as $({\bf
s}_{i}\cdot{\bf s}_{j})^{2}$ in the expansion of
$\left(g^{ab}\frac{\partial{\bf s}}{\partial x^{a}}\cdot\frac{\partial{\bf
s}}{\partial x^{b}}\right)^{2}$ are used.

In these discrete calculations, we extract an intensive part of the energies,
the so-called position- and direction-dependent coupling constants

$\displaystyle\Gamma_{ij}=\hat{\Gamma}^{-1}\gamma_{ij},\quad\hat{\Gamma}=\bar{n}\,\langle\gamma^{\rm{iso}}\rangle,\qquad({\rm
FMI,DMI}),$ (22)
$\displaystyle\begin{split}&\gamma_{12}=\frac{v_{12}}{v_{13}v_{14}}+\frac{v_{21}}{v_{23}v_{24}},\quad\gamma_{13}=\frac{v_{13}}{v_{12}v_{14}}+\frac{v_{31}}{v_{32}v_{34}},\\\
&\gamma_{14}=\frac{v_{14}}{v_{13}v_{14}}+\frac{v_{41}}{v_{42}v_{43}},\quad\gamma_{23}=\frac{v_{23}}{v_{21}v_{24}}+\frac{v_{32}}{v_{31}v_{34}},\\\
&\gamma_{24}=\frac{v_{24}}{v_{21}v_{23}}+\frac{v_{42}}{v_{41}v_{43}},\quad\gamma_{34}=\frac{v_{34}}{v_{31}v_{32}}+\frac{v_{43}}{v_{41}v_{42}}\end{split}$
(23)

for $H_{\rm FM}$ and $H_{\rm DM}$, and

$\displaystyle\Omega_{ij}=\hat{\Omega}^{-1}\omega_{ij},\quad\hat{\Omega}=\bar{n}\,\langle\omega^{\rm{iso}}\rangle,\qquad({\rm
MEC}),$ (24)
$\displaystyle\begin{split}&\omega_{12}=\frac{v_{12}^{3}}{v_{13}v_{14}}+\frac{v_{21}^{3}}{v_{23}v_{24}},\quad\omega_{13}=\frac{v_{13}^{3}}{v_{12}v_{14}}+\frac{v_{31}^{3}}{v_{32}v_{34}},\\\
&\omega_{14}=\frac{v_{14}^{3}}{v_{13}v_{14}}+\frac{v_{41}^{3}}{v_{42}v_{43}},\quad\omega_{23}=\frac{v_{23}^{3}}{v_{21}v_{24}}+\frac{v_{32}^{3}}{v_{31}v_{34}},\\\
&\omega_{24}=\frac{v_{24}^{3}}{v_{21}v_{23}}+\frac{v_{42}^{3}}{v_{41}v_{43}},\quad\omega_{34}=\frac{v_{34}^{3}}{v_{31}v_{32}}+\frac{v_{43}^{3}}{v_{41}v_{42}}\end{split}$
(25)

for $H_{\rm ME}$. The difference between $\gamma_{ij}$ and $\omega_{ij}$ comes
from the quadratic form of the integrand of $H_{\rm ME}$ in Eq. (21).

Here, we describe only $\Gamma_{ij}$ for simplicity; the definitions of these
two are the same except for the expressions of $\gamma_{ij}$ and $\omega_{ij}$
and their $v_{ij}$ in Eqs. (15) and (18). We should note that $\hat{\Gamma}$
in $\Gamma_{ij}$ is an irrelevant normalization constant, and a non-trivial
interaction comes from $\gamma_{ij}$ in Eq. (23). The position- and direction-
dependence of interaction is effectively introduced as an intensive part of
energy in Hamiltonian via $\Gamma_{ij}$ in the FG modeling.

In Eq. (22), $\bar{n}(=\\!\frac{\sum_{\Delta(ij)}1}{\sum_{ij}1})$ denotes the
mean value of $n_{ij}\\!=\\!\sum_{\Delta(ij)}1$, which is the total number of
tetrahedrons $\Delta$ sharing bond $ij$ (Fig. 2(c)), and
$\langle\gamma^{\rm{iso}}\rangle\\!=\\!\langle\frac{\sum_{\Delta}\sum_{ij(\Delta)}\gamma_{ij}}{\sum_{\Delta}\sum_{ij(\Delta)}1}\rangle$
denotes the ensemble average calculated with 1000 isotropic configurations of
$\vec{\tau}$. Note that the random configurations for
$\langle\gamma^{\rm{iso}}\rangle$ in $\Omega_{ij}$ include a randomly
distributed $\zeta$ on the bonds. The reason for introducing the
$v_{0}$-dependent factor $\langle\gamma^{\rm{iso}}\rangle$ is that
$\gamma_{ij}$ varies depending on $v_{0}$ in Eqs. (15) and (18), and the
parameters $\lambda$, $D$ and $\alpha$ should be changed according to any
variations in $v_{0}$ if $\langle\gamma^{\rm{iso}}\rangle$ is not included.

We briefly show the surface effect for GC. There are two possible sources;
static and dynamical, for the surface effect. The first one comes from the
lattice structure and the other from a dynamical effect in FG modeling. The
first one, which is the static part, is shown using the coupling constant in
Eq. (22) as follows. It is noteworthy that the lattice average
$\overline{\Gamma}\\!=\\!\frac{\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}}{\sum_{\Delta}\sum_{ij(\Delta)}}$
of $\Gamma_{ij}$ for an isotropic configuration of $\vec{\tau}$ satisfy
$\overline{\Gamma}(=\\!\overline{\Gamma}_{ij})\\!\simeq\\!\bar{n}^{-1}$.
Indeed, let $\gamma_{ij}$ be a coefficient obtained with an isotropic
configuration of $\vec{\tau}$ and $\overline{\gamma}$ be the lattice average
of $\gamma_{ij}$. Then, we have that
$\overline{\Gamma}\\!=\\!\bar{n}^{-1}\langle\gamma^{\rm{iso}}\rangle^{-1}\overline{\gamma}\\!=\\!\frac{1}{\bar{n}\langle\gamma^{\rm{iso}}\rangle}\frac{\sum_{\Delta}\sum_{ij(\Delta)}\gamma_{ij}}{\sum_{\Delta}\sum_{ij(\Delta)}1}\\!\simeq\\!\bar{n}^{-1}$
because the lattice average
$\overline{\gamma}(\\!=\\!\frac{\sum_{\Delta}\sum_{ij(\Delta)}\gamma_{ij}}{\sum_{\Delta}\sum_{ij(\Delta)}1})$
for isotropic $\vec{\tau}$ is close to $\langle\gamma^{\rm{iso}}\rangle$.
Using this $\overline{\Gamma}$, the Hamiltonian of DMI in Eq. (2) for
isotropic $\vec{\tau}$ can be written as

$\displaystyle\begin{split}H_{\rm
DM}=&\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j}\\\
\to&\sum_{ij}(\bar{n}^{-1}\sum_{\Delta(ij)}1)\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j}=\sum_{ij}(n_{ij}/\bar{n})\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j},\end{split}$ (26)

where $\Gamma_{ij}$ is replaced by $\overline{\Gamma}\\!\simeq\\!\bar{n}^{-1}$
and the summation convention
$\sum_{ij}\sum_{\Delta(ij)}\\!=\\!\sum_{\Delta}\sum_{ij(\Delta)}$ is used. The
replacement of $\Gamma_{ij}$ with a constant $\overline{\Gamma}$ is to neglect
all dynamical effects implemented in FG modeling and to leave only
contribution from the lattice structure. This form of $H_{\rm
DM}\\!=\\!\sum_{ij}(n_{ij}/\bar{n})\vec{e}_{ij}\cdot{\bf s}_{i}\times{\bf
s}_{j}$ includes the weight $n_{ij}/\bar{n}$, which modifies the interaction
$\vec{e}_{ij}\cdot{\bf s}_{i}\times{\bf s}_{j}$ to be dependent on the
position of bond $ij$ whether it is on the surface or inside. The number
$n_{ij}\\!=\\!\sum_{\Delta(ij)}1$ is smaller (larger) than the mean value
$\bar{n}$ on the surface (bulk) (Table 2 and Fig. 2(c)). Thus, we confirm that
this weight difference is one possible origin of the surface effect in the
small strain region $f\\!\to\\!0$, where $\vec{\tau}$ is isotropic. This GC
effect at $f\\!=\\!0$ is confirmed in Section III.1.

When an anisotropic configuration of $\vec{\tau}$ (and $\zeta$) is reflected
in $\gamma_{ij}$, the weight difference is further enlarged in the region of
$f\\!>\\!0$ in the incomplete skyrmion phase (Fig. 9(a)). This enlargement
caused by anisotropic $\vec{\tau}$ is the dynamical part of the surface effect
for GC. Here, we should emphasize that the surface effect is the difference in
the interaction weight in $H_{\rm DM}$; effective DMI coupling on the side
surface is considerably smaller than the bulk one (Fig. 9(a)).

##  Appendix C Decomposition of interaction coefficient and the corresponding
energies

Figure 13:  (a) Polar coordinate position of vertex $i$. The coordinate origin
is at the center of nanodot. (b) Spin vector ${\bf s}_{i}$ is decomposed into
the radial $\vec{e}_{i}^{\;r}$, tangential $\vec{e}_{i}^{\;\theta}$ and
$\vec{e}^{\;z}$ directions at vertex $i$, where ${\bf s}_{i}^{\mu}\\!=\\!({\bf
s}_{i}\\!\cdot\\!\vec{e}_{i}^{\;\mu})\vec{e}_{i}^{\;\mu},(\mu\\!=\\!r,\theta,z)$.
(c) The unit vector $\vec{e}_{ij}$ from vertices $i$ to $j$ can be decomposed
into the $\vec{e}^{\;r}$, $\vec{e}^{\;\theta}$ and $\vec{e}^{\;z}$ components,
where $\vec{e}^{\;r}$ and $\vec{e}^{\;\theta}$ are given by Eq. (27) and
$\vec{e}_{ij}^{\;{\mu}}\\!=\\!(\vec{e}_{ij}\\!\cdot\\!\vec{e}^{\;\mu})\vec{e}^{\;\mu},(\mu\\!=\\!r,\theta,z)$.

Let $\vec{e}_{ij}$ be the unit vector from vertices $i$ to $j$ and
$\vec{e}_{i}^{\;r}$, $\vec{e}_{i}^{\;\theta}$ and $\vec{e}^{\;z}$ be the unit
vectors along the $r$, $\theta$ and $z$ directions at the vertex $i$ as
illustrated in Figs. 13(a),(b), where the coordinate origin is at the center
of nanodot and $\vec{e}_{i}^{\;z}$ is independent of $i$ and written as
$\vec{e}^{\;z}$. Here, we define unit vectors

$\displaystyle\vec{e}^{\;\mu}=\left(\vec{e}_{i}^{\;\mu}+\vec{e}_{j}^{\;\mu}\right)/\|\vec{e}_{i}^{\;\mu}+\vec{e}_{j}^{\;\mu}\
\|,\quad(\mu=r,\theta)$ (27)

to represent the $r$ and $\theta$ directions at the center of bond $ij$ (Fig.
13(c)). Then, we have a decomposition of $\vec{e}_{ij}$ such that
$\vec{e}_{ij}\\!=\\!(\vec{e}_{ij}\cdot\vec{e}^{\;r})\vec{e}^{\;r}\\!+\\!(\vec{e}_{ij}\cdot\vec{e}^{\;\theta})\vec{e}^{\;\theta}\\!+\\!(\vec{e}_{ij}\cdot\vec{e}^{\;z})\vec{e}^{\;z}$.
Since $\Gamma_{ij}$ is defined on bond $ij$ along $\vec{e}_{ij}$,
$\Gamma_{ij}$ is naturally direction-dependent. To extract this direction
dependence, we use the vector $\Gamma_{ij}\vec{e}_{ij}$, which is decomposed
as

$\displaystyle\begin{split}&\Gamma_{ij}\vec{e}_{ij}=\Gamma_{ij}^{r}\vec{e}^{\;r}+\Gamma_{ij}^{\theta}\vec{e}^{\;\theta}+\Gamma_{ij}^{z}\vec{e}^{\;z},\\\
&\Gamma_{ij}^{r}(\vec{\tau})=\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;r}),\;\Gamma_{ij}^{\theta}(\vec{\tau})=\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;\theta}),\;\Gamma_{ij}^{z}(\vec{\tau})=\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;z}).\end{split}$
(28)

This decomposition of $\Gamma_{ij}\vec{e}_{ij}$ corresponds to the
decomposition of DMI vector (Appendix D) and is considered a physically
meaningful decomposition. Using these ”microscopic” direction-dependent
factors $\Gamma_{ij}^{\mu}$, the ”macroscopic” direction-dependent
coefficients are naturally defined by

$\displaystyle\Gamma^{\mu}=\frac{1}{N_{B}}\sum_{ij}\sum_{\Delta(ij)}|\Gamma_{ij}^{\mu}(\vec{\tau})|,\quad(\mu=r,\theta,z).$
(29)

The meaning of $\Gamma^{\mu}$ is as follows: If
$\Gamma^{r}\\!>\\!\Gamma^{\theta,z}$ for example, we understand from this
inequality that $\Gamma_{ij}$ along bond $ij$ parallel or almost parallel to
$\pm\vec{e}^{\;r}$ is larger than $\Gamma_{ij}$ along bond $ij$ parallel or
almost parallel to $\pm\vec{e}^{\;\theta,z}$, because the direction of
$\vec{e}_{ij}$ is uniformly distributed.

Using these $\Gamma^{\mu}$, the corresponding $H_{\rm DM}$ can also be
decomposed into direction-dependent energies:

$\displaystyle\begin{split}&H_{\rm
DM}^{r}=\frac{\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf
s}_{j}}{\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{r}|},\\\
&H_{\rm
DM}^{\theta}=\frac{\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;\theta})^{2}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf
s}_{j}}{\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{\theta}|},\\\
&H_{\rm
DM}^{z}=\frac{\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})(\vec{e}_{ij}\cdot\vec{e}^{\;z})^{2}\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf
s}_{j}}{\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{z}|},\end{split}$
(30)

where $N_{B}\\!=\\!\sum_{ij}1$ is the total number of bonds. Using these
expressions, we have

$\displaystyle\begin{split}H_{\rm
DM}=&\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j}\\\
=&\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})\left[(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}+(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}+(\vec{e}_{ij}\cdot\vec{e}^{\;r})^{2}\right]\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j}\\\
=&\left(\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{r}(\vec{\tau})|\right)H_{\rm
DM}^{r}+\left(\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{\theta}(\vec{\tau})|\right)H_{\rm
DM}^{\theta}\\\
&+\left(\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{z}(\vec{\tau})|\right)H_{\rm
DM}^{z},\end{split}$ (31)

and therefore, the lattice average is given by

$\displaystyle H_{\rm DM}=\Gamma^{r}H_{\rm DM}^{r}+\Gamma^{\theta}H_{\rm
DM}^{\theta}+\Gamma^{z}H_{\rm DM}^{z}.$ (32)

Eq. (32) is satisfied only in the lattice average. However, for the ensemble
average $\langle*\rangle$ calculated via MC simulations, it is reasonable to
assume $\Gamma^{\mu}\\!\simeq\\!\langle\Gamma^{\mu}\rangle$ and $H_{\rm
DM}^{\mu}\\!\simeq\\!\langle H_{\rm DM}^{\mu}\rangle$ in the equilibrium
configurations, and hence, we use this decomposition of $H_{\rm DM}$ for the
simulation data to evaluate the direction-dependent effective coefficients and
energies. We should note that $\Gamma_{ij}$ in $H_{\rm FM}$ and $\Omega_{ij}$
in $H_{\rm ME}$ in Eq. (2) are defined on bond $ij$ and can also be decomposed
into $\mu$ components in the same procedure described in Eq. (28). The
directional dependence of $H_{\rm FM}$ and $H_{\rm ME}$ can also be evaluated
with the same procedure for $H_{\rm DM}$.

For the calculation of $H_{\rm DM}^{\mu}({\rm side})$ on the side surface
(Figs. 9(d),(e)), the sum
$\sum_{ij}\sum_{\Delta(ij)}(=\\!\sum_{\Delta}\sum_{ij(\Delta)})$ in the
expressions of Eq. (30) should be limited to the bonds $ij$ on the side
surface. In this case, $N_{B}$ is the total number of bonds on the side
surface.

Figure 14:  (a) Position-dependent Finsler length $v$ vs. $f$ for the cosine
type $v_{ij}\\!=\\!|\cos\theta_{ij}|\\!+\\!v_{0}$ ($v_{0}\\!=\\!0.1$) for DMI
(and FMI) obtained on the $N\\!=\\!5430$ lattice, where
$\cos\theta_{ij}\\!=\\!\vec{\tau}\cdot\vec{e}_{ij}$ (Fig. 12(b)) and (b) the
direction-dependent ($\Leftrightarrow(r,\theta,z)$-dependent) Finsler length
$v^{\mu}$ vs. $f$.

Finally, in this appendix, we plot the unit Finsler length for DMI in Fig.
14(a) calculated using the lattice average

$\displaystyle\begin{split}v=\frac{1}{N_{B}}\dot{\sum}_{ij}\frac{v_{ij}+v_{ji}}{2}\end{split}$
(33)

on the side surface (side), the upper/lower disk (disk), and the inside
(bulk). The components $v^{\mu},(\mu\\!=\\!r,\theta,z)$ (Fig. 14(b)) are
obtained with the absolute components of the vector
$\vec{v}_{ij}(=\\!\frac{v_{ij}\\!+\\!v_{ji}}{2}\vec{e}_{ij}$) such that

$\displaystyle\begin{split}v^{\mu}=\frac{1}{N_{B}}\dot{\sum}_{ij}\frac{v_{ij}+v_{ji}}{2}|e_{ij}^{\mu}|,\end{split}$
(34)

where $N_{B}=\dot{\sum}_{ij}1$. The meaning of $\dot{\sum}_{ij}$, which is the
same as that in Eq. (7), is that the sum over bonds ${\sum}_{ij}$ is limited
to the disk, or side, or bulk. The index symmetrization of $v_{ij}$ for $v$
and $v^{\mu}$ is necessary owing to $v_{ij}\\!\not=\\!v_{ji}$, in contrast
with the case of $\Gamma_{ij}$, where $\Gamma_{ij}\\!=\\!\Gamma_{ji}$ (Eq.
(22)). Without the symmetrization, ${\sum}_{ij}v_{ij}$ depends on the labeling
of vertex $i$ and becomes ill-defined because the symbol $\sum_{ij}$ implies
the sum over pairs of vertices $i$ and $j$ connected by bond under the
condition $i\\!>\\!j$. The sum over bonds
$(1/N_{B}){\sum}_{ij}\frac{v_{ij}+v_{ji}}{2}$ is the same as the sum over
bonds $(1/2N_{B}){\sum}_{ij}v_{ij}$ without the condition $i\\!>\\!j$. Finsler
length $v$ is direction independent, whereas $v^{\mu}$ depends on the
direction $\mu(=\\!r,\theta,z)$.

We find from Fig. 14(a) that $v({\rm bulk})$ is almost constant and $v({\rm
side})$ and $v({\rm disk})$ depend on $f$ for positive $f$. In the incomplete
skyrmion phase, $v({\rm side})\\!\ll\\!v({\rm bulk})$, which is considered to
be the origin of the surface effect in the FG modeling. From Fig. 14(b), we
find that $v^{r}({\rm bulk})$ ($v^{\theta,z}({\rm bulk})$) increases
(decrease) with increasing $f$, and hence, the variation of $v^{r}({\rm
bulk})$ is opposite to those of $v^{\theta,z}({\rm bulk})$. The observation
that $v^{r}({\rm bulk})$ increases for $f\\!>\\!0$ indicates that $v_{ij}$ on
bond $ij$ in the direction $\vec{e}_{ij}$ almost parallel to
$\vec{e}_{ij}^{\;r}$ is larger than that on bonds $ij$ of $\vec{e}_{ij}$
almost perpendicular to $\vec{e}_{ij}^{\;r}$. It is noteworthy that this
behavior of $v^{\mu}$ is reflected in $\Gamma({\rm bulk})$ and
$\Gamma^{\mu}({\rm bulk})$ in Figs. 9(a),(c), though $\Gamma_{ij}$ is a
rational function of $v_{ij}$ in Eqs. (22) and (23). We should note that in
the isotropic case at $f\\!\to\\!0$, $v_{ij}$ is at random and independent of
$v_{ji}$, whereas $v_{ij}\\!\simeq\\!v_{ji}$ in the anisotropic case $|f|>0$
where $\tau$ is uniformly aligned. We find no position dependence of $v$ on
$f$ at $f\\!=\\!0$. The directional independence of $v$ implies that the
original isotropic model ($\Leftrightarrow\\!\Gamma_{ij}\\!=\\!{\rm constant}$
) is restored in the limit of $f\\!\to\\!0$, where $\vec{\tau}$ randomly
distributes and effectively has no influence on ${\bf s}$.

##  Appendix D DMI vector

Effective coupling constants
$\sum_{\Delta(ij)}\Gamma_{ij}^{\mu}(\vec{\tau}),\;(\mu=r,\theta,z)$ are
understood to be $\mu$ components of a DMI vector. It is interesting to
observe that DMI vector is composed of primitive objects in FG modeling.
First, $H_{\rm DM}$ can be written as

$\displaystyle H_{\rm DM}=\sum_{ij}{\bf D}_{ij}\cdot{\bf
s}_{i}\\!\times\\!{\bf s}_{j},$ (35)

where ${\bf D}_{ij}$ is a DMI vector and $\sum_{ij}$ denotes the sum over
bonds $ij$. Second, $H_{\rm DM}$ in Eq. (2) can be expressed by the sum over
bonds:

$\displaystyle H_{\rm
DM}=\sum_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\vec{e}_{ij}\cdot{\bf
s}_{i}\times{\bf s}_{j},$ (36)

where $\sum_{\Delta(ij)}$ denotes the sum over tetrahedrons sharing bond $ij$.
The expression in Eq. (36) is obtained from $H_{\rm DM}$ in Eq. (2) by using
the identity $\sum_{\Delta}\sum_{ij(\Delta)}\\!=\\!\sum_{ij}\sum_{\Delta(ij)}$
(Fig. 2 and the caption). Therefore, comparing Eqs. (35) and (36), we obtain

$\displaystyle\begin{split}{\bf
D}_{ij}=&\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\vec{e}_{ij}\\\
=&\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\left[(\vec{e}_{ij}\cdot\vec{e}^{\;r})\vec{e}^{\;r}\\!+\\!(\vec{e}_{ij}\cdot\vec{e}^{\;\theta})\vec{e}^{\;\theta}\\!+\\!(\vec{e}_{ij}\cdot\vec{e}^{\;z})\vec{e}^{\;z}\right],\end{split}$
(37)

and the $\mu$ component

$\displaystyle{D}_{ij}^{\mu}={\bf
D}_{ij}\\!\cdot\\!\vec{e}^{\;\mu}=\left(\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})\right)\vec{e}_{ij}\cdot\vec{e}^{\;\mu},\;(\mu=r,\theta,z),$
(38)

where $\vec{e}^{\;\mu}$ represents $\vec{e}^{\;r}$, $\vec{e}^{\;\theta}$ and
$\vec{e}^{\;z}$ (see Eq. (27) and Fig. 13(c)). Since
$\vec{e}_{ij}\cdot\vec{e}^{\;\mu}$ in Eq. (38) is independent of the sum
$\sum_{\Delta(ij)}$ and $\Gamma_{ij}(\vec{\tau})\\!>\\!0$, we find
$|{D}_{ij}^{\mu}|=\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|$.
If we write the lattice average of $|{D}_{ij}^{\mu}|$ by $D^{\mu}$ such that
$D^{\mu}\\!=\\!\frac{1}{N_{B}}\sum_{ij}|D_{ij}^{\mu}(\vec{\tau})|\\!=\\!\frac{1}{N_{B}}\sum_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|$,
then
$\Gamma^{\mu}(=\\!\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}|\Gamma_{ij}^{\mu}(\vec{\tau})|)$
is identified with $D^{\mu}$:

$\displaystyle\begin{split}{D}^{\mu}=&\frac{1}{N_{B}}\sum_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|\\\
=&\frac{1}{N_{B}}\sum_{\Delta}\sum_{ij(\Delta)}\Gamma_{ij}(\vec{\tau})|\vec{e}_{ij}\cdot\vec{e}^{\;\mu}|\\\
=&\Gamma^{\mu},\;(\mu=r,\theta,z).\end{split}$ (39)

Using these relations, Eq. (32) can also be written as

$\displaystyle H_{\rm DM}=D^{r}H_{\rm DM}^{r}+D^{\theta}H_{\rm
DM}^{\theta}+D^{z}H_{\rm DM}^{z}.$ (40)

Therefore, the extensive parts $H_{\rm DM}^{\mu},(\mu\\!=\\!r,\theta,z)$ are
considered to be physically meaningful direction-dependent energies, expressed
to be determined as a response to the intensive parts $D^{\mu}$, which are
also dynamically determined.

The symbol $\Gamma$ is used for the lattice average defined by

$\displaystyle\Gamma=\frac{1}{N_{B}}\sum_{ij}\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau}),$
(41)

which corresponds to the lattice average of $D_{ij}\\!=\\!\|{\bf
D}_{ij}\|\\!=\\!\sum_{\Delta(ij)}\Gamma_{ij}(\vec{\tau})$.

## References

  * (1) N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka and R. Wiesendanger, Writing and deleting single magnetic skyrmions, Science 341 (6146), 636-639 (2013), https://doi.org/10.1126/science.1240573.
  * (2) A. Fert, N. Reyren and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews 2, 17031 (2017), https://doi.org/10.1038/natrevmats.2017.31.
  * (3) X. Zhang, Y. Zhou, K. M. Song, T.-E. Park, J. Xia, M. Ezawa, X. Liu, W. Zhao, G. Zhao and S. Woo, Skyrmion-electronics: writing, deleting, reading and processing magnetic skyrmions toward spintronic applications, J. Phys.: Condens. Matter 32 143001 (2020), https://doi.org/10.1088/1361-648X/ab5488.
  * (4) Y. Tokura and N. Kanazawa, Magnetic Skyrmion Materials, Chem. Rev. 121, 5, pp. 2857–2897 (2021), https://doi.org/10.1021/acs.chemrev.0c00297.
  * (5) B$\ddot{\rm o}$rge G$\ddot{\rm o}$bel, I. Mertig, O. A. Tretiakov, Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles, Phys. Rep. 895 1-28, (2021), https://doi.org/10.1016/j.physrep.2020.10.001.
  * (6) K. Wang, V. Bheemarasetty, J. Duan, S. Zhou and G. Xiao, Fundamental physics and applications of skyrmions: A review, J. Mag. Mag. Mat. 563, 169905 (2022), https://doi.org/10.1016/j.jmmm.2022.169905.
  * (7) Sheng Li, Xuewen Wang, Theo Rasing, Magnetic skyrmions: Basic properties and potential applications, Interdisciplinary Materials 563, 169905 (2022), https://doi.org/10.1002/idm2.12072.
  * (8) M. Uchida, Y. Onose,Y. Matsui and Y. Tokura, Real-Space Observation of Helical Spin Order, Science 311, pp.359-361 (2006), https://doi.org/10.1126/science.11206.
  * (9) S. M$\ddot{\rm u}$hlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, P. B$\ddot{\rm o}$ni, Skyrmion Lattice in a Chiral Magnet, Science, 123, 915 (2009), https://doi.org/10.1126/science.116676.
  * (10) X. 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), https://doi.org/10.1038/nature09124.
  * (11) A. N. Bogdanov, and U. K. R${\ddot{\rm o}}{\rm\beta}$ler, Chiral Symmetry Breaking in Magnetic Thin Films and Multilayers, Phys. Rev. Lett., 87, 037203 (2001), https://doi.org/10.1103/PhysRevLett.87.037203.
  * (12) A. B. Butenko, A. A. Leonov, U. K. R${\ddot{\rm o}}{\rm\beta}$ler, and A. N. Bogdanov, Stabilization of skyrmion textures by uniaxial distortions in noncentrosymmetric cubic helimagnets, Phys. Rev. B 82, 052403 (2010), https://doi.org/10.1103/PhysRevB.82.052403.
  * (13) S. Rohart and A. Thiaville, Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction, Phys. Rev. B 88, 184422 (2013), https://doi.org/10.1103/PhysRevB.88.184422.
  * (14) H. Du, R. Che, L. Kong X. Zhao, C. Jin, C. Wang, J. Yang, W. Ning, R. Li, C. Jin, X. Chen, J. Zang, Y. Zhang, and M. Tian, Edge-mediated skyrmion chain and its collective dynamics in a confined geometry, Nature Comm. 6, 8504 (2015), https://doi.org/10.1038/ncomms9504.
  * (15) T. Matsumoto, Y.-G. So, Y. Kohno, Y. Ikuhara, and N. Shibata, Stable Magnetic Skyrmion States at Room Temperature Confined to Corrals of Artificial Surface Pits Fabricated by a Focused Electron Beam, Nano Lett. 18, 754-762 (2018), https://doi.org/10.1021/acs.nanolett.7b03967.
  * (16) C. Jin, Zi-An Li, A. Kov${\rm\acute{a}}$cs, J. Caron, F. Zheng, F. N. Rybakov, N. S. Kiselev, H. Du, S. Bl${\rm\ddot{u}}$gel, M. Tian, Y. Z., M. Farle, and Rafal E. Dunin-Borkowski, Control of morphology and formation of highly geometrically confined magnetic skyrmions, Nature Comm. 8, 15569 (2017), https://doi.org/10.1038/ncomms15569.
  * (17) Z. Hou, Q. Zhang, G. Xu, S. Zhang, C. Gong, B. Ding, H. Li, F. Xu, Y. Yao, E. Liu, G. Wu, X. Zhang, and W. Wa, Manipulating the Topology of Nanoscale Skyrmion Bubbles by Spatially Geometric Confinement, ACS Nano. 13, 922-929 (2019), https://doi.org/10.1021/acsnano.8b09689.
  * (18) Pin Ho, Anthony K.C. Tan, S. Goolaup, A.L. Gonzalez Oyarce, M. Raju, L.S. Huang, Anjan Soumyanarayanan, and C. Panagopoulos, Geometrically Tailored Skyrmions at Zero Magnetic Field in Multilayered Nanostructures, Phys. Rev. Appl. 11, 024064 (2019), https://doi.org/10.1103/PhysRevApplied.11.024064.
  * (19) S. Seki, Y. Okamura, K. Shibata, R. Takagi, N. D. Khanh, F. Kagawa, T. Arima, and Y. Tokura, Stabilization of magnetic skyrmions by uniaxial tensile strain, Phys.Rev. B 96, 220404(R) (2017), https://doi.org/10.1103/PhysRevB.96.220404.
  * (20) Y. Wang, L. Wang, J. Xia, Z. Lai, G. T. X. Zhang, Z. Hou, X. Gao, W. Mi, C. Feng, M. Zeng, G. Zhou, G. Yu, G. Wu, Y. Zhou, W. Wang, X. Zhang, and J. Liu, Electric-field-driven non-volatile multi-state switching of individual skyrmions in a multiferroic heterostructure, Nature Comm. 11, 3577 (2020), https://doi.org/10.1038/s41467-020-17354-7.
  * (21) H. Li, C. A. Akosa, P. Yan, Y. Wang, and Z. Cheng, Stabilization of Skyrmions in a Nanodisk Without an External Magnetic Field , Phys. Rev. Applied. 13, 034046 (2020), https://doi.org/10.1038/s41467-020-17354-7.
  * (22) A. O. Leonov, U. K. R$\ddot{\rm o}$ßler, and M. Mostovoy, Target-skyrmions and skyrmion clusters in nanowires of chiral magnets, EPJ Web of Conferences, 75, 05002 (2014), https://doi.org/10.1051/jepjconf/20147505002.
  * (23) H. Du, W. Ning, M. Tian and Y. Zhang, Field-driven evolution of chiral spin textures in a thin helimagnet nanodisk, Phys. Rev. B 87, 014401 (2013);, https://doi.org/10.1103/PhysRevB.87.014401.
  * (24) F, Zheng, H. Li, S. Wang, D. Song, C. Jin, W. Wei, A. Kov$\acute{\rm a}$cs, J. Zang, M. Tian, Y. Zhang, H. Du and R. E. Dunin-Borkowski, Direct Imaging of a Zero-Field Target Skyrmion and Its Polarity Switch in a Chiral Magnetic Nanodisk, Phys. Rev. Lett. 119, 197205 (2017), https://doi.org/10.1103/PhysRevLett.119.197205.
  * (25) N. Kent, R. Streubel, C. -H. Lambert, A. Ceballos, S.-G. Je, S. Dhuey, Mi-Y. Im, F. B$\ddot{\rm u}$uttner, F. Hellman, S. Salahuddin and P. Fischer, Generation and stability of structurally imprinted target skyrmions in magnetic multilayers, Appl. Phys. Lett. 115, 112404 (2019), https://doi.org/10.1063/1.5099991.
  * (26) M. Beg, R. Carey, W. Wang, D. Cort$\acute{\rm e}$s-Ortu$\tilde{\rm n}$o, M. Vousden, M.-A. Bisotti, M. Albert, D. Chernyshenko, O. Hovorka, R. L. Stamps and H. Fangohr, Ground state search, hysteretic behaviour, and reversal mechanism of skyrmionic textures in confined helimagnetic nanostructures, Scientific Reports 5, 17137 (2015), https://doi.org/10.1038/srep17137.
  * (27) M. Beg, M. Albert, M. -A. Bisotti, D. Cort$\acute{\rm e}$s-Ortu$\tilde{\rm n}$o, W. Wang, R. Carey, M. Vousden, O. Hovorka, C. Ciccarelli, C. S. Spencer, C. H. Marrows and H. Fangohr, Dynamics of skyrmionic states in confined helimagnetic nanostructures, Phys. Rev. B 95, 014433 (2017), https://doi.org/10.1103/PhysRevB.95.014433.
  * (28) N. Mehmood, J. Wang, C. Zhang, Z. Zeng, J. Wang, and Q. Liu, Magnetic skyrmion shape manipulation by perpendicular magnetic anisotropy excitation within geometrically confined nanostructures, J. Mag. Mag. Mat. 545, 168775 (2022), https://doi.org/10.1016/j.jmmm.2021.168775.
  * (29) K. Niitsu, Y. Liu, A. C. Booth, X. Yu, N. Mathur, M. J. Stolt, D. Shindo, S. Jin, J. Zang, N. Nagaosa, and Y. Tokura, Geometrically stabilized skyrmionic vortex in FeGe tetrahedral nanoparticles, Nature Materials 21, 305-310 (2022), https://doi.org/10.1038/s41563-021-01186-x.
  * (30) G. Diguet, B. Ducharne, S. El Hog, F. Kato, H. Koibuchi, T. Uchimoto, and H. T. Diep, Finsler Geometry Modeling and Monte Carlo Study on Geometrically Confined skyrmions in Nanodot, Proc. ICFD2023, http://arxiv.org/abs/2308.16199.
  * (31) R. Friedberg and H.-C. Ren, Field theory on a computationally constructed random lattice, Nucl. Phys. B 235, 310-320 (1984), https://doi.org/10.1016/0550-3213(84)90501-7.
  * (32) N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953), https://doi.org/10.1063/1.1699114.
  * (33) D.P. Landau, Finite-size behavior of the Ising square lattice, Phys. Rev. B 13, 2997 (1976), https://doi.org/10.1103/PhysRevB.13.2997.
  * (34) G. Diguet, B. Ducharne, S. El Hog, F. Kato, H. Koibuchi, T. Uchimoto, and H.T. Diep, Monte Carlo studies of skyrmion stabilization under geometric confinement and uniaxial strain, J. Mag. Mag. Mat. 579, 170819 (2023), https://doi.org/10.1016/j.jmmm.2023.170819.
  * (35) D. Andrienko, Introduction to liquid crystals, J. Mol. Liq., 267, pp. 520-541 (2018), https://doi.org/10.1016/j.molliq.2018.01.175.
  * (36) T. Shinjo T. Okuno R. Hassdorf, K. Shigeto and T. Ono, Magnetic Vortex Core Observation in Circular Dots of Permalloy, Science 289, 930 (2000), https://doi.org/10.1126/science.289.5481.930.
  * (37) T. Okuno, K. Shigeto, T. Ono, K. Mibu and T. Shinjo, MFM study of magnetic vortex cores in circular permalloy dots: behavior in external field, J. Mag. Mag. Mat. 240, pp. 1-6 (2002), https://doi.org/10.1016/S0304-8853(01)00708-9.