# Particle-continuum-medium duality of skyrmions

X. R. Wang<EMAIL_ADDRESS>Physics Department, The Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong  HKUST Shenzhen
Research Institute, Shenzhen 518057, China  William Mong Institute of Nano
Science and Technology, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  X. C. Hu  Physics Department, The
Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong  HKUST Shenzhen Research Institute, Shenzhen 518057, China

###### Abstract

Topological solitons are crucial to many branches of physics, such as models
of fundamental particles in quantum field theory, information carriers in
nonlinear optics, and elementary entities in quantum and classical
computations. Chiral magnetic materials are a fertile ground for studying
solitons. In the past a few years, a huge number of all kinds of topologically
protected localized magnetic solitons have been found. The number is so large,
and a proper organization and classification is necessary for their future
developments. Here we show that many topological magnetic solitons can be
understood from the duality of particle and elastic continuum-medium nature of
skyrmions. In contrast to the common belief that a skyrmion is an elementary
particle that is indivisible, skyrmions behave like both particle and
continuum media that can be tore apart to bury other objects, reminiscing
particle-wave duality in quantum mechanics. Skyrmions, like indivisible
particles, can be building blocks for cascade skyrmion bags and target
skyrmions. They can also act as bags and glues to hold one or more skyrmions
together. The principles and rules for stable composite skyrmions are
explained and presented, revealing their rich and interesting physics.

##  I Introduction

Magnetic skyrmions have attracted much attention in recent years for their
academic interest and potential applications in information technology [1, 2,
3, 4, 5, 7, 6, 8]. Various aspects of magnetic skyrmions have been extensively
and intensively studied [9, 12, 10, 11, 13, 14, 15, 16, 14, 18, 19, 33, 34,
20, 23, 24, 25, 26, 27, 28, 29, 30, 31, 17, 32, 21, 35, 36, 22, 37, 38, 39,
40, 41, 42, 43], including characterization, generation, and manipulations.
Recently, it was realized that only skyrmions with positive formation energy
are intrinsically circular [44, 45, 46, 47]. The natural morphology of
skyrmions of negative formation energy are stripes of well-defined width. In
this case, condensed stripy skyrmions, ranging from irregular maze to
periodically arranged helical states, is the preferred thermodynamic
equilibrium state [47]. Randomly or orderly arranged condensed stripe
skyrmions transform into skyrmion crystals (SkXs) smoothly and continuously as
stripe width increases via material parameter engineering [48] or as the
skyrmion density increases through thermal agitations and/or the assistances
of external forces [46, 47].

With all these advances in skyrmion physics, there are still many unsolved
problems. For example, other than different elementary skyrmions mentioned
above, in the past eight years, many localized magnetic structures,
topologically different from the elementary skyrmions and with arbitrary
integer skyrmion numbers, were found [49, 50, 51, 52, 53, 54, 55, 56, 57, 58,
59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]. Skyrmionium, a spin
texture of one skyrmion inside another larger skyrmion with zero net skyrmion
number, have been observed both numerically and experimentally [49, 50, 51,
52, 53, 54, 55, 56, 57, 58]. A skyrmionium can also be embedded inside another
even larger skyrmion to form a target skyrmion of 3 layers with net skyrmion
number 1. In fact, this cascade process can continue to form target skyrmions
of any cascade levels with skyrmion number 0 or 1 for even and odd numbers of
layers, respectively. Target skyrmions cannot be continuously transformed into
a single ferromagnetic domain structure or an elementary skyrmion. Also, so-
called skyrmion bags, spin structures of arbitrary number of skyrmions inside
another larger skyrmion, have also been observed [49, 59, 60, 61, 62, 63, 64].
Their skyrmion numbers can be arbitrary integers $Q$ and they cannot be
continuously deformed into $Q$ spatially separated elementary skyrmions.
Similar to the target skyrmions, one can put one or many skyrmion bags into
another even bigger skyrmion to form a cascade skyrmion bag of 3 layers. The
process can also continue to have various cascade skyrmion bags of arbitrary
layers. Although composite skyrmions have been found in various simulations
and experiments, there is no clear understanding of their existence
conditions. They were found through trial and error or hunch so far.

The increasing number of newly discovered topologically non-trivial localized
spin structures in chiral magnetic films reminisces the discovery of a huge
number of elementary particles in the first half of last century that resulted
in the quark model of elementary particles. It calls for a “quark model” to
properly organize these newly discovered topological magnetic structures. In
this paper, we consider a perpendicularly magnetized chiral magnetic film with
the Dzyaloshinskii–Moriya interaction (DMI) and normal ferromagnetic exchange
interaction characterized by parameter $D$ and $A$, respectively, as well as
magnetic anisotropy measured by $K$. We show that the stability and properties
of various types of composite skyrmions is very sensitive to a parameter
defined as $\kappa\equiv(\pi D)^{2}/(16AK)$. Near $\kappa=1$ that separates
isolated circular skyrmions from condensed stripe skyrmions [45], the
skyrmions show strong particle-continuum-medium duality. In contrast to the
general belief that a skyrmion is a fundamental object indivisible. A skyrmion
behaves sometimes like a continuum medium that can be tore apart to embed one
or more than one skyrmions or composited skyrmions although it often appears
as a particle. In another word, one or many circular and stripe skyrmions can
be hold or be glued by another larger skyrmion, which acts as a glue or a bag,
to form a skyrmionium or skyrmion bag/cluster. The embedding can continue
layer by layer to form all kinds of target skyrmions and cascade skyrmion
bags.

Specifically, following interesting results were obtained in the absence of
the external magnetic field for $\kappa\leq 1$ where isolated circular
skyrmions are metastable and $\kappa>1$ where condensed stripe skyrmions are
stable. 1) Topologies of stable static magnetic textures are fully determined
by $\kappa$, and $A/D$ defines their length scale. 2) The maximal number of
cascade layers in target skyrmions for $\kappa\leq 1$ increases monotonically
from 1 at $\kappa=0.7$, to 2 at $\kappa=0.8$, and to 4 at $\kappa=0.94$, and
eventually diverges at $\kappa=1$ for an infinite large system. 3) An unstable
target skyrmion or skyrmion bag for $\kappa<1$ can be stabilized by inserting
enough skyrmions inside the innermost bag and/or the next innermost bag. The
maximal number of skyrmions needed to stabilize a skyrmion bag decreases from
a larger value at small $\kappa$, say 4 at $\kappa=0.7$, to the lowest
possible 1 when $\kappa$ approaches 1. This number increases also as the
number of cascade layers increases. Furthermore, skyrmion size inside the
innermost bag increases with the number of skyrmions inside the bag, and
approaches the size of isolated elementary skyrmions. 4) For $\kappa>1$, the
number of target stripe skyrmions and cascade stripe skyrmion bags can be any
number as long as the average space-size occupied by each skyrmion is larger
than twice of elementary stripe width. Furthermore, stripe width and stripe
spin profile of any target stripe skyrmions and cascade stripe skyrmion bags
are the same as those of elementary stripe skyrmions. 5) When the average
skyrmion-skyrmion distance inside a stripe skyrmion bag is comparable to the
elementary stripe width for $\kappa>1$, stripe skyrmions in the innermost
layer become disk-like objects and form an SkX. The paper is organized as
follows. The model and methodology is present in the next section. We prove
there that stable/metastable spin structures are fully determined by $\kappa$
with the fundamental length scale of $4A/(\pi D)$. Section III is the results,
and discussion and conclusion are given in Sec. IV.

##  II Model and Methodology

To demonstrate particle-continuum-medium duality of skyrmions, we consider a
thin chiral magnetic film of thickness $d$ in the xy-plane. The magnetic
energy of a spin structure $\vec{m}$ is

$\displaystyle
E=d\iint\\{A|\nabla\vec{m}|^{2}+D[(\vec{m}\cdot\nabla)m_{z}-m_{z}\nabla\cdot\vec{m}]+$
(1) $\displaystyle K_{\rm u}(1-m_{z}^{2})-\frac{1}{2}\mu_{0}M_{\rm
s}\vec{H}_{\rm d}\cdot\vec{m}+\mu_{0}HM_{\rm s}(1-m_{z})\\}\,\mathrm{d}S,$

where $K_{\rm u}$, $H$, $M_{\rm s}$, $\vec{H}_{\rm d}$ and $\mu_{0}$ are the
magneto-crystalline anisotropy, perpendicular magnetic field, the saturation
magnetization, the demagnetizing field and the vacuum permeability,
respectively. The energy of ferromagnetic state $m_{z}=1$ is chosen as the
energy reference of $E=0$. For an ultra thin film, demagnetization effect can
be included in the effective anisotropy $K=K_{\rm u}-\mu_{0}M_{\rm s}^{2}/2$.
This is a good approximation when the film thickness $d$ is much smaller than
the exchange length [9].

One can recast Eq. (1) to reveal how model parameters affect stable/metastable
magnetic structures,

$\displaystyle E$
$\displaystyle=d\int\int\\{A|\nabla\mathbf{m}|^{2}+D[m_{z}\nabla\cdot\mathbf{m}-(\mathbf{m}\cdot\nabla)m_{z}]$
(2) $\displaystyle\quad+K(1-m_{z}^{2})+\mu_{0}HM_{\rm
s}(1-m_{z})\\}\mathrm{d}x\mathrm{d}y$ $\displaystyle=d\kappa
K\int\int\\{|L\nabla\mathbf{m}|^{2}+\frac{4L}{\pi}[m_{z}\nabla\cdot\mathbf{m}-(\mathbf{m}\cdot\nabla)m_{z}]$
$\displaystyle\quad+\frac{1}{\kappa}(1-m_{z}^{2})+\frac{2}{\kappa^{\prime}}(1-m_{z})\\}\mathrm{d}x\mathrm{d}y$
$\displaystyle=dA\int\int\\{|\nabla\mathbf{m}|^{2}+\frac{4}{\pi}[m_{z}\nabla\cdot\mathbf{m}-(\mathbf{m}\cdot\nabla)m_{z}]$
$\displaystyle\quad+\frac{1}{\kappa}(1-m_{z}^{2})+\frac{2}{\kappa^{\prime}}(1-m_{z})\\}\mathrm{d}x\mathrm{d}y,$

where $\kappa=\frac{\pi^{2}D^{2}}{16AK}$,
$\kappa^{\prime}=\frac{\pi^{2}D^{2}}{8A\mu_{0}M_{s}H}$, and $L=4A/(\pi D)$.
The stable/metastable spin structures minimize the last integral of above
equation [47] where $x$ and $y$ is in the units of $L$. In the absence of
external magnetic field ($1/\kappa^{\prime}=0$) considered here, all
stable/metastable magnetic structures satisfy following dimensionless
equation,

$\nabla^{2}\vec{m}+\frac{4}{\pi}\left[(\nabla\cdot\vec{m})\hat{z}-\nabla
m_{z}\right]+\frac{1}{\kappa}m_{z}\hat{z}=0.$ (3)

In another word, stable/metastable magnetic structures are fully determined by
$\kappa$, not by $A$, $D$, $M_{s}$, and $K$ separately. $L$ is the length
scale. This fact greatly simplifies our study of skyrmion static properties,
and all stable magnetic structures in general. It allows us to compare
skyrmions in different chiral magnets in terms of $\kappa$. Varying $\kappa$
over all possible values corresponds to studying all possible magnetic films.
Isolated circular skyrmions are metastable and their energies are $E=8\pi
Ad\sqrt{1-\kappa}$ when $H=0$ and $\kappa\leq 1$ [9]. Obviously, the formation
energy of an isolated skyrmion goes to zero as $\kappa$ approaches 1. Earlier
studies [44, 45, 46, 47] showed that $\kappa=1$ separates isolated circular
skyrmions from condensed stripe skyrmions. The natural skyrmion morphologies
are stripes in order to take the advantages of negative skyrmion formation
energy when $\kappa>1$ [46].

Spin dynamics in a magnetic field is governed by the Landau-Lifshitz-Gilbert
(LLG) equation,

$\frac{\partial\vec{m}}{\partial t}=-\gamma\vec{m}\times\vec{H}_{\rm
eff}+\alpha\vec{m}\times\frac{\partial\vec{m}}{\partial t},$ (4)

where $\gamma$ and $\alpha$ are respectively gyromagnetic ratio and Gilbert
damping constant. $\vec{H}_{\rm eff}=\frac{2A}{\mu_{0}M_{\rm
s}}\nabla^{2}\vec{m}+\frac{2K_{\rm u}}{\mu_{0}M_{\rm
s}}m_{z}\hat{z}+H\hat{z}+\vec{H}_{\rm d}+\vec{H}_{\rm DM}+\vec{h}$ is the
effective field which includes the exchange field, the anisotropy field, the
external magnetic field along $\hat{z}$, the demagnetizing field, the DMI
field $\vec{H}_{\rm DM}$, and the thermal fluctuating field $\vec{h}$ (in the
case of finite temperature), respectively [72]. In the absence of an energy
source, the steady state solutions of Eq. (4) are stable/metastable spin
textures. We will use MuMax3 [72] to obtain all kinds of composite skyrmions.
The initial configurations are important in this endeavour. According to our
previous study [44], each nucleation domain in a chiral magnetic film, in
which skyrmions are stable/metastable, develops into a skyrmion. Thus, various
domains embedded inside other domains can develop into various composite
skyrmions, as demonstrated below.

