# Mechanisms for Magnetic Skyrmion Catalysis and Topological Superconductivity

Yun-Peng Huang<EMAIL_ADDRESS>Beijing National Laboratory for
Condensed Matter Physics, and Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  Panagiotis Kotetes<EMAIL_ADDRESS>CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China

###### Abstract

We propose an alternative route to stabilize magnetic skyrmion textures which
does not require Dzyaloshinkii-Moriya interaction, magnetic anisotropy, or an
external Zeeman field. Instead, it solely relies on the emergence of flux in
the system’s ground state. We discuss scenarios that lead to a nonzero flux,
and identify the magnetic skyrmion ground states which become accessible in
its presence. Moreover, we explore the chiral superconductors obtained for the
surface states of a topological crystalline insulator when two types of
magnetic skyrmion crystals coexist with a pairing gap. Our work opens
perspectives for engineering topological superconductivity in a minimal
fashion, and promises to unearth functional topological materials and devices
which may be more compatible with electrostatic control than the currently
explored skyrmion-Majorana platforms.

##  I Introduction

Magnetic textures are advantageous ingredients for engineering Majorana
quasiparticles [1, 2, 3, 4, 5, 6, 7, 8, 9], since they simultaneously break
time reversal symmetry (TRS) and induce Rashba spin-orbit coupling (SOC) [10,
11, 12]. One finds two distinct ways through which magnetic textures give rise
to these excitations. The first relies on magnetic textures crystals (MTCs),
which are magnetic textures that repeat periodically in space. The interplay
of superconductivity with MTCs of the helix and skyrmion [13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25] types gives rise to effective $p\,$-wave
topological superconductors (TSCs) in 1D [26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37] and 2D [30, 31, 38, 39, 40, 41, 42], respectively. Recent scanning
tunneling microscopy (STM) measurements in helical Fe chains deposited on top
of Re(0001) hint towards the presence of edge Majorana zero modes (MZMs) [43],
while the coexistence of magnetic skyrmions and superconductivity has also
been demonstrated in 2D Fe/Ir magnets on Re(0001) [44]. In the second
scenario, MZMs bind to magnetic skyrmion defects [45, 46, 47, 48, 49, 50, 51],
which are imposed in a chiral ferromagnet [52, 53, 54, 55] coupled to a
conventional superconductor. Important experimental advancements along these
lines were recently made in [IrFeCoPt]/Nb hybrids [56]. In particular,
magnetic skyrmions and (anti)vortices composites were experimentally realized
[57], thus setting the stage for the observation of MZMs in these systems in
the near future.

The above exciting experimental successes towards novel MZM platforms should
be also scrutinized with respect to their degree of functionality and their
potential role as the hardware for a topological quantum computer [58, 59], in
which, a large number of MZMs would need to be simultaneously manipulated and
braided [60]. On one hand, Fe/Ir/Re systems typically contain a large number
of domains and inhomogeneities which may render them prone to disorder, while
they are not amenable to electrostatic tuning. Instead, these systems can be
tailored using STM tips, whose number however is difficult to scale up. This
obstacle may lead to limitations when it comes to braiding multiple MZMs. In
contrast, [IrFeCoPt]/Nb hybrids appear to feature a higher degree of
functionality, since it is envisaged that MZMs will be manipulated using
magnetic skyrmion racetracks [61]. However, this platform presents a caveat of
different nature. Such a system harbors both vortices and antivortices, which
are induced by the applied magnetic field and the magnetic skyrmions,
respectively. The magnetic field is required to assist the nucleation of
skyrmions [55], which appear in a restricted window of the parameter space
spanned by the strengths of the field, the Dzyaloshinkii-Moriya interaction
(DMI) [62, 63], and the magnetic anisotropy. Therefore, here it is the
mechanism for skyrmions itself that sets constraints on the controllability
and reliability of the system as a potential quantum computing platform, since
the coexistence of different types of vortices leads to complex dynamics that
can hinder the implementation of braiding and thus give rise to an inherent
source of noise and decoherence.

In this Manuscript, we argue that these experimental hurdles can be surpassed
by considering routes for realizing magnetic skyrmions which are free from
DMI, magnetic anisotropy, and a Zeeman field. As we demonstrate, it is
sufficient to violate TRS in the system by means of a nonzero “flux”
$\vartheta$, which acts as a source for the skyrmion charge $C$. Here, the
flux enters as a Lagrange multiplier which shifts the intensive free energy of
the system $F$ according to $F\mapsto F-\vartheta{\cal C}$, where:

$\displaystyle{\cal C}=\frac{1}{4\pi}\int
d\bm{r}\phantom{.}\bm{M}(\bm{r})\cdot\big{[}\partial_{x}\bm{M}(\bm{r})\times\partial_{y}\bm{M}(\bm{r})\big{]}\,.$
(1)

In the above, $\bm{M}(\bm{r})$ corresponds to the magnetization profile of the
magnetic skyrmion crystal or defect. Since ${\cal C}$ and $C$ become nonzero
simultaneously, the term $-\vartheta{\cal C}$ promotes the formation of
magnetic skyrmions. For this reason, we term this mechanism as magnetic
skyrmion catalysis (MSC).

Notably, a number of the above preconditions are known to be fulfilled in
quantum Hall systems. Indeed, magnetic skyrmions have already been
theoretically predicted [64, 65] and experimentally observed [66] in
ferromagnetic quantum wells more than two decades ago. Here, however, we are
not interested in situations where a strong quantizing magnetic field is
present. Instead, we focus on weakly doped quasi-2D itinerant magnets with
energy bands which feature a nonzero Berry curvature [67]. Even more, as we
show here, MSC is also accessible in systems which contain band touching
points and thus a singular Berry curvature. In the latter case, we find that
chiral fluctuations promote magnetic skyrmions either in a spontaneous
fashion, or, due to the orbital coupling to an external magnetic field.
Notably, in this second scenario, MSC is mediated by a phenomenon akin to the
magnetic catalysis known from high-energy physics [68], but without the
formation of Landau levels.

The magnets of interest are also assumed to feature a pairing gap, which is
either inherent due to spontaneous Cooper pair formation, or, is induced by
means of proximity to a conventional superconductor [69, 70]. In either
situation, we show that the synergy of a magnetic skyrmion crystal and
superconductivity renders the system an effective $p+ip$ superconductor which
harbors chiral Majorana edge modes [5, 6]. As a consequence, such systems
further set the stage for the appearance of MZMs at point-like defects, such
as, in the cores of vortices introduced in the phase of the superconducting
[2, 3] or the magnetic texture [71] fields.

Our manuscript is organized as follows. In Sec. II we derive the expression of
the flux $\vartheta$ which is the coefficient of the cubic term of the free
energy which appears in the presence of the TRS violation. In Sec. III we
study the behavior of $\vartheta$ for two standard models, which describe a
Chern insulator and a generalized Dirac-type of model. Section IV proposes
possible routes to violate TRS and demonstrates that fluctuations promote the
emergence of MSC. In Sec. V we identify all the possible magnetic ground
states that dictate itinerant magnets with tetragonal symmetry when spin-
rotational symmetry is intact but TRS is broken. Afterwards, we show in Sec.
VI that MSC can be realized in a Dirac model with a winding of two units. Then
we move on and discuss in Sec. VII the implications that MSC can have on
engineering topological superconductivity. This work concludes in Sec. VIII
with a summary and a discussion of possible experimental platforms that appear
prominent to harbor MSC. Finally, further explanations and details regarding
our calculations are given in Appendices A-E.

##  II Cubic term in the Landau expansion and flux

We proceed with the exposition of the MSC. We consider that the system is free
of any kind of SOC and other sources of magnetism, hence, also preserving spin
rotational invariance. In addition, we assume that the magnetic instability is
driven by Fermi surface nesting, which in turn implies that the spin
susceptibility $\chi(\bm{q})$ peaks at the star of the nesting vector
$\bm{Q}$.

The band structure results from a matrix Hamiltonian
$\hat{{h}}(\bm{k})\mathds{1}_{\sigma}$ where $\mathds{1}_{\sigma}$ is the unit
matrix in spin space and $\bm{k}$ the wave vector. We develop a Landau theory
for the magnetic order parameter $\bm{M}(\bm{r})$, which is expressed in
energy units, and results from the mean-field decoupling of a Hubbard
interaction with strength $U$. The magnetic order parameter couples to the
electrons through an exchange term $\propto\bm{M}(\bm{r})\cdot\bm{\sigma}$,
where $\bm{\sigma}$ define the spin Pauli matrices.

Under the above conditions, the single-particle Hamiltonian in the magnetic
phase takes the following form:

$\displaystyle\hat{{h}}_{\bm{M}}(\bm{k}+\bm{q},\bm{k})=\hat{{h}}(\bm{k})\mathds{1}_{\sigma}(2\pi)^{2}\delta(\bm{q})+\bm{M}(\bm{q})\cdot\bm{\sigma}\,,$
(2)

where we used the Fourier transform:

$\displaystyle\bm{M}(\bm{q})=\int
d\bm{r}\phantom{.}e^{-i\bm{q}\cdot\bm{r}}\bm{M}(\bm{r})\,,$ (3)

and the Dirac delta function $\delta(\bm{q})$. The Hamiltonian is written in
the formalism of second quantization as:

$\displaystyle{h}_{\bm{M}}=\int\frac{d\bm{q}}{(2\pi)^{2}}\int
d\bm{k}\phantom{.}\bm{\psi}^{{\dagger}}(\bm{k}+\bm{q})\hat{{h}}_{\bm{M}}(\bm{k}+\bm{q},\bm{k})\bm{\psi}(\bm{k})\,,$
(4)

where
$\bm{\psi}^{\dagger}(\bm{k})=\big{(}\psi_{\uparrow}^{\dagger}(\bm{k}),\,\psi_{\downarrow}^{\dagger}(\bm{k})\big{)}$.
Here, $\psi_{\uparrow,\downarrow}^{\dagger}(\bm{k})$ creates an electron of
spin projection $\uparrow,\downarrow$ and wave vector $\bm{k}$.

To derive a Landau expansion in powers of $\bm{M}(\bm{q})$, we employ the
Matsubara formalism and integrate out the electronic degrees of freedom in a
standard fashion [72]. The lowest term in this expansion is the quadratic one
since, in contrast to the situation taking place in quantum Hall systems [64,
65, 66], here we assume that there is no net ferromagnetic moment in the
ground state when $\bm{M}(\bm{r})$ is absent. The expansion coefficient at
quadratic order $\propto|\bm{M}(\bm{q})|^{2}$ is given by $2/U-\chi(\bm{q})$
[38], and determines the leading magnetic instabilities.

TRS violation allows for a cubic term in the expansion that contains the term
${\cal
L}(\bm{q},\bm{p})=\bm{M}(-\bm{q}-\bm{p})\cdot\big{[}\bm{M}(\bm{q})\times\bm{M}(\bm{p})\big{]}$,
which relates to the skyrmion charge density. Similar terms have been
previously discussed in a phenomenological fashion for chiral spin liquids
[73, 74, 75, 76], and have been also shown to be inducible by a magnetic field
[77]. However, little attention has been paid so far to the emergence of this
term as a result of the topological properties of the band structure and its
interplay with fluctuations.

In the following, we provide an analytical expression for the coefficient
$\vartheta$ of the respective cubic free energy term, which is written in the
compact form:

$\displaystyle{F}^{(3)}=\frac{1}{4\pi}\int\frac{d\bm{q}}{(2\pi)^{2}}\int\frac{d\bm{p}}{(2\pi)^{2}}\phantom{.}\vartheta(\bm{q},\bm{p})\phantom{.}{\cal
L}(\bm{q},\bm{p})\,.$

In the above we inserted the flux density:

$\displaystyle\vartheta(\bm{q},\bm{p})$ $\displaystyle=$ $\displaystyle
i\int\frac{d\bm{k}}{3\pi}\phantom{.}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\Big{\\{}{\rm
tr}\big{[}\hat{{G}}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\hat{{G}}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{q})$
(5)
$\displaystyle\qquad\qquad\hat{{G}}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{q}-\bm{p})\big{]}-\bm{q}\leftrightarrow\bm{p}\Big{\\}}\,,$

which is expressed in terms of the Matsubara Green function in the nonmagnetic
phase:

$\displaystyle\hat{G}^{-1}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})={\color[rgb]{0,0,0}i\omega_{\nu}}-\hat{h}(\bm{k})\,.$

${\color[rgb]{0,0,0}\omega_{\nu}}$ denote fermionic Matsubara frequencies, $T$
is the temperature in energy units, and ${\rm tr}$ stands for the trace
operation in the matrix space spanned by $\hat{h}(\bm{k})$.

The cubic order term ${F}^{(3)}$ takes the form $-\vartheta{\cal C}$ when the
system is in the vicinity of a ferromagnetic instability. In this case, we
have $\bm{q},\bm{p}\approx\bm{0}$ and thus we expand Eq. (5) at lowest order
in $\bm{q}$ and $\bm{p}$, which leads to:

$\displaystyle\vartheta=2i\varepsilon_{znm}\int\frac{d\bm{k}}{3\pi}\phantom{.}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}{\rm
tr}\left[\hat{{G}}(k)\frac{\partial\hat{{G}}(k)}{\partial
k_{n}}\frac{\partial\hat{{G}}(k)}{\partial k_{m}}\right]$ (6)

where we set $k=({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})$ and introduced the
totally antisymmetric Levi-Civita symbol $\varepsilon_{znm}$, with
$n,m=\\{x,y\\}$. Repeated index summation is implied. Straightforward
calculations lead to one of the key results of this work:

$\displaystyle\vartheta=\sum_{s=0,1,2}\frac{4}{(s+1)!s!}\sum_{\alpha}\int\frac{d\bm{k}}{\pi}\phantom{.}i\varepsilon_{znm}\left<\partial_{k_{n}}\bm{u}_{\alpha}(\bm{k})\right|\big{[}\mathds{1}-\hat{{\cal
P}}_{\alpha}(\bm{k})\big{]}\big{[}\varepsilon_{\alpha}(\bm{k})-\hat{h}(\bm{k})\big{]}^{s-2}\left|\partial_{k_{m}}\bm{u}_{\alpha}(\bm{k})\right>\partial_{\mu}^{s}f[\varepsilon_{\alpha}(\bm{k})-\mu],$
(7)

where $\left|\bm{u}_{\alpha}(\bm{k})\right>$ are the eigenvectors of
$\hat{h}(\bm{k})$ with corresponding energy dispersions
$\varepsilon_{\alpha}(\bm{k})$ and projectors $\hat{{\cal
P}}_{\alpha}(\bm{k})=\left|\bm{u}_{\alpha}(\bm{k})\right>\left<\bm{u}_{\alpha}(\bm{k})\right|$.
$f(\epsilon-\mu)$ is the Fermi-Dirac distribution defined in the grand
canonical ensemble and evaluated at an energy $\epsilon$ for a chemical
potential $\mu$.

We remark that, while in the above $\vartheta$ is given by a sum of
contributions arising from each band separately, it still originates from
interband-only transitions, see App. A. In fact, this is reflected in the
presence of the projectors $\mathds{1}-\hat{{\cal P}}_{\alpha}(\bm{k})$. Note
also that for $s=2$, the term inbetween the eigenstates can be replaced by the
identity operator, and thus leads to the Berry curvature [67]

$\displaystyle\Omega_{\alpha}(\bm{k})=i\varepsilon_{znm}\big{<}\partial_{k_{n}}\bm{u}_{\alpha}(\bm{k})\big{|}\partial_{k_{m}}\bm{u}_{\alpha}(\bm{k})\big{>}\,.$
(8)

of the $\alpha$-th band. The above reflects the topological origin of the flux
and Eq. (7) reveals how $\vartheta$ relates to the band structure properties.

##  III Flux in two-band models

We evaluate $\vartheta$ for a two band model
$\hat{{h}}(\bm{k})=\bm{d}(\bm{k})\cdot\bm{\kappa}$, where $\bm{\kappa}$ denote
Pauli matrices that may relate to valley or similar degrees of freedom. We use
the energy dispersions $\varepsilon_{\pm}(\bm{k})=\pm\varepsilon(\bm{k})$ with
$\varepsilon(\bm{k})=|\bm{d}(\bm{k})|$, and obtain:

$\displaystyle\vartheta=\sum_{s}^{0,1,2}\sum_{\alpha}^{{\color[rgb]{0,0,0}\pm
1}}\int\frac{d\bm{k}}{\pi}\phantom{.}\frac{(-2)^{s}\alpha^{1+s}\Omega(\bm{k})}{(s+1)!s![\varepsilon(\bm{k})]^{2-s}}\phantom{.}\partial_{\mu}^{s}f[-\alpha\varepsilon(\bm{k})-\mu]$
(9)

where the quantity:

$\displaystyle\Omega(\bm{k})=\frac{1}{2}\phantom{.}\hat{\bm{d}}(\bm{k})\cdot\Big{[}\partial_{k_{x}}\hat{\bm{d}}(\bm{k})\times\partial_{k_{y}}\hat{\bm{d}}(\bm{k})\Big{]}\,$
(10)

denotes the Berry curvature of the valence band [67], with
$\hat{\bm{d}}(\bm{k})=\bm{d}(\bm{k})/|\bm{d}(\bm{k})|$. In the following, we
explore the behaviour of $\vartheta$ for two kinds of models. First, the case
of an extended “Dirac model” which exhibits a single $n$-th order band
touching point and, second, the Qi-Wu-Zhang model [78] which describes a Chern
insulator or metal.

###  III.1 Case of extended Dirac model

We first obtain $\vartheta$ for the generalized “Dirac model”:

$\displaystyle\bm{d}(\bm{k})=\big{(}\upsilon_{1}k^{\ell}\sin(\ell\phi),m,\upsilon_{2}k^{\ell}\cos(\ell\phi)\big{)},$
(11)

where we used polar coordinates $k_{x}=k\cos\phi$ and $k_{y}=k\sin\phi$. This
model yields an $\ell$-th order band touching point for $m=0$. At zero
temperature and equal “velocities” $\upsilon_{1,2}=\upsilon$, we find the
closed form result:

$\displaystyle\vartheta_{T=0}(m)=\frac{\ell}{3}\frac{{\rm
sgn}(m)}{m^{2}}\Big{[}\Theta\big{(}|m|-|\mu|\big{)}-|m|\delta\big{(}|m|-|\mu|\big{)}\Big{]},$
(12)

which is manifestly independent of $\upsilon$. $\Theta(\epsilon)$ denotes the
Heaviside unit step function at energy $\epsilon$. We observe that the flux
term is proportional to the vorticity of the band touching point. Quite
remarkably, we find that given these conditions, the flux is nonzero only when
the chemical potential lies within the band gap, or, when it exactly touches
the band gap edge, which yields a resonance condition arising from the delta
function.

This result appears at first sight discouraging, since the MSC seems not to be
accessible for itinerant magnets in which magnetism is driven by the presence
of a Fermi surface which exhibits a substantial degree of nesting. However,
the observation that the flux peaks for chemical potential values $|\mu|$
close to $|m|$, implies that introducing a nonzero temperature can unlock MSC,
since the temperature effectively broadens the energy levels.

Figure 1: Behavior of the coefficient $\vartheta$ as a function of temperature
$T$ and mass/gap parameter $m$. Results for the model of the generalized model
of Sec. III.1 are shown in (a), while calculations for the Chern insulator
model of Sec. III.2 are shown in (b). $\vartheta$ mostly exhibits a
qualitative similar behavior for the two models. However, a crucial difference
is that at $T=0$, $\vartheta$ is precisely zero for $|m|<|\mu|$ in the model
of Eq. (11), while such a restriction is not present for the Chern insulator
of Eq. (14) as shown in the inset of (b). For the calculations we used
$\mu=1$.

Indeed this is confirmed by the numerical evaluation of $\vartheta$ for
nonzero temperatures shown in Fig. 1(a). Therefore, when $T$ is comparable to
$|m|-|\mu|$, the system fulfills the resonance condition and leads to nonzero
flux. In the low temperature regime, the flux can be approximated by the
following form:

$\displaystyle\vartheta_{\rm
res}\approx-\frac{\ell}{3\pi}\frac{m}{\mu^{2}}\frac{T}{\big{(}|\mu|-|m|\big{)}^{2}+T^{2}}\,.$
(13)

The above nontrivial dependence of $\vartheta$ on the chemical potential
further implies that the magnetic skyrmions contribute to the electric charge
[64] since $\partial\vartheta/\partial\mu\neq 0$. However, at exactly zero
temperature the magnetic skyrmions do not carry electric charge since, in
contrast to Ref. 64, our starting point is a paramagnetic ground state and on
top of that, no SOC is present [79].

Besides the inclusion of a nonzero temperature, at this stage, it is equally
important to comment on the effects of velocity anisotropy. We find that the
system is practically insensitive to changing the Fermi surface properties.

###  III.2 Case of Chern insulator

In this pararaph we briefly discuss results for a two band Chern insulator
model described by the $\bm{d}(\bm{k})$ vector:

$\displaystyle\bm{d}(\bm{k})=\big{(}\sin k_{x},\sin k_{y},\cos k_{x}+\cos
k_{y}-m\big{)}$ (14)

with $\bm{k}\in[-\pi,\pi]^{2}$. The above model supports a topologically
nontrivial phases with a Chern number of a single unit $\propto{\rm sgn}(m)$
in the interval $|m|<2$.

The expression for $\vartheta$ is once again given by Eq. (9). Since obtaining
analytical results is tedious, we restrict to a numerical investigation that
reveals the salient features that characterize the behaviour of this system.
Our calculations shown in Fig. 1(b) confirm that a nonzero $\vartheta$ is
obtainable even for $|\mu|>|m|$, albeit with a relatively small magnitude.
Nevertheless, the coefficient can be further boosted by considering a low
nonzero temperature.

##  IV Mass induction mechanisms

So far, the origin of the flux $\vartheta$ and the mass term $m$ in Eq. (11)
have remained unspecified. The sole requirement to generate these, is that the
energy bands see a nonzero Berry curvature. Therefore, prominent candidates
appear to be narrow gap magnetic semiconductors. Preferably, these should be
of a low electronic density in order to be tunable via electrostatic gating.

MSC also arises in weakly doped semimetals which contain higher order band
touching points as in Eq. (11) for $m=0$. However, this is possible only in
the presence of fluctuations in the channel which gives rise to $m$. Indeed,
the presence of such chiral fluctuations enables the system to spontaneoussly
induce a nonzero $m$, since this allows for MSC to take place and, in turn,
minimize the free energy by stabilizing a magnetic skyrmion crystal ground
state. In the remainder of this subsection, we focus on $\ell$-th order band
touching points as in Eq. (11).

We make the above scenario plausible by adopting a phenomenological approach
along the lines of prior works discussing fluctuations mechanisms in different
contexts [80, 81, 82]. We focus on the mass sector of the free energy, which
takes the form:

$\displaystyle
F_{m}=\left(\frac{1}{V_{m}}-\chi_{m}\right)\frac{m^{2}}{2}-\vartheta(m){\cal
C}-{\cal M}_{z}(m)B_{z},$

where $\chi_{m}$ denotes the respective “mass-mass” susceptibility. $B_{z}$
denotes an out-of-plane magnetic field, which is sufficiently weak to consider
the Zeeman effect fully negligible. Thus, $B_{z}$ couples to the electrons via
the orbital effects. The latter is here introduced by means of the zero
temperature orbital magnetization ${\cal M}_{z}$ [67, 83], which is a function
of $m$, and is obtained from the expression:

$\displaystyle{\cal
M}_{z}(m)=\frac{\mu}{\Phi_{0}}\sum_{\alpha}\int\frac{d\bm{k}}{2\pi}\phantom{.}\Omega_{\alpha}(\bm{k})f\big{[}\varepsilon_{\alpha}(\bm{k})-\mu\big{]}\,,$

where we introduced the magnetic flux quantum $\Phi_{0}$. Note that here we
neglected the contribution of ${\cal C}$ to the orbital magnetization, since
${\cal C}$ is examined at the level of an instability analysis. For a two band
model of the form Eq. (11), the orbital magnetization reads:

$\displaystyle{\cal M}_{z}(m)=\frac{\ell
m\mu}{\Phi_{0}}\left[\frac{\Theta\big{(}|m|-|\mu|\big{)}}{|m|}+\frac{\Theta\big{(}|\mu|-|m|\big{)}}{|\mu|}\right]\,,$
(15)

which already includes a factor of $2$ due to spin.

In the remainder we assume low temperatures, as well as a value $|\mu|>|m|$
but for a chemical potential slightly above the conduction band bottom. These
conditions allow us to approximately write ${\cal M}_{z}(m)\approx{\rm
sgn}(\mu)\ell m/\Phi_{0}$. The same conditions allow us to obtain an analogous
expression for $\vartheta(m)$. Specifically, we assume a thermally activated
MSC and employ Eq. (13). After further simplifications, we consider
$\vartheta\approx-\ell m/3\pi\mu^{2}T$.

We now move on with identifying $m$. The interaction potential $V_{m}$ is non-
negative and thus attractive in the channel of $m$. Here, $m$ is considered
not to emerge, which is reflected in the relation $1/V_{m}-\chi_{m}>0$ assumed
to hold throughout. However, the situation changes when the second and/or the
third terms are considered. The value of $m$ induced by ${\cal C}$ and/or
$B_{z}$, that we here term $m_{\rm ind}$ is found by extremizing the free
energy, i.e., $dF_{m}/dm=0$, which yields:

$\displaystyle m_{\rm ind}=\frac{\ell V_{m}}{1-V_{m}\chi_{m}}\left[\frac{{\cal
C}}{3\pi\mu^{2}T}-{\rm sgn}(\mu)\frac{B_{z}}{\Phi_{0}}\right].$

We observe that even for ${\cal C}=0$, a nonzero $m$ becomes readily induced,
due to the applied magnetic field. This phenomenon is a manifestation of the
weak-field analog of magnetic catalysis [68, 84], and was previously discussed
for field-induced chiral superconductors [85] and density waves [86, 87, 80].

For the value $m_{\rm ind}$ the energy of the system is reduced. This is
because a nonzero $m$ emerges in spite of the system being detuned from an
instability in the mass channel. This can be more transparently demonstrated
by introducing the field $\tilde{m}$, which is conjugate to $m$ in the
statistical mechanics sense. Following Ref. [88], we find:

$\displaystyle\tilde{m}=\chi_{m}m+\ell\left[\frac{{\cal C}}{3\pi\mu^{2}T}-{\rm
sgn}(\mu)\frac{B_{z}}{\Phi_{0}}\right].$

The next step is to perform a Legendre transform which expresses the free
energy in terms of $\tilde{m}$. We thus find:

$\displaystyle F_{\tilde{m}}=\frac{\big{\\{}\tilde{m}-\ell\big{[}{\cal
C}/3\pi\mu^{2}T-{\rm
sgn}(\mu)B_{z}/\Phi_{0}\big{]}\big{\\}}^{2}}{2\chi_{m}}-V_{m}\frac{\tilde{m}^{2}}{2}\,.$

By plugging in the above the induced value of $\tilde{m}$ due to ${\cal C}$
and/or $B_{z}$, we find that the energy is reduced by the amount:

$\displaystyle\delta E=-\frac{\ell^{2}}{2}\frac{\big{[}{\cal
C}/3\pi\mu^{2}T-{\rm sgn}(\mu)B_{z}/\Phi_{0}\big{]}^{2}}{1/V_{m}-\chi_{m}}\,.$
(16)

The above also reveals that chiral fluctuations mediate a coupling between
${\cal C}$ and $B_{z}$, with the latter acting now as a source of magnetic
skyrmion density.

Conclusively, we find that chiral fluctuations favor a nonzero ${\cal C}$, as
they open up an additional route to reduce the free energy. Whether a nonzero
${\cal C}$ finally appears, will be decided by the remaining magnetic
contributions to the free energy.

##  V Magnetic ground states in the presence of flux

We now proceed with discussing aspects of the magnetic instability and phase
diagram for generic antiferromagnets. We assume a system with square lattice
symmetry which exhibits tendency to develop magnetic ground states at the star
of the ordering wave vectors $\pm\bm{Q}_{1,2}$, with $\bm{Q}_{1}=Q(1,0)$ and
$\bm{Q}_{2}=Q(0,1)$. In the presence of the nonzero flux, the cubic term also
enforces ordering at wave vectors
$\bm{Q}_{\pm}=\bm{Q}_{1}\pm\bm{Q}_{2}$.111Such a multi-$\bm{Q}$ magnetic
ground state can be alternatively stabilized by an external field instead of
flux, as it was recently experimentally observed in Ref. 89. Notably,
$\bm{Q}_{\pm}$ do not belong to the star $\pm\bm{Q}_{1,2}$. Hence, if
$\bm{Q}_{1,2}$ are to define the leading magnetic instability, the
susceptibility at $\bm{Q}_{\pm}$ is not expected to show an equally sharp
peak. Therefore, the respective magnetic orders
$\bm{M}_{\bm{Q}_{\pm}}\equiv\bm{M}_{\pm}$ would not undergo an instability
when considered alone222Note that in crystal lattices with trigonal and
hexagonal symmetry, the cubic term involves three magnetic order parameters
with all entering the magnetic skyrmion charge and determining the leading
instability..

Let us now discuss the above rather general implementation for a square
lattice in more detail. We focus on the relevant order parameter subspace and
write:

$\displaystyle
F_{\bm{M}_{1,2,\pm}}=F_{\bm{M}_{1,2}}+\bar{\alpha}\left(|\bm{M}_{+}|^{2}+|\bm{M}_{-}|^{2}\right)$
$\displaystyle-\vartheta_{1,2}\big{[}\bm{M}_{+}^{*}\cdot\big{(}\bm{M}_{1}\times\bm{M}_{2}\big{)}-\bm{M}_{-}^{*}\cdot\big{(}\bm{M}_{1}\times\bm{M}_{2}^{*}\big{)}+{\rm
c.c.}\big{]},$ (17)

where $\vartheta_{1,2}$ corresponds to $\vartheta(\bm{q},\bm{p})$ after being
evaluated for all permutations of the momenta
$\bm{q},\bm{p}=\big{\\{}\bm{Q}_{1},\bm{Q}_{2},\bm{Q}_{\pm}\big{\\}}$. In the
above, we introduced:

$\displaystyle F_{\bm{M}_{1,2}}$ $\displaystyle=$
$\displaystyle\alpha\big{(}|\bm{M}_{1}|^{2}+|\bm{M}_{2}|^{2}\big{)}+\frac{\tilde{\beta}}{2}\big{(}|\bm{M}_{1}|^{2}+|\bm{M}_{2}|^{2}\big{)}^{2}$
(18) $\displaystyle+$
$\displaystyle\frac{\beta-\tilde{\beta}}{2}\big{(}|\bm{M}_{1}^{2}|^{2}+|\bm{M}_{2}^{2}|^{2}\big{)}+(g-\tilde{\beta})|\bm{M}_{1}|^{2}|\bm{M}_{2}|^{2}$
$\displaystyle+$
$\displaystyle\frac{\tilde{g}}{2}\big{(}|\bm{M}_{1}\cdot\bm{M}_{2}|^{2}+|\bm{M}_{1}\cdot\bm{M}_{2}^{\ast}|^{2}\big{)}\,.$

Phases | Magnetic Order Parameters: $\bm{M}_{1}$, $\bm{M}_{2}$, $\bm{M}_{\pm}$ | Skyrmion Charge in a Magnetic Unit Cell
---|---|---
Sk-SVC |  $\bm{M}_{1}=M\big{(}1,0,0\big{)}$, $\bm{M}_{2}=M\big{(}0,1,0\big{)}$, $\bm{M}_{\pm}=\pm\tilde{M}\big{(}0,0,1\big{)}$ | $|C|=2,\qquad\tilde{M}\neq 0$
Sk-MS$\|$MH |  $\bm{M}_{1}=M\cos\eta\big{(}0,0,1\big{)}$, $\bm{M}_{2}=M\sin\eta\big{(}i\sin\lambda,0,\cos\lambda\big{)}$, $\bm{M}_{\pm}=\tilde{M}\big{(}0,i\sin\lambda,0\big{)}$ | $|C|=2,\qquad\tilde{M}\neq 0$
Sk-SWC4 | $\bm{M}_{1}=M\big{(}i\cos\lambda,0,\sin\lambda\big{)}$, $\bm{M}_{2}=M\big{(}0,i\cos\lambda,\sin\lambda\big{)}$ | $\begin{cases}\big{|}C\big{|}=1&|M|>|\tilde{M}\cot\lambda|\\\ \big{|}C\big{|}=2&|M|<|\tilde{M}\cot\lambda|\end{cases}$
$\bm{M}_{\pm}=\tilde{M}\big{(}i\sin\lambda,\pm
i\sin\lambda,\cos\lambda\big{)}$
Sk-SWC2 | $\bm{M}_{1}=M\cos\eta\big{(}1,0,i\big{)}$, $\bm{M}_{2}=\sqrt{2}M\sin\eta\big{(}\cos\lambda,i\sin\lambda,0\big{)}$ | $\begin{cases}\big{|}C\big{|}=1&\big{|}\tilde{M}\big{|}<|M||\cos\eta+\sqrt{2}\sin\eta\cos\lambda|\\\ \big{|}C\big{|}=2&\big{|}\tilde{M}\big{|}>|M||\cos\eta+\sqrt{2}\sin\eta\cos\lambda|\end{cases}$
$\bm{M}_{\pm}=\tilde{M}\big{(}1,\pm i\cot\lambda,i\big{)}/2$

Table 1: The four magnetic leading instabilities of Eq. (17) which possess a
nonzero skyrmion charge in the presence of a nonzero $\gamma$. The four
magnetic phases build upon the nonskyrmion analog phases which are obtained
for $\gamma=0$ and are solutions of Eq. (18). The magnetic phases Sk-MS$||$MH
and Sk-SWC2 are found only as local minima (metastable) solutions.

All the magnetic ground states of $F_{\bm{M}_{1,2}}$, have been analytically
and precisely identified in Ref. [90]. Notably, skyrmion type of crystals do
not belong to the thermodynamically stable magnetic ground states of
$F_{\bm{M}_{1,2}}$. However, among the accessible ground states of the free
energy in Eq. (18) one finds noncollinear and noncoplanar phases which can be
converted into skyrmionic textures crystals. These phases correspond to the
so-called spin vortex (SVC) [91] and spin whirl (SWC4) [90] crystal phases,
and are going to be of relevance in this work.

In contrast to Eq. (18), the free energy of Eq. (17) in principle allows
skyrmion phases by virtue of the cubic term which appears due to the TRS
violation. Indeed, we verify this by obtaining an analytical solution also
here, in spite of the two extra terms. In fact, an exact solution becomes
possible here because the order parameters $\bm{M}_{\pm}$ are essentially
driven by the magnetic instability at $\bm{Q}_{1,2}$, and $\bar{\alpha}$ is a
non-negative number. These two conditions allow us to integrate out
$\bm{M}_{\pm}$. The Euler-Lagrange equations of motion yield the constraints:

$\displaystyle\bm{M}_{+}=+\frac{\vartheta_{1,2}}{\bar{\alpha}}\bm{M}_{1}\times\bm{M}_{2}\phantom{{\dagger}}\,{\rm
and}\,\phantom{{\dagger}}\bm{M}_{-}=-\frac{\vartheta_{1,2}}{\bar{\alpha}}\bm{M}_{1}\times\bm{M}_{2}^{*}.$
(19)

Plugging the above back into Eq. (17), we find an effective free energy
depending only on $\bm{M}_{1,2}$:

$\displaystyle{\cal
F}_{\bm{M}_{1,2}}=F_{\bm{M}_{1,2}}-\frac{\gamma}{2}\left(|\bm{M}_{1}\times\bm{M}_{2}|^{2}+|\bm{M}_{1}\times\bm{M}_{2}^{*}|^{2}\right)\,,$

where we introduced the variable

$\displaystyle\gamma=\frac{2\vartheta_{1,2}^{2}}{\bar{\alpha}}\geq 0\,.$ (20)

One observes that the emergence of MSC in square lattice systems relies on a
fluctuations type of mechanism, which is mediated now by the noncritical order
parameters $\bm{M}_{\pm}$. This is a peculiarity of the given crystal symmetry
group, since it does not support ordering at three or more wave vectors which
add up to the null vector. In addition, it is worth noting that although
${\cal F}_{\bm{M}_{1,2}}$ alone respects TRS, magnetic skyrmion crystal ground
states are now possible thanks to the additional presence of $\bm{M}_{\pm}$.

To show this, we minimize ${\cal F}_{\bm{M}_{1,2}}$ with respect
$\bm{M}_{1,2}$. This task is carried out in App. B. Remarkably, we find that
the types of accessible magnetic ground states obtained in Ref. [90] persist
even after switching on $\gamma$. Moreover, the type of magnetic ground state
which minimizes ${\cal F}_{\bm{M}_{1,2}}$, can be obtained from the magnetic
phase diagram of $F_{\bm{M}_{1,2}}$, by employing the renormalized
coefficients:

$\displaystyle g\mapsto g-\gamma\quad{\rm
and}\quad\tilde{g}\mapsto\tilde{g}+\gamma\,.$ (21)

A proof for the above mapping is given in App. B. With $\bm{M}_{1,2}$ at hand,
we obtain $\bm{M}_{\pm}$ using Eq. (19) and, after a Fourier transform, we
determine $\bm{M}(\bm{r})$. Appendix C presents the order parameters and
respective magnetization profiles of all the double-$\bm{Q}$ magnetic ground
states of the free energy in Eq. (17).

Table 1 discusses the properties of the four magnetic ground states of
$F_{\bm{M}_{1,2,\pm}}$, which feature a nonzero skyrmion charge for
$\gamma\neq 0$. These are related to the following four respective ground
states of Eq. (18), i.e., the SVC, the SWC4, as well as the MS$||$MH and the
SWC2. The difference between SWC4 and SWC2 is that the former respects and the
latter violates fourfold rotational symmetry. Note that SWC2,4 generate
noncoplanar magnetic profiles, while the SVC and MS$||$MH collinear. Evenmore,
we remark that only the SVC and SWC4 constitute thermodynamically stable
magnetic ground states of the free energy in Eq. (18). As a result, also Sk-
MS$||$MH and Sk-SWC2 have a metastable character.

We find that starting from a SVC (SWC4) ground state for $\gamma=0$ with
$\bm{M}_{\pm}=\bm{0}$, can induce a Sk-SVC (Sk-SWC4) phase with $|C|=2$
($|C|=1,2$) when $\gamma$ is switched on. This is because $\bm{M}_{\pm}$ are
now nonzero and modify $\bm{M}(\bm{r})$.

Figure 2: Trajectories of the leading instabilities upon varying the mass term
$m$ which is added to the model of Eq. (22). The dashed line indicates the
trajectory for $\gamma=0$, i.e., without including it for the minimization of
the free energy. Instead, the solid line is obtained after including the
effects of a nonzero $\gamma$ which stems from a nonzero $m$. The two panels
show two general Fermi surface shapes which are obtained when varying the
Schrödinger masses $m_{1,2}$. For the Fermi surface in panel (a) [(b)], we
have $m_{2}=0.3m_{1}$ [$m_{2}=0.7m_{1}$]. In all cases, the chemical potential
is accordingly tuned so that the Fermi surface exhibits nesting with roughly
the same $\bm{Q}_{1,2}$ throughout each phase diagram. When $m=0$, the CSDW
(SVC) phase constitutes the leading instability in the (a) [(b)] case.
Remarkably, for a nonzero $\gamma$ the original trajectory is pushed towards
the Sk-SVC phase which now replaces the SVC phase which is obtained only for
$m=0$. The axes of the phase diagram are given in terms of the variables
$G=\big{(}g-\gamma-\tilde{\beta}\big{)}/|\beta-\tilde{\beta}|$ and
$\tilde{G}=\big{(}\tilde{g}+\gamma\big{)}/|\beta-\tilde{\beta}|$. The values
of the coefficients employed to draw the phase diagram are given in App. D.

##  VI Magnetism on the surface of a topological crystalline insulator

In this section we demonstrate the mechanism of MSC using as our starting
point a concrete 2D model:

$\displaystyle\hat{h}_{0}(\bm{k})=\frac{k_{x}k_{y}}{m_{1}}\kappa_{1}+\frac{k_{x}^{2}-k_{y}^{2}}{2m_{2}}\kappa_{3}\,,$
(22)

which leads to a single quadratic band touching point. Among other possible
applications, this model also describes the protected surface states of a bulk
topological crystalline insulator (TCI) [92].

To explore the emergence of MSC for the present system, we phenomenologically
introduce a mass term $m\kappa_{2}$ to the Hamiltonian of Eq. (22), thus
obtaining the model of Eq. (11) for $\ell=2$. The mass term can be induced
through the mechanisms discussed in Sec. IV and its origin is not further
specified in the following. We further assume that the chemical potential lies
energetically outside the gap opened by the mass term, therefore leading to a
Fermi surface. For instance, in the panels (a) and (b) of Fig. 2 we show two
distinct type of Fermi surfaces which generally emerge for the TCI surface
states. These are obtained by means of varying the chemical potential and/or
the Schrödinger masses $m_{1,2}$. For suitable parameter values, the
conduction band is dictated by well-nested Fermi segments. In the presence of
a Hubbard interaction which is diagonal in valley $\kappa$ space, nesting
promotes magnetic ordering at $\bm{Q}_{1,2}$.

To identify the accessible magnetic ground states for this system, we evaluate
the Landau coefficients of the free energy in Eq. (17). We follow a standard
approach that we detail in App. D. We find that a Fermi surface of the type
shown in Fig. 2(a), favors the stabilization of the collinear so-called
charge-spin density wave (CSDW) phase [91]. In contrast, a Fermi surface of
type shown in Fig. 2(b) promotes the establishment of the SVC phase. Such a
tendency is in accordance with previous results for models of Fe-based systems
[93]. There, it was shown that Fermi surfaces and interactions which are
featureless in valley space, such as in Fig. 2(b), promote the SVC over the
CSDW phase.

These results hold both in the presence and absence of $m$. However, in the
presence of the latter additional phenomena take place. When the system is
originally in the SVC phase, the noncritical $\bm{M}_{\pm}$ components are
also generated. Hence, a nonzero $m$ converts the SVC phase into the Sk-SVC
phase with $|C|=2$. In stark contrast, the MSC does not manifest itself in the
same fashion when the starting point is the CSDW phase. This is because the
order parameters $\bm{M}_{1,2}$, which miminize the free energy in Eq. (18),
are parallel in spin space. Nonetheless, for a system with CSDW instability at
$m=0$, switching on $m$ can allow for a first order transition to the Sk-SVC
phase with $|C|=2$.

Numerical results confirming the above are shown in Fig. 2. The values for the
related coefficients of the free energy in Eq. (17) are given in App. D.

##  VII Topological superconductivity

In this section, we continue with the investigation of the types of
topological phases which become accessible in the event of coexistence of the
Sk-SVC and Sk-SWC4 ground states with spin-singlet $s\,$-wave
superconductivity. Here we have in mind situations of hybrid systems where,
e.g., the surface of a bulk TCI in proximity a conventional superconductor.
Nonetheless, our results can also find applicability to individual material
candidates, which exhibit the microscopic coexistence of magnetism and
superconductivity. Note that, while in this work we restrict to conventional
pairing terms, unconventional pairing gaps also lead to topological phases in
conjunction with noncollinear and noncoplanar magnetism.

Figure 3: (a) Example of a nested Fermi surface for the model of Eq. (22) with
$m_{1}=1$, $m_{2}=0.7$, $m=0.038$ and $\mu=0.082$. We find that for these
parameter values the leading magnetic instability is of the Sk-SVC type. (b)
Fermi surface after downfolding to the magnetic Brillouin zone (MBZ) for the
nesting vectors $\bm{Q}_{1,2}$. Note that the ${\rm\Gamma}$ point of the
original Brillouin zone becomes the M point in the MBZ. On the other hand, the
points X and Y of the MBZ originate from points of the Fermi surface
experiencing single-$\bm{Q}$ nesting with $\bm{Q}_{1,2}$, respectively.
Finally, the ${\rm\Gamma}$ point of the MBZ (orange dot) stems from a point
which lies energetically away from zero, and experiences nesting with
$\bm{Q}_{\pm}$.

A detailed exploration of the types of topological phases under various
configurations of magnetic texture crystals and (un)conventional
superconductivity has been presented in Ref. [42]. There, it was found that
the symmetries dictating the magnetic profile and the superconducting gap are
pivotal for stabilizing fully-gapped and nodal TSCs, which support a variety
of strong, weak and crystalline phases. In the same work, the topological
properties of a SWC4 phase in coexistence with conventional superconductivity
were investigated. For a single band model, the symmetries dictating the SWC4
magnetization profile impose that only nodal topological superconductivity is
accessible in such a case. In the presence of an edge, these nodal TSCs harbor
so-called bidirectional Majorana modes, which are dispersive modes that do not
have a fixed sign for their group velocity. The topological mechanism that
leads to these modes has a one-dimensional origin and the modes are protected
by weak $\mathbb{Z}_{2}$ topological invariants [42] which are linked to band
inversions at the X and Y points of the magnetic Brillouin zone (MBZ).

Transitions to fully-gapped topological phases become possible in the SWC4
phase only after certain types of perturbations are added [42]. As we show in
the following paragraphs, a nonzero $\gamma$ which stabilizes the Sk-SWC4
leads to a gap in the spectrum and engineers a chiral TSC harboring a number
of $|C|=1,2$ branches of chiral Majorana modes per given edge. A similar
analysis for the Sk-SVC phase leads to a chiral TSC with $|C|=2$ number of
chiral Majorana modes per edge. Remarkably, the number of the arising chiral
edge modes is given by the respective skyrmion charge of the magnetic phase.
These results qualitatively hold for any single band model that respects the
square lattice symmetries and features the same type of nesting.

In the next paragraphs, we confirm the above mentioned topological phases
stemming from the Sk-SVC and Sk-SWC4 magnetic orders for the model of Eq.
(22), with the remaining parameter values shown in Fig. 3. Note that, while
the Sk-SWC4 magnetic order does not emerge as the leading magnetic instability
for this model, the results of our investigation have a general character and
can be applicable to systems that actually harbor this magnetic ground state.

###  VII.1 Bogoliubov - de Gennes formalism

The description of topological superconductivity from magnetic texture
crystals requires folding the energy dispersions of the Hamiltonian in Eq.
(22) down to the MBZ [42]. For a continuum model, this process generates a
band structure with an infinite number of bands. Nonetheless, only the bands
crossing the Fermi level are important for inferring the topological
properties of the systems under consideration. This holds as long as the
magnetic and pairing gaps are sufficiently weak to allow us to restrict to
only the bands contributing to the Fermi surface. In the remainder, we assume
that this assumption holds.

Within this framework, the respective many-particle Hamiltonian takes the form
$H_{\rm BdG}=\frac{1}{2}\int d\bm{k}\phantom{.}\bm{{\cal
X}}^{\dagger}(\bm{k})\hat{H}_{\rm BdG}(\bm{k})\bm{{\cal X}}(\bm{k})$, where we
introduced the low energy Bogoliubov - de Gennes (BdG) Hamiltonian:

$\displaystyle\hat{H}_{\rm
BdG}(\bm{k})=\left(\begin{array}[]{cccc}\hat{H}_{\rm
sc}(\bm{k}-\bm{q}_{+})&\bm{M}_{1}\cdot\bm{\sigma}&\bm{M}_{2}\cdot\bm{\sigma}&\bm{M}_{+}\cdot\bm{\sigma}\\\
\bm{M}_{1}^{*}\cdot\bm{\sigma}&\hat{H}_{\rm
sc}(\bm{k}+\bm{q}_{-})&\bm{M}_{-}^{*}\cdot\bm{\sigma}&\bm{M}_{2}\cdot\bm{\sigma}\\\
\bm{M}_{2}^{*}\cdot\bm{\sigma}&\bm{M}_{-}\cdot\bm{\sigma}&\hat{H}_{\rm
sc}(\bm{k}-\bm{q}_{-})&\bm{M}_{1}\cdot\bm{\sigma}\\\
\bm{M}_{+}^{*}\cdot\bm{\sigma}&\bm{M}_{2}^{*}\cdot\bm{\sigma}&\bm{M}_{1}^{*}\cdot\bm{\sigma}&\hat{H}_{\rm
sc}(\bm{k}+\bm{q}_{+})\end{array}\right)\,,$ (27)

which is defined in terms of the spinor:

$\displaystyle\bm{{\cal
X}}^{\dagger}(\bm{k})=\Big{(}\bm{\Psi}^{\dagger}(\bm{k}-\bm{q}_{+}),\,\bm{\Psi}^{\dagger}(\bm{k}+\bm{q}_{-}),\,\bm{\Psi}^{\dagger}(\bm{k}-\bm{q}_{-}),\,\bm{\Psi}^{\dagger}(\bm{k}+\bm{q}_{+})\Big{)},\,$
(28)

with $\bm{q}_{\pm}=\bm{Q}_{\pm}/2$. $\bm{{\cal X}}(\bm{k})$ is in turn defined
in terms of the subspinor:

$\displaystyle\bm{\Psi}^{\dagger}(\bm{k})=\Big{(}\bm{\chi}_{\uparrow}^{\dagger}(\bm{k}),\,\,\bm{\chi}_{\downarrow}^{\dagger}(\bm{k}),\,\,\bm{\chi}_{\downarrow}(-\bm{k}),\,\,-\bm{\chi}_{\uparrow}(-\bm{k})\Big{)}\,,$
(29)

where
$\bm{\chi}_{\uparrow,\downarrow}^{\dagger}(\bm{k})/\bm{\chi}_{\uparrow,\downarrow}(\bm{k})$
are creation/annihilation operators of electrons with $\uparrow,\downarrow$
spin projection. Furthermore, for a fixed spin, $\bm{\chi}$ consists of two
components which correspond to the valley degree of freedom $\kappa$ which
leads to the valence and conduction bands of the TCI surface states of the
Hamiltonian in Eq. (22). With no loss generality, the magnetic order is
assumed to be diagonal in valley space.

In the general presence of a mass term $m$, the Hamiltonian $\hat{H}_{\rm
sc}(\bm{k})$ describes the arising superconductor in the nonmagnetic phase,
and is given as:

$\displaystyle\hat{H}_{\rm
sc}(\bm{k})=\big{[}\hat{h}_{0}(\bm{k})-\mu\big{]}\tau_{3}+m\kappa_{2}+\Delta\tau_{1}\,,$
(30)

where we assumed that the surface states feel a uniform pairing gap $\Delta$
induced by means of proximity to a conventional superconductor. The
Hamiltonian above is expressed in terms of the Nambu Pauli matrices
$\bm{\tau}$. The latter are further supplemented by the respective identity
matrix $\mathds{1}_{\tau}$. Note that we omit writing unit matrices and
Kronecker product symbols throughout.

The Hamiltonian $\hat{H}_{\rm sc}(\bm{k})$ belongs to symmetry class C since
it is dictated by a charge conjugation symmetry
$\hat{\Xi}_{C}^{\dagger}\hat{H}_{\rm sc}(\bm{k})\hat{\Xi}_{C}=-\hat{H}_{\rm
sc}(-\bm{k})$, which is effected by the operator
$\hat{\Xi}_{C}=\tau_{2}\hat{{\cal K}}$. $\hat{{\cal K}}$ defines complex
conjugation. Superconductors in class C exhibit topologically nontrivial
properties in 2D which are classified by a $\mathbb{Z}$ topological index [94,
95, 96, 97, 98]. In the present situation, the emergence of topologically
nontrivial behaviour is linked to the band inversion at $\bm{k}=\bm{0}$, which
takes place when the criterion $|m|=\sqrt{\mu^{2}+\Delta^{2}}$ is satisfied.

Figure 4: (a) Comparison of the bare band structures ($\Delta=0$) for the full
and approximate models of Eq. (30) (solid black line) and (31) (red dashed
line), after downfolding to the magnetic Brillouin zone for the nesting
vectors $\bm{Q}_{1,2}$. For the model of Eq. (30) we used the same parameter
values employed in Fig. (3). The same values of $m_{1,2}$ were also considered
for the model in Eq. (31). However, in the latter case, we reduced the
chemical potential by $\approx 11\%$, in order to recover the same nesting
properties and be in a position to compare the two band structures. Note that
the full model yields twice as many bands, due to the additional presence of
the valence band. As we show, the two band structures agree quite well near
the Fermi level. As the energy is increased, one finds that there is a
discrepancy in the vicinity of the M point, which corresponds to the
$\bm{k}=\bm{0}$ point in the unfolded scheme. Therefore, as long as the
pairing and magnetic gaps to be added are safely smaller than $|\mu|-|m|$, the
approximate model describes quite accurately the band structure. (b) Here we
depict the bare band structure obtained only from the model of Eq. (31) after
downfolding to its respective MBZ. (c) In this panel we consider the effect
due to the addition of a conventional superconducting gap on the band
structure of panel (b). In all panels the energy scales are given in units of
the chemical potential value $\mu=0.082$.

###  VII.2 Projection onto the conduction band of the TCI surface states