Figure 1: Typical structures and profiles of composite circular skyrmions for
$\kappa=0.9$ (a1)-(a3) and composite stripe skyrmions for $\kappa=4$
(b1)-(b3). One skyrmion inside another one (a1) and (b1) from the same initial
configuration shown in the inset of (a1); two skyrmions inside another larger
one (a2) and (b2) from the initial configuration shown in the inset of (a2);
two composite skyrmions of (a2) or (b2) inside another larger skyrmion (a3)
and (b3) from the initial configuration shown in the inset of (a2). (c) Time
evolution of topological charges for structures (a1) (the black solid line);
(a2) (the red solid line); (a3) (the blue solid line); (b1) (the cyan dash
line); (b2) (the green dash line); and (b3) (the red dash line). $t$ is in the
logarithmic scale. Skyrmion charge density $\rho$ is encoded by colours (the
blue for positive and the red for negative) while the gray-scale encodes
$m_{z}$. (d) Spin profiles $m_{z}(x)$ along the green lines labelled by ⓝ in
(b1)-(b3). The symbols are numerical data and the solid lines are the
theoretical fits to
$\cos\Theta(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=25.10$ nm and $w=7.63$nm.

##  III Results

###  III.1 Target skyrmions and cascade skyrmion bags at $\kappa=0.9$ and
$\kappa=4$

Stable/metastable composite skyrmions are expected to more favourably appear
in chiral magnetic films with $\kappa$ slightly smaller than 1 or $\kappa>1$
because creating an individual skyrmion gains energy for $\kappa>1$ and costs
nearly zero energy near $\kappa=1^{-}$. To demonstrate that this expectation
is indeed true, we choose $A=4\,\rm{pJ\,m^{-1}}$, $D=1\,\rm mJ\,m^{-2}$, and
$M_{s}=0.15\,\rm MA\,m^{-1}$, values in the range of $A=2\sim 16\,\rm pJ\,\rm
m^{-1}$, $D=0.68\sim 4\,\rm mJ\,\rm m^{-2}$, $K_{u}=0.2\sim 2.5\,\rm MJ\,\rm
m^{-3}$ and $M_{s}=0.65\sim 1.1\,\rm MA\,\rm m^{-1}$ for chiral magnets of
$\rm Ir/\rm Fe/\rm Co/\rm Pt$ multilayer and $\rm PdFe/\rm Ir$ bilayers [8,
29, 24, 73]. $K_{u}$ is used to simulate samples with different $\kappa$ since
in reality $K_{u}$ is sensitive to chemical compositions, fabrication process,
structures, and the temperature. $K_{u}=0.185\,\rm MJ\,m^{-3}$ and
$K_{u}=0.0526\,\rm MJ\,m^{-3}$ give $\kappa=0.9$ and $\kappa=4$, respectively.
As explained in the previous section, the stable/metastable structures for the
same $\kappa$ are invariant when they are scaled by $4A/(\pi D)$. Figure 1
shows the typical structures of composite skyrmions in a sample of $400\,$nm
$\times 400\,$nm $\times 0.5\,$nm in size with periodic boundary conditions
along both $x$ and $y$ directions. Throughout this study, skyrmion structures
are encoded in the gray-scale for $m_{z}$ and by colors for skyrmion charge
density $\rho$ defined as
$\rho=\mathbf{m}\cdot(\partial_{x}\mathbf{m}\times\partial_{y}\mathbf{m})/(4\pi)$.

Figure 1(a1) is a stable skyrmionium, starting from an initial configuration
of a $30\,\rm nm\times 30\,\rm nm\times 0.5\,\rm nm$ $(m_{z}=-1)$-domain
inside another $(m_{z}=1)$-domain of $100\,\rm nm\times 100\,\rm nm\times
0.5\,\rm nm$ at the center in the background of $m_{z}=-1$ as shown in the
inset (the red-line square). It forms shortly a stable skyrmionium of skyrmion
number (charge) 0 with an inner core of topological charge -1 and an outer
shell of topological charge 1, separated by the white-color circle on which
$\rho=0$. The topological charge remains 0 shortly after $10\,$ps as shown by
the black solid curve in Fig. 1(c). Figure 1(a2) is a stable composite
topological skyrmion with two skyrmions inside another larger skyrmion,
starting from an initial configuration (shown in the inset) of two $30$ nm
$\times$ $30\,\rm nm\times 0.5\,\rm nm$ $(m_{z}=-1)$-domains inside another
$(m_{z}=1)$-domain of $100\,\rm nm\times 100\,\rm nm\times 0.5\,\rm nm$ in the
background of $m_{z}=-1$. Two skyrmions enclosed by the white-color contour of
$\rho=0$ have skyrmion number -1 each while the total charge outside of the
contour is 1. The red solid curve in Fig. 1(c) is the time evolution of the
total topological charge $Q(t)$, $Q(t)=-1$ shortly after $10\,$ps evolution.
Similarly, Fig. 1(a3) is a stable cascade skyrmion bag of 3 cascade layers:
two composite skyrmions shown in Fig. 1(a2) inside another larger skyrmion.
This cascade skyrmion bag comes from an initial configuration consisting of
two composite domains as that in the inset of Fig. 1(a2) inside another larger
domain. As expected, this composite skyrmion has a skyrmion number of
$4\times(-1)+2\times 1+(-1)=-3$ as shown by the blue curve in Fig. 1(c).
Interestingly, $Q(t)$ exhibits a peak with a positive value larger than 1
before reaching its final negative value and its cause is unknown yet.

In order to investigate whether the above results are also true for stripe
skyrmions when skyrmion formation energy is negative, we repeat the above
simulations for film with $\kappa=4$. Figures 1(b1)-(b3) show three stable
structures of composite stripe skyrmions: stripe skyrmionium of topological
charge 0 or one stripe skyrmion inside another larger stripe skyrmion (b1); a
composite stripe skyrmion bag of topological charge -1 with two stripe
skyrmions inside another larger stripe skyrmions (b2); a composite cascade
stripe skyrmion bag of topological charge -3 with two composite stripe
skyrmion bags as shown in Fig. 1(b2) inside another stripe skyrmions (b3). As
a direct comparison to the case of $\kappa=0.9$, the initial configurations
are exactly the same as their corresponding counterparts in Figs. 1(a1)-(a3).
The dash curves in Fig. 1(c) are the time evolution of topological charge of
different structures. Different from the case of $\kappa<1$, $Q(t)$ reaches
monotonically to its final stable values. Also, different from the case of
$\kappa<1$ where the skyrmion sizes in different layers are different, the
stripe width in the current case are the same. Even more interestingly, spin
profile are the same as that for elementary stripe skyrmion [45],
$\Theta(x)=2\arctan[\frac{\sinh(L/w)}{\sinh(x/w)}]$ for $m_{z}>0$ and
$\Theta(x)=2\arctan[\frac{\sinh(x/w)}{\sinh(L/w)}]$ for $m_{z}<0$,
respectively, with $|x|\leq L/2$. $\Theta$ is the polar angle of the
magnetization at $x$ and $x=0$ is the center of a stripe where $m_{z}=\pm 1$.
The symbols in Fig. 1(d) are numerical data along the green line labelled by
green ⓝ in Figs. 1(b1)-(b3). The solid curve is the fit of spin profiles to
$\cos\Theta(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=25.10\,$nm and $w=7.63\,$nm. All data from different
stripes falling onto the same curve demonstrates that stripes, building blocks
of the pattern, are identical.

###  III.2 Maximal layer number of target skyrmions for $\kappa\leq 1$

Table 1: The value of material parameters $K$, $K_{u}$ and $\kappa$ for 8 different films.  Films | $K_{u}(\,\rm MJ\,m^{-3}\,)$ | $K(\,\rm MJ\,m^{-3}\,)$ | $\kappa$
---|---|---|---
1 | 0.2344 | 0.2203 | 0.7
2 | 0.2197 | 0.2056 | 0.75
3 | 0.2068 | 0.1927 | 0.8
4 | 0.1955 | 0.1814 | 0.85
5 | 0.1854 | 0.1713 | 0.9
6 | 0.1817 | 0.1676 | 0.92
7 | 0.1781 | 0.1640 | 0.94
8 | 0.1747 | 0.1606 | 0.96
Figure 2: (a) $\kappa$-dependence of the maximal number of layers of target
skyrmions. The inset shows how skyrmion number of unstable composite spin
structure varies with time $t$ when $\kappa=0.7$. The final stable spin
texture is an elementary skyrmion of $Q=1$. (b) The initial configuration of a
cascaded domain-inside-domain structure of $n=11$. (c)-(g) The final
metastable target skyrmions from the initial configuration of (b) for various
$\kappa$: Skyrmionium for $\kappa=0.8$ (c) and target skyrmions of 3 layers
for $\kappa=0.9$ (d), of 4 layers for $\kappa=0.92$ (e), of 5 layers for
$\kappa=0.94$ (f), and of 8 layers for $\kappa=0.96$ (f). Color bar denotes
skyrmion charge density and the gray-scale encodes $m_{z}$.  Figure 3: Various
stable cascade skyrmion bags of $n=2,3,4,$ and 5 cascade layers for
$\kappa=0.7,0.8,$ and 0.9. (a1)-(a4) Stable cascade skyrmion bags of $n=2,$ 3,
4, and 5 layers with minimal number of $N_{c}=4$ (a1), 18 (a2), 38 (a3), and
66 (a4) skyrmions in the innermost bag for $\kappa=0.7$. (a5) A stable
skyrmion bag of $N=64\gg N_{c}$ and $n=2$ for $\kappa=0.7$. (b1)-(b4) Stable
cascade skyrmion bags of $n=3$, 4, 5, and 6 layers with minimal number of
$N_{c}=4$ (b1), 10 (b2), 17 (b3), and 28 (b4) skyrmions in the innermost bag
for $\kappa=0.8$. (b5) A stable cascade skyrmion bag of $N=64\gg N_{c}$ and
$n=3$ for $\kappa=0.8$. (c1)-(c4) Stable cascade skyrmion bags of $n=4,$ 5, 6
and 7 layers with minimal number of $N_{c}=2$ (c1), 3 (c2), 5 (c3), and 7 (c4)
in the innermost bag for $\kappa=0.9$. The insets are the initial
configurations. (c5) A stable cascade skyrmion bag of $n=4$ and $N=64\gg
N_{c}$ skyrmions in the innermost bag for $\kappa=0.9$. (d) $N_{c}$ vs. n.
Dots are numerical data for $\kappa=0.7$ (the black); 0.8 (the red) and 0.9
(the blue). The solid lines are the fits of data to
$N_{c}=a(\kappa)(n-1)^{2}$, with $a(0.7)=4.2$, $a(0.8)=1.1$, and $a(0.9)=0.2$.
(e) Skyrmion size in the innermost bag as a function of the number of
skyrmions $N$ in the innermost bag for $\kappa=0.7$ and $n=2$ (the black full
line); $\kappa=0.7$ and $n=3$ (the black dash line; $\kappa=0.7$ and $n=4$
(the black point line); $\kappa=0.8$ and $n=3$ (the red full line);
$\kappa=0.8$ and $n=4$ (the red dash line); $\kappa=0.8$ and $n=5$ (the red
point line); $\kappa=0.9$ and $n=4$ (the blue full line); $\kappa=0.9$ and
$n=5$ (the blue dash line); $\kappa=0.9$ and $n=6$ (the blue point line) and
elementary skyrmion size for $\kappa=0.7$ (the black short dash line);
$\kappa=0.8$ (the red short dash line); $\kappa=0.9$ (the blue short dash
line). Color bar denotes skyrmion charge density and the gray-scale encodes
$m_{z}$.

When $\kappa<1$, the skyrmion formation energy of an elementary skyrmion is
$E=8\pi A\sqrt{1-\kappa}$ [9], and skyrmion size and skyrmion wall width are
$R=\pi D/(4K\sqrt{1-\kappa})$ and $w=\pi D/(4K)$, respectively. The energy
cost is mainly from the skyrmion wall region. The size and wall regions in a
composite skyrmion is obviously larger than that in an elementary skyrmion,
resulting in an increase of energy-cost. Thus, stability of various composite
skyrmions should be sensitive to how far $\kappa$ is away from 1 at which the
energy cost for creating a wall is zero. Although an isolated skyrmion is
expected to be metastable when $\kappa\leq 1$ as long as the skyrmion size is
much bigger than the lattice constant $a$, $R>a$ or $\kappa>1-[\pi
D/(4Ka)]^{2}$, in order to validate the continuum description of the films by
Eqs. (1), (2), and (3), the maximal number of cascade layers of target
skyrmions should change with $\kappa$. In order to test this hypothesis, we
consider 8 samples with various $K_{u}$ listed in Tab. 1 while all other model
parameters are unchanged. The corresponding $\kappa$ values are 0.7, 0.75,
0.8, 0.9 0.92, 0.94 and 0.96. We started with a cascade domain-inside-domain
structure as shown in Fig. 2(b), consisting of 6 white $(m_{z}=1)$-domains and
6 gray $(m_{z}=-1)$-domains. They were segregated by 11 concentric squares of
side lengths of
$750,\,700,\,650,\,600,\,550,\,500,\,450,\,400,\,350,\,300,\,250\,$nm at the
film center. The sample size is of $2400\,\rm nm\times 2400\,\rm nm\times
0.5\,\rm nm$, large enough so that the final stable spin structures do not
depend on the sample size. The mesh size is $1\,$nm $\times 1\,$nm $\times
0.5\,$nm in the simulations. For $\kappa=0.7$, the initial configuration
quickly evolves into an unstable shell structure of skyrmions number $Q=1$ in
less than $10\,$ps as shown in the inset of Fig. 2(a) in the first jump around
$t=0$ from $Q=0$ to 1 in $Q(t)$-$t$ curve. Simulations show that the shell
structure shrinks its size, and the innermost spin structure disappears one by
one from the center. As revealed in $Q(t)$-$t$ curve, $Q$ jumps back and forth
between 1 and 0 each time when one shell disappears. The final metastable spin
structure is an elementary skyrmion of $Q=1$, indicating that target skyrmions
and skyrmionium are not stable when $\kappa=0.7$.

For $\kappa=0.8$, the final stable spin structure is a skyrmionium as shown in
Fig. 2(c). Figures 2(d)-(g) are the final stable target skyrmions for
$\kappa=0.9$, 0.92, 0.94, and 0.96 with maximal layer number 3 (d), 4 (e), 5
(f), and 8 (g), respectively. The $\kappa$-dependence of the maximal number of
cascade layers in the final stable/metastable target skyrmions is plotted in
Fig. 2(a), showing the number increases to a very large value as $\kappa$
approaches 1.

###  III.3 Stabilization and properties of cascade skyrmion bags for
$\kappa\leq 1$

Figure 4: Various composite skyrmion of $\kappa=0.7$. (a1)-(a3) Stable cascade
skyrmions of $N=7$ and minimal $N_{1,c}=32$ (a1); $N=10$ and $N_{1,c}=18$
(a2); $N=16$ and $N_{1,c}=2$ (a3). (a4) A stable cascade skyrmion of $N=7$ and
$N_{1}=80$ which is much larger than the critical value of 32. (b1)-(b3)
Stable cascade skyrmion bags of $n=4,$ 5, and 6 layers and $N=16$ and
$N_{1,c}=44,104,$ and 200. (b4) A stable cascade skyrmions of $n=4$, $N=10$,
and $N_{1}=128$ which is much larger than the critical value of 44. (c1)-(c3)
Four stable skyrmion superstructures of $n=3,$ 4, and 5 layers with minimal 2
(c1), 7 (c2), 12 (c3) skyrmion bags in the innermost layers. Each skyrmion bag
contains $N=16$ elementary skyrmions. (c4) A stable skyrmion superstructure of
3 layers with 25 skyrmion bags of $N=16$. Color bar denotes skyrmion charge
density and the gray-scale encodes $m_{z}$.

Figure 2 shows the maximal number of layers in stable target skyrmions being 1
for $\kappa=0.7$, i.e. even a skyrmionium of $Q=0$ is not stable. The energy
of a skyrmion bag with many skyrmions in the innermost layer is higher than
that of a skyrmionium. Thus, one may expect that a spin bag is more unstable
than a skyrmionium for $\kappa=0.7$. Surprisingly and counter-intuitively,
skyrmion bags with $N$ skyrmions inside another larger skyrmion are stable as
long as $N$ is larger than a critical value of $N_{c}=4$ as shown in Fig.
3(a1) for $N=N_{c}=4$ and in Fig. 3(a5) for $N=64$. Their total net
topological charges are respectively $Q=-3$ and $Q=-63$ obtained from the
initial configurations of 4 and 64 $(m_{z}=-1)$-domains inside a
$(m_{z}=1)$-domain as shown in the insets. Similarly, the cascade skyrmion
bags of $n$ layers can also be stabilized by putting $N$ skyrmions in the
innermost bag when $N$ is larger than a critical value $N_{c}(n)$ which
depends on layer number of cascade skyrmion bag. Figures 3(a2)-(a4) are the
stable cascade skyrmion bags of $n=3,$ 4 and 5 with corresponding $N_{c}=18$,
38 and 66 when $\kappa=0.7$ respectively. These cascade skyrmion bags are
obtained from the initial configurations shown in the insets. The net
topological charges of these cascade skyrmion bags are $Q=-18$, -37 and -66.
For $\kappa=0.8$, Figs. 3(b1)-(b4) show cascade skyrmion bags of $n=3$, 4, 5
and 6 with corresponding $N_{c}=4$, 10, 17 and 28. Figure 3(b5) shows a stable
cascade skyrmion of 3 layers with $N=64$ which is larger than $N_{c}=4$. For
$\kappa=0.9$, Figs. 3(c1)-(c4) show cascade skyrmion bags of $n=4$, 5, 6 and 7
with corresponding $N_{c}=2$, 3, 5 and 7. Figure 3(c5) shows a stable cascade
skyrmion of 3 layers with $N=64$ which is larger than $N_{c}=2$.

Interestingly, we find that $N_{c}$ is proportion to $(n-1)^{2}$ for all our
$\kappa$’s. Figure 3(d) shows the relation between $N_{c}$ and $n$ for
$\kappa=0.7$ (black dots); 0.8(red dots) and 0.9 (blue dots). Numerical data
(symbols) fit well to $N_{c}=4.2(n-1)^{2}$ (the black line);
$N_{c}=1.1(n-1)^{2}$ (the red line); $N_{c}=0.2(n-1)^{2}$ (the blue line).
More interestingly, the skyrmion size $R$ of skyrmions in the innermost bag is
smaller than that of elementary skyrmions, $R_{0}=\pi D/(4K\sqrt{1-\kappa})$,
but it increases with $N$ as shown in Fig. 3(e).

In order to explain why a cascade skyrmion bag at a lower $\kappa$, say
$\kappa=0.7$, needs a minimal number of skyrmions in the innermost bag to
stabilize itself, one may suspect that innermost bag size needs to be large
enough so that geometrical restriction would not annihilates the skyrmions in
the bag. If this is correct, one would expect that adding more skyrmions to
the next innermost bag tends to destabilize a cascade skyrmion bag because
skyrmions in the next innermost bag tends to reduce innermost bag size. And
minimal number of skyrmions in the innermost bag of a stable cascade skyrmions
increases with number of skyrmions in the next innermost bag. Strangely, this
seemingly reasonable conjecture turns out to be incorrect. Contrary to this
conjecture, we find the opposite results. As shown in Fig. 4(a1) for a stable
cascade skyrmion bag of 3 layers ($n=3$) at $\kappa=0.7$, the cascade skyrmion
bag have 7 skyrmions in the innermost bag and 32 skyrmions in the next
innermost bag. In Fig. 3(a2) we have seen that the minimal number of skyrmions
in the innermost bag required to stabilize the cascade skyrmion is $N_{c}=18$
if no extra skyrmions are in the next innermost bag. In fact, all cascade
skyrmion bags with 7 skyrmions in the innermost bag are stable as long as the
number of skyrmions in the next innermost bag is larger than 32. Figure 4(a4)
shows a cascade skyrmion bag of 3 layers with 7 skyrmions in the innermost bag
and 80 skyrmions in the next innermost bag. The minimal number of skyrmions in
the next innermost bag decreases as the number of skyrmions in the innermost
bag increases. Equivalently, the minimal number of skyrmions in the innermost
bag decreases with the increase of the number of skyrmions in the next
innermost bag. Figures 4(a2)-(a3) are two stable cascade skyrmion bags of 3
layers with 10 and 16 skyrmions in the innermost bag and corresponding minimal
18 and 2 skyrmions in the next innermost bag, respectively.

Above results are very sensitive to layer number $n$ in cascade skyrmion bags.
If one keeps the number of skyrmions in the innermost bag of a cascade
skyrmion bag below its minimal value such that the composite skyrmion by
itself is unstable, the minimal number of skyrmions in the next innermost bag
$N_{1,c}$(n) required to stabilize the cascade skyrmion bag increases with the
cascade level. Figures 4(b1)-(b3) are stable cascade skyrmion bags of 4, 5,
and 6 and layers with 16 skyrmions in the innermost bag and minimal 44, 104
and 200 skyrmions in the next innermost bag, respectively. Figure 4(b4) is a
stable cascade skyrmion of 4 layers of $N=16$ and $N_{1}=128$ which is much
larger than the corresponding critical value $N_{1,c}=44$.

Instead of inserting elementary skyrmions into one particular cascade layer,
one can also put skyrmion bags into it to obtain a more complicated skyrmion
superstructure. Figures 4(c1)-(c3) are three such stable super skyrmion bags
for $\kappa=0.7$ where skyrmioniun and other target skyrmions are unstable.
Different from Fig. 3(a1) where adding 4 elementary skyrmions into an unstable
skyrmionium can obtain a stable skyrmion bag, Fig. 4(c1) is a stable skyrmion
superstructure when we insert two skyrmion bags of $N=16$ into an unstable
skyrmionium. Of course, this skyrmion superstructure has $n=3$ cascade layers
and the total skyrmion number is $2\times(16-1)+1=31$. 2 is the minimal number
of bags needed to stabilise the structure. All skyrmion superstructures with
more than 2 skyrmion bags in an unstable skyrmionium are stable, as shown in
Fig. 4(c4) with 25 skyrmion bags of $N=16$. The minimal number of bags
required to stabilize an unstable target skyrmion of more layers increases
with the target skyrmion layers. Figures 4(c2) and (c3) are the final stable
skyrmion superstructures when minimal 7 and 12 skyrmion bag of $N=16$ are
inserted into a target skyrmion of 3 (c2) and 4 (c3) layers, respectively. One
notes that different sample sizes are used for these two cases in order to
make sure that the results are not affected by the boundaries.

###  III.4 Resilience of stripe width of composite skyrmions for $\kappa>1$

Figure 5: Stable target stripe skyrmions of 4, 6, and 8 layers for $A/D=4\,$nm
and $\kappa=1.1$ (a1)-(a3) and $\kappa=4$ (b1)-(b3), respectively. The sample
sizes are $600\,\rm nm\times 600\,\rm nm\times 0.5\,nm$ and the average layer
width is larger than the elementary strip width. (a1) and (b1), (a2) and (b2),
and (a3) and (b3) start respectively from the same initial configurations
shown in the inset of (a1) (a2), and (a3). (c)-(d) Spin profile of $m_{z}(x)$
along the green lines and orange line labelled by ⓝ in (a1)-(b3). The symbols
are numerical data and the solid lines are the theoretical fits to
$\cos\Theta(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=32.50\,$nm and $w=5.23\,$nm (c); and $L=25.10\,$nm
and $w=7.63\,$nm (d). Color bar denotes skyrmion charge density and the gray-
scale encodes $m_{z}$.  Figure 6: Target stripe skyrmions with the maximal
numbers of cascade layers $n_{max}$ for $\kappa=1.1$ and 4. The maximal
numbers increase with sample sizes and are slightly larger than the ratio of
sample size and the elementary stripe width. (a1)-(a3) Stable target skyrmions
of $n_{max}=4$, 8 and 11 in $200\,\rm nm\times 200\,\rm nm\times 0.5\,nm$,
$400\,\rm nm\times 400\,\rm nm\times 0.5\,nm$ and $600\,\rm nm\times 600\,\rm
nm\times 0.5\,nm$ samples for $\kappa=1.1$. (b1)-(b3) Stable target skyrmions
of $n_{max}=5$, 8 and 12 in $200\,\rm nm\times 400\,\rm nm\times 0.5\,nm$,
$400\,\rm nm\times 800\,\rm nm\times 0.5\,nm$ and $600\,\rm nm\times 1200\,\rm
nm$ samples for $\kappa=1.1$. (c1)-(c3) Stable target skyrmions of
$n_{max}=5$, 9 and 13 in $200\,\rm nm\times 200\,\rm nm\times 0.5\,nm$,
$400\,\rm nm\times 400\,\rm nm\times 0.5\,nm$ and $600\,\rm nm\times 600\,\rm
nm\times 0.5\,nm$ samples for $\kappa=4$. (d1)-(d3) Stable target skyrmions of
$n_{max}=5$, 10 and 14 in $200\,\rm nm\times 400\,\rm nm\times 0.5\,nm$,
$400\,\rm nm\times 800\,\rm nm\times 0.5\,nm$ and $600\,\rm nm\times 1200\,\rm
nm$ for $\kappa=4$. Insets of (a1)-(d3) are initial configurations. (e)-(f)
Spin profiles of stripes. The symbols are numerical data and the solid lines
are the fits to
$\cos\Theta(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=32.50\,$nm and $w=5.23\,$nm (e); and $L=25.10\,$nm
and $w=7.63\,$nm (f). Color bar denotes skyrmion charge density and the gray-
scale encodes $m_{z}$. Figure 7: Morphology and structures of various stable
cascade stripe skyrmion bags for $\kappa=4$. (a1)-(a3) Stripe skyrmion bags of
$n=2$ with $N=4$ (a1), 16 (a2) and 100 (a3) skyrmions in the innermost bag.
(b1)-(b3) Cascade stripe skyrmion bags of $n=3$ (b1), 4 (b2) and 5 (b3) with
$N=16$ skyrmions in the innermost bag. (c1)-(c3) Cascade stripe skyrmion bags
of $n=3$ with $N=4$ skyrmions in the innermost bag and $4$ skyrmions (c1), 16
skyrmions (c2), and 36 skyrmions (c3) in the next innermost bag. Insets are
the initial configurations of these cascade stripe skyrmion bags. Color bar
denotes skyrmion charge density and the gray-scale encodes $m_{z}$.

Skyrmions are stripes with negative formation energy when $\kappa>1$, and
fundamental stripe width is $a(\kappa)2\pi A/D$, where $a(\kappa)$ quickly and
monotonically decreases from $+\infty$ at $\kappa=1^{+}$ to 1 at $\kappa\gg 1$
[44, 45]. Thus, low energy states of a chiral magnetic film are condensed
stripe skyrmions with the optimal stripe width. Different from the isolated
skyrmions that cost energy to create skyrmion walls, stripes at low skyrmion
density are flexible and irregular in order to fill up the film as shown in
Figs. 1(b1)-(b3). Their widths should be resilient to cascade structures no
matter whether they are in a single stripe skyrmion, or a target stripe
skyrmion, or in a cascade stripe skyrmion bags, as long as the distance
between two cascade layers are larger than the fundamental stripe width of
$a(\kappa)2\pi A/D$. Figure 5 shows stable target stripe skyrmions of $n=4$, 6
and 8 layers in a sample of $600\,\rm nm\times 600\,\rm nm\times 0.5\,\rm nm$
for $\kappa=1.1$ (a1)-(a3) and $\kappa=4$ (b1)-(b3). They are generated from
the initial configurations shown in the insets of (a1)-(a3) for $n=4$, 6, and
8, respectively. The reasons for using these initial configurations are the
same for those in Fig. 1. The average available layer space length of the
target skyrmions are $600\,\rm nm/(2\times 4)=75\,\rm nm$ for Figs. 5(a1) and
5(b1), $50\,\rm nm$ for Figs. 5(a2) and 5(b2), and $37\,\rm nm$ for Figs.
5(a3) and 5(b3), larger than the elementary stripe width of $L=32.5\,\rm nm$
for $\kappa=1.1$ and $L=25.10\,\rm nm$ for $\kappa=4.$ Spin profiles of
stripes in all layers are the same as that for elementary stripes as shown in
Figs. 5(c)-(d). The symbols are numerical data along the green lines labelled
by the green ⓝ in Figs. 5(a1)-(a3) for (c) and along the orange lines labelled
by the orange ⓝ in Figs. 5(b1)-(b3) for (d). The solid curves are the fit to
$\cos\Theta(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L/2-L/2)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L/2-L/2)/w\right]}$
for $L/2<x<L/2+L$ with $L=32.50\,$nm and $w=5.2\,$nm [Fig. 5(c)], and with
$L=25.10\,$nm and $w=7.6\,$nm [Fig. 5(d)].

Different from the case of $\kappa<1$, the number of layers of a stable target
stripe skyrmion can be any value in the thermodynamic limit, and has a maximal
value for a finite sample. The maximal number is slightly larger than the half
of the ratio between the sample size and the elementary stripe width. Figures
6(a1)-(a3) are three stable target stripe skyrmions for $\kappa=1.1$ with
maximal layers of $n_{max}=4$, 8 and 11 when the sample sizes are $200\,\rm
nm\times 200\,\rm nm\times 0.5\,nm$, $400\,\rm nm\times 400\,\rm nm\times
0.5\,nm$ and $600\,\rm nm\times 600\,\rm nm\times 0.5\,nm$ for $\kappa=1.1$,
respectively. They are from the initial configurations shown in the insets of
corresponding figures with more layers than their final stable values. Similar
to the evolutions of those stable circular target skyrmions in Fig. 2, the
extra layers shrink and disappear one by one from the sample center. $n_{max}$
is very close to $200/(2\times 32.5)=3.15$; $400/(2\times 32.5)=6.3$; and
$600/(2\times 32.5)=9.45$, corresponding to average available space size of
$25\,\rm nm$ (a1), $25\,\rm nm$ (a2), and $\approx 27\,\rm nm$ (a3). The
general features do not depend on the sample shape. As shown in Figs.
6(b1)-(b3) for rectangular samples of $200\,\rm nm\times 400\,\rm nm\times
0.5\,nm$, $400\,\rm nm\times 800\,\rm nm\times 0.5\,nm$ and $600\,\rm nm\times
1200\,\rm nm\times 0.5\,nm$ with exactly the same model parameters of
$\kappa=1.1$ as those in Figs. 6(a1)-(a3). $n_{max}=5$ (b1), 8 (b2), and 12
(b3) are mainly determined by the shorter length, and are similar to those in
(a1)-(a3). Remarkably, the skyrmion in the innermost layer becomes a circular
skyrmion of size $21\,$nm (a1), $10.3\,$nm (a2) and $14.5\,$nm (a3), much
smaller than the elementary stripe width of $L=32.5\,$nm. However, stripe
width and spin profile in the outer layers are the same as the elementary
stripes as shown in Fig. 6(e) with the symbols from the numerical simulations
along the green lines labelled by the green ⓝ in Fig. 6(b1)-(b3) and the solid
curves from the fit to
$\cos\Theta(x)=-\frac{\sinh^{2}(L_{1}/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=32.50\,$nm and $w=5.2\,$nm.

To show that this feature is general for all $\kappa>1$, Figs. 6(c1)-(c3) are
three target stripe skyrmions with $n_{max}=5$, 9 and 13 in $200\,\rm nm\times
200\,\rm nm\times 0.5\,nm$, $400\,\rm nm\times 400\,\rm nm\times 0.5\,nm$ and
$600\,\rm nm\times 600\,\rm nm\times 0.5\,nm$ samples for $\kappa=4$,
respectively. Since the elementary stripe width is reduced to $L=25.1\,$nm for
$\kappa=4$, $n_{max}$ increases slightly and is again slightly larger than the
ratio of sample size to $2L$ that is $\simeq 4$ (c1), 8 (c2) and 12 (c3).
Figures 6(d1)-(d3) are the stable target stripe skyrmions with the maximal
layers of $n_{max}=5$, 10 and 14 for rectangular samples of $200\,\rm nm\times
400\,\rm nm\times 0.5\,nm$, $400\,\rm nm\times 800\,\rm nm\times 0.5\,nm$ and
$600\,\rm nm\times 1200\,\rm nm\times 0.5\,nm$, respectively, with exactly the
same model parameters as those in Figs. 6(c1)-(c3). Again, the skyrmion in the
innermost layer is a circular object of $5.73\,$nm (d1), $3.76\,$nm (d2) and
$6.44\,$nm (d3) in diameter that is smaller than the elementary stripe width
of $L=25.1\,$nm. Stripe width and spin profile in the outer layers are the
same as the elementary stripes as shown in Fig. 6(f) with the symbols from the
numerical simulations along the orange lines labelled by the orange ⓝ in Figs.
6(d1)-(d3) and the solid curves from the fit to
$\cos\Theta(x)=-\frac{\sinh^{2}(L_{1}/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$\cos\Theta(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}\left[(x-L)/w\right]}{\sinh^{2}(L/2w)+\sinh^{2}\left[(x-L)/w\right]}$
for $L/2<x<3L/2$ with $L=25.1\,$nm and $w=7.63\,$nm.

Figures 5 and 6 show clearly that stripes in all kinds of composite stripe
skyrmions keep the properties of elementary stripe skyrmions, such as their
widths and spin profiles, at low skyrmion density such that the distance
between two neighbouring layers is larger than the elementary width.
Naturally, one interesting question is what will happen to the structure and
morphology of composite skyrmions when more skyrmions are added into one
particular layer or several layers. From the fact that a collection of
elementary stripe skyrmions form a SkX when the skyrmion-skyrmion distance is
comparable to elementary stripe width [46, 47], one expects that highly
irregularly arranged flexible stripes of well-defined width as shown in Fig.
7(a1) deform their shapes to circular objects when the mean skyrmion-skyrmion
distance inside one particular bag is reduced, and they eventually form SkXs
in the layer(s) when the distance is around $L=a(\kappa)2\pi A/D$. Figure 7
shows morphology change of cascade stripe skyrmions bags of various layers $n$
and $\kappa=4$ when the number of skyrmions in the innermost and next
innermost bags increase from a low skyrmion density to a higher skyrmion
density in a sample of $400\,\rm nm\times 400\,\rm nm\times 0.5\,nm$. Figures
7(a1)-(a3) are the stable cascade stripe skyrmion bags of $n=2$ with $N=4$
(a1), 16 (a2), and 100 (a3) syrmions in the innermost bag. Skyrmions in the
innermost bag are highly irregular stripes when the skyrmion density is low,
as shown in Fig. 7(a1). It becomes a mixture of 4 circular skyrmions and 12
short straight stripe skyrmions when the average skyrmion-skyrmion distance is
around $\sqrt{(400-2\times 25.10)^{2}/16}/2\approx 44\,\rm nm\simeq 1.8L$ as
shown in Fig. 7(a2). 100 skyrmions in the innermost bag for a nice triangular
crystal as shown in Fig. 7(a3) when the average skyrmion-skyrmion distance is
around $0.7L$. Figures 7(b1)-(b3) are the stable cascade stripe skyrmion bags
of $n=3$ (b1), 4 (b2) and 5 (b3) with $N=16$ skyrmions in the innermost bag.
The bag size decreases as $n$ increases such that the effective skyrmion
density in the bag increases. The 16 skyrmions form a mixture of circular and
stripe skyrmions (b1), a SkX-like structure (b2), and an almost perfect SkX
(b3). Figures 7(c1)-(c3) are cascade stripe skyrmion bags of $n=3$ with $N=4$
skyrmions in the innermost bag and 4 skyrmions (c1), 16 skyrmions (c2), and 36
skyrmions (c3) in the next innermost bag. When skyrmion densities in both bags
are low, 4 skyrmions in both bags are irregular stripes as shown in (c1).
Skyrmions in both bags are circular when the average skyrmion-skyrmion
distance is around elementary stripe width as shown in (c2). Skyrmions in both
bags form triangular lattice structures when the effective skyrmion densities
are high enough as shown in (c3).

##  IV Discussion and conclusion

It should be pointed out that many composite skyrmions discussed here were
reported in previous micromagnetic simulations and in experiments. For
example, skyrmionium were observed in space-coiling meta-structure [55],
target skyrmions were observed in FeGe nanodisks [53], in space-coiling meta-
structure [55], and in $\mathrm{Pd/Fe/Ir}$ hexagon systems[62], more exotic
cascaded skyrmion bags were observed in liquid crystals [59], in
$\mathrm{Fe_{3-x}GeTe_{2}}$ [70], and in chiral nematic liquid crystals [71].
Here we propose a coherent way to organize these exotic skyrmion-like
structures.

Composite skyrmions reveal fundamental nature of skyrmion duality. The duality
allows all kinds of superstructures from elementary skyrmions, very similar to
using bricks to construct different buildings in daily life. These
superstructures provide countless new information carriers with distinct
topologies and make new device concepts and new designs possible. For example,
topology-based memory and gates should have high stability, high error
tolerance and low noise [64, 65, 66, 67, 68]. Topologically protected particle
can also find its applications in electro-optic, microfluidic, and
nanoparticle transport [59, 69, 70, 71]. Charge and spin in solid state are
coupled. Electron transport can be affected by the magnetic texture either
through the emerged magnetic and electric fields, or simply through spin-orbit
interactions. Thus, it should be highly interesting to identify those
localized composite skyrmions which can be easily made on-demand, high
controllability, and easily distinguish one from the others, especial
electronically.

Composite skyrmions can be used for magnetic data storage due to high
stability and topological protection. They can also be used as magnetic
sensors due to their sensitivity to external magnetic fields[69]. Another
possible application of composite skyrmions is in energy-efficient computing
applications, such as low-power spintronic devices and neuromorphic computing
due to their varieties in morphology huge number of accessible states such
that they can emulate the behavior of neurons[74]. One may use composite
skyrmions to generate spin currents, which can be used to control
magnetization.

Although we did not consider magnetic field effect in this study, external
magnetic fields can surely modify magnetic anisotropy that, in turn, can
affect stability and properties of composite skyrmions[70]. Magnetic fields
are surely an effective control knob. It should be very important and
interesting to understand in detail how a magnetic field affects composite
skyrmions, not only in terms of practical applications but also fundamentals,

In conclusion, topology and static properties of stable/metastable magnetic
textures of chiral magnets in the absence of magnetic field are fully
determined by parameter $\kappa$, which measures relative strength of chiral
interaction to the ferromagnetic exchange interaction and anisotropy, not by
individual material parameters. $L=4A/(\pi D)$ is the fundamental length scale
of all possible magnetic structures in a chiral magnet. Skyrmions have
particle-continuum-medium duality that results in various interesting
composite skyrmions such as target skyrmions and skyrmion bags. The
complexities and the maximal number of cascade layers of a skyrmion
superstructure with a given topology is sensitive to value of $\kappa$ for
$\kappa\leq 1$ which support isolated circular skyrmions. For chiral magnetic
films of $\kappa>1$, where condensed stripe skyrmions is the thermodynamic
equilibrium state at finite temperature, cascade stripe skyrmions maintain the
width of elementary stripe skyrmions as long as the space between two cascade
layers is larger than the width of elementary stripe skyrmions. When the space
between two cascade layers is less than the fundamental stripe width, the
cascade stripe skyrmions tend to be unstable. In the case of inserting many
skyrmions inside one particular layer, whether in the innermost layer or any
other cascade layer, such that the average skyrmion-skyrmion distance in
around the fundamental stripe width, skyrmions because circular and form a
triangular SkX. Our findings open the possibilities of new concepts and new
designs in applications.

###### Acknowledgements.

This work is supported by the National Key Research and Development Program of
China (No. 2020YFA0309600), the NSFC Grant (No. 11974296), and Hong Kong RGC
Grants (No. 16301619, 16300522, and 16302321).

##  Appendix A Evidences of structure determination by $\kappa$

The putative assertion from Eq. (3) is that (meta)stable spin structures are
completely determined by $\kappa$ up to a proper scaling is a sweet
surprising. In order to further confirm this assertion, we use a film of
$400\,$nm$\times 400\,$nm$\times 0.5\,$nm described by Eq. (1) to simulate
various metastable structures with the same $\kappa$ and different $A/D$. We
set $\kappa=4$ that supports condensed stripe skyrmions and with 10 nucleation
domain randomly distributed in the film. Initially, the model parameters are
$A=1\,\rm pJ\,m^{-1}$, $D=1\,\rm mJ\,m^{-2}$, $K_{u}=0.154\,\rm MJ\,m^{-3}$,
$M_{s}=0.58\,\rm MA\,m^{-1}$. Figure 8(a) below shows the morphology of the
final stable structure with 10 ramified stripe skyrmions. Then we gradually
increase $A/D$ from $1\,$nm to $5\,$nm and decrease $K_{u}$ to keep the
$\kappa=4$ unchanged, and Figs. 8(b)-(e) are the final stable structures for
$A/D=2\,$nm (b); $3\,$nm (c); $4\,$nm (d); $5\,$nm (e). Although these
structures are visually different, showing multiply solutions of Eq. (3), they
are in fact similar. Fig. 8(f) plot $M_{z}(x)$ along the green lines labelled
by ①-⑤ in Figs. 8(a)-(e). The spin profiles fall on the same curve when $x$ is
plotted in the units of $A/D$. The results are exactly what the assertion is
expected.

Figure 8: (a)-(e) The structure of final state when A/D=1 nm (a); 2 nm (b); 3
nm (c); 4 nm (d); 5 nm (e). The inset of (a) shows the structure of the
initial state. (f) The Spin profiles along the green lines labelled by ①-⑤ in
(a)-(e). The symbols are numerical data and the solid lines are the
theoretical fits to
$M_{z}(x)=\frac{\sinh^{2}(L/w)-\sinh^{2}(x/2w)}{\sinh^{2}(L/w)+\sinh^{2}(x/2w)}$
for $-L/2\textless x\textless L/2$ and
$M_{z}(x)=\frac{\sinh^{2}\left[(x-L)/2w\right]-\sinh^{2}(L/w)}{\sinh^{2}(L/w)+\sinh^{2}\left[(x-L)/2w\right]}$
for $L/2\textless x\textless 3L/2$.

## References

  * [1] A. N. Bogdanov and U. K. Rößler, Chiral Symmetry Breaking in Magnetic Thin Films and Multilayers, Phys. Rev. Lett. 87, 037203 (2001).
  * [2] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, 797 (2006).
  * [3] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfeiderer, A. Rosch, A. Neubauer, G. Georgii, and P. Böni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009).
  * [4] 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 (2010).
  * [5] X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet $\mathrm{FeGe}$, Nat. Mater. 10, 106 (2011).
  * [6] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nat. Nanotech. 8, 899 (2013).
  * [7] A. Fert, V. Cros, and J. Sampaio, Skyrmions on the track, Nat. Nanotech. 8, 152 (2013); J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures, ibid. 8, 839 (2013).
  * [8] 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, 636 (2013).
  * [9] X. S. Wang, H. Y. Yuan, and X. R. Wang, A theory on skyrmion size, Commun. Phys. 1, 31 (2018).
  * [10] J. Li, A. Tan, K. W. Moon, A. Doran, M. A. Marcus, A. T. Young, E. Arenholz, S. Ma, R. F. Yang, C. Hwang et al., Tailoring the topology of an artificial magnetic skyrmion, Nat. Commun. 5, 4704 (2014).
  * [11] H. Y. Yuan and X. R. Wang, Skyrmion Creation and Manipulation by Nano-Second Current Pulses, Sci. Rep. 6, 22638 (2016).
  * [12] Y. Zhou and M. A. Ezawa, A reversible conversion between a skyrmion and a domain-wall pair in a junction geometry, Nat. Commun. 5, 4652 (2014).
  * [13] W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen et al., Blowing magnetic skyrmion bubbles, Science 349, 283 (2015).
  * [14] H. Du, R. Che, L. Kong, X. Zhao, C. Jin, C. Wang, J. Yang, W. Ning, R. Li, C. Jin et al., Edge-mediated skyrmion chain and its collective dynamics in a confined geometry, Nat. Commun. 6, 8504 (2015).
  * [15] Y. Onose, Y. Okamura, S. Seki, S. Ishiwata, and Y. Tokura, Observation of magnetic excitations of skyrmion crystal in a helimagnetic insulator $\mathrm{Cu}_{2}\mathrm{OSeO}_{3}$, Phys. Rev. Lett. 109, 037603 (2012).
  * [16] H. S. Park, X. Z. Yu, S. Aizawa, T. Tanikaki, T. Akashi, Y. Takahashi, T. Matsuda, N. Kanazawa, Y. Onose, D. Shindo et al., Observation of the magnetic flux and three-dimensional structure of skyrmion lattices by electron holography, Nat. Nanotech. 9, 337 (2014).
  * [17] X. Gong, H. Y. Yuan, and X. R. Wang, Current-driven skyrmion motion in granular films, Phys. Rev. B 101, 064421 (2020).
  * [18] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and Stefan Blügel, Spontaneous atomic-scale Magnetic skyrmion lattice in two dimensions, Nat. Phys. 7, 713 (2011).
  * [19] S. Krause and R. Wiesendanger, Spintronics: Skyrmionics gets hot, Nat. Mater. 15, 493 (2016).
  * [20] P. Dürrenfeld, Y. Xu, J. Åkerman, and Y. Zhou, Controlled skyrmion nucleation in extended magnetic layers using a nanocontact geometry, Phys. Rev. B 96, 054430 (2017).
  * [21] H. Y. Yuan, X. S. Wang, Man-Hong Yung, and X. R. Wang, Wiggling skyrmion propagation under parametric pumping, Phys. Rev. B 99, 014428 (2019).
  * [22] S. Rohart and A. Thiaville, Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction, Phys. Rev. B 88, 184422 (2013).
  * [23] M. N. Wilson, A. B. Butenko, A. N. Bogdanov, and T. L. Monchesky, Chiral skyrmions in cubic helimagnet films: The role of uniaxial anisotropy, Phys. Rev. B 89, 094411 (2014).
  * [24] N. Romming, A. Kubetzka, C. Hanneken, K. von Bergmann, and R. Wiesendanger, Field-Dependent Size and Shape of Single Magnetic Skyrmions, Phys. Rev. Lett. 114, 177203 (2015).
  * [25] N. Vidal-Silva, A. Riveros, and J. Escrig, Stability of Néel skyrmions in ultra-thin nanodots considering Dzyaloshinskii-Moriya and dipolar interactions, J. Magn. Magn. Mater. 443, 116 (2017).
  * [26] 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).
  * [27] H.-B. Braun, Fluctuations and instabilities of ferromagnetic domain-wall pairs in an external magnetic field, Phys. Rev. B 50, 16485 (1994).
  * [28] A. Siemens, Y. Zhang, J. Hagemeister, E. Y. Vedmedenko, and R. Wiesendanger, Minimal radius of magnetic skyrmions: statics and dynamics, New. J. Phys. 18, 045021 (2016).
  * [29] E. Simon, K. Palotás, L. Rózsa, L. Udvardi, and L. Szunyogh, Formation of magnetic skyrmions with tunable properties in $\mathrm{PdFe}$ bilayer deposited on $\mathrm{Ir(111)}$, Phys. Rev. B 90, 094410 (2014).
  * [30] S. Jaiswal , K. Litzius, I. Lemesh, F. Büttner, S. Finizio, J. Raabe, M. Weigand, K. Lee, J. Langer, B. Ocker et al., Investigation of the Dzyaloshinskii-Moriya interaction and room temperature skyrmions in $\mathrm{W/CoFeB/MgO}$ thin films and microwires, Appl. Phys. Lett. 111, 022409 (2017).
  * [31] E. A. Karhu, U. K. Rößler, A. N. Bogdanov, S. Kahwaji, B. J. Kirby, H. Fritzsche, M. D. Robertson, C. F. Majkrzak, and T. L. Monchesky, Chiral modulations and reorientation effects in $\mathrm{MnSi}$ thin films, Phys. Rev. B 85, 094429 (2012).
  * [32] K. Litzius, I. Lemesh, B. Krüger, P. Bassirian, L. Caretta, K. Richter, F. Büttner, K. Sato, O. A. Tretiakov, J. Förster et al., Skyrmion Hall effect revealed by direct time-resolved X-ray microscopy, Nat. Phys. 13, 170 (2017).
  * [33] S. Woo, K. Litzius, B. Krüger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand et al., Observation of room-temperature magnetic skyrmions and their current-driven dynamics in ultrathin metallic ferromagnets, Nat. Mater. 15, 501 (2016).
  * [34] J. Iwasaki, M. Mochizuki, and N. Nagaosa, Universal current-velocity relation of skyrmion motion in chiral magnets, Nat. Commun. 4, 1463 (2013).
  * [35] C. Reichhardt, D. Ray, and C. J. Olson Reichhardt, Collective Transport Properties of Driven Skyrmions with Random Disorder, Phys. Rev. Lett. 114, 217202 (2015).
  * [36] C. Reichhardt and C. J. Olson Reichhardt, Noise fluctuations and drive dependence of the skyrmion Hall effect in disordered systems, New J. Phys. 18, 095005 (2016).
  * [37] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures, Nat. Nanotech. 8, 839 (2013).
  * [38] S.-Z. Lin, C. Reichhardt, C. D. Batista, and A. Saxena, Particle model for skyrmions in metallic chiral magnets: Dynamics, pinning, and creep, Phys. Rev. B 87, 214419 (2013).
  * [39] W. Koshibae and N. Nagaosa, Theory of current-driven skyrmions in disordered magnets, Sci. Rep. 8, 6238 (2018).
  * [40] R. Juge, S.-G. Je, D. de Souza Chaves, S. Pizzini, L. D. Buda-Prejbeanu, L. Aballe, M. Foerster, A. Locatelli, T. O. Mentes, A. Sala et al., Magnetic skyrmions in confined geometries: Effect of the magnetic field and the disorder, J. Magn. Magn. Mater. 455, 3 (2018).
  * [41] J.-V. Kim and M. -W. Yoo, Current-driven skyrmion dynamics in disordered films, Appl. Phys. Lett. 110, 132404 (2017).
  * [42] S. Hoshino and N. Nagaosa, Theory of the magnetic skyrmion glass, Phys. Rev. B 97, 024413 (2018).
  * [43] S. Woo, K. M. Song, X. Zhang, Y, Zhou, M. Ezawa, X. Liu, S. Finizio, J. Raabe, N. J. Lee, S.-I. Kim et al., Current-driven dynamics and inhibition of the skyrmion Hall effect of ferrimagnetic skyrmions in $\mathrm{GdFeCo}$ films, Nat. Commun. 9, 959 (2018).
  * [44] X. R. Wang, X. C. Hu, and H. T. Wu, Stripe skyrmions and skyrmion crystals, Commun. Phys. 4, 142 (2021).
  * [45] H. T. Wu, X. C. Hu, K. Y. Jing, and X. R. Wang, Size and profile of skyrmions in skyrmion crystals, Commun. Phys. 4, 210 (2021).
  * [46] H. T. Wu, X. C. Hu, and X. R. Wang, Nematic and smectic stripe phases and stripe-SkX transformations, Science China 65, 247512 (2022).
  * [47] X. C. Hu, H. T. Wu, and X. R. Wang, A theory of skyrmion crystal formation, Nanoscale 14, 7516 (2022).
  * [48] X. R. Wang, X. C. Hu, and Z. Z. Sun, Topological Equivalence of Stripy States and Skyrmion Crystals, Nano Letters (in press), (2023).
  * [49] A. O. Leonov, U. K. Rößler, and M. Mostovoy, Target-skyrmions and skyrmion clusters in nanowires of chiral magnets, EDP Sciences 75, 05002 (2014).
  * [50] X. Zhang, J. Xia, Y. Zhou, D. Wang, X. Liu, W. Zhao, and M. Ezawa, Control and manipulation of a magnetic skyrmionium in nanostructures, Phys. Rev. B 94, 094420 (2016).
  * [51] S. Zhang, F. Kronast, G. van der Laan, and T. Hesjedal, Real-space observation of skyrmionium in a ferromagnet-magnetic topological insulator heterostructure, Nano Lett. 18, 1057 (2018).
  * [52] M. M. M. Bhukta, M. Abhilash, P. Gajanan, M. Sougata, B. S. Braj, and B. Subhankar, A novel chiral spin texture: Antiferromagnetic Skyrmionium, arXiv 1810.08262, (2018).
  * [53] F. Zheng, H. Li, S. Wang, D. Song, C. Jin, W. Wei, A. Kovács, J. Zang, M. Tian, Y. Zhang et al., Direct imaging of a zero-field target skyrmion and its polarity switch in a chiral magnetic nanodisk, Phys. Rev. Lett. 119, 197205 (2017).
  * [54] Y. Zhang, C. Xu, P. Chen, Y. Nahas, S. Prokhorenko, and L. Bellaiche, Emergence of skyrmionium in a two-dimensional $\mathrm{CrGe(Se,Te)}_{3}$ Janus monolayer, Phys. Rev. B 102, 241107 (2020).
  * [55] Z. L. Deng, T. Shi, A. Krasnok, X. P. Li, and A. Alú , Observation of localized magnetic plasmon skyrmions, Nat. Commun. 13, 1 (2022).
  * [56] A. G. Kolesnikov, M. E. Stebliy, A. S. Samardak, and A. V. Ognev, Skyrmionium–high velocity without the skyrmion Hall effect, Sci. Rep. 8, 1 (2018).
  * [57] L. Shen, X. Li, Y. Zhao, J. Xia, G. P. Zhao, and Y. Zhou, Current-induced dynamics of the antiferromagnetic skyrmion and skyrmionium, Phys. Rev. Appl. 12, 064033 (2019).
  * [58] S. Li, J. Xia, X. Zhang, M. Ezawa, W. Kang, X. X. Liu, Y. Zhou, and W. S. Zhao, Dynamics of a magnetic skyrmionium driven by spin waves, Appl. Phys. Lett. 112, 142404 (2018).
  * [59] D. Foster, C. Kind, P. J. Ackerman, J. -S. B. Tai, M. R. Dennis, and I. I. Smalyukh, Two-dimensional skyrmion bags in liquid crystals and ferromagnets, Nat. Phys. 15, 655 (2019).
  * [60] J. Tang, L. Kong, Y. Wu, W. Wang, Y. Chen, Y. Wang, J. Li, Y. Soh, Y. Xiong, M. Tian et al., Target bubbles in $\mathrm{Fe_{3}Sn_{2}}$ nanodisks at zero magnetic field, ACS Nano 14, 10986 (2020).
  * [61] X. Shen, R. Zhao, L. Ji, C. Hu, W. Ren, W. Chen, Y. Li, J. Zhang, X. Zhang, and X. Dong, Breaking the topological protection of target skyrmions by the excitation of spin wave modes under microwave magnetic field, J. Magn. Magn. Mater. 541, 168521 (2022).
  * [62] D. Cortés-Ortuño, N. Romming, M. Beg, K. v. Bergmann, A. Kubetzka, O. Hovorka, H. Fangohr, and R. Wiesendanger, Nanoscale magnetic skyrmions and target states in confined geometries, Phys. Rev. B 99, 214408 (2019).
  * [63] C. Kind, S. Friedemann, and D. Read, Existence and stability of skyrmion bags in thin magnetic films, Appl. Phys. Lett. 116, 022413 (2020).
  * [64] J. Tang, Y. D. Wu, W. W. Wang, L. Y. Kong, B. Y. Lv, W. S. Wei, J. D. Zang, M. L. Tian, and H. F. Du, Magnetic skyrmion bundles and their current-driven dynamics, Nature Nanotech. 16, 1086 (2021).
  * [65] C. Kind and D. Foster, Magnetic skyrmion binning, Phys. Rev. B 103, L100413 (2023).
  * [66] Z. Zhang, M. Xua, G. Jiang, J. Zhang, D. Meng, W. Chen, Y. Chen, and C. Hu, High-density racetrack memory based on magnetic skyrmion bags controlled by voltage gates, J. Appl. Phys. 132, 113901 (2022).
  * [67] A. O. Leonov, Clusters and chains in bulk and thin-layered cubic helimagnets, Phys. Rev. B 105, 094404 (2022).
  * [68] R. Chen and Y. Li, Voltage-Controlled Skyrmionic Interconnect with Multiple Magnetic Information Carriers, ACS Appl. Mater. Interfaces 14, 30420 (2022).
  * [69] T. K. Bhowmick, A. De, and R. K. Lake, High figure of merit magneto-optics from interfacial skyrmions on topological insulator, Phys. Rev. B 98, 024424 (2018).
  * [70] L. Powalla, M. T. Birch, K. Litzius, S. Wintz, F. S. Yasin, L. A. Turnbull, F. Schulz, D. A. Mayoh, G. Balakrishnan, M. Weigand et al., Seeding and Emergence of Composite Skyrmions in a van der Waals Magnet, Adv. Mater. 2208930 (2023).
  * [71] Y. Shen and I. Dierking, Electrically driven formation and dynamics of skyrmionic solitons in chiral nematics, Phys. Rev. Appl. 15, 054023 (2023).
  * [72] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia- Sanchez, and F. B. V. Waeyenberge, The design and verification of MuMax3, AIP Adv. 4, 107133 (2014).
  * [73] A. Soumyanarayanan, M. Raju, A. L. Gonzalez Oyarce, A. K. C. Tan, M.-Y. Im, A. P. Petrović, P. Ho, K. H. Khoo, M. Tran, C. K. Gan, F. Ernult et al., Tunable room-temperature magnetic skyrmions in $\mathrm{Ir/Fe/Co/Pt}$ multilayers, Nature Mater. 16, 898 (2017).
  * [74] K. M. Song, J.-S. Jeong, B. Pan, X. Zhang, J. Xia, S. Cha, T.-E. Park, K. Kim, S. Finizio, J. Raabe et al., Skyrmion-based artificial synapses for neuromorphic computing, Nat. Electron. 3, 148 (2020).


$\tanh$-mark where the correction is only approximate, the relative error does
increase to around 1% on scales larger than $25\leavevmode\nobreak\
\,h^{-1}\,{\rm Mpc}$. Although, at very large scales the shot-noise correction
does not greatly impact the uncertainty since the overall contribution of shot
noise on those scales is marginal. It is evident from the figure that most of
the uncertainty from the correction is introduced on scales below
$25\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. Here we are within the smoothing
radii and shot noise is expected to be the strongest. It is important to note
that this does not capture the effect on the relative error between not
applying a correction at all and applying a correction. We can only conclude
that while we might be able to accurately recover the true signal, the step of
applying the correction via the fitting introduces additional uncertainties
that increase towards smaller scales. This uncertainty is expected to be
smaller the higher $\bar{N}$ is, since the fitting process should be less
prone to uncertainties. Intuitively, this means that the points that we need
to fit as depicted in Figure 4 and 6 are distributed closer to $1/\bar{N}=0$
and hence the extrapolation to the $y$-axis is more robust.

Figure 20: Relative error of the weighted correlation function in Covmos. The
mark is set to the $\tanh$ with parameters fixed to $(a,b)=(0.6,-0.5)$ for the
solid line and the toy model for the dashed line. Green colours refer to the
correction using 30 realisations of depletions and blue colours stand for the
undepleted case. The relative error is computed as the mean standard deviation
over 5 realisations divided by the mean.

We have briefly touched upon the impact of the MAS kernel on our methodology
in Section 6.5 with Figure 11 and it is crucial to assess this in more detail.
One can show that the shot-noise correction is smaller the higher the order of
the MAS is, which can be understood as a larger smoothing scale as higher-
order MAS kernels are more extended in configuration space. In Figure 21 and
22 we present an illustration of how the smoothing scale as well as the shape
of the MAS affect the behaviour of shot noise, particularly on the amplitude
of the different powers of $1/\bar{N}$. We show the result of the standard
fitting procedure for obtaining the shot-noise-corrected signal and fits do
include the undepleted catalogue, which is not accessible in real data. This
figure illustrates the contributions from different powers of $1/\bar{N}$ in
the shot-noise polynomial. Having included the undepleted case enables an
accurate estimate of the shot-noise behaviour and identifying any breakdown at
a given polynomial order. As we can already seen from the fits in Figure 21,
once an NGP MAS is employed the low orders for the polynomials do not seem to
be sufficient anymore to describe the data. This is further underlined by
Figure 22 showing an increase in amplitude for the polynomial coefficients
when lowering the MAS order. This means that the higher the smoothing scale
the less contributions are coming from higher-order shot-noise expressions,
and the better we can fit the dependency with a low-order polynomial.
Intuitively this makes sense if we hypothetically increase the smoothing scale
to infinity. In that case, the obtained density field would be a constant in
space and therefore all galaxies will have the same weight. In that scenario,
the weighted correlation function will reduce to the unweighted correlation
function that is only affected by shot noise at zero-lag. Hence in our
measurements there would be no contamination and the polynomial fits would
just be a constant. Such a trend can also be seen to some extent in Figure 21,
where the curve becomes more and more linear and converges to a vertical line
the higher the order of the MAS.

(a)

(b)

Figure 21: Impact of the smoothing scale on the polynomial behaviour of shot
noise. The mark is fixed to the White mark with $(\rho_{*},p)=(4.0,10.0)$.
Left and right panels show the fit of $w_{f}$ and $\bar{m}_{f}$ for one
realisation of Covmos, respectively. Depleted measurements are obtained via
taking the mean of 30 realisations of depletions. Linestyles refer to
different MAS used for obtaining the density field. The error of the
undepleted case and 1.7% depletion are computed as 10% of the smallest
uncertainty of the other depletions.

(a)

(b)

Figure 22: Coefficient amplitudes of the polynomial fits corresponding to
Figure 21.

Using two different types of catalogues for gauging the shot-noise correction
has its limitations, which we discuss in the following. The Covmos catalogues
were absolutely necessary in the first place in order to have an almost noise-
free signal to test the method. An inherent difference, which has to be kept
in mind when interpreting our findings, is the fact that Covmos realisations
are not simply an upscaled version of the ELEPHANT suite, rather a completely
different set of catalogues. First and foremost, Covmos catalogues consist of
dark matter particles instead of galaxies and second, the redshift and assumed
cosmology are different compared to ELEPHANT. Most importantly, the Covmos
catalogues are not extracted from full N-body simulations. This means that the
features and the general shape of the weighted correlation functions look
different between Covmos and ELEPHANT simulations. Nevertheless, the
functional form of the Taylor expansion of the mark stays the same and having
access to realisations with very large number densities of points served the
purpose of validating the method for estimating the corrected signal. Care has
to be taken, however, in the choice of the mark to comply with the convergence
criterion of the Taylor expansion as described in Section 6.4, which is
universally valid for both sets of catalogues. By construction, the ELEPHANT
suite is limited by the the small number of galaxies in the catalogues. This
results in the shot-noise correction to introduce larger uncertainties than
would do in larger catalogues. Having only twice as many objects in the
simulations would already half the distance to extrapolate from the undepleted
case to $1/\bar{N}=0$ on a linear scale. The ELEPHANT simulations are
therefore a particularly difficult case and we expect better results if
applied to state-of-the-art N-body simulations with higher densities.
Nevertheless, as long as the mark is chosen appropriately, a robust recovery
of the true signal is possible. For future analysis we recommend to use
simulations with higher resolution in order to mitigate the impact of shot
noise.

We conclude the discussion with a brief summary of the main points raised in
this section as well as Section 6 to be considered when applying our
methodology to correct for shot noise in marked correlation functions:

  * –

Enough realisations of depletion should be used in order to get converged
depleted point sets. In our setting 30 depletions appeared to be enough.

  * –

The polynomial order should be chosen accordingly such that the shot noise
dependency is well modelled without overfitting.

  * –

A higher-order MAS (e.g. TSC/PCS) is preferred to reduce possible bias due to
shot noise at low scales.

  * –

In order to apply our methodology to correct for shot noise, a Taylor
expansion of the mark function in powers of $\delta_{Rf}$ has to exist.

  * –

The convergence radius $B$ of the Taylor series should be larger than
$\sqrt{\langle\delta_{Rf}^{2}\rangle_{\rho}}$ in order to satisfy the
assumption of a convergent Taylor expansion.

  * –

In addition to the criterion, as given in the last bullet point, also the
fraction of points inside the convergence radius can be checked, which is more
agnostic about the actual PDF of $\delta_{Rf}$ in the catalogues.

  * –

A Taylor expansion with non-zero and decreasing coefficients both for odd and
even powers of $\delta_{Rf}$ are preferred in order to allow for a robust
estimation of shot noise.

  * –

We do not recommend to use marks which have a mean mark very close to zero
after correcting for shot noise as measurements get very unstable and
uncertainties diverge.

##  9 Conclusions

In this work, we have studied marked correlation functions in the context of
detecting MG and how discreteness effects from estimating marks on a finite
point set propagate into the measurement of marked correlation functions. We
utilised the Covmos realisations (Baratta et al., 2023) that have a
particularly high density of points, making them the most suitable for this
purpose, and of ELEPHANT simulations with HOD galaxies (Alam et al., 2021) to
investigate the discriminating power of marked correlation functions between
MG and GR. The latter is comprised of several realisations of GR as well as of
$f(R)$ and nDGP gravity theories. These are two particularly interesting
modifications to GR to be studied with marked correlation functions because
they exhibit screening mechanisms making the fifth force dependent on the
environment. We proposed several marks based on large-scale environments using
the T-Web formalism as well as local density. This includes marks that creates
anti-correlation between galaxies in different environments or between
galaxies in low- and high-density regions.

For the first time we undertook a thorough investigation of a possible bias
due to shot noise in marked correlation functions. We showed on a toy model
that the effect of shot noise can be treated analytically by computation of a
small amount of terms. We were able to recover the signal
$w_{f}(\mathbf{r})=(1+W_{f}(\mathbf{r}))\bar{m}_{f}$ from the undepleted
catalogue to sub-percent precision. For general marks, under the assumption
that the mark function can be Taylor-expanded in powers of the density
contrast, we showed that an analytic treatment is hopeless due to the
necessity of computing an infinity of higher-order correlators. Instead, we
developed a methodology for estimating the shot-noise-corrected signal from
measuring the weighted correlated function and mean mark in catalogues
depleted to different densities. This is possible by noting the resummation
behaviour of shot-noise contributions to
$w_{f}(\mathbf{r})=(1+W_{f}(\mathbf{r}))\bar{m}_{f}^{2}$ and $\bar{m}_{f}$ as
a power series of the reciprocal of the mean number of points per grid cell
$1/\bar{N}$. By applying our method to the $\tanh$-mark in the Covmos
realisations, we were able to recover an unbiased signal of
$1+W_{f}(\mathbf{r})$ within 5% accuracy for all tested scales. We proceeded
with an extension of the formalism in redshift space, where we found the same
method to be applicable. Furthermore, we derived a measurable criterion based
on the work of Philcox et al. (2020) to assess the validity of assuming a
Taylor expansion of the mark and provide guidelines for the application of our
methodology for shot-noise correction. We found effects of shot noise mostly
on scales below 20-30$\,h^{-1}\,{\rm Mpc}$ when using the White mark, which
might be important for the modelling of marked correlation functions, although
the impact on the relative difference between GR and MG appears to be mild.
Moreover, we found that the NGP MAS to give biased results due to higher-order
terms in the $1/\bar{N}$ series being non-negligible. This makes the NGP MAS a
sub-optimal choice compared to higher-order schemes, such as PCS.

Equipped with a robust method to recover the true signal in measured marked
correlation function, we tested the performance of the previously proposed
marks on the ELEPHANT simulations. Concerning marks based on the local
density, we did not find the White mark to be particularly powerful on large
scales. Only on the very lowest scales, below $20\leavevmode\nobreak\
\,h^{-1}\,{\rm Mpc}$, we reported significant differences. In redshift space,
the situation changes slightly for the N1 model where we found differences in
the quadrupole at $s>20\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. We found
that the novel $\tanh$-mark that we introduced is very effective. It allows
significant differences for $f(R)$-gravity with $\log(|f_{R0}|)=-6$ compared
to GR, and uniquely up to scales of 80$\,h^{-1}\,{\rm Mpc}$. At those scales,
modelling the weighted correlation function is more tractable; making this
mark an excellent candidate to test for deviations from GR in real surveys.
Furthermore, we found promising results when using the tidal torque as a mark,
with significant differences up to scales of $r\approx 60\leavevmode\nobreak\
\,h^{-1}\,{\rm Mpc}$. The use of anti-correlation together with large-scale
environments, as in the $\textrm{Void}_{\textrm{AC}}$\- and
$\textrm{Wall}_{\textrm{AC}}$-mark, exhibits significant discriminating power
for several MG theories on moderate scales. However, these findings have to be
tested with high-density simulations to assess if they are biased by
discreteness effects, as our correction method cannot be applied
straightforwardly to those types of mark.

In summary, this work demonstrates that correcting for shot noise in marked
correlation functions is of paramount importance to measure unbiased
amplitudes without being plagued by shot noise, and in turn to be able to
distinguish MG from GR. This is also particularly important for the modelling
of the weighted correlation function in the same way as it is for the power
spectrum. Generally, we found shot noise to have the strongest impact on small
to intermediate scales. Marks that incorporates an anti-correlation between
objects in high- and low-density regions by switching signs in the mark are
found to be the most effective for distinguishing between GR from MG, also
when using scales beyond $20\leavevmode\nobreak\ \,h^{-1}\,{\rm Mpc}$. In the
future, extending the concept of anti-correlation in weighted correlation
functions to different marks could alleviate the current constraint due to the
convergence radius of the Taylor expansion as is the case for the
$\tanh$-mark. In general, a consolidation of the $\tanh$-mark performances on
improved MG and GR simulations with higher densities would be desired.
Moreover, a thorough investigation of a kind of model-independent shot-noise
effect on general marked correlation functions would enable an appropriate
correction for marks based on large-scale environments or the tidal
torque/field, which we showed to be interesting candidates. A future study
should assess the capability of the Lagrangian perturbation theory model to
capture the behaviour of our novel mark based on the local density on
intermediate scales. The $\tanh$-mark could then serve as an optimal choice
for a weighted clustering analysis in current and future galaxy surveys since
no accurate modelling of small scales is required. Having a relevant model of
marked correlation functions together with a high-performance mark should add
a powerful observable to find the needle in the haystack of gravity theories.

###### Acknowledgements.

The project leading to this publication has received funding from the
Excellence Initiative of Aix-Marseille University - A*MIDEX, a French
”Investissements d’Avenir” programme (AMX-19-IET-008 - IPhU). M. Kärcher is
funded by the Excellence Initiative of Aix-Marseille University - A*MIDEX, a
French ”Investissements d’Avenir” programme (AMX-19-IET-008 - IPhU). This
research made use of matplotlib, a Python library for publication quality
graphics (Hunter, 2007).

## References

  * Alam et al. (2015) Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12
  * Alam et al. (2021) Alam, S., Arnold, C., Aviles, A., et al. 2021, J. Cosmology Astropart. Phys., 2021, 050
  * Alam et al. (2019) Alam, S., Zu, Y., Peacock, J. A., & Mandelbaum, R. 2019, MNRAS, 483, 4501
  * Armijo et al. (2024a) Armijo, J., Baugh, C. M., Norberg, P., & Padilla, N. D. 2024a, MNRAS, 529, 2866
  * Armijo et al. (2024b) Armijo, J., Baugh, C. M., Norberg, P., & Padilla, N. D. 2024b, MNRAS, 528, 6631
  * Armijo et al. (2018) Armijo, J., Cai, Y.-C., Padilla, N., Li, B., & Peacock, J. A. 2018, MNRAS, 478, 3627
  * Aubert et al. (2022) Aubert, M., Cousinou, M.-C., Escoffier, S., et al. 2022, MNRAS, 513, 186
  * Aviles et al. (2020) Aviles, A., Koyama, K., Cervantes-Cota, J. L., Winther, H. A., & Li, B. 2020, J. Cosmology Astropart. Phys., 2020, 006
  * Baratta et al. (2023) Baratta, P., Bel, J., Gouyou Beauchamps, S., & Carbone, C. 2023, A&A, 673, A1
  * Baratta et al. (2020) Baratta, P., Bel, J., Plaszczynski, S., & Ealet, A. 2020, A&A, 633, A26
  * Barreira et al. (2015) Barreira, A., Bose, S., & Li, B. 2015, J. Cosmology Astropart. Phys., 2015, 059
  * Battye et al. (2018) Battye, R. A., Bolliet, B., & Pace, F. 2018, Phys. Rev. D, 97, 104070
  * Bautista et al. (2021) Bautista, J. E., Paviot, R., Vargas Magaña, M., et al. 2021, MNRAS, 500, 736
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Beisbart & Kerscher (2000) Beisbart, C. & Kerscher, M. 2000, ApJ, 545, 6
  * Bertotti et al. (2003) Bertotti, B., Iess, L., & Tortora, P. 2003, Nature, 425, 374
  * Beutler et al. (2012) Beutler, F., Blake, C., Colless, M., et al. 2012, MNRAS, 423, 3430
  * Blake et al. (2020) Blake, C., Amon, A., Asgari, M., et al. 2020, A&A, 642, A158
  * Blake et al. (2011) Blake, C., Brough, S., Colless, M., et al. 2011, MNRAS, 415, 2876
  * Bonnaire et al. (2022) Bonnaire, T., Aghanim, N., Kuruvilla, J., & Decelle, A. 2022, A&A, 661, A146
  * Bose et al. (2015) Bose, S., Hellwing, W. A., & Li, B. 2015, J. Cosmology Astropart. Phys., 2015, 034
  * Brax et al. (2022) Brax, P., Davis, A.-C., & Elder, B. 2022, Phys. Rev. D, 106, 044040
  * Castorina et al. (2015) Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, J. Cosmology Astropart. Phys., 2015, 043
  * Cautun et al. (2013) Cautun, M., van de Weygaert, R., & Jones, B. J. T. 2013, MNRAS, 429, 1286
  * Cautun et al. (2014) Cautun, M., van de Weygaert, R., Jones, B. J. T., & Frenk, C. S. 2014, MNRAS, 441, 2923
  * Chan & Blot (2017) Chan, K. C. & Blot, L. 2017, Phys. Rev. D, 96, 023528
  * Chan et al. (2012) Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509
  * Chaniotis & Poulikakos (2004) Chaniotis, A. K. & Poulikakos, D. 2004, Journal of Computational Physics, 197, 253
  * Clifton et al. (2012) Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep, 513, 1
  * Cognola et al. (2008) Cognola, G., Elizalde, E., Nojiri, S., et al. 2008, Phys. Rev. D, 77, 046009
  * Damour & Polyakov (1994) Damour, T. & Polyakov, A. M. 1994, Nuclear Physics B, 423, 532
  * De Felice & Tsujikawa (2010) De Felice, A. & Tsujikawa, S. 2010, Living Reviews in Relativity, 13, 3
  * de la Torre et al. (2013) de la Torre, S., Guzzo, L., Peacock, J. A., et al. 2013, A&A, 557, A54
  * de la Torre et al. (2017) de la Torre, S., Jullo, E., Giocoli, C., et al. 2017, A&A, 608, A44
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036
  * Desmond & Ferreira (2020) Desmond, H. & Ferreira, P. G. 2020, Phys. Rev. D, 102, 104060
  * Dvali et al. (2000) Dvali, G., Gabadadze, G., & Porrati, M. 2000, Physics Letters B, 485, 208
  * Euclid Collaboration et al. (2024) Euclid Collaboration, Mellier, Y., Abdurro’uf, et al. 2024, arXiv e-prints, arXiv:2405.13491
  * Falck et al. (2012) Falck, B. L., Neyrinck, M. C., & Szalay, A. S. 2012, ApJ, 754, 126
  * Forero-Romero et al. (2009) Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., & Yepes, G. 2009, MNRAS, 396, 1815
  * Guth (1981) Guth, A. H. 1981, Phys. Rev. D, 23, 347
  * Guzzo et al. (2008) Guzzo, L., Pierleoni, M., Meneux, B., et al. 2008, Nature, 451, 541
  * Hamaus et al. (2022) Hamaus, N., Aubert, M., Pisani, A., et al. 2022, A&A, 658, A20
  * Heavens & Peacock (1988) Heavens, A. & Peacock, J. 1988, MNRAS, 232, 339
  * Hernández-Aguayo et al. (2018) Hernández-Aguayo, C., Baugh, C. M., & Li, B. 2018, MNRAS, 479, 4824
  * Hinshaw et al. (2013) Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19
  * Hinterbichler & Khoury (2010) Hinterbichler, K. & Khoury, J. 2010, Phys. Rev. Lett., 104, 231301
  * Hu & Sawicki (2007) Hu, W. & Sawicki, I. 2007, Phys. Rev. D, 76, 064004
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Ishak (2019) Ishak, M. 2019, Living Reviews in Relativity, 22, 1
  * Jullo et al. (2019) Jullo, E., de la Torre, S., Cousinou, M. C., et al. 2019, A&A, 627, A137
  * Khoury & Weltman (2004a) Khoury, J. & Weltman, A. 2004a, Phys. Rev. D, 69, 044026
  * Khoury & Weltman (2004b) Khoury, J. & Weltman, A. 2004b, Phys. Rev. Lett., 93, 171104
  * Koyama & Silva (2007) Koyama, K. & Silva, F. P. 2007, Phys. Rev. D, 75, 084040
  * Landy & Szalay (1993) Landy, S. D. & Szalay, A. S. 1993, ApJ, 412, 64
  * Layzer (1956) Layzer, D. 1956, AJ, 61, 383
  * Libeskind et al. (2018) Libeskind, N. I., van de Weygaert, R., Cautun, M., et al. 2018, MNRAS, 473, 1195
  * Liu et al. (2021) Liu, R., Valogiannis, G., Battaglia, N., & Bean, R. 2021, Phys. Rev. D, 104, 103519
  * Llinares & McCullagh (2017) Llinares, C. & McCullagh, N. 2017, MNRAS, 472, L80
  * Lombriser (2014) Lombriser, L. 2014, Annalen der Physik, 264, 259
  * Lombriser et al. (2009) Lombriser, L., Hu, W., Fang, W., & Seljak, U. 2009, Phys. Rev. D, 80, 063536
  * Lombriser et al. (2015) Lombriser, L., Simpson, F., & Mead, A. 2015, Phys. Rev. Lett., 114, 251101
  * Manera et al. (2013) Manera, M., Scoccimarro, R., Percival, W. J., et al. 2013, MNRAS, 428, 1036
  * Martin (2012) Martin, J. 2012, Comptes Rendus Physique, 13, 566
  * Massara et al. (2021) Massara, E., Villaescusa-Navarro, F., Ho, S., Dalal, N., & Spergel, D. N. 2021, Phys. Rev. Lett., 126, 011301
  * Neyrinck (2008) Neyrinck, M. C. 2008, MNRAS, 386, 2101
  * Nicolis et al. (2009) Nicolis, A., Rattazzi, R., & Trincherini, E. 2009, Phys. Rev. D, 79, 064036
  * Paillas et al. (2021) Paillas, E., Cai, Y.-C., Padilla, N., & Sánchez, A. G. 2021, MNRAS, 505, 5731
  * Peebles & Hauser (1974) Peebles, P. J. E. & Hauser, M. G. 1974, ApJS, 28, 19
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565
  * Philcox et al. (2021) Philcox, O. H. E., Aviles, A., & Massara, E. 2021, J. Cosmology Astropart. Phys., 2021, 038
  * Philcox et al. (2020) Philcox, O. H. E., Massara, E., & Spergel, D. N. 2020, Phys. Rev. D, 102, 043516
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A14
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Reyes et al. (2010) Reyes, R., Mandelbaum, R., Seljak, U., et al. 2010, Nature, 464, 256
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009
  * Riess et al. (2022) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7
  * Satpathy et al. (2019) Satpathy, S., A C Croft, R., Ho, S., & Li, B. 2019, MNRAS, 484, 2148
  * Schaap & van de Weygaert (2000) Schaap, W. E. & van de Weygaert, R. 2000, A&A, 363, L29
  * Schmidt (2009) Schmidt, F. 2009, Phys. Rev. D, 80, 123003
  * Sefusatti et al. (2016) Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624
  * Sheth (2005) Sheth, R. K. 2005, MNRAS, 364, 796
  * Simpson et al. (2013) Simpson, F., Heavens, A. F., & Heymans, C. 2013, Phys. Rev. D, 88, 083510
  * Simpson et al. (2011) Simpson, F., James, J. B., Heavens, A. F., & Heymans, C. 2011, Phys. Rev. Lett., 107, 271301
  * Sinha & Garrison (2019) Sinha, M. & Garrison, L. 2019, in Software Challenges to Exascale Computing, ed. A. Majumdar & R. Arora (Singapore: Springer Singapore), 3–20
  * Sinha & Garrison (2020) Sinha, M. & Garrison, L. H. 2020, MNRAS, 491, 3022
  * Sotiriou & Faraoni (2010) Sotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451
  * Sousbie (2011) Sousbie, T. 2011, MNRAS, 414, 350
  * Tröster et al. (2020) Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10
  * Tsujikawa (2010) Tsujikawa, S. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, ed. G. Wolschin, Vol. 800, 99–145
  * Vainshtein (1972) Vainshtein, A. I. 1972, Physics Letters B, 39, 393
  * Valogiannis & Bean (2018) Valogiannis, G. & Bean, R. 2018, Phys. Rev. D, 97, 023535
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2
  * White (2016) White, M. 2016, J. Cosmology Astropart. Phys., 2016, 057
  * White & Padmanabhan (2009) White, M. & Padmanabhan, N. 2009, MNRAS, 395, 2381
  * Williams et al. (2004) Williams, J. G., Turyshev, S. G., & Boggs, D. H. 2004, Phys. Rev. Lett., 93, 261101
  * Williams et al. (2012) Williams, J. G., Turyshev, S. G., & Boggs, D. H. 2012, Classical and Quantum Gravity, 29, 184004
  * Xiao et al. (2022) Xiao, X., Yang, Y., Luo, X., et al. 2022, MNRAS, 513, 595
  * Yang et al. (2020) Yang, Y., Miao, H., Ma, Q., et al. 2020, ApJ, 900, 6
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84
  * Zhang et al. (2007) Zhang, P., Liguori, M., Bean, R., & Dodelson, S. 2007, Phys. Rev. Lett., 99, 141302
  * Zheng et al. (2007) Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760

##  Appendix A Marks and cross-correlation

We assume a population of $N$ galaxies that can be split up into a
$1$-population and a $2$-population with respective numbers $N_{1}$ and
$N_{2}$, summing up to $N=N_{1}+N_{2}$. This split could have been done by
using the $\textrm{Void}_{\textrm{AC}}$ mark where we assign a mark of -1 to
all galaxies residing in voids and +1 otherwise. We assume boxes with periodic
boundary conditions, as in the main text, and $RR_{n}$ denotes the normalised
$RR$ counts such that $RR_{n}=RR/(N_{R}(N_{R}-1))$. Next we define the
correlation functions for each population of galaxies as well as the cross-
correlation among the two sub-populations to be

$\begin{split}\xi&=\frac{DD}{N(N-1)RR_{n}}-1,\\\
\xi_{11}&=\frac{D_{1}D_{1}}{N_{1}(N_{1}-1)RR_{n}}-1,\\\
\xi_{22}&=\frac{D_{2}D_{2}}{N_{2}(N_{2}-1)RR_{n}}-1,\\\
\textrm{and}\quad\xi_{12}&=\frac{D_{1}D_{2}}{N_{1}N_{2}RR_{n}}-1.\end{split}$
(77)

We note here that when we cross-correlate we do not assume double counting
hence the normalisation by the total number of possible pairs is only
$N_{1}N_{2}$. Terms of the form $D_{i}D_{j}$ refer to unnormalised counts of
pairs between the $i$\- and $j$-population.

The total double counted pair counts can be split up into contributions like

$DD=D_{1}D_{1}+D_{2}D_{2}+2D_{1}D_{2},$ (78)

where the factor of 2 is necessary as the cross-counts are not double counted
on their own. This allows for the following split of the total correlation
function

$\begin{split}\xi=&\frac{D_{1}D_{1}}{N(N-1)RR_{n}}+\frac{D_{2}D_{2}}{N(N-1)RR_{n}}\\\
&+2\frac{D_{1}D_{2}}{N(N-1)RR_{n}}-1\\\
=&\frac{N_{1}(N_{1}-1)}{N(N-1)}\frac{D_{1}D_{1}}{N_{1}(N_{1}-1)RR_{n}}+\frac{N_{2}(N_{2}-1)}{N(N-1)}\frac{D_{2}D_{2}}{N_{2}(N_{2}-1)RR_{n}}\\\
&+2\frac{N_{1}N_{2}}{N(N-1)}\frac{D_{1}D_{2}}{N_{1}N_{2}RR_{n}}-1\\\
=&f_{11}(\xi_{11}+1)+f_{22}(\xi_{22}+1)+2f_{12}(\xi_{12}+1)-1,\end{split}$
(79)

where we defined

$\begin{split}f_{11}&=\frac{N_{1}(N_{1}-1)}{N(N-1)},\\\
f_{22}&=\frac{N_{2}(N_{2}-1)}{N(N-1)},\\\ \textrm{and}\quad
f_{12}&=\frac{N_{1}N_{2}}{N(N-1)}.\end{split}$ (80)

By realising that $f_{11}+f_{12}+2f_{12}=1$ we can finally write

$\xi=f_{11}\,\xi_{11}+f_{22}\,\xi_{22}+2f_{12}\,\xi_{12}.$ (81)

Completely analogous can be proceeded if we consider weighted correlation
functions where the mark can only take two values, like in the
$\textrm{Void}_{\textrm{AC}}$ mark. We define the individual weighted
correlation functions to be

$\begin{split}W&=\frac{WW}{(\sum m_{i})^{2}-\sum
m_{i}^{2}}\frac{1}{RR_{n}}-1,\\\
W_{11}&=\frac{W_{1}W_{1}}{(\sum_{1}m_{i})^{2}-\sum_{1}m_{i}^{2}}\frac{1}{RR_{n}}-1,\\\
W_{22}&=\frac{W_{2}W_{2}}{(\sum_{2}m_{i})^{2}-\sum_{2}m_{i}^{2}}\frac{1}{RR_{n}}-1,\\\
\textrm{and}\quad
W_{12}&=\frac{W_{1}W_{2}}{\sum_{1}m_{i}\sum_{2}m_{i}}\frac{1}{RR_{n}}-1,\end{split}$
(82)

for which we used the shorthand notation $\sum=\sum_{1}+\sum_{2}$ to indicate
sums over weights belonging solely to galaxies from either population 1 or
population 2. The total sum $WW$ over products of weights can be split up into
contributions from $W_{1}W_{1}$, $W_{2}W_{2}$ and $W_{1}W_{2}$ in an analogous
way as the $DD$ counts in the unweighted case. Defining prefactors as

$\begin{split}f^{W}_{11}&=\frac{(\sum_{1}m_{i})^{2}-\sum_{1}m_{i}^{2}}{(\sum
m_{i})^{2}-\sum m_{i}^{2}},\\\
f^{W}_{22}&=\frac{(\sum_{2}m_{i})^{2}-\sum_{2}m_{i}^{2}}{(\sum m_{i})^{2}-\sum
m_{i}^{2}},\\\ \textrm{and}\quad
f^{W}_{12}&=\frac{\sum_{1}m_{i}\sum_{2}m_{i}}{(\sum m_{i})^{2}-\sum
m_{i}^{2}},\end{split}$ (83)

and also noting that $(\sum w_{i})^{2}-\sum
w_{i}^{2}=(\sum_{1}w_{i})^{2}-\sum_{1}w_{i}^{2}+(\sum_{2}w_{i})^{2}-\sum_{2}w_{i}^{2}+2\sum_{1}w_{i}\sum_{2}w_{i}$
we arrive at

$W=f^{W}_{11}\,W_{11}+f^{W}_{22}\,W_{22}+2f^{W}_{12}\,W_{12}.$ (84)

This shows that for specific marks the weighted correlation function can be
split up into a sum of auto-correlations and a cross-correlation. This can be
generalised if the mark e.g. takes three or more different values and the
result will include contributions from all possible auto- and cross-
correlations.

One particularly interesting case is the $\textrm{Void}_{\textrm{AC}}$ mark
where the two values the mark can take is simply $-1$ and $+1$. Let us assume
that the 1-population has $-1$ as a mark and the 2-population has $+1$. First
of all we realise that in that case $W_{1}W_{1}=D_{1}D_{1}$ as well as
$W_{2}W_{2}=D_{2}D_{2}$ because the sum of pair-product weights will simply be
a sum of 1’s. Furthermore $W_{1}W_{2}=-D_{1}D_{2}$ as the product of two
weights will always be -1 for pairs in the cross-correlation. The
normalisation also simplifies yielding
$(\sum_{1}w_{i})^{2}-\sum_{1}w_{i}^{2}=N_{1}(N_{1}-1)$ and analogous for the
2-population. For the cross-correlation we get
$\sum_{1}w_{i}\sum_{2}w_{i}=-N_{1}N_{2}$. Hence the individual auto- and
cross-correlations are the same between the weighted and unweighted case
$W_{11}=\xi_{11}$, $W_{22}=\xi_{22}$ and $W_{12}=\xi_{12}$. This implies that
the total weighted correlation function has the same individual contributions
as the unweighted one but with different prefactors

$W=f^{W}_{11}\,\xi_{11}+f^{W}_{22}\,\xi_{22}+2f^{W}_{12}\,\xi{12},$ (85)

since the prefactors simplify to

$\begin{split}f^{W}_{11}&=\frac{N_{1}(N_{1}-1)}{(N_{1}-N_{2})^{2}-N},\\\
f^{W}_{22}&=\frac{N_{2}(N_{2}-1)}{(N_{1}-N_{2})^{2}-N},\\\ \textrm{and}\quad
f^{W}_{12}&=-\frac{N_{1}N_{2}}{(N_{1}-N_{2})^{2}-N}.\end{split}$ (86)

If we define the constant $\mathcal{C}=((N_{1}+N_{2})^{2}-N)/(N(N-1))$ then we
can write the marked correlation function for the
$\textrm{Void}_{\textrm{AC}}$-mark as

$\mathcal{M}=\frac{1+\frac{f_{11}}{\mathcal{C}}\,\xi_{11}+\frac{f_{22}}{\mathcal{C}}\,\xi_{22}-2\frac{f_{12}}{\mathcal{C}}\,\xi_{12}}{1+f_{11}\,\xi_{11}+f_{22}\,\xi_{22}-2f_{12}\,\xi_{12}}.$
(87)

This illustrates the fact that in this case the marked correlation function is
nothing else than a specific combination of unweighted auto and cross-
correlations.

##  Appendix B Convolution of the density contrast

In this section we show how an additional convolution of the already smoothed
density contrast can be treated simply as a single convolution with a higher-
order smoothing kernel. Let us start with

$\delta_{RR}(\mathbf{x})=\frac{1}{a^{6}}\int_{\mathbf{x}^{\prime\prime}}\left[\int_{\mathbf{x}^{\prime}}F\left(\frac{\mathbf{x}^{\prime\prime}-\mathbf{x}^{\prime}}{a}\right)\delta(\mathbf{x}^{\prime})\,\text{d}^{3}x^{\prime}\right]G\left(\frac{\mathbf{x}-\mathbf{x}^{\prime\prime}}{a}\right)\,\text{d}^{3}x^{\prime\prime}$
(88)

where we can identify the smoothed density contrast inside the square brackets
and two smoothing kernels $F$ and $G$. This can be rewritten as an integral
over $\mathbf{x}^{\prime\prime}$ involving only the two kernels with a
coordinate transformation $\mathbf{y}=\mathbf{x}-\mathbf{x}^{\prime\prime}$

$\begin{split}\delta_{RR}(\mathbf{x})&=\frac{1}{a^{6}}\int_{\mathbf{x}^{\prime}}\left[\int_{\mathbf{x}^{\prime\prime}}F\left(\frac{\mathbf{x}-\mathbf{x}^{\prime}-\mathbf{y}}{a}\right)G\left(\frac{\mathbf{y}}{a}\right)\,\text{d}^{3}y\right]\delta(\mathbf{x}^{\prime})\,\text{d}^{3}x^{\prime}\\\
&=\frac{1}{a^{3}}\int_{\mathbf{x}^{\prime}}\delta(\mathbf{x}^{\prime})H\left(\frac{\mathbf{x}-\mathbf{x}^{\prime}}{a}\right)\,\text{d}^{3}x^{\prime},\end{split}$
(89)

where we identified in the first equality that after the coordinate change the
two kernels are convolved in the variable $\mathbf{y}$ resulting in a new
kernel $H$ at location $\mathbf{x}-\mathbf{x}^{\prime}$. Hence the density
field is only convolved once with the $H$-kernel that is the convolution of
both $G$ and $F$. Now, if the two kernels are a PCS and NGP kernel,
respectively, then the $H$-kernel would be a quartic kernel.