The above topological properties are however not relevant for the cases of
interest. Here, the Fermi level crosses the conduction band and the desired
topological properties arise from energies near the Fermi level. Hence, under
the assumption that the gaps induced by the magnetic and pairing terms are
much smaller than the energy difference $|\mu|-|m|$, it is eligible to project
onto the conduction band of TCI surface states, and even fully discard the
presence of $m$ in the resulting Hamiltonian since it leads to a minor
modification of the Fermi surface shape.

The legitimacy of this approach is confirmed by comparing the resulting band
structure in the MBZ of the full model and the model which only restricts to
the conduction band and $m\kappa_{2}$ is dropped. As one observes in Fig. 4
for $\Delta=0$, the main difference between the full and projected model
arises near the M point of the MBZ, which corresponds to the $\bm{k}=\bm{0}$
point of the unfolded $\bm{k}$ space. Nevertheless, under the weak-coupling
assumption for the pairing and magnetic gaps, no band inversion takes place at
M, which lies sufficiently high in energy. Therefore, using the projected
model allows to qualitatively capture all the features which arise from the
interplay of superconductivity and the magnetic skyrmion texture crystal near
the Fermi level.

Given the above observations, in the following we approximate $\hat{H}_{\rm
sc}(\bm{k})$ of Eq. (30), according to:

$\displaystyle\hat{H}_{\rm
sc}(\bm{k})\approx\varepsilon(\bm{k})\tau_{3}+\Delta\tau_{1}$ (31)

which involves the energy dispersion:

$\displaystyle\varepsilon(\bm{k})=\sqrt{\big{(}k_{x}k_{y}/m_{1}\big{)}^{2}+\big{[}\big{(}k_{x}^{2}-k_{y}^{2}\big{)}/2m_{2}\big{]}^{2}}-\mu\,,$
(32)

where we also incorporated the chemical potential $\mu$ in it for convenience.
Note, that the approximate model in Eq. (31) sees new nesting vectors, which
have slightly modified lengths compared to the ones of the full model. Since
our analysis has mainly a qualitative character, we perform our upcoming study
in the MBZ defined for the nesting vectors of Eq. (31). Figure 4(c) depicts
the low-energy band structure in the MBZ obtained from the approximate
Hamiltonian in Eq. (31).

To this end, we remark that dropping $m\kappa_{2}$ from the Hamiltonian in
order to facilitate the exploration of TSC phases stemming solely from the
conduction band, does not at all imply that $m$ is fully discarded. The
presence of $m$ is still captured in an indirect fashion. This is because we
account for the additional emergence of the magnetic order components
$\bm{M}_{\pm}$ which appear only when $\vartheta$ is present. We remind the
reader that the latter is induced by a nonzero $m$ via the occurence of MSC.

###  VII.3 Compact Representation of the BdG Hamiltonian

The Hamiltonian in Eq. (31), as well as the one in Eq. 27, can be both written
in a compact manner, which exposes their symmetries in a more transparent way.
This is achieved by first introducing the Pauli matrices $\lambda_{1,2,3}$ and
$\rho_{1,2,3}$ related to foldings in the $k_{y}$ and $k_{x}$ directions,
respectively. These matrices correspondingly act in
$\left\\{\bm{k},\bm{k}+\bm{Q}_{2}\right\\}$ and
$\left\\{\bm{k},\bm{k}+\bm{Q}_{1}\right\\}$ spaces. Subsequently we consider
the unitary transformation:

$\displaystyle\hat{H}_{\rm BdG}^{\prime}(\bm{k})=\hat{{\cal
U}}_{\lambda}^{\dagger}\hat{{\cal U}}_{\rho}^{\dagger}\hat{H}_{\rm
BdG}(\bm{k})\hat{{\cal U}}_{\rho}\hat{{\cal U}}_{\lambda},$ (33)

where $\hat{\cal U}_{\zeta}=(\zeta_{3}+\zeta_{2})/\sqrt{2}$. Using the above,
the downfolded Hamiltonian of Eq. (31) takes the form:

$\displaystyle\hat{H}_{{\rm BdG;sc}}^{\prime}(\bm{k})$ $\displaystyle=$
$\displaystyle\varepsilon_{s}(\bm{k})\tau_{3}+\varepsilon_{x}(\bm{k})\rho_{2}\tau_{3}+\varepsilon_{y}(\bm{k})\lambda_{2}\tau_{3}$
(34) $\displaystyle+$
$\displaystyle\varepsilon_{xy}(\bm{k})\lambda_{2}\rho_{2}\tau_{3}+\Delta\tau_{1}\,,$

where we introduced the linear combinations:

$\displaystyle\varepsilon_{s}(\bm{k})$ $\displaystyle=$
$\displaystyle\frac{\varepsilon(\bm{k}-\bm{q}_{+})+\varepsilon(\bm{k}+\bm{q}_{-})+\varepsilon(\bm{k}-\bm{q}_{-})+\varepsilon(\bm{k}+\bm{q}_{+})}{4}\,$
$\displaystyle\varepsilon_{x}(\bm{k})$ $\displaystyle=$
$\displaystyle\frac{\varepsilon(\bm{k}-\bm{q}_{+})-\varepsilon(\bm{k}+\bm{q}_{-})+\varepsilon(\bm{k}-\bm{q}_{-})-\varepsilon(\bm{k}+\bm{q}_{+})}{4}\,$
$\displaystyle\varepsilon_{y}(\bm{k})$ $\displaystyle=$
$\displaystyle\frac{\varepsilon(\bm{k}-\bm{q}_{+})+\varepsilon(\bm{k}+\bm{q}_{-})-\varepsilon(\bm{k}-\bm{q}_{-})-\varepsilon(\bm{k}+\bm{q}_{+})}{4}\,$
$\displaystyle\varepsilon_{xy}(\bm{k})$ $\displaystyle=$
$\displaystyle\frac{\varepsilon(\bm{k}-\bm{q}_{+})-\varepsilon(\bm{k}+\bm{q}_{-})-\varepsilon(\bm{k}-\bm{q}_{-})+\varepsilon(\bm{k}+\bm{q}_{+})}{4}.$

The indices appearing above are chosen after atomic orbitals, in order to
reflect the behavior of the above functions under the inversions
$k_{x}\mapsto-k_{x}$ and/or $k_{y}\mapsto-k_{y}$.

Figure 5: Energy band structure for a Sk-SVC coexisting with a conventional
pairing gap $\Delta=0.25$, given the bare energy dispersion shown in Fig.
4(b). The latter is obtained using the parameters for $m_{1,2}$ and $\mu$
discussed in Fig. 3. Red (black) color denotes dispersions in the presence
(absence) of the $\bm{M}_{\pm}$ magnetic order parameters on top of a SVC
phase. (a)-(d) and (e)-(h) show the full (top row) and a zoom-in (bottom row)
of the band structure, respectively. The Sk-SVC leads to a fully-gapped
spectrum with a Hamiltonian which belongs to class D$\oplus$D due to a unitary
symmetry. The spectra of the two blocks are identical and allow for phases
with an even number of chiral Majorana modes per edge. In all panels the
energy scales are given in units of the chemical potential value $\mu=0.082$.
Finally, we assumed the ratio $\tilde{M}=M/4$.

The Hamiltonian in Eq. (34) is invariant under spin rotations and commutes
with $\lambda_{2}\rho_{2}$. Moreover, it is invariant under complex
conjugation $\hat{{\cal K}}$ and usual time reversal operation $\hat{{\cal
T}}=i\sigma_{y}\hat{{\cal K}}$, while it possesses a chiral symmetry effected
by $\hat{\Pi}=\tau_{2}$ and a charge conjugation symmetry
$\hat{\Xi}_{D}=\tau_{2}\sigma_{y}\hat{{\cal K}}$. Note that additional
antiunitary symmetries can be defined by virtue of the unitary symmetry [30].
As we discuss in the next section, a subset of these symmetries persist
inspite of the addition of magnetism.

###  VII.4 Topological phases: Sk-SVC ground state

For a Sk-SVC phase, the magnetic order parameters take the form
$\bm{M}_{1}=M(1,0,0)$, $\bm{M}_{2}=M(0,1,0)$, and
$\bm{M}_{\pm}=\pm\tilde{M}(0,0,1)$. This leads to the BdG Hamiltonian in the
new frame:

$\displaystyle\hat{H}_{\rm BdG}^{\prime}(\bm{k})=\hat{H}_{{\rm
BdG;sc}}^{\prime}(\bm{k})-M\rho_{1}\sigma_{x}-M\lambda_{1}\sigma_{y}-\tilde{M}\lambda_{3}\rho_{3}\sigma_{z}\,.$
(35)

The above Hamiltonian preserves the earlier found antiunitary charge
conjugation symmetry effected by the operator
$\hat{\Xi}_{D}=\tau_{2}\sigma_{y}\hat{{\cal K}}$, whose presence ensures the
emergence of Majorana excitations. Moreover, Eq. (35) features a unitary
symmetry generated by the operator ${\cal O}=\lambda_{2}\rho_{2}\sigma_{z}$.
Its presence allows us to block diagonalize the Hamiltonian by means of an
additional unitary transformation which renders ${\cal O}$ diagonal. We thus
consider the unitary transformation effected by ${\cal S}=\big{(}{\cal
O}+\sigma_{x})/\sqrt{2}$ and obtain the two-block-diagonal Hamiltonian:

$\displaystyle\hat{H}_{{\rm BdG}}^{\prime\prime}(\bm{k})=\hat{H}_{{\rm
BdG;sc}}^{\prime}(\bm{k})-M\rho_{1}\sigma_{x}+M\lambda_{3}\rho_{2}+\tilde{M}\lambda_{1}\rho_{1}\sigma_{x}\,.$
(36)

An additional unitary transformation effected by the operator ${\cal
S}^{\prime}=\big{(}\rho_{3}+\rho_{1}\sigma_{x})/\sqrt{2}$ decomposes the above
Hamiltonian into the two identical blocks:

$\displaystyle\hat{h}_{{\rm BdG},\sigma}(\bm{k})=\hat{h}_{{\rm
BdG;sc}}^{\prime}(\bm{k})-M\rho_{3}-M\lambda_{3}\rho_{2}+\tilde{M}\lambda_{1}\rho_{3}\,,$
(37)

Figure 6: Energy dispersions for a Sk-SWC4 phase coexisting with a
conventional pairing gap $\Delta=0.25$. given the bare energy dispersion shown
in Fig. 4(b). The latter is obtained using the parameters for $m_{1,2}$ and
$\mu$ discussed in Fig. 3. Red (black) color denotes dispersions in the
presence (absence) of the $\bm{M}_{\pm}$ magnetic order parameters on top of a
SVC phase. (a)-(d) and (e)-(h) show the full (top row) and a zoom-in (bottom
row) of the band structure, respectively. The Sk-SWC4 leads to a fully-gapped
spectrum with the respective Hamiltonian belonging to class D. Note that, in
contrast to the Sk-SVC case, here there is no additional twofold degeneracy.
For the calculations we considered the Sk-SWC4 structure shown in Table 1,
with $\lambda=\pi/4$ and $\tilde{M}=M/4$. Finally, note that in all panels the
energy scales are given in units of the chemical potential value $\mu=0.082$.

which are labelled by the quantum number $\sigma$, that behaves as spin in the
rotated frame. $\hat{h}_{{\rm BdG;sc}}^{\prime}(\bm{k})$ is identical to
$\hat{H}_{{\rm BdG;sc}}^{\prime}(\bm{k})$ but restricted to a single $\sigma$
block. The appearance of two identical blocks can be viewed as an emergent
SO(2) spin rotational invariance in the new frame.

In Fig. 5, we present the energy bands resulting from the Hamiltonian in Eq.
(35). We indeed confirm that the Sk-SVC phase leads to a fully-gapped twofold
degenerate spectrum. The gap is opened thanks to the presence of
$\bm{M}_{\pm}$ which are induced due the violation of TRS.

Gap closings and reopenings are also accessible, and these happen only at the
${\rm\Gamma}$ point. The related band inversions stabilize chiral TSC phases.
Indeed, we find that each one of these two blocks is dictated by a charge
conjugation symmetry generated by $\hat{\Xi}=\rho_{2}\tau_{2}\hat{{\cal K}}$,
which renders each Hamiltonian in class D, and allows for chiral TSCs with an
integer topological invariant ${\cal N}$. Since the two blocks are identical
we in fact have ${\cal N}\in 2\mathbb{Z}$, and an even number of chiral
Majorana modes are expected to appear per edge.

To demonstrate the emergence of chiral TSC phases, we focus on one of the two
blocks and rewrite it as:

$\displaystyle\hat{h}_{{\rm BdG},\sigma}(\bm{k})=\hat{h}_{{\rm
BdG;sc}}^{\prime}(\bm{k})-\sqrt{2}Me^{i\pi\lambda_{3}\rho_{1}/4}\rho_{3}+\tilde{M}\lambda_{1}\rho_{3}\,.$

As a next step, we perform a unitary transformation which transfers the phase
factor $e^{i\pi\lambda_{3}\rho_{1}/4}\rho_{3}$ from the magnetic part of the
Hamiltonian to $\hat{h}_{{\rm BdG;sc}}^{\prime}(\bm{k})$. We proceed by
focusing near the ${\rm\Gamma}$ point, which undergoes band inversions. We
thus expand $\hat{h}_{{\rm BdG;sc}}^{\prime}(\bm{k})$ up to linear order in
$k_{x}$ and $k_{y}$ about $\bm{k}=\bm{0}$, and obtain the expression:

$\displaystyle\hat{h}_{{\rm BdG},\sigma}^{\prime}(\bm{k})$
$\displaystyle\approx$
$\displaystyle\sqrt{2}e^{i\pi\lambda_{3}\rho_{1}/4}\big{(}\upsilon_{x}k_{x}\rho_{2}+\upsilon_{y}k_{y}\lambda_{2}\big{)}\tau_{3}-\mu_{\rm\Gamma}\tau_{3}$
(38)
$\displaystyle+\Delta\tau_{1}-\big{(}\sqrt{2}M-\tilde{M}\lambda_{1}\big{)}\rho_{3}\,,$

where
$\sqrt{2}\upsilon_{x,y}=\partial_{\bm{k}_{x,y}}\varepsilon_{x,y}(\bm{k})|_{\bm{k}=\bm{0}}$.
$\mu_{\rm\Gamma}$ defines the effective chemical potential at ${\rm\Gamma}$.
Depending on the values of $M$ and $\tilde{M}$, the band inversion at
${\rm\Gamma}$ occurs for an eigenstate with fixed projection of $\lambda_{1}$,
i.e. $\pm 1$. Therefore, it is eligible to project the above Hamiltonian onto,
say, the $\lambda_{1}=-1$ sector, and find:

$\displaystyle\hat{h}_{{\rm BdG},\sigma,\lambda_{1}=-1}^{\prime}(\bm{k})$
$\displaystyle=$
$\displaystyle\big{(}\upsilon_{x}k_{x}\rho_{2}+\upsilon_{y}k_{y}\rho_{1}\big{)}\tau_{3}-\mu_{\rm\Gamma}\tau_{3}$
(39)
$\displaystyle+\Delta\tau_{1}-\big{(}\sqrt{2}M+\tilde{M}\big{)}\rho_{3}\,.$

The topologically nontrivial properties of the above Hamiltonian are known
from previous works [99, 100], and arise from phase transition mediated by the
band inversion at $\bm{k}=\bm{0}$. Therefore, each class D Hamiltonian block
of interest gives rise to a single branch of chiral Majorana modes at a line
defect [98].

Concluding with the analysis of this system, it is important to comment on the
stability of the here-obtained chiral TSC phases with ${\cal N}=2$, since
these were derived under the presence of the unitary symmetry generated by
${\cal O}$. Hence, it is particularly interesting to examine the fate of these
phases in the presence of perturbations such as charge inhomogeneities and
Zeeman fields, which violate ${\cal O}$. Such a study is carried out in App. E
and we find that there indeed exists a certain window of stability for these
phases against ${\cal O}$-symmetry-breaking perturbations.

###  VII.5 Topological phases: Sk-SWC4 ground state

The analysis of the Sk-SWC4 can be immediately inferred from the previous
investigations of topological superconductivity stemming from the SWC4 phase.
Specifically, Ref. [42] showed that only nodal TSCs become accessible, and
support Majorana bidirectional modes. In fact, the nodes in the spectrum are
not removable due to an emergent Kramers degeneracy arising at the
${\rm\Gamma}$ point of the band structure. As it was discussed there, lifting
this degeneracy, as for instance by an out-of-plane field, permitted the
appearance of chiral TSCs with a single chiral Majorana mode branch per edge.

Interestingly, the addition of the magnetization components $\bm{M}_{\pm}$,
renders the spectrum fully gapped and achieves the removal of this Kramers
degeneracy at the ${\rm\Gamma}$. Therefore, a band inversion is now possible,
and allows the system to enter the topologically nontrivial phase. Given the
fact that the Hamiltonian resides in class D, a single band inversion at
${\rm\Gamma}$ converts the system into a chiral TSC with ${\cal N}=1$ and a
single chiral Majorana branch per line defect.

The above statements are backed by our numerical calculations. Representative
results are presented in Fig. 6, where we show the pivotal character of the
extra magnetic terms $\propto\bm{M}_{\pm}$ in splitting the degeneracy at
${\rm\Gamma}$, which appears in a system with only the SWC4 magnetic phase.

##  VIII Conclusions and Outlook

We bring forward an alternative route to controllably stabilize magnetic
skyrmion ground states. Our mechanism that we here term as magnetic skyrmion
catalysis (MSC) is free from the requirement of Dzyaloshinkii-Moriya
interaction, magnetic anisotropy and the application of an external Zeeman
field. In contrast, it solely relies on time-reversal symmetry (TRS)
violation, and the presence of flux $\vartheta$ in the ground state. As such,
it opens perspectives for functional topological superconducting platforms
that may be amenable to electrostatic or other types of control.

For our exploration, we adopt a Landau-type of approach, and consider that the
magnetization field $\bm{M}(\bm{r})$ is a weak perturbation to a paramagnetic
ground state. Due to the violation of TRS, the Landau free energy contains a
term proportional to the quantity $\int
d\bm{r}\bm{M}(\bm{r})\cdot\big{[}\partial_{x}\bm{M}(\bm{r})\times\partial_{y}\bm{M}(\bm{r})\big{]}$
which is crucial for the catalysis of magnetic skyrmions. In our study the
spatial variation of both the modulus $|\bm{M}(\bm{r})|$ and the orientation
$\bm{M}(\bm{r})/|\bm{M}(\bm{r})|$ of the magnetization are important for the
stabilization of the skyrmion magnetic ground state. This has to be contrasted
with the approach adopted in a recent study [101], which considers the
stabilization of smooth skyrmion textures originating from a ferromagnetic
ground state. While also the authors of Ref. 101 independently discuss the
possibility of promoting magnetic skyrmions due to Dirac points in the band
structure and TRS violation, the underlying mechanism stabilizing skyrmions
differs from ours.

Under the above assumptions, we provide a closed-form analytical expression
for the flux $\vartheta$, which acts as a source of magnetic skyrmion charge.
Using this expression, we determine the value of this coefficient for two
typical models characterized by a nonzero Berry curvature. In these cases, we
study the behavior of the intrinsic flux and show that it becomes enhanced
when the temperature is tuned to the energy difference between the Dirac mass
$m$ and the chemical potential $\mu$. Since $m$ acts as a source of
$\vartheta$ and, thus, is responsible for the MSC, we further discuss possible
scenarios that allow stabilizing a nonzero $m$. An appealing possibility is to
generate $m$ spontaneoussly via an attractive interaction in this channel.
Nonetheless, even away from criticality, $m$ and in turn a magnetic skyrmion
crystal ground state can be favored, solely by means of fluctuations.
Alternatively, in the presence of fluctuations, a nonzero $m$ can be imposed
by means of the orbital coupling to an applied magnetic field. The above
results also reveal that skyrmion crystal ground states and topological phases
generated via the MSC are compatible with the additional presence of an
external orbital magnetic field and thus, in turn, with the appearance of
superconducting vortices. Hence, new promising possibilities open up for
pinning Majorana zero modes using skyrmion-vortex excitations, in analogy to
the ones currently explored [57]. Within the MSC framework, this can become
possible by means of locally controlling flux in space, and generating
isolated magnetic skyrmion excitations.

In addition, we obtain the full set of possible magnetic ground states which
become accessible for itinerant magnets with tetragonal symmetry, and preserve
the full spin-rotational group, but violate TRS. This is achieved by
restricting to the principal harmonics entering the magnetization profile. We
find that there are two thermodynamically stable phases which acquire a
nonzero skyrmion charge by means of MSC. These are the Sk-SVC and Sk-SWC4 and
emerge as the skyrmion variants of the spin-vortex (SVC) and spin-whirl (SWC4)
crystal phases which are obtained in the presence of TRS. We demonstrate the
emergence of the Sk-SVC phase for a concrete extended Dirac model, which gives
rise to a quadratic band crossing in the massless case, and can describe the
protected surface states of a topological crystalline insulator.

Apart from the magnetic phase diagram, we also investigate the topological
phases obtained for the above extended Dirac model when the two magnetic
skyrmion crystal ground states mentioned above, coexist with a conventional
pairing gap. This can be for instance induced to the system either by
proximity to a conventional superconductor, or, spontaneously due to
interactions. We find that the Sk-SVC leads to a chiral topological
superconductor (TSC) with a topological invariant ${\cal N}=2$, which implies
that chiral Majorana modes can be trapped at line defects. Instead, the Sk-
SWC4 allows for phases with both ${\cal N}=1,2$, depending on the details of
the magnetic order parameter. Therefore, our results open alternative
perspectives for engineering topological superconductivity using the mechanism
of MSC, which is in principle less demanding than existing schemes for
inducing magnetic skyrmion phases.

Concerning the experimental realization of MSC, we note that the SVC phase has
been already found in CaKFe4As4 [102], while the SWC4 has been theoretically
predicted [90] for hole-doped BaFe2As2. Hence, Fe-based systems with
coexisting magnetism and superconductivity [103, 104, 105, 106, 107, 108]
appear promising to exhibit intrinsic chiral TSCs once flux emerges. Another
category of potential intrinsic TSCs is the recently discovered family of
Kagome superconductors [109, 110, 111, 112]. These are known to exhibit both
superconductivity and a flux phase [113, 114, 115] which can lead to a nonzero
Berry curvature in the energy bands [116, 117]. Therefore, the possible
additional emergence of magnetism promises to enable the MSC and allow in turn
for chiral TSC phases.

Apart from materials, further potential candidates include hybrids structures
based on topological crystalline insulators [92] or graphene systems [118,
119]. These are prominent for manipulating MZMs, since they can be in
principle electrostatically gated. Quite remarkably, one notes that twisted
bilayer graphene appears to provide all the necessary ingredients for MSC to
take place, since it harbors superconductivity [120], magnetism [121], Chern
phases [122, 123], while it has been also theoretically predicted to harbor
magnetic skyrmions due to long range Coulomb interactions in its quantum Hall
ferromagnetic phases [124].

## Acknowledgements

We wish to thank Achim Rosch for his valuable feedback on our work and for
further helpful discussions. In addition, PK acknowledges funding from the
National Natural Science Foundation of China (Grant No. 12250610194).

##  Appendix A Additional details regarding the calculation of $\vartheta$

The expansion of Eq. (5) for $\bm{q},\bm{p}\approx\bm{0}$, yields:

$\displaystyle\vartheta(\bm{q},\bm{p})=i\int\frac{d\bm{k}}{3\pi}\phantom{.}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\bigg{\\{}{\rm
tr}\left[\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{x}}\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{y}}\right]-\frac{\partial}{\partial
k_{x}}\leftrightarrow\frac{\partial}{\partial
k_{y}}\bigg{\\}}\big{(}q_{x}p_{y}-q_{y}p_{x}\big{)}\,.$

Since the integrand of $F^{(3)}$ is antisymmetric in
$\bm{p}\leftrightarrow\bm{q}$, the factor $q_{x}p_{y}-q_{y}p_{x}$ can be
replaced by $2q_{x}p_{y}$. Therefore, near a ferromagnetic instability the
cubic term reads $F^{(3)}=-{\vartheta}{\cal C}$ where:

$\displaystyle\vartheta=2i\int\frac{d\bm{k}}{3\pi}\phantom{.}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\left\\{{\rm
tr}\left[\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{x}}\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{y}}\right]-\frac{\partial}{\partial
k_{x}}\leftrightarrow\frac{\partial}{\partial k_{y}}\right\\}\,.$

Using the relation for a general Green function
$\hat{G}\hat{G}^{-1}=1\Rightarrow\partial\hat{G}=-\hat{G}(\partial\hat{G}^{-1})\hat{G}=\hat{G}(\partial\hat{{\cal
H}})\hat{G}$ we obtain:

$\displaystyle{\rm
tr}\left[\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{x}}\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{y}}\right]={\rm
tr}\left[\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial
k_{x}}\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial
k_{y}}\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\right].$

We now proceed by expressing the Green functions with the help of projectors.
At this point we consider that the system is described by a band structure
with energy dispersions $\varepsilon_{\alpha}(\bm{k})$. For calculating
${\vartheta}$ we will additionaly introduce a chemical potential $\mu$ and
work in the grand-canonical ensemble. Under these conditions we obtain:

$\displaystyle{\rm
tr}\left[\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{x}}\frac{\partial\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})}{\partial
k_{y}}\right]=\sum_{\alpha,\beta}{\rm tr}\left[\frac{\hat{{\cal
P}}_{\alpha}(\bm{k})}{[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\alpha}(\bm{k})]^{3}}\frac{\partial\hat{h}(\bm{k})}{\partial
k_{x}}\frac{\hat{{\cal
P}}_{\beta}(\bm{k})}{[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\beta}(\bm{k})]^{2}}\frac{\partial\hat{h}(\bm{k})}{\partial
k_{y}}\right]\,.$

Notably, in the above sum only the $\alpha\neq\beta$ contributions will
survive due to the antisymmetry under the exchange of the derivatives.
Therefore we have:

$\displaystyle{\vartheta}=2i\int\frac{d\bm{k}}{3\pi}\sum_{\alpha\neq\beta}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\frac{1}{[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\alpha}(\bm{k})]^{3}}\frac{1}{[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\beta}(\bm{k})]^{2}}\left\\{{\rm
tr}\left[\hat{{\cal
P}}_{\alpha}(\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial k_{x}}\hat{{\cal
P}}_{\beta}(\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial
k_{y}}\right]-\frac{\partial}{\partial
k_{x}}\leftrightarrow\frac{\partial}{\partial k_{y}}\right\\}$

and after symmetrizing, we find:

$\displaystyle{\vartheta}=i\int\frac{d\bm{k}}{3\pi}\sum_{\alpha\neq\beta}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\frac{\varepsilon_{\alpha}(\bm{k})-\varepsilon_{\beta}(\bm{k})}{\big{\\{}[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\alpha}(\bm{k})][{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\beta}(\bm{k})]\big{\\}}^{3}}\left\\{{\rm
tr}\left[\hat{{\cal
P}}_{\alpha}(\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial k_{x}}\hat{{\cal
P}}_{\beta}(\bm{k})\frac{\partial\hat{h}(\bm{k})}{\partial
k_{y}}\right]-\frac{\partial}{\partial
k_{x}}\leftrightarrow\frac{\partial}{\partial k_{y}}\right\\}\,.$

We can further simplicify the above expression using relations satisfied by
the eigenstates of the system $\left|\bm{u}_{\alpha}(\bm{k})\right>$. We
specifically obtain:

$\displaystyle{\vartheta}=i\int\frac{d\bm{k}}{3\pi}\sum_{\alpha\neq\beta}T\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}}}\frac{[\varepsilon_{\alpha}(\bm{k})-\varepsilon_{\beta}(\bm{k})]^{3}}{\big{\\{}[{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\alpha}(\bm{k})][{\color[rgb]{0,0,0}i\omega_{\nu}}+\mu-\varepsilon_{\beta}(\bm{k})]\big{\\}}^{3}}\Big{[}\big{<}\partial_{k_{x}}\bm{u}_{\alpha}(\bm{k})\big{|}\bm{u}_{\beta}(\bm{k})\big{>}\big{<}\bm{u}_{\beta}(\bm{k})\big{|}\partial_{k_{y}}\bm{u}_{\alpha}(\bm{k})\big{>}-\partial_{k_{x}}\leftrightarrow\partial_{k_{y}}\Big{]}.$

After the Matsubara summation we obtain

$\displaystyle{\vartheta}$ $\displaystyle=$ $\displaystyle
i\int\frac{d\bm{k}}{6\pi}\sum_{\alpha\neq\beta}\big{<}\partial_{k_{x}}\bm{u}_{\alpha}(\bm{k})\big{|}\bm{u}_{\beta}(\bm{k})\big{>}\Bigg{\\{}12\frac{f[\varepsilon_{\alpha}(\bm{k})-\mu]-f[\varepsilon_{\beta}(\bm{k})-\mu]}{\big{[}\varepsilon_{\alpha}(\bm{k})-\varepsilon_{\beta}(\bm{k})\big{]}^{2}}-6\frac{f^{\prime}[\varepsilon_{\alpha}(\bm{k})-\mu]+f^{\prime}[\varepsilon_{\beta}(\bm{k})-\mu]}{\varepsilon_{\alpha}(\bm{k})-\varepsilon_{\beta}(\bm{k})}$
$\displaystyle\qquad\qquad+f^{\prime\prime}[\varepsilon_{\alpha}(\bm{k})-\mu]-f^{\prime\prime}[\varepsilon_{\beta}(\bm{k})-\mu]\Bigg{\\}}\big{<}\bm{u}_{\beta}(\bm{k})\big{|}\partial_{k_{y}}\bm{u}_{\alpha}(\bm{k})\big{>}-\partial_{k_{x}}\leftrightarrow\partial_{k_{y}}\,.$

To obtain the result of Eq. (7), we first exchange the band indices
$\alpha\leftrightarrow\beta$ in the terms containing the Fermi-Dirac
distributions of the $\beta$ bands. This gives rise to an overall factor of
$2$ to the terms containing the Fermi-Dirac distribution of the $\alpha$
bands. Afterwards, we eliminate the dependence on the $\beta$ bands by
expressing $\varepsilon_{\beta}(\bm{k})$ in the operator form
$\hat{h}(\bm{k})$. One finally arrives at Eq. (7) is finally by expressing the
sum containing the various orders of derivatives of the Fermi-Dirac
distribution in a compact fashion.

##  Appendix B Minimization of the magnetic free energy for a nonzero flux

In this appendix we minimize the free energy in Eq. (17) with respect
$\bm{M}_{1,2}$ and identify also the constrained $\bm{M}_{\pm}$. Our analysis
follows the spirit of Ref. [90]. At quadratic order of the expansion of the
free energy for $\bm{M}_{1,2}$, a degeneracy and enhanced symmetry emerges.
Hence, this allows us to adopt the parametrization
$\bm{M}_{1}=M\cos\eta\phantom{.}\hat{\bm{n}}_{1}$ and
$\bm{M}_{2}=M\sin\eta\phantom{.}\hat{\bm{n}}_{2}$, with
$|\hat{\bm{n}}_{1,2}|^{2}=1$ and $\eta\in[0,\pi/2]$. Thus the Landau
functional reads:

$\displaystyle{\cal F}_{\bm{M}_{1,2}}=\alpha
M^{2}+\frac{\tilde{\beta}}{2}M^{4}+\frac{\beta-\tilde{\beta}}{2}\big{(}\cos^{4}\eta|\hat{\bm{n}}_{1}^{2}|^{2}+\sin^{4}\eta|\hat{\bm{n}}_{2}^{2}|^{2}\big{)}M^{4}$
$\displaystyle\phantom{{\dagger}}\qquad+\left[g-\tilde{\beta}+\tilde{g}\phantom{.}\frac{|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}|^{2}+|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}^{*}|^{2}}{2}-\gamma\phantom{.}\frac{|\hat{\bm{n}}_{1}\times\hat{\bm{n}}_{2}|^{2}+|\hat{\bm{n}}_{1}\times\hat{\bm{n}}_{2}^{*}|^{2}}{2}\right]\sin^{2}(2\eta)\frac{M^{4}}{4}\,.\qquad$
(40)

The presence of a nonzero $\gamma$ solely affects the double-$\bm{Q}$ phases.
There we focus on these. In order to extremize the free energy functional in
the case of double-$\bm{Q}$ phases ($\eta\neq 0,\pi/2$) it is convenient first
to write the explicit form of the two complex spin-vectors as
$\hat{\bm{n}}_{s}=(|a_{s}|e^{i\zeta_{s}},|b_{s}|e^{i\xi_{s}},|c_{s}|e^{i\omega_{s}})$
with $s=1,2$. By virtue of translational invariance we can arbitrarily choose
the overall phase factor of each complex vector. This allows us to set
$\omega_{s}=0$ ($s=1,2$). We can now proceed by further decomposing the
complex vectors $\hat{\bm{n}}_{s}$ ($s=1,2$) into real and imaginary parts, as
in the previous paragraph. After taking into account translational invariance
and by also exploiting the SO(3) spin-rotational invariance allowing us to fix
the spin-orientation of one of the four real vectors
($\hat{\bm{n}}_{s,\operatorname{Re}}\,,\hat{\bm{n}}_{s,\operatorname{Im}}$) we
can write:

$\displaystyle\hat{\bm{n}}_{1,\operatorname{Re}}$ $\displaystyle=$
$\displaystyle\cos\lambda_{1}\left(0\,,0\,,1\right)\,,\qquad$
$\displaystyle\hat{\bm{n}}_{1,\operatorname{Im}}$ $\displaystyle=$
$\displaystyle\sin\lambda_{1}\left(\cos\chi\,,\sin\chi\,,0\right)\,,\qquad$
$\displaystyle\hat{\bm{n}}_{2,\operatorname{Re}}$ $\displaystyle=$
$\displaystyle\cos\lambda_{2}\left(\sin\rho\cos\omega\,,\sin\rho\sin\omega\,,\cos\rho\right)\,,\qquad$
$\displaystyle\hat{\bm{n}}_{2,\operatorname{Im}}$ $\displaystyle=$
$\displaystyle\sin\lambda_{2}\left(\cos\phi\,,\sin\phi\,,0\right)\,.$

For compactness, we now define the coefficients:

$\displaystyle
P=\frac{|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}|^{2}+|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}^{*}|^{2}}{2}$
$\displaystyle=$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}\,,$
$\displaystyle
K=\frac{|\hat{\bm{n}}_{1}\times\hat{\bm{n}}_{2}|^{2}+|\hat{\bm{n}}_{1}\times\hat{\bm{n}}_{2}^{*}|^{2}}{2}$
$\displaystyle=$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\times\hat{\bm{n}}_{2,\operatorname{Re}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Im}}\times\hat{\bm{n}}_{2,\operatorname{Im}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Re}}\times\hat{\bm{n}}_{2,\operatorname{Im}})^{2}+(\hat{\bm{n}}_{1,\operatorname{Im}}\times\hat{\bm{n}}_{2,\operatorname{Re}})^{2}\,.$

Let us now calculate a number of terms that appear in the free energy

$\displaystyle|\hat{\bm{n}}_{1}^{2}|^{2}$ $\displaystyle=$
$\displaystyle\cos^{2}(2\lambda_{1})\,,$
$\displaystyle|\hat{\bm{n}}_{2}^{2}|^{2}$ $\displaystyle=$
$\displaystyle\cos^{2}(2\lambda_{2})+\sin^{2}(2\lambda_{2})\sin^{2}\rho\cos^{2}\delta_{2}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}$
$\displaystyle=$
$\displaystyle\cos^{2}\lambda_{1}\cos^{2}\lambda_{2}\cos^{2}\rho\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}$
$\displaystyle=$
$\displaystyle\sin^{2}\lambda_{1}\sin^{2}\lambda_{2}\cos^{2}\delta_{1}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}$
$\displaystyle=$ $\displaystyle 0\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}$
$\displaystyle=$
$\displaystyle\sin^{2}\lambda_{1}\cos^{2}\lambda_{2}\sin^{2}\rho\cos^{2}\delta_{3}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\times\hat{\bm{n}}_{2,\operatorname{Re}})^{2}$
$\displaystyle=$
$\displaystyle\cos^{2}\lambda_{1}\cos^{2}\lambda_{2}-(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Im}}\times\hat{\bm{n}}_{2,\operatorname{Im}})^{2}$
$\displaystyle=$
$\displaystyle\sin^{2}\lambda_{1}\sin^{2}\lambda_{2}-(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Re}}\times\hat{\bm{n}}_{2,\operatorname{Im}})^{2}$
$\displaystyle=$
$\displaystyle\cos^{2}\lambda_{1}\sin^{2}\lambda_{2}-(\hat{\bm{n}}_{1,\operatorname{Re}}\cdot\hat{\bm{n}}_{2,\operatorname{Im}})^{2}\,,$
$\displaystyle(\hat{\bm{n}}_{1,\operatorname{Im}}\times\hat{\bm{n}}_{2,\operatorname{Re}})^{2}$
$\displaystyle=$
$\displaystyle\sin^{2}\lambda_{1}\cos^{2}\lambda_{2}-(\hat{\bm{n}}_{1,\operatorname{Im}}\cdot\hat{\bm{n}}_{2,\operatorname{Re}})^{2}\,,$

where we set $\delta_{1}=\phi-\chi$, $\delta_{2}=\omega-\phi$ and
$\delta_{3}=\chi-\omega$. From the above we find the relation $K=1-P$. This
implies that one can go back to the start and rewrite the Landau functional
as:

$\displaystyle{\cal F}_{\bm{M}_{1,2}}$ $\displaystyle=$ $\displaystyle\alpha
M^{2}+\frac{\tilde{\beta}}{2}M^{4}+\frac{\beta-\tilde{\beta}}{2}\big{(}\cos^{4}\eta|\hat{\bm{n}}_{1}^{2}|^{2}+\sin^{4}\eta|\hat{\bm{n}}_{2}^{2}|^{2}\big{)}M^{4}$
(41)
$\displaystyle+\left[g-\gamma-\tilde{\beta}+\left(\tilde{g}+\gamma\right)\frac{|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}|^{2}+|\hat{\bm{n}}_{1}\cdot\hat{\bm{n}}_{2}^{*}|^{2}}{2}\right]\sin^{2}(2\eta)\frac{M^{4}}{4}\,.$

Thus, we find that the presence of a nonzero $\chi$, renormalizes the
following parameters at $\chi=0$ according to the replacement $g\mapsto
g-\gamma$ and $\tilde{g}\mapsto\tilde{g}+\gamma$.

##  Appendix C Magnetic ground states in the presence of nonzero flux

Based on the above results, we find that all the magnetic ground states
obtained in the presence of a nonzero flux, can be obtained using the ground
states for zero flux discussed in Ref. [90]. First of all, we note that the
single-$\bm{Q}$ phase, i.e., the magnetic stripe and helix do not induce
nonzero $\bm{M}_{\pm}$. Therefore, their expression remain as previous. From
the double-$\bm{Q}$ phases, we also find that neither the CSDW induces the two
additional magnetization components, since it is a collinear phase. The
expression for the remaining six magnetic ground states are presented in Table
2.

Phases |  $\bm{M}_{1}\big{(}\bm{r}\big{)}$, $\bm{M}_{2}\big{(}\bm{r}\big{)}$, $\bm{M}_{\pm}\big{(}\bm{r}\big{)}$
---|---
Sk-SVC | $\bm{M}_{1}=M\big{(}\cos\big{(}\bm{Q}_{1}\cdot\bm{r}\big{)},0,0\big{)}$, $\bm{M}_{2}=M\big{(}0,\cos\big{(}\bm{Q}_{2}\cdot\bm{r}\big{)},0\big{)}$, $\bm{M}_{\pm}=\pm\tilde{M}\cos\big{(}\bm{Q}_{\pm}\cdot\bm{r}\big{)}\big{(}0,0,1\big{)}$
Sk-MS$\|$MH |  $\bm{M}_{1}=\left(0,0,M\cos\eta\cos\left(\bm{Q}_{1}\cdot\bm{r}\right)\right)$, $\bm{M}_{2}=M\sin\eta\big{(}\sin\lambda\sin\big{(}\bm{Q}_{2}\cdot\bm{r}\big{)},0,\cos\lambda\cos\big{(}\bm{Q}_{2}\cdot\bm{r}\big{)}\big{)}$
$\bm{M}_{\pm}=\big{(}0,\tilde{M}\sin\lambda\sin\eta\cos\eta\sin\big{(}\bm{Q}_{\pm}\cdot\bm{r}\big{)},0\big{)}$
$\mathrm{Sk-SWC_{4}}$ |  $\bm{M}_{1}=M\left(\cos\lambda\sin\left(\bm{Q}_{1}\cdot\bm{r}\right),0,\sin\lambda\cos\left(\bm{Q}_{1}\cdot\bm{r}\right)\right)$ , $\bm{M}_{2}=M\left(0,\cos\lambda\sin\left(\bm{Q}_{2}\cdot\bm{r}\right),\sin\lambda\cos\left(\bm{Q}_{2}\cdot\bm{r}\right)\right)$
$\bm{M}_{\pm}=\tilde{M}\left(\sin\lambda\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right),\pm\sin\lambda\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right),\cos\lambda\cos\left(\bm{Q}_{\pm}\cdot\bm{r}\right)\right)$
$\mathrm{Sk-SWC_{2}}$ |  $\bm{M}_{1}=M\cos\eta\left(\cos\left(\bm{Q}_{1}\cdot\bm{r}\right),0,\sin\left(\bm{Q}_{1}\cdot\bm{r}\right)\right)$, $\bm{M}_{2}=\sqrt{2}M\sin\eta\left(\cos\lambda\cos\left(\bm{Q}_{2}\cdot\bm{r}\right),\sin\lambda\sin\left(\bm{Q}_{2}\cdot\bm{r}\right),0\right)$
$\bm{M}_{\pm}=\tilde{M}\left(\cos\left(\bm{Q}_{\pm}\cdot\bm{r}\right),\pm\cot\lambda\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right),\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right)\right)/2$
$\widetilde{\mathrm{MS\perp MH}}$ |  $\bm{M}_{1}=\left(0,0,M\cos\eta\cos\left(\bm{Q}_{1}\cdot\bm{r}\right)\right)$, $\bm{M}_{2}=M\sin\eta\left(\sin\left(\bm{Q}_{2}\cdot\bm{r}\right),\cos\left(\bm{Q}_{2}\cdot\bm{r}\right),0\right)/\sqrt{2}$
$\bm{M}_{\pm}=\tilde{M}\sin\eta\cos\eta\left(\mp\cos\left(\bm{Q}_{\pm}\cdot\bm{r}\right),\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right),0\right)$
$\widetilde{\mathrm{DPMH}}$ |  $\bm{M}_{1}=M\left(\sin\left(\bm{Q}_{1}\cdot\bm{r}\right),0,\cos\left(\bm{Q}_{1}\cdot\bm{r}\right)\right)/\sqrt{2}$ , $\bm{M}_{2}=M\left(\sin\left(\bm{Q}_{2}\cdot\bm{r}\right),0,\cos\left(\bm{Q}_{2}\cdot\bm{r}\right)\right)/\sqrt{2}$
$\bm{M}_{\pm}=\tilde{M}\left(0,\left(i\mp
i\right)\sin\left(\bm{Q}_{\pm}\cdot\bm{r}\right),0\right)$

Table 2: The coordinate-space magnetization profile $\bm{M}(\bm{r})$, for the
six double-$\bm{Q}$ magnetic phases of Ref. [90], which become modified in the
presence of a nonzero flux in the free energy of Eq. (17).

##  Appendix D Landau free energy coefficients and phase diagram details

We can explicitly calculate the coefficients for the Landau free energy Eq.
(17) given the concrete model we consider in Eq. (22) after phenomenologically
adding a mass term $m$. The Hamiltonian of the system can be written as
$\hat{h}(\bm{k})=\bm{d}(\bm{k})\cdot\bm{\kappa}$ with the help of the Pauli
matrix vector $\bm{\kappa}$ and
$\bm{d}(\bm{k})=\left(k_{x}k_{y}/m_{1},m,\big{(}k_{x}^{2}-k_{y}^{2}\big{)}/2m_{2}\right)$.
We introduce the bare Matsubara Green function (for convenience we employ a
slighly modified notation below):

$\hat{G}\left({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}\right)=\frac{1}{{\color[rgb]{0,0,0}i\omega_{\nu}}-\hat{h}(\bm{k})}=\frac{1}{2}\underset{\alpha=\pm}{\sum}\frac{\mathds{1}+\alpha\hat{\bm{d}}(\bm{k})\cdot\bm{\kappa}}{{\color[rgb]{0,0,0}i\omega_{\nu}}-\varepsilon_{\alpha}(\bm{k})}=\underset{\alpha=\pm}{\sum}\frac{\hat{\hat{\mathcal{P}}}_{\alpha}(\bm{k})}{{\color[rgb]{0,0,0}i\omega_{\nu}}-\varepsilon_{\alpha}(\bm{k})}\\\
=G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{\alpha}\hat{\mathcal{P}}_{\alpha}(\bm{k}).$
(42)

Here,
$\hat{\mathcal{P}}_{\alpha}(\bm{k})=\big{[}\mathds{1}+\alpha\hat{\bm{d}}(\bm{k})\cdot\bm{\kappa}\big{)}]/2$
is the projector representing the valley degree of freedom,
$\hat{\bm{d}}(\bm{k})\equiv\bm{d}(\bm{k})/\left|\bm{d}(\bm{k})\right|$, and
$G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{\alpha}=\left({\color[rgb]{0,0,0}i\omega_{\nu}}-\varepsilon_{\alpha}(\bm{k})\right)^{-1}$.
It is helpful to define the following quantities:

$A_{\bm{k}_{1},\alpha_{1};\bm{k}_{2},\alpha_{2};\cdots\cdots;\bm{k}_{n}\alpha_{n}}^{\left(n\right)}=\mathrm{tr}\big{\\{}\hat{\mathcal{P}}_{\alpha_{1}}\left(\bm{k}_{1}\right)\hat{\mathcal{P}}_{\alpha_{2}}\left(\bm{k}_{2}\right)\cdots\cdots\hat{\mathcal{P}}_{\alpha_{n}}\left(\bm{k}_{n}\right)\big{\\}}.$
(43)

The coefficients of the Landau free energy are given by taking the derivatives
of it with respect to the order parameters $\bm{M}_{\bm{q}}$. The free energy
functional is expressed in terms of the perturbation in operator form
$\hat{V}=\bm{M}_{\bm{q}}\cdot\bm{\sigma}$:

$F=F_{0}-T\,\mathrm{Tr}\ln\left(1-\hat{G}\hat{V}\right)={F}_{0}+T\,\sum_{n=1}^{\infty}\frac{1}{n}\mathrm{Tr}\big{(}\hat{G}\hat{V}\big{)}^{n}\equiv
F_{0}+\sum_{n=1}^{\infty}F^{(n)}\,.$ (44)

### The Second Order Terms

The second order contribution to the free energy $F^{\left(2\right)}$ is
written as:

${F}^{\left(2\right)}=\alpha\left(\left|\bm{M}_{1}\right|^{2}+\left|\bm{M}_{2}\right|^{2}\right)+\bar{\alpha}\left(|\bm{M}_{+}|^{2}+|\bm{M}_{-}|^{2}\right)\,.$
(45)

The coefficients $\alpha$ and $\bar{\alpha}$ contain the spin susceptibility,
which is given by the expression:

$\displaystyle\chi(\bm{q})=-\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{\alpha,\beta=\pm}\,A_{\bm{k},\alpha;\bm{k}+\bm{q},\beta}^{\left(2\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{\alpha}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{q}}^{\beta}$
(46)

with the projector
$A_{\bm{k},\alpha;\bm{k}^{\prime},\beta}^{(2)}=\big{(}1+\alpha\beta\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime}}\big{)}/2$.
Considering a Hubbard interaction with strength $U$, yields the coefficients
$\alpha=2/U-\chi(\bm{Q}_{1,2}),$ and $\bar{\alpha}=2/U-\chi(\bm{Q}_{\pm})$.
Specifically, the interaction strength $U$ is set to meet the Stoner criterion
$2=U\chi(\bm{Q}_{1,2})$, and thus we have:
$\bar{\alpha}=\chi(\bm{Q}_{1,2})-\chi(\bm{Q}_{\pm})$.

### The Third Order Terms

The third-order free energy is written as:

$\displaystyle F^{(3)}=$
$\displaystyle\frac{T}{3N}\sum_{\bm{k}}\mathrm{Tr}\left\\{\hat{G}\left({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}\right)\bm{M}_{\bm{k}-\bm{k}^{\prime}}\cdot\bm{\sigma}\hat{G}\left({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}^{\prime}\right)\bm{M}_{\bm{k}^{\prime}-\bm{k}^{\prime\prime}}\cdot\bm{\sigma}\hat{G}\left({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}^{\prime\prime}\right)\bm{M}_{\bm{k}^{\prime\prime}-\bm{k}}\cdot\bm{\sigma}\right\\}$
$\displaystyle=$
$\displaystyle-\vartheta_{1,2}\left[\bm{M}^{*}_{+}\cdot\left(\bm{M}_{1}\times\bm{M}_{2}\right)-\bm{M}^{*}_{-}\cdot\left(\bm{M}_{1}\times\bm{M}^{*}_{2}\right)+\mathrm{c.c}\right].$
(47)

We can start with calculating the projector:

$\displaystyle
A_{\bm{k},s;\bm{k}^{\prime},s^{\prime};\bm{k}^{\prime\prime},s^{\prime\prime}}^{\left(3\right)}=$
$\displaystyle\mathrm{Tr}\left\\{\hat{\mathcal{P}}_{s}(\bm{k})\hat{\mathcal{P}}_{s^{\prime}}\left(\bm{k}^{\prime}\right)\hat{\mathcal{P}}_{s^{\prime\prime}}\left(\bm{k}^{\prime\prime}\right)\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\left(1+ss^{\prime}\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime}}+ss^{\prime\prime}\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime}}+s^{\prime}s^{\prime\prime}\hat{\bm{d}}_{\bm{k}^{\prime}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime}}+i\epsilon_{ij\ell}ss^{\prime}s^{\prime\prime}\hat{\bm{d}}_{\bm{k}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}^{\ell}\right).$
(48)

For term of $\bm{M}_{+}^{*}\cdot\left(\bm{M}_{1}\times\bm{M}_{2}\right)$, the
expression of the coefficient is generated by considering the functional
derivatives:

$\displaystyle\vartheta_{1,2}$ $\displaystyle=$
$\displaystyle-\frac{2iT}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\underset{s,s^{\prime},s^{\prime\prime}=\pm}{\sum}\Bigg{(}A_{\bm{k},s;\bm{k}-\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1},s^{\prime\prime}}^{\left(3\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{Q}_{2}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime\prime}}$
(49)
$\displaystyle\qquad\qquad\qquad\qquad\phantom{{\dagger}}-A_{\bm{k},s;\bm{k}-\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime}}^{\left(3\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime}}\Bigg{)}$

where we note that:

$\displaystyle
A_{\bm{k},s;\bm{k}^{\prime},s^{\prime};\bm{k}^{\prime\prime},s^{\prime\prime}}^{\left(3\right)}=$
$\displaystyle
A_{\bm{k}^{\prime},s^{\prime};\bm{k}^{\prime\prime},s^{\prime\prime};\bm{k},s;}^{\left(3\right)}=A_{\bm{k}^{\prime\prime},s^{\prime\prime};\bm{k},s;\bm{k}^{\prime},s^{\prime}}^{\left(3\right)}.$
(50)

If TRS holds for the system, the projector $A^{\left(3\right)}$ is real and
the coeffcient $\vartheta_{1,2}$ is zero. When the TRS of the system breaks,
the triple product term of $\hat{\bm{d}}_{\bm{k}}$ in the
$A_{\bm{k},s;\bm{k}^{\prime},s^{\prime};\bm{k}^{\prime\prime},s^{\prime\prime}}^{\left(3\right)}$
becomes nonzero and leads to an imaginary factor. Notice that here the
imaginary parts of
$A_{\bm{k},s;\bm{k}-\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1},s^{\prime\prime}}^{\left(3\right)}$
and
$A_{\bm{k},s;\bm{k}-\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime}}^{\left(3\right)}$
have the opposite sign when integrated, so the expression of the coefficients
can be simplified as:

$\vartheta_{1,2}=\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\underset{s,s^{\prime},s^{\prime\prime}=\pm}{\sum}2\operatorname{Im}\left(A_{\bm{k},s;\bm{k}-\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1},s^{\prime\prime}}^{\left(3\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{Q}_{2}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime\prime}}\right).$
(51)

### The Fourth Order Terms

Similarly, we can write down the fourth order free energy of the system:

$\displaystyle F^{\left(4\right)}$ $\displaystyle=$
$\displaystyle\frac{T}{4N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\mathrm{Tr}\Big{\\{}\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k})\bm{M}_{\bm{k}-\bm{k}^{\prime}}\cdot\bm{\sigma}\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}^{\prime})\bm{M}_{\bm{k}^{\prime}-\bm{k}^{\prime\prime}}\cdot\bm{\sigma}\hat{G}({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}^{\prime\prime})\bm{M}_{\bm{k}^{\prime\prime}-\bm{k}^{\prime\prime\prime}}\bm{\sigma}\hat{G}\left({\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}^{\prime\prime\prime}\right)\bm{M}_{\bm{k}^{\prime\prime\prime}-\bm{k}}\cdot\bm{\sigma}\Big{\\}}$
(52) $\displaystyle=$
$\displaystyle\frac{\tilde{\beta}}{2}\left(\left|\bm{M}_{1}\right|^{2}+\left|\bm{M}_{2}\right|^{2}\right)^{2}+\frac{\beta-\tilde{\beta}}{2}\left(\left|\bm{M}^{2}_{1}\right|^{2}+\left|\bm{M}^{2}_{2}\right|^{2}\right)+\left(g-\tilde{\beta}\right)\left|\bm{M}_{1}\right|^{2}\left|\bm{M}_{2}\right|^{2}+\frac{\tilde{g}}{2}\left(\left|\bm{M}_{1}\cdot\bm{M}_{2}\right|^{2}+\left|\bm{M}_{1}\cdot\bm{M}^{*}_{2}\right|^{2}\right).\quad\qquad$

Once again we start with the calculation of the trace of the projectors:

$\displaystyle
A_{\bm{k},s;\bm{k}^{\prime},s^{\prime};\bm{k}^{\prime\prime},s^{\prime\prime};\bm{k}^{\prime\prime\prime},s^{\prime\prime\prime}}^{\left(4\right)}=\mathrm{Tr}\left\\{\hat{\mathcal{P}}_{s}(\bm{k})\hat{\mathcal{P}}_{s^{\prime}}\left(\bm{k}^{\prime}\right)\hat{\mathcal{P}}_{s^{\prime\prime}}\left(\bm{k}^{\prime\prime}\right)\hat{\mathcal{P}}_{s^{\prime\prime\prime}}\left(\bm{k}^{\prime\prime\prime}\right)\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{8}\left(1+ss^{\prime}\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime}}+ss^{\prime\prime}\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime}}+ss^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}+s^{\prime}s^{\prime\prime}\hat{\bm{d}}_{\bm{k}^{\prime}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime}}\right.+s^{\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}^{\prime}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}+s^{\prime\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}\cdot\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}$
$\displaystyle+i\epsilon_{ijk}ss^{\prime}s^{\prime\prime}\hat{\bm{d}}_{\bm{k}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}^{k}+i\epsilon_{ijk}ss^{\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}^{k}+i\epsilon_{ijk}ss^{\prime\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}^{k}+i\epsilon_{ijk}s^{\prime}s^{\prime\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}^{\prime}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}^{k}$
$\displaystyle+\left.\left(\delta_{ij}\delta_{kl}-\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)ss^{\prime}s^{\prime\prime}s^{\prime\prime\prime}\hat{\bm{d}}_{\bm{k}}^{i}\hat{\bm{d}}_{\bm{k}^{\prime}}^{j}\hat{\bm{d}}_{\bm{k}^{\prime\prime}}^{k}\hat{\bm{d}}_{\bm{k}^{\prime\prime\prime}}^{l}\right)\,.$
(53)

The fourth order free energy contains two type of terms. The one only involves
one of the two order parameters $\bm{M}_{1,2}$, while the other involves both.
We find the following expressions:

$\displaystyle\frac{\tilde{\beta}}{2}$ $\displaystyle=$
$\displaystyle\frac{8T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}^{\pm}\left(A_{\bm{k},s;\bm{k}+\bm{Q}_{1\left(2\right)},s^{\prime};\bm{k},s^{\prime\prime};\bm{k}+\bm{Q}_{1\left(2\right)},s^{\prime\prime\prime}}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1\left(2\right)}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1\left(2\right)}}^{s^{\prime\prime\prime}}\right)\,,$
$\displaystyle\frac{\beta-\tilde{\beta}}{2}$ $\displaystyle=$
$\displaystyle\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}^{\pm}\Bigg{\\{}4A_{\bm{k}+\bm{Q}_{1\left(2\right)},s;\bm{k},s^{\prime};\bm{k}-\bm{Q}_{1\left(2\right)},s^{\prime\prime};\bm{k},s^{\prime\prime\prime}}^{(4)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1\left(2\right)}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}-\bm{Q}_{1\left(2\right)}}^{s^{\prime\prime\prime}}$
$\displaystyle\qquad\quad\qquad\quad-2A_{\bm{k},s;\bm{k}+\bm{Q}_{1\left(2\right)},s^{\prime};\bm{k},s^{\prime\prime};\bm{k}+\bm{Q}_{1\left(2\right)},s^{\prime\prime\prime}}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1\left(2\right)}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1\left(2\right)}}^{s^{\prime\prime\prime}}\Bigg{\\}}\,,$
$\displaystyle g$ $\displaystyle=$
$\displaystyle\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}^{\pm}\Bigg{\\{}4A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k},s^{\prime\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime\prime}}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}$
$\displaystyle\qquad\qquad-
A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{2},d}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{a}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}$
$\displaystyle\qquad\qquad-
A_{\bm{k},s;\bm{k}+\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{1},d}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{a}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime\prime\prime}}\Bigg{\\}},\quad\qquad$
$\displaystyle\frac{\tilde{g}}{2}$ $\displaystyle=$
$\displaystyle\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}^{\pm}\Bigg{\\{}A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime\prime}}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}$
$\displaystyle\quad\qquad+A_{\bm{k},s;\bm{k}+\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{1},s^{\prime\prime\prime}}^{\left(4\right)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime\prime\prime}}\Bigg{\\}}.$

By using the following relation:

$\displaystyle\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}=\pm}A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime\prime}}^{\left(4\right)}+A_{\bm{k},s;\bm{k}+\bm{Q}_{2},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{1},s^{\prime\prime\prime}}^{\left(4\right)}$
$\displaystyle\qquad\qquad\qquad\qquad=2\operatorname{Re}\left(A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q_{2},s^{\prime\prime\prime}}}^{\left(4\right)}\right)$
(54)

the above expressions are simplified as follows:

$\displaystyle g$ $\displaystyle=$
$\displaystyle\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\sum_{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}}^{\pm}\Bigg{\\{}4A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k},s^{\prime\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime\prime}}^{(4)}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}$
$\displaystyle\quad-2\operatorname{Re}\left(A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{2},d}^{(4)}\right)G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{a}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}\Bigg{\\}},$
$\displaystyle\frac{\tilde{g}}{2}$ $\displaystyle=$
$\displaystyle\frac{2T}{N}\sum_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}\underset{s,s^{\prime},s^{\prime\prime},s^{\prime\prime\prime}=\pm}{\sum}2\operatorname{Re}\left(A_{\bm{k},s;\bm{k}+\bm{Q}_{1},s^{\prime};\bm{k}+\bm{Q}_{1}+\bm{Q}_{2},s^{\prime\prime};\bm{k}+\bm{Q}_{2},s^{\prime\prime\prime}}^{\left(4\right)}\right)G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}}^{s}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}}^{s^{\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{1}+\bm{Q}_{2}}^{s^{\prime\prime}}G_{{\color[rgb]{0,0,0}i\omega_{\nu}},\bm{k}+\bm{Q}_{2}}^{s^{\prime\prime\prime}}.$

The numerical calculation of the above coefficient is performed with a cut-off
of $\left|k_{x}\right|,\left|k_{x}\right|\leq 1.5$ and a mesh of $N=100^{2}$
points in $\bm{k}$ space. Figure 7 collects our results for the magnetic phase
diagrams discussed in the main text.

Figure 7: (a) The table shows the data employed to obtain Fig. 2(a) using the
expressions defined in App. D. Each column of the table corresponds, one-to-
one from left to right, to a pair of points in Fig. 2 which are obtained for
by discarding and including, respectively, the contribution of $\gamma$ which
is induced by a nonzero $m$. Our calculations are carried out for the
temperature value $T=1/30$ and the chemical potential is $\mu=0.17$ when
$m=0$. The Schrödinger masses are set to be $m_{1}=1$ and $m_{2}=0.3$. For all
our calculations, the chemical potential is suitably adjusted throughout to
fix the electron density of the system. Thus the shape of the Fermi surface
stays practically unchanged and the corresponding nesting vectors are pinned
around $\bm{Q}_{1}\simeq\big{(}0.7,0\big{)}$ and
$\bm{Q}_{2}\simeq\big{(}0,0.7\big{)}$ with
$\chi\big{(}\bm{Q}_{1,2}\big{)}\simeq 0.98$. To ensure nesting at the above
nesting vectors, the chemical potential is suitably adjusted throughout. (b)
The table shows the data of the corresponding point from right to left in Fig.
2(b), upon varying $m$. The calculation is obtained when the temperature is
$T=1/70$ and the chemical potential is $\mu=1/14$ when $m=0$. The Schrödinger
massses are set to be $m_{1}=1$ and $m_{2}=0.7$. The chemical potential is
tuned to fix the shape of the Fermi surface, and thus the corresponding
nesting vectors are pinned around $\bm{Q}_{1}\sim\big{(}0.65,0\big{)}$ and
$\bm{Q}_{2}\sim\big{(}0,0.65\big{)}$.

##  Appendix E Additional calculations on the robustness of the Sk-SVC chiral
topological superconductor

This appendix presents additional numerical calculations in order to
investigate the stability of the chiral TSC with ${\cal N}=2$ which emerges
for the case of a Sk-SVC phase. Figure 8 (9) discusses the stability against
the addition of an external field (of charge density wave terms).

Figure 8: Band structure for a Sk-SVC topological superconductor in the
additional presence of a Zeeman field $\bm{B}$. The field enters the BdG
Hamiltonian through the term $\bm{B}\cdot\bm{\sigma}$. Blue (red) color shows
the dispersions in the presence (absence) of the field. (a)-(d) and (e)-(h)
show two different cases of field configurations. For the top [bottom] row we
have $\bm{B}=0.01(1,\sqrt{2},\sqrt{5})$ [$\bm{B}=0.1(1,\sqrt{2},\sqrt{5})$].
From the above results, we find that the twofold degeneracy becomes generally
split but persists at the ${\rm\Gamma}$ point. Hence, the band inversion takes
place simultaneously for the two class D Hamiltonian blocks, until the
magnetic field strength becomes sufficiently strong to fully split the bands.
For the remaining parameters we used the values in Fig. 5. Figure 9: Band
structure for a Sk-SVC topological superconductor in the additional presence
of a combination of charge density wave terms with parameters
$W_{\bm{Q}_{1,2,\pm}}\equiv W_{1,2,\pm}$. Each charge density wave component
enters the BdG Hamiltonian in a similar fashion to
$\bm{M}_{1,2,\pm}\cdot\bm{\sigma}$, albeit that the respective Hamiltonian
term is diagonal in spin space and contains a $\tau_{3}$ matrix instead of
$\mathds{1}_{\tau}$. We consider the configuration where $W_{1}=We^{-i\pi/5}$,
$W_{2}=We^{-i\pi/4}$, $W_{+}=We^{-i\pi/3}$ and $W_{-}=We^{-i\pi/7}$. Black
(red) color shows the dispersions in the presence (absence) of the field.
(a)-(d) and (e)-(h) show two different cases for the strength of the charge
density waves. For the top (bottom) row $W=0.001$ ($W=0.1$). From the above
results, we find that the twofold degeneracy persists at all point for weak
values of $W$. However, also here, when $W$ becomes sufficiently strong, it
splits the twofold degeneracy in the entire MBZ. For the remaining parameters
we used the values considered in Fig. 5.

## References

  * [1] G. E. Volovik, The Universe in a Helium Droplet, Oxford University Press, New York, USA, ISBN 9780199564842, (2003).
  * [2] N. Read and D. Green, Paired States of Fermions in Two dimensions with Breaking of Parity and Time-Reversal Symmetries, and the Fractional Quantum Hall Effect, Phys. Rev. B 61, 10267 (2000). 10.1103/PhysRevB.61.10267.
  * [3] G. E. Volovik, Fermion zero modes on vortices in chiral superconductors, JETP Lett. 70, 609 (1999). 10.1134/1.568223.
  * [4] A. Y. Kitaev, Unpaired Majorana Fermions in Quantum Wires, Phys. Usp. 44, 131 (2001). 10.1070/1063-7869/44/10S/S29.
  * [5] M. Z. Hasan and C. L. Kane, Colloquium: Topological Insulators, Rev. Mod. Phys. 82, 3045 (2010). 10.1103/RevModPhys.82.3045.
  * [6] X.-L. Qi and S.-C. Zhang, Topological Insulators and Superconductors, Rev. Mod. Phys. 83, 1057 (2011). 10.1103/RevModPhys.83.1057.
  * [7] K. Flensberg, F. von Oppen, and A. Stern, Engineered platforms for topological superconductivity and Majorana zero modes, Nat. Rev. Mater. 6, 944 (2021). 10.1038/s41578-021-00336-6.
  * [8] A. O. Zlotnikov, M. S. Shustin, and A. D. Fedoseev, Aspects of topological superconductivity in 2D systems: noncollinear magnetism, skyrmions, and higher-order topology, J. Supercond. Nov. Magn. 34, 3053 (2021), (2021). 10.1007/s10948-021-06029-z.
  * [9] U. Güngördü and A. A. Kovalev, Majorana bound states with chiral magnetic textures, J. Appl. Phys. 132, 041101 (2022). 10.1063/5.0097008.
  * [10] M. Mostovoy, Ferroelectricity in Spiral Magnets, Phys. Rev. Lett. 96, 067601 (2006). 10.1103/PhysRevLett.96.067601.
  * [11] B. Braunecker, G. I. Japaridze, J. Klinovaja, and D. Loss, Spin-selective Peierls transition in interacting 1D conductors with spin-orbit interaction, Phys. Rev. B 82, 045127 (2010). 10.1103/PhysRevB.82.045127.
  * [12] T.-P. Choy, J. M. Edge, A. R. Akhmerov, and C. W. J. Beenakker, Majorana Fermions Emerging from Magnetic Nanoparticles on a Superconductor without Spin-Orbit Coupling, Phys. Rev. B 84, 195442 (2011). 10.1103/PhysRevB.84.195442.
  * [13] U. K. Rössler, A. N. Bogdanov, and C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, pages 797 (2006). 10.1038/nature05056.
  * [14] I. Martin and C. D. Batista, Itinerant Electron-Driven Chiral Magnetic Ordering and Spontaneous Quantum Hall Effect in Triangular Lattice Models, Phys. Rev. Lett. 101, 156402 (2008). 10.1103/PhysRevLett.101.156402.
  * [15] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion Lattice in a Chiral Magnet, Science 323, 915 (2009). 10.1126/science.1166767.
  * [16] 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). 10.1038/nature09124.
  * [17] 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 FeGe, Nat. Mater. 10, 106 (2011). 10.1038/nmat2916.
  * [18] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions, Nat. Phys. 7, 713 (2011). 10.1038/nphys2045.
  * [19] T. Okubo, S. Chung, and H. Kawamura, Multiple-q states and skyrmion lattice of the triangular-lattice Heisenberg antiferromagnet under magnetic fields, Phys. Rev. Lett. 108, 017206 (2012). 10.1103/PhysRevLett.108.017206.
  * [20] A. O. Leonov and M. Mostovoy, Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet, Nat. Commun. 6, 8275 (2015). 10.1038/ncomms9275.
  * [21] R. Ozawa, S. Hayami, and Y. Motome, Zero-Field Skyrmions with a High Topological Number in Itinerant Magnets, Phys. Rev. Lett. 118, 147205 (2017). 10.1103/PhysRevLett.118.147205.
  * [22] M. Hervé, B. Dupé, R. Lopes, M. Böttcher, M. D. Martins, T. Balashov, L. Gerhard, J. Sinova, and W. Wulfhekel, Stabilizing spin spirals and isolated skyrmions at low magnetic field exploiting vanishing magnetic anisotropy, Nat. Commun. 9, 2198 (2018). 10.1038/s41467-018-04680-0.
  * [23] J. Spethmann, S. Meyer, K. von Bergmann, R. Wiesendanger, S. Heinze, and A. Kubetzka, Discovery of Magnetic Single- and Triple-$\bm{q}$ States in Mn/Re(0001), Phys. Rev. Lett. 124, 227203 (2020). 10.1103/PhysRevLett.124.227203.
  * [24] S. Hayami and Y. Motome, Topological spin crystals by itinerant frustration, J. Phys.: Condens. Matter 33, 443001 (2021). 10.1088/1361-648X/ac1a30.
  * [25] S. Hayami, Rectangular and square skyrmion crystals on a centrosymmetric square lattice with easy-axis anisotropy, Phys. Rev. B 105, 174437 (2022). 10.1103/PhysRevB.105.174437.
  * [26] M. Kjaergaard, K. Wölms, and K. Flensberg, Majorana Fermions in Superconducting Nanowires without Spin-Orbit Coupling Phys. Rev. B 85, 020503(R) (2012). 10.1103/PhysRevB.85.020503.
  * [27] I. Martin and A. F. Morpurgo, Majorana fermions in superconducting helical magnets, Phys. Rev. B 85, 144505 (2012). 10.1103/PhysRevB.85.144505.
  * [28] J. Klinovaja and D. Loss, Giant Spin-Orbit Interaction Due to Rotating Magnetic Fields in Graphene Nanoribbons, Phys. Rev. X 3, 011008 (2013). 10.1103/PhysRevX.3.011008.
  * [29] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, and A. Yazdani, Proposal for Realizing Majorana Fermions in Chains of Magnetic Atoms on a Superconductor, Phys. Rev. B 88, 020407(R) (2013). 10.1103/PhysRevB.88.020407.
  * [30] P. Kotetes, Classification of Engineered Topological Superconductors, New J. Phys. 15, 105027 (2013). 10.1088/1367-2630/15/10/105027.
  * [31] S. Nakosai, Y. Tanaka, and N. Nagaosa, Two-Dimensional P-Wave Superconducting States with Magnetic Moments on a Conventional S-Wave Superconductor, Phys. Rev. B 88, 180503(R) (2013). 10.1103/PhysRevB.88.180503.
  * [32] B. Braunecker and P. Simon, Interplay between Classical Magnetic Moments and Superconductivity in Quantum One-Dimensional Conductors: Toward a Self-Sustained Topological Majorana Phase, Phys. Rev. Lett. 111, 147202 (2013). 10.1103/PhysRevLett.111.147202.
  * [33] J. Klinovaja, P. Stano, A. Yazdani, and D. Loss, Topological Superconductivity and Majorana Fermions in RKKY Systems, Phys. Rev. Lett. 111, 186805 (2013). 10.1103/PhysRevLett.111.186805.
  * [34] M. M. Vazifeh, and M. Franz, Self-Organized Topological State with Majorana Fermions, Phys. Rev. Lett. 111, 206802 (2013). 10.1103/PhysRevLett.111.206802.
  * [35] F. Pientka, L. I. Glazman, and F. von Oppen, Topological Superconducting Phase in Helical Shiba Chains, Phys. Rev. B 88, 155420 (2013). 10.1103/PhysRevB.88.155420.
  * [36] K. Pöyhönen, A. Westström, J. Röntynen, and T. Ojanen, Majorana States in Helical Shiba Chains and Ladders, Phys. Rev. B 89, 115109 (2014). 10.1103/PhysRevB.89.115109.
  * [37] F. Pientka, L. I. Glazman, and F. von Oppen, Unconventional topological phase transitions in helical Shiba chains, Phys. Rev. B 89, 180505(R) (2014). 10.1103/PhysRevB.89.180505.
  * [38] D. Mendler, P. Kotetes, and G. Schön, Magnetic order on a topological insulator surface with warping and proximity-induced superconductivity, Phys. Rev. B 91, 155405 (2015). 10.1103/PhysRevB.91.155405.
  * [39] W. Chen and A. P. Schnyder, Majorana edge states in superconductor-noncollinear magnet interfaces, Phys. Rev. B 92, 214502 (2015). 10.1103/PhysRevB.92.214502.
  * [40] E. Mascot, J. Bedow, M. Graham, S. Rachel, and D. K. Morr, Topological Superconductivity in Skyrmion Lattices, Npj Quantum Mater. 6, 6 (2021). 10.1038/s41535-020-00299-x.
  * [41] J. Bedow, E. Mascot, T. Posske, G. S. Uhrig, R. Wiesendanger, S. Rachel, and D. K. Morr, Topological superconductivity induced by a triple-Q magnetic structure, Phys. Rev. B 102, 180504(R) (2020). 10.1103/PhysRevB.102.180504.
  * [42] D. Steffensen, M. H. Christensen, B. M. Andersen, and P. Kotetes, Topological Superconductivity Induced by Magnetic Texture Crystals, Phys. Rev. Research 4, 013225 (2022). 10.1103/PhysRevResearch.4.013225.
  * [43] H. Kim, A. Palacio-Morales, T. Posske, L. Rózsa, K. Palotás, L. Szunyogh, M. Thorwart, and R. Wiesendanger, Toward Tailoring Majorana Bound States in Artificially Constructed Magnetic Atom Chains on Elemental Superconductors, Sci. Adv. 4, eaar5251 (2018). 10.1126/sciadv.aar5251.
  * [44] A. Kubetzka, J. M. Bürger, R. Wiesendanger, and K. von Bergmann, Towards skyrmion superconductor hybrid systems, Phys. Rev. Materials 4, 081401(R) (2020). 10.1103/PhysRevMaterials.4.081401.
  * [45] S. S. Pershoguba, S. Nakosai, and A. V. Balatsky, Skyrmion-induced bound states in a superconductor, Phys. Rev. B 94, 064513 (2016). 10.1103/PhysRevB.94.064513.
  * [46] G. Yang, P. Stano, J. Klinovaja, and D. Loss, Majorana bound states in magnetic skyrmions, Phys. Rev. B 93, 224505 (2016). 10.1103/PhysRevB.93.224505.
  * [47] K. M. D. Hals, M. Schecter, and M. S. Rudner, Composite Topological Excitations in Ferromagnet-Superconductor Heterostructures, Phys. Rev. Lett. 117, 017001 (2016). 10.1103/PhysRevLett.117.017001.
  * [48] U. Güngördü, S. Sandhoefner, and A. A. Kovalev, Stabilization and control of Majorana bound states with elongated skyrmions, Phys. Rev. B 97, 115136 (2018). 10.1103/PhysRevB.97.115136.
  * [49] S. Rex, I. V. Gornyi, and A. D. Mirlin, Majorana bound states in magnetic skyrmions imposed onto a superconductor, Phys. Rev. B 100, 064504 (2019). 10.1103/PhysRevB.100.064504.
  * [50] M. Garnier, A. Mesaros, and P. Simon, Topological superconductivity with orbital effects in magnetic skyrmion based heterostructures, arXiv:1909.12671 (2019). 10.48550/arXiv.1909.12671.
  * [51] M. Garnier, A. Mesaros, and P. Simon, Topological superconductivity with deformable magnetic skyrmions, Commun. Phys. 2, 126 (2019). 10.1038/s42005-019-0226-5.
  * [52] N. Bogdanov and D. A. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets. Sov. Phys. JETP 68, 101 (1989).
  * [53] A. Bogdanov and A. Hubert, Thermodynamically stable magnetic vortex states in magnetic crystals, JMMM 138, 255 (1994). 10.1016/0304-8853(94)90046-9.
  * [54] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions. Nature Nanotechnol. 8, 899 (2013). 10.1038/nnano.2013.243.
  * [55] A. N. Bogdanov and C. Panagopoulos, Physical foundations and basic properties of magnetic skyrmions, Nat. Rev. Phys. 2, 492 (2020). 10.1038/s42254-020-0203-7.
  * [56] 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, and C. Panagopoulos, Tunable room-temperature magnetic skyrmions in Ir/Fe/Co/Pt multilayers, Nat. Mater. 16, 898 (2017). 10.1038/nmat4934.
  * [57] A. P. Petrović, M. Raju, X. Y. Tee, A. Louat, I. Maggio-Aprile, R. M. Menezes, M. J. Wyszyński, N. K. Duong, M. Reznikov, Ch. Renner, M. V. Milo$\check{{\rm s}}$ević, and C. Panagopoulos, Skyrmion-(Anti)Vortex Coupling in a Chiral Magnet-Superconductor Heterostructure, Phys. Rev. Lett. 126, 117205 (2021). 10.1103/PhysRevLett.126.117205.
  * [58] A. Y. Kitaev, Fault-Tolerant Quantum Computation by Anyons, Ann. Phys. 303, 2 (2003). 10.1016/S0003-4916(02)00018-0.
  * [59] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian Anyons and Topological Quantum Computation, Rev. Mod. Phys. 80, 1083 (2008). 10.1103/RevModPhys.80.1083.
  * [60] D. A. Ivanov, Non-Abelian Statistics of Half-Quantum Vortices in P-Wave Superconductors, Phys. Rev. Lett. 86, 268 (2001). 10.1103/PhysRevLett.86.268.
  * [61] S. A. Díaz, J. Klinovaja, D. Loss, and S. Hoffman, Majorana bound states induced by antiferromagnetic skyrmion textures, Phys. Rev. B 104, 214501 (2021). 10.1103/PhysRevB.104.214501.
  * [62] I. Dzyaloshinskii, A Thermodynamic Theory of Weak Ferromagnetism of Antiferromagnetics, J. Phys. Chem. Solids 4, 241 (1958). 10.1016/0022-3697(58)90076-3.
  * [63] T. Moriya, Anisotropic Superexchange Interaction and Weak Ferromagnetism, Phys. Rev. 120, 91 (1960). 10.1103/PhysRev.120.91.
  * [64] S. L. Sondhi, A. Karlhede, S. A. Kivelson, and E. H. Rezayi, Skyrmions and the crossover from the integer to fractional quantum Hall effect at small Zeeman energies, Phys. Rev. B 47, 16419 (1993). 10.1103/PhysRevB.47.16419.
  * [65] H. A. Fertig, L. Brey, R. C${\rm\hat{o}}$té, and A. H. MacDonald, Charged spin-texture excitations and the Hartree-Fock approximation in the quantum Hall effect, Phys. Rev. B 50, 11018 (1994). 10.1103/PhysRevB.50.11018.
  * [66] S. E. Barrett, G. Dabbagh, L. N. Pfeiffer, K. W. West, and R. Tycko, Optically Pumped NMR Evidence for Finite-Size Skyrmions in GaAs Quantum Wells near Landau Level Filling $\nu=1$, Phys. Rev. Lett. 74, 5112 (1995). 10.1103/PhysRevLett.74.5112.
  * [67] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010). 10.1103/RevModPhys.82.1959.
  * [68] V. P. Gusynin, V. A. Miransky, and I. A. Shovkovy, Catalysis of Dynamical Flavor Symmetry Breaking by a Magnetic Field in 2+1 Dimensions, Phys. Rev. Lett. 73, 3499 (1994). 10.1103/PhysRevLett.73.3499.
  * [69] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Robustness of Majorana fermions in proximity-induced superconductors, Phys. Rev. B 82, 094522 (2010). 10.1103/PhysRevB.82.094522.
  * [70] A. C. Potter and P. A. Lee, Engineering a p+ip superconductor: Comparison of topological insulator and Rashba spin-orbit-coupled materials, Phys. Rev. B 83, 184520 (2011). 10.1103/PhysRevB.83.184520.
  * [71] D. Steffensen, B. M. Andersen, and P. Kotetes, Trapping Majorana zero modes in vortices of magnetic texture crystals coupled to nodal superconductors, Phys. Rev. B 104, 174502 (2021). 10.1103/PhysRevB.104.174502.
  * [72] H. Bruus and K. Flensberg, Many-Body Quantum Theory in Condensed Matter Physics, Oxford University Press, New York, USA, ISBN 9780198566335 (2004).
  * [73] C. Hickey, L. Cincio, Z. Papić, and A. Paramekanti, Emergence of chiral spin liquids via quantum melting of noncoplanar magnetic orders, Phys. Rev. B 96, 115115 (2017). 10.1103/PhysRevB.96.115115.
  * [74] Y.-F. Jiang and H.-C. Jiang, Topological Superconductivity in the Doped Chiral Spin Liquid on the Triangular Lattice, Phys. Rev. Lett. 125, 157002 (2020). 10.1103/PhysRevLett.125.157002.
  * [75] Y. Huang and D. N. Sheng, Topological Chiral and Nematic Superconductivity by Doping Mott Insulators on Triangular Lattice, Phys. Rev. X 12, 031009 (2022). 10.1103/PhysRevX.12.031009.
  * [76] A. Bose, A. Haldar, E. S. Sørensen, and A. Paramekanti, Chiral Broken Symmetry Descendants of the Kagomé Lattice Chiral Spin Liquid, arXiv:2204.10329 (2022). 10.48550/arXiv.2204.10329.
  * [77] D. Sen and R. Chitra, Large-U limit of a Hubbard model in a magnetic field: Chiral spin interactions and paramagnetism, Phys. Rev. B 51, 1922 (1995). 10.1103/PhysRevB.51.1922.
  * [78] X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Topological quantization of the spin Hall effect in two-dimensional paramagnetic semiconductors, Phys. Rev. B 74, 085308 (2006). 10.1103/PhysRevB.74.085308.
  * [79] F. Freimuth, R. Bamler, Y. Mokrousov, and A. Rosch, Phase-space Berry phases in chiral magnets: Dzyaloshinskii-Moriya interaction and the charge of skyrmions, Phys. Rev. B 88, 214409 (2013). 10.1103/PhysRevB.88.214409.
  * [80] P. Kotetes, A. Aperis, G. Varelogiannis Magnetic-field-induced chiral hidden order in URu2Si2, Philos. Mag. 94, 3789 (2014). 10.1080/14786435.2014.909614.
  * [81] R. M. Fernandes and A. J. Millis, Nematicity as a Probe of Superconducting Pairing in Iron-Based Superconductors, Phys. Rev. Lett. 111, 127001 (2013). 10.1103/PhysRevLett.111.127001.
  * [82] G. Livanas, A. Aperis, P. Kotetes, and G. Varelogiannis, Nematicity from mixed $S_{\pm}+d_{x^{2}-y^{2}}$ states in iron-based superconductors, Phys. Rev. B 91, 104502 (2015). 10.1103/PhysRevB.91.104502.
  * [83] J. Shi, G. Vignale, D. Xiao, and Q. Niu, Quantum Theory of Orbital Magnetization and Its Generalization to Interacting Systems, Phys. Rev. Lett. 99, 197202 (2007). 10.1103/PhysRevLett.99.197202.
  * [84] P. Kotetes, Topological Insulators, Morgan & Claypool, San Rafael, USA, ISBN 9781643276540 (2019). 10.1088/978-1-68174-517-6.
  * [85] R. B. Laughlin, Magnetic Induction of $d_{x^{2}-y^{2}}+id_{xy}$ Order in High-Tc Superconductors, Phys. Rev. Lett. 80, 5188 (1998). 10.1103/PhysRevLett.80.5188.
  * [86] J.-X. Zhu and A. V. Balatsky, Field induced $d_{x^{2}-y^{2}}+id_{xy}$ state in d-density-wave metals, Phys. Rev. B 65, 132502 (2002). 10.1103/PhysRevB.65.132502.
  * [87] P. Kotetes and G. Varelogiannis, Chirality Induced Tilted-Hill Giant Nernst Signal, Phys. Rev. Lett. 104, 106404 (2010). 10.1103/PhysRevLett.104.106404.
  * [88] J. W. Negele and H. Orland, Quantum Many-Particle Systems, CRC Press, Boca Raton, USA, ISBN 9780429497926 (1998). 10.1201/9780429497926.
  * [89] R. Takagi, N. Matsuyama, V. Ukleev, L. Yu, J. S. White, S. Francoual, J. R. L. Mardegan, S. Hayami, H. Saito, K. Kaneko, K. Ohishi, Y. ${\rm\bar{O}}$nuki, T. Arima, Y. Tokura, T. Nakajima, and S. Seki, Square and rhombic lattices of magnetic skyrmions in a centrosymmetric binary compound, Nat. Commun. 13, 1472 (2022). 10.1038/s41467-022-29131-9.
  * [90] M. H. Christensen, B. M. Andersen, and P. Kotetes, Unravelling incommensurate magnetism and its emergence in iron-based superconductors, Phys. Rev. X 8, 041022 (2018). 10.1103/PhysRevX.8.041022.
  * [91] J. Lorenzana, G. Seibold, C. Ortix, and M. Grilli, Competing Orders in FeAs Layers, Phys. Rev. Lett. 101, 186402 (2008). 10.1103/PhysRevLett.101.186402.
  * [92] L. Fu, Topological Crystalline Insulators, Phys. Rev. Lett. 106, 106802 (2011). 10.1103/PhysRevLett.106.106802.
  * [93] M. H. Christensen, D. D. Scherer, P. Kotetes, and B. M. Andersen, Role of multiorbital effects in the magnetic phase diagram of iron pnictides, Phys. Rev. B 96, 014523 (2017). 10.1103/PhysRevB.96.014523.
  * [94] A. Altland and M. R. Zirnbauer, Nonstandard Symmetry Classes in Mesoscopic Normal-Superconducting Hybrid Structures, Phys. Rev. B 55, 1142 (1997). 10.1103/PhysRevB.55.1142.
  * [95] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008). 10.1103/PhysRevB.78.195125.
  * [96] A. Kitaev, Periodic table for topological insulators and superconductors, AIP Conf. Proc. 1134, 22 (2009). 10.1063/1.3149495.
  * [97] S. Ryu, A. P. Schnyder, A. Furusaki and A. W. W. Ludwig, Topological Insulators and Superconductors: Ten-Fold Way and Dimensional Hierarchy, New J. Phys. 12, 065010 (2010). 10.1088/1367-2630/12/6/065010.
  * [98] J. C. Y. Teo and C. L. Kane, Topological Defects and Gapless Modes in Insulators and Superconductors, Phys. Rev. B 82, 115120 (2010). 10.1103/PhysRevB.82.115120.
  * [99] S. Fujimoto, Topological order and non-Abelian statistics in noncentrosymmetric s-wave superconductors, Phys. Rev. B 77, 220501(R) (2008). 10.1103/PhysRevB.77.220501.
  * [100] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Generic New Platform for Topological Quantum Computation Using Semiconductor Heterostructures, Phys. Rev. Lett. 104, 040502 (2010). 10.1103/PhysRevLett.104.040502.
  * [101] Z. Dong and L. Levitov, Chiral Stoner magnetism in Dirac bands, arXiv:2208.02051 (2022). 10.48550/arXiv.2208.02051.
  * [102] W. R. Meier, Q.-P. Ding, A. Kreyssig, S. L. Bud’ko, A. Sapkota, K. Kothapalli, V. Borisov, R. Valentí, C. D. Batista, P. P. Orth, R. M. Fernandes, A. I. Goldman, Y. Furukawa, A. E. Böhmer, and P. C. Canfield, Hedgehog spin-vortex crystal stabilized in a hole-doped iron-based superconductor, npj Quant. Mat. 3, 5 (2018). 10.1038/s41535-017-0076-x.
  * [103] N. Ni, M. E. Tillman, J.-Q. Yan, A. Kracher, S. T. Hannahs, S. L. Bud’ko, and P. C. Canfield. Effects of Co substitution on thermodynamic and transport properties and anisotropic ${H}_{c2}$ in $\text{Ba}{({\text{Fe}}_{1-x}{\text{Co}}_{x})}_{2}{\text{As}}_{2}$ single crystals. Phys. Rev. B 78, 214515 (2008). 10.1103/PhysRevB.78.214515.