# Magnetic interactions in AB-stacked kagome lattices: magnetic structure,
symmetry, and duality

A. Zelenskiy  Department of Physics and Atmospheric Science, Dalhousie
University, Halifax, Nova Scotia, Canada B3H 3J5  T. L. Monchesky  Department
of Physics and Atmospheric Science, Dalhousie University, Halifax, Nova
Scotia, Canada B3H 3J5  M. L. Plumer  Department of Physics and Atmospheric
Science, Dalhousie University, Halifax, Nova Scotia, Canada B3H 3J5
Department of Physics and Physical Oceanography, Memorial University of
Newfoundland, St. John’s, Newfoundland, A1B 3X7, Canada  B. W. Southern
Department of Physics and Astronomy, University of Manitoba, Winnipeg,
Manitoba, Canada R3T 2N2

###### Abstract

We present the results of an extensive study of the phase diagram and spin
wave excitations for a general spin model on a hexagonal AB-stacked kagome
system. The boundaries of the magnetic phases are determined via a combination
of numerical (Monte Carlo) and analytical (Luttinger-Tisza) methods. We also
determine the phase coexistence regions by considering the instabilities in
the spin wave spectra. Depending on the strength of the spin-orbit coupling
(SOC), some spin and lattice rotations become decoupled, leading to
considerably larger symmetry groups than typical magnetic groups. Thus, we
provide a detailed symmetry description of the magnetic Hamiltonian with
negligible, weak, and intermediate strength of SOC. The spin symmetry in these
three cases has a strong effect on the splittings observed in the spin
excitation spectra. We further identify a number of self-duality
transformations that map the Hamiltonian onto itself. These transformations
describe the symmetry of the parameter space and provide an exact mapping
between the properties of different magnetic orders and lead to accidental
degeneracies. Finally, we discuss the physical relevance of our findings in
the context of $\mathrm{Mn}_{3}X$ compounds.

###### pacs:

Valid PACS appear here

##  I Introduction

Magnetic materials on lattices comprised of equilateral triangles continue to
attract attention from both the experimental and theoretical viewpoint due to
the richness of physical properties provided by geometric frustration Sadoc
and Mosseri (1999); Diep (2005); Moessner and Ramirez (2006); Lacroix  _et
al._ (2011). Frustrated antiferromagnets are often characterized by non-
collinear spin configurations, since the topology of the lattice forbids a
conventional Néel order. Among the many realizations of frustrated systems,
kagome magnets with antiferromagnetic nearest-neighbor (nn) interactions have
become a staple example of systems with macroscopic degeneracy in the ground
state, giving rise to high sensitivity of various symmetry-breaking
perturbations Chalker  _et al._ (1992); Harris  _et al._ (1992); Huse and
Rutenberg (1992); Balents  _et al._ (2002); Fu  _et al._ (2015); Fujihala  _et
al._ (2020).

Recently, a family of magnetic materials with AB-stacked kagome layer
structure and a general formula $\mathrm{Mn}_{3}X$ have been experimentally
shown to host the anomalous Hall effect (AHE) as well as anomalous Nernst
effect (ANE) Chen  _et al._ (2014); Nakatsuji  _et al._ (2015); Kiyohara  _et
al._ (2016); Nayak  _et al._ (2016); Ikhlas  _et al._ (2017); Hong  _et al._
(2020). These discoveries prompted recent theoretical and experimental studies
of the magnetic properties in $\mathrm{Mn}_{3}X$ magnets Liu and Balents
(2017); Park  _et al._ (2018); Li  _et al._ (2019); Reichlova  _et al._
(2019); Soh  _et al._ (2020); Chen  _et al._ (2020); Zelenskiy  _et al._
(2021a).

Over the last decade, magnets with strong SOC have been under an intense
investigation motivated by an ongoing search for unconventional magnetic
phases. On the one hand, these include strongly correlated disordered states
with large degeneracy in the ground states, such as various types of spin
liquids Kitaev (2006); Ran  _et al._ (2007); Balents (2010); Yan  _et al._
(2011); Messio  _et al._ (2012); Bauer  _et al._ (2014); Savary and Balents
(2016); Takagi  _et al._ (2019). On the other hand, however, ordered magnetic
textures such as skyrmion lattices Bogdanov and Hubert (1994); Rößler  _et
al._ (2006), and multi-$\mathbf{Q}$ structures consisting of linear
superpositions of non-colinear spin density waves Hayami and Motome (2021);
Leonov and Mostovoy (2015), have attracted increasing interest in the
literature. These non-trivial magnetic orders often have nonzero scalar spin
chirality, which serves as a source of the emergent electromagnetic fields,
within the Berry phase formalism Schulz  _et al._ (2012); Everschor-Sitte and
Sitte (2014), giving rise to important transport properties, such as
topological Hall effect (THE) Neubauer  _et al._ (2009); Bruno  _et al._
(2004); Everschor-Sitte and Sitte (2014); Kurumaji  _et al._ (2019); He  _et
al._ (2022) and spin Hall effect (SHE) Hirsch (1999); Chen and Byrnes (2019).
More recently, magnetic frustration was identified as one of the stabilizing
factors for multi-$\mathbf{Q}$ spin configurations, leading to magnetic orders
beyond those typically observed in chiral ferromagnets Hayami and Motome
(2021).

In the case of $\mathrm{Mn}_{3}X$ compounds, both the THE and SHE have been
observed experimentally, and studies have also established that both the AHE,
and ANE, as well as the magnetic structure are strongly anisotropic, implying
that the spin couplings beyond isotropic exchange are crucial for
understanding the magnetic properties of these systems Kiyohara  _et al._
(2016); Nayak  _et al._ (2016); Zhang  _et al._ (2017). Previous studies have
also established that different types of anisotropic interactions compete with
each other leading to additional frustration Liu and Balents (2017); Soh  _et
al._ (2020); Chen  _et al._ (2020); Zelenskiy  _et al._ (2021a). These facts
motivate a systematic study of the magnetic ground state properties in an
extended parameter space.

In our previous work in Ref. Zelenskiy  _et al._ (2021a) (hereafter referred
to as Ref. I), we derived a magnetic model for these magnetic compounds using
general symmetry principles. Apart from the typical exchange couplings, the
symmetry-allowed terms consist of various anisotropic couplings, including
Dzyaloshinskii-Moriya (DM), bond-dependent anisotropic exchange, and single-
ion anisotropy (SIA). These anisotropic interactions arise from the coupling
of spins to the underlying lattice via SOC. In Ref. I, we have studied the
ground state properties and spin wave excitations of this model relevant to
the magnetism of $\mathrm{Mn}_{3}X$ systems. In particular, this work provided
a detailed analysis of the interplay between various anisotropic interactions
and their effect on the static and dynamic properties of $\mathrm{Mn}_{3}X$
compounds. We determined that both SIA and bond-dependent anisotropy compete
with the DM interactions, leading to an induced magnetic moment and an
excitation spectrum with broken six-fold symmetry. However, we also found that
these properties are extremely sensitive to both the signs and relative
magnitudes of the anisotropic interactions.

In the present work, we investigate semi-classical ground state properties of
hexagonal AB-stacked kagome systems for an extended range of magnetic
interactions using a combination of analytical Luttinger-Tisza (LT) and
numerical Monte-Carlo (MC) techniques.

Figure 1: (a) A sketch of the AB-stacked kagome lattice. The conventional unit
cell defined by the lattice vectors $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ is
shaded in blue, while the hexagon formed by the four sublattices is shaded in
pink. (b) The three types of nn and nnn interactions appearing in (a). Solid
black lines indicate the nn in-plane interactions, while the dashed red and
green lines represent the nn and nnn out-of-plane interactions respectively
($\bm{\mathcal{A}}_{i}=\\{J_{i},D_{i},A_{i}^{(z)},A_{i}^{(xy)}\\}$). (c) An
enlarged diagram of the unit cell convention used in this work.
$\mathbf{S}_{i}$ label the spins on the six sublattices. Spins 1, 2, and 3
reside in layer A, while spins 4, 5, and 6 are in layer B.

To tackle the large parameter space of the magnetic model, we group the
magnetic interactions based on the effective spin symmetry imposed by the
relative strength of the SOC. We determine three SOC symmetry regimes and
present the corresponding group structure, along with the relevant irreducible
representations (irreps). Furthermore, for each of these symmetry regimes, we
identify a set of self-duality transformations that reduce the number of
independent points in the parameter space. We find that in the weak SOC limit
the magnetic Hamiltonian has the largest number of dualities which comprise a
group with non-Abelian structure. The numerical and analytical calculations
reveal a variety of magnetic phases, including single-$\mathbf{Q}$,
multi-$\mathbf{Q}$, as well as more complicated structures with delocalized
structure factors in the Brillouin zone. We parameterize the magnetic phases
and study the elementary spin wave excitations. Finally, we analyze the
effects of bond-dependent exchange and SIA on the magnetic phases stabilized
by exchange and DM interactions, and discuss the implications for the case of
$\mathrm{Mn}_{3}X$ systems.

The rest of this paper is organized as follows. In Sec. II we introduce the
magnetic model and briefly outline the analytical and numerical methods. In
Sec. III we identify the connection between the strength of the SOC and the
effective symmetry of the magnetic system. Next, in Sec. IV, the self-duality
transformations are introduced and derived for each SOC symmetry case. The
magnetic ground state phase diagrams for models with exchange and DM
interactions are presented in Sec. V.

The distinct types of magnetic order are described in Sec. VI, and the
corresponding spin-wave excitation spectra are given in Sec. VII. The effects
of anisotropic interactions are discussed in Sec. VIII. Finally, Section IX is
devoted to concluding remarks and a summary of the results.

##  II Model and Methods

###  II.1  Model

The spin Hamiltonian for $\mathrm{Mn}_{3}\mathrm{X}$-type AB-stacked kagome
lattice systems has been derived in Ref. I from symmetry principles and
contains four different types of interactions:

$\displaystyle\mathcal{H}=\mathcal{H}_{J}+\mathcal{H}_{D}+\mathcal{H}_{A}+\mathcal{H}_{K}$
(1)
$\displaystyle\mathcal{H}_{J}=\frac{1}{2}\sum_{\mathbf{r}\mathbf{r^{\prime}}}\sum_{ij}J_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{S}_{i}(\mathbf{r})\cdot\mathbf{S}_{j}(\mathbf{r^{\prime}})$
$\displaystyle\mathcal{H}_{D}=\frac{1}{2}\sum_{\mathbf{r}\mathbf{r^{\prime}}}\sum_{ij}D_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{\hat{z}}\cdot\left(\mathbf{S}_{i}(\mathbf{r})\times\mathbf{S}_{j}(\mathbf{r^{\prime}})\right)$
$\displaystyle\mathcal{H}_{A}=\frac{1}{2}\sum_{\mathbf{r}\mathbf{r^{\prime}}}\sum_{ij}\sum_{\alpha}A_{ij\alpha}(\mathbf{r}-\mathbf{r^{\prime}})\left(\mathbf{\hat{n}}_{i\alpha}\cdot\mathbf{S}_{i}(\mathbf{r})\right)\left(\mathbf{\hat{n}}_{j\alpha}\cdot\mathbf{S}_{j}(\mathbf{r^{\prime}})\right),$
$\displaystyle\mathcal{H}_{K}=\sum_{\mathbf{r}}\sum_{i}\sum_{\alpha}K_{\alpha}\left(\mathbf{\hat{n}}_{i\alpha}\cdot\mathbf{S}_{i}(\mathbf{r})\right)^{2}.$

Here, $\mathcal{H}_{J}$ is the isotropic Heisenberg exchange,
$\mathcal{H}_{D}$ is the DM interaction, $\mathcal{H}_{K}$ is the SIA, and
$\mathcal{H}_{A}$ is the symmetric anisotropic exchange interaction. Sum
indices $\mathbf{r}$, $\mathbf{r^{\prime}}$ label unit cells,
$i,j\in\\{1,...,6\\}$ label atoms in each unit cell (Fig. 1 (c)), and
$\alpha\in\\{x,y,z\\}$ labels the spin vector components. Vectors
$\mathbf{n}_{i\alpha}$ represent local anisotropy axes and can be written as

$\mathbf{\hat{n}}_{ix}=\begin{bmatrix}\cos{\alpha_{i}}\\\ \sin{\alpha_{i}}\\\
0\end{bmatrix},\mathbf{\hat{n}}_{iy}=\begin{bmatrix}-\sin{\alpha_{i}}\\\
\cos{\alpha_{i}}\\\ 0\end{bmatrix},\mathbf{\hat{n}}_{iz}=\begin{bmatrix}0\\\
0\\\ 1\end{bmatrix},$ (2)

where $\alpha_{i}$ give the angle of the anisotropy axes with respect to the
global $x$-direction. In this paper, we will restrict our attention to in- and
out-of-plane nearest-neighbor (nn) and out-of-plane next-nearest-neighbor
(nnn) interactions shown in Fig. 1 (a). The interaction labels are based on
$\mathrm{Mn}_{3}X$ bond distances, with index 1 labeling out-of-plane nn
interactions, 2 and 3 labeling the in-plane nn, and 4 and 5 labeling the out-
of-plane nnn. In this work, we will ignore the breathing anisotropy Chen  _et
al._ (2020), in order to simplify the analysis, in which case
$\bm{\mathcal{A}}_{2}=\bm{\mathcal{A}}_{3}$, and
$\bm{\mathcal{A}}_{4}=\bm{\mathcal{A}}_{5}$ (Fig. 1).

It has recently been shown that an effective spin Hamiltonians of the form (1)
can be derived through perturbation theory from lattice Kondo model with SOC
Hayami and Motome (2018); Akagi  _et al._ (2012); Ghosh  _et al._ (2016). The
isotropic exchange terms in this case correspond to the Ruderman-Kittel-
Kasuya-Yosida (RKKY) interactions while the remaining anisotropic spin
interactions originate from the SOC. However, while the DM interactions depend
linearly on the strength of the SOC, the SIA and anisotropic exchange yield
quadratic dependence Moriya (1960); Hayami and Motome (2018). Therefore,
within the perturbation theory, we generally expect the magnitudes of the
latter two types of the anisotropic interactions to be smaller than that of
the DM interaction.

It is sometimes useful to write the Hamiltonian in the general quadratic form:

$\mathcal{H}=\frac{1}{2}\sum_{\mathbf{r}\mathbf{r}^{\prime}}\sum_{ij}\mathbf{S}^{T}_{i}(\mathbf{r})\bm{\mathcal{A}}_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{S}_{j}(\mathbf{r}^{\prime}),$
(3)

where the coupling matrix $\bm{\mathcal{A}}_{ij}(\bm{\delta})$ between two
inequivalent spins is given by

$\bm{\mathcal{A}}_{ij}(\bm{\delta})=\begin{bmatrix}V_{ij}^{+}(\bm{\delta})&W_{ij}^{+}(\bm{\delta})&0\\\
W_{ij}^{-}(\bm{\delta})&V_{ij}^{-}(\bm{\delta})&0\\\
0&0&X_{ij}(\bm{\delta})\end{bmatrix},$ (4)

with

$\displaystyle V_{ij}^{\pm}(\bm{\delta})$
$\displaystyle=J_{ij}(\bm{\delta})\pm
A^{(xy)}_{ij}(\bm{\delta})\cos(\bar{\alpha}_{ij}),$ $\displaystyle
W_{ij}^{\pm}(\bm{\delta})$ $\displaystyle=\pm
D_{ij}(\bm{\delta})+A^{(xy)}_{ij}(\bm{\delta})\sin(\bar{\alpha}_{ij}),$
$\displaystyle X_{ij}(\bm{\delta})$
$\displaystyle=J_{ij}(\bm{\delta})+A^{(z)}_{ij}(\bm{\delta}).$

Here, $\bm{\delta}=\mathbf{r}-\mathbf{r}^{\prime}$,
$\bar{\alpha}_{ij}=\alpha_{i}+\alpha_{j}$, $A_{ij}^{(xy)}$ and $A_{ij}^{(z)}$
are the two types of anisotropic exchange parameters allowed by symmetry. The
part of the coupling matrix that corresponds to the single-ion anisotropy is
given by

$\bm{\mathcal{A}}_{ii}(0)=K^{+}\bm{\mathbb{I}}+\begin{bmatrix}K^{-}\cos
2\alpha_{i}&K^{-}\sin 2\alpha_{i}&0\\\ K^{-}\sin 2\alpha_{i}&-K^{-}\cos
2\alpha_{i}&0\\\ 0&0&K_{Z}\end{bmatrix},$ (5)

where $\bm{\mathbb{I}}$ is a $3\times 3$ identity matrix, $K^{+}=K_{x}+K_{y}$,
$K^{-}=K_{x}-K_{y}$, and $K_{Z}=2K_{z}-K^{+}$. Note that the first term in
this expression is just a constant energy shift, and therefore $K^{+}$ can be
ignored in further calculations.

###  II.2 Luttinger-Tisza method

To provide an initial characterization of the classical ground states, we
define lattice Fourier transforms of the spin vectors as

$\mathbf{S}_{i}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum_{\mathbf{q}}\mathbf{S}_{i}(\mathbf{q})e^{-i\mathbf{q}\cdot\mathbf{r}},$
(6)

where $\mathbf{q}$ are the wave-vectors restricted to the first Brillouin
zone, and $\mathbf{S}(\mathbf{q})$ are the Fourier amplitudes. From Eq. (1)
and (3), the total energy of the system is

$\mathcal{H}=\frac{1}{2}\sum_{\mathbf{q}}\sum_{ij}\mathbf{S}^{T}_{i}(\mathbf{q})\bm{\mathcal{A}}_{ij}(\mathbf{q})\mathbf{S}_{j}(-\mathbf{q}),$
(7)

where the Fourier transform of the magnetic interactions are given by

$\bm{\mathcal{A}}_{ij}(\mathbf{q})=\sum_{\bm{\delta}}\bm{\mathcal{A}}_{ij}(\bm{\delta})e^{-i\mathbf{q}\cdot\bm{\delta}}.$
(8)

The true ground state is calculated by minimizing (7), subject to local
normalization constraints $|\mathbf{S}_{i}(\mathbf{r})|=1$ for all spins in
the system. This strong constraint significantly complicates the problem and
often makes it impossible to solve. Instead, the Luttinger-Tisza (LT) method
Luttinger and Tisza (1946) replaces the local constraints by a global
constraint, whereby the sum of all spin magnitudes is set to be equal to the
number of spins. This simplification allows one to recast the energy
minimization in the form of an eigenvalue problem,

$\sum_{j}\bm{\mathcal{A}}_{ij}(\mathbf{q})\mathbf{S}_{j}(\mathbf{q})=\varepsilon(\mathbf{q})\mathbf{S}_{i}(\mathbf{q})$
(9)

Under ideal circumstances, the smallest eigenvalue
$\varepsilon_{\text{LT}}(\mathbf{q})$ and the corresponding eigenvector give
the ground state of the magnetic system. However, this method often produces
unphysical solutions for systems with strong anisotropic interactions and
multiple sublattices Zaliznyak and Zhitomirsky (1995); Maksimov  _et al._
(2019). As a result, we used the LT method to determine the approximate
locations of the phase boundaries, as well as the lower bounds on the ground
state energies to guide the numerical calculations.

###  II.3 Monte Carlo

To identify the ground state magnetic configurations, we utilized classical
Monte Carlo (MC) simulations, which were carried out using standard local
heat-bath updates. To resolve the individual phases, we have used system sizes
ranging from $6^{3}$ to $24^{3}$ unit cells (1296 to 82944 spins respectively)
and between $10^{4}$ and $10^{6}$ MC steps. In each simulation, the
temperature is reduced down to $T\approx 10^{-6}$, to ensure energy
convergence. All phases presented in this work have $Q_{z}=0$, meaning that
every unit cell along the $c$-axis has exactly the same magnetic structure.
This fact allowed us to determine the phase boundaries using smaller system
sizes (between $6\times 6\times 2$ and $18\times 18\times 2$ unit cells). The
MC data is Fourier transformed to obtain the spin structure factor

$S(\mathbf{q})=\frac{1}{6}\sum_{ij}\big{\langle}\mathbf{S}_{i}(\mathbf{q})\cdot\mathbf{S}_{j}(-\mathbf{q})\big{\rangle}e^{i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})},$
(10)

where $\mathbf{r}_{i}$ are the positions of atoms inside of the unit cell.

###  II.4 Spin waves

The dynamics of the magnetic system can be described by the Landau-Lifshitz
equation

$\frac{\text{d}\mathbf{S}_{i}(\mathbf{r},t)}{\text{d}t}=\mathbf{H}_{i}(\mathbf{r},t)\times\mathbf{S}_{i}(\mathbf{r},t),$
(11)

where $\mathbf{H}_{i}(\mathbf{r})$ is the effective field at each site

$\mathbf{H}_{i}(\mathbf{r},t)=\sum_{\bm{\delta}}\sum_{j}\bm{\mathcal{A}}_{ij}(\bm{\delta})\mathbf{S}_{j}(\mathbf{r}+\bm{\delta},t).$
(12)

In this work, we will look for solutions of the linearized form of Eq. (11)
which correspond to the low-energy spin wave excitations. This is done by
first changing into a local coordinate system where the local $z$-components
are aligned with the ground state spin configuration. Note that as long as
this coordinate transformation is described by a local rotation
$\mathbf{U}_{i}(\mathbf{r})$, the dynamic evolution of the local spin
components can be written in the same form as Eq. (11), replacing
$\mathbf{S}_{i}(\mathbf{r},t)$ and $\mathbf{H}_{i}(\mathbf{r},t)$ with
$\widetilde{\mathbf{S}}_{i}(\mathbf{r},t)=\mathbf{U}_{i}(\mathbf{r})\mathbf{S}_{i}(\mathbf{r},t)$
and
$\widetilde{\mathbf{H}}_{i}(\mathbf{r},t)=\mathbf{U}_{i}(\mathbf{r})\mathbf{H}_{i}(\mathbf{r},t)$
respectively.

##  III Spin symmetry

When studying the properties of a magnetic model, it is important to take a
proper account of the symmetries that leave the Hamiltonian invariant. Apart
from the space group transformations imposed by the underlying crystalline
lattice, magnetic systems also include symmetries associated with spin
rotations and reflections. The most common spin symmetry is the time reversal,
$T$, which must be broken in order to establish magnetic order. A combination
of crystallographic group operations with the time reversal operator leads to
magnetic point and space groups. The addition of this one operator
significantly extends the number of symmetrically distinct systems: in three
dimensions there are 1651 magnetic space groups compared to 230 “regular”
crystal space groups Glazer and Burns (2013); Bradley and Cracknell (2010). In
materials with large SOC, the spins are typically pinned to the lattice,
meaning that for each lattice transformation there is a corresponding spin
transformation. The symmetry group of the Hamiltonian is then the paramagnetic
group which is a direct product

$\mathcal{G}_{\text{SOC}}=\mathcal{G}_{L}\otimes Z_{2}^{(T)},$ (13)

where $\mathcal{G}_{L}$ is the space group and $Z_{2}^{(T)}=\\{E,T\\}$.
However, in the limit of decoupled spin and orbital degrees of freedom one has
a magnetic system with isotropic Heisenberg exchange interactions which are
invariant under all global spin rotations. In fact, as was noted by Brinkman
and Elliott, there are many instances of magnetic systems where the spins are
at least partially decoupled from the lattice Brinkman  _et al._ (1966);
Brinkman and Elliott (1966); Brinkman (1967). The symmetry of these systems is
then described by the spin point and space groups Litvin and Opechowski
(1974); Litvin (1977) which are typically much larger than the corresponding
magnetic groups. Although the effects of the extended spin symmetry operations
have previously been considered almost exclusively in the context of either
isotropic spins or single-ion anisotropy, their importance in the intermediate
cases, which include DM and anisotropic exchange interactions, remains
relatively underrepresented and has only began gaining interest in recent
years Corticelli  _et al._ (2022); Liu  _et al._ (2022). For more information
about the spin space groups we refer the reader to reference Corticelli  _et
al._ (2022) which provides an excellent review of the subject. In compounds
with $\mathrm{Mn}_{3}X$-like structure, depending on the strength of the SOC,
one can identify three distinct cases for magnetic models depicted in Fig. 2.
Each case corresponds to a different group of spin symmetries.

In this section, we present the symmetry analysis of these three cases by
deriving the corresponding spin groups and determining the resulting irreps
for $\mathbf{Q}=0$. Note that the analysis presented here is for classical
magnetic moments (rotations in SO(3)), but could be readily extended to
quantum spin operators (rotations in SU(2)). The details of some derivations
can be found in the [Supplemental Material].

###  III.1 Decoupled case

Figure 2: A diagram illustrating the effects of SOC on the symmetry of the
spin Hamiltonian. In $\mathrm{Mn}_{3}X$ systems, we identify three SOC limits
resulting in distinct symmetry groups: decoupled, weak SOC, and intermediate
SOC. These correspond to, in the same order, the isotropic, XY-anisotropy, and
Ising anisotropy (top panel). The total spin group of the Hamiltonian and the
corresponding irreps are indicated by the colored blocks. The positive and
negative $z$-components of spins are indicated by the filled circles and
pluses respectively, and the blue and red colors of the spin components are
used to indicate parallel and antiparallel out-of-plane nnn (labels “g” and
“u” respectively). There are four irreps in the decoupled limit, which are
first split into ten irreps in the weak SOC limit and then further split into
twelve irreps in the intermediate SOC limit. The irreps with the out-of-plane
spin components (bottom row) are unchanged going from the weak to intermediate
SOC limit, and therefore have the same labels. The remaining planar irreps are
labeled using two labels corresponding to the weak and intermediate limits.

As mentioned previously, when the spin and orbital degrees of freedom are
completely decoupled, the magnetic interactions correspond to the isotropic
exchange, $\mathcal{H}=\mathcal{H}_{J}$, meaning that the spin Hamiltonian is
invariant with respect to pure space transformations (lattice site
permutations) and global spin rotations by an arbitrary angle. The crystal
symmetries pertaining to this work form a space group P6${}_{3}/mmc$, which we
denote as $\mathrm{D}_{6h}^{(L)}$ (indicating also the point group), and the
spin rotations combined with the time reversal symmetry form a group
$\text{SO}^{(S)}(3)\otimes Z_{2}^{(T)}$. Since the elements in these two
groups commute, the full group of the isotropic exchange Hamiltonian is simply

$\mathcal{G}_{J}=\mathrm{D}_{6h}^{(L)}\otimes\mathrm{SO}^{(S)}(3)\otimes
Z_{2}^{(T)}.$ (14)

As a result, the irreps of $\mathcal{G}_{J}$ are obtained by taking a direct
product of the irreps of spin and lattice degrees of freedom. The former is
related to the set of spherical harmonics with $l=1$ and has dimension 3,
while the latter depends on the periodicity of the magnetic structure. For
$\mathbf{Q}=0$, the symmetries of $\mathrm{D}_{6h}^{(L)}$ reduce to the point
group, and the irrep decomposition becomes $A_{1g}\oplus E_{2g}\oplus
B_{2u}\oplus E_{1u}$. Therefore, there are four irreps in the decomposition of
magnetic states: two of dimension 3 ($T_{g}$, $T_{u}$), and two of dimension 6
($Q_{g}$, $Q_{u}$), as illustrated in Fig. 2. The two triplets correspond to
the collinear configurations, while the two six-dimensional irreps are related
to the 120 degree configurations. The labels “g” and “u” are used to indicate
the parity of the irreps under the spatial inversion, as per usual group
theory notation.

$\mathcal{G}_{D}$ | $C_{1}(\phi)$ | $C_{6}(\phi)$ | $C_{3}(\phi)$ | $C_{2}(\phi)$ | $C_{2}^{\prime}$ | $C_{2}^{\prime\prime}$ | $I(\phi)$ | $IC_{6}(\phi)$ | $IC_{3}(\phi)$ | $IC_{2}(\phi)$ | $IC_{2}^{\prime}$ | $IC_{2}^{\prime\prime}$
---|---|---|---|---|---|---|---|---|---|---|---|---
$E_{g}^{(nm)}$ | $c^{(0)}_{n}$ | $c^{(m)}_{n}$ | $c^{(2m)}_{n}$ | $c^{(3m)}_{n}$ | 0 | 0 | $c^{(0)}_{n}$ | $c^{(m)}_{n}$ | $c^{(2m)}_{n}$ | $c^{(3m)}_{n}$ | 0 | 0
$E_{u}^{(nm)}$ | $c^{(0)}_{n}$ | $c^{(m)}_{n}$ | $c^{(2m)}_{n}$ | $c^{(3m)}_{n}$ | 0 | 0 | -$c^{(0)}_{n}$ | -$c^{(m)}_{n}$ | -$c^{(2m)}_{n}$ | -$c^{(3m)}_{n}$ | 0 | 0

Table 1: Characters of the two-dimensional irreps of group $\mathcal{G}_{D}$
that do not also appear in group $\mathrm{D}_{6h}\otimes Z_{2}^{(T)}$. For
simplicity, the elements with time reversal operator are not included in the
table. Here, the classes $C_{k}(\phi)$ have two elements,
$\\{C_{z}(\pm\phi,\pm\phi_{k})\\}$, which corresponds to a rotation around
$z$-axis of spin by an angle $\phi$ and lattice by angle
$\phi_{k}=\frac{2\pi}{k}$. $C_{2}^{\prime}$ and $C_{2}^{\prime\prime}$ include
simultaneous 180 degree rotations of lattice (same as in the $\mathrm{D}_{6h}$
point group) and spin with axes in the $xy$-plane. The second half of classes
is obtained by combining the first half with the spatial inversion operator.
The non-zero characters are defined as
$c_{n}^{(m)}=2\cos\left(n\phi+\frac{m\pi}{3}\right)$, where
$n\in\mathbb{Z}\setminus\\{0\\}$, $m\in\mathbb{Z}_{6}$.

###  III.2 Weak coupling case

As discussed in Sec. II, when the SOC is small but non-zero, it is reasonable
to assume that the DM interactions are the dominant type of anisotropy in the
system, giving $\mathcal{H}\approx\mathcal{H}_{J}+\mathcal{H}_{D}$. The
addition of DM coupling to the spin Hamiltonian significantly complicates the
symmetry analysis of the model. However, it is still possible to determine the
structure of the corresponding spin group as well as all irreps. We also note
that the SIA and bond-dependent interactions that couple the $z$-components of
spins do not change the symmetry of this spin group.

It can be shown that a spin rotation applied to a DM coupling between two
spins gives

$\displaystyle
D\mathbf{\hat{z}}\cdot[\mathbf{S}_{i}^{\prime}(\mathbf{r})\times\mathbf{S}_{j}^{\prime}(\mathbf{r}^{\prime})]$
$\displaystyle=D\mathbf{\hat{z}}\cdot[(\mathbf{M}\mathbf{S}_{i}(\mathbf{r}))\times(\mathbf{M}\mathbf{S}_{j}(\mathbf{r}^{\prime}))]$
$\displaystyle=D(\mathbf{M}^{T}\mathbf{\hat{z}})\cdot[\mathbf{S}_{i}(\mathbf{r})\times\mathbf{S}_{j}(\mathbf{r}^{\prime})],$
(15)

where $\mathbf{M}$ is the rotation matrix. Therefore, all rotations that leave
$\mathbf{\hat{z}}$ invariant belong to the symmetry group of DMI. This
constitutes a group of axial rotations in spin space $\mathrm{SO}^{(S)}(2)$,
implying XY-anisotropy. The complications arise from the fact that not all
lattice permutations leave the DM interaction invariant. In particular,
transformations that include $C_{2}$ rotations around axes parallel to the
kagome layers, and the corresponding reflections reverse the direction of the
bonds, which flips the sign of the DM vector. In order for these to become
proper symmetry operators, they must be combined with the corresponding spin
rotations/reflections. Since the spin and lattice operations are now coupled
and do not necessarily commute with each other, we can no longer write the
total group of the Hamiltonian as a direct product of spin and lattice
symmetry groups. Nevertheless, the group of DM coupling can then be written as
a semidirect product

$\mathcal{G}_{D}=\left(\mathrm{C}_{6h}^{(L)}\otimes\mathrm{SO}^{(S)}(2)\right)\ltimes\mathrm{C}_{2}^{SL}\otimes
Z_{2}^{(T)},$ (16)

where $\mathrm{C}_{6h}^{(L)}$ is the group of lattice symmetries that leave DM
interaction invariant, and $\mathrm{C}_{2}^{(SL)}$ has one non-trivial element
$C_{2}^{S}C_{2}^{L}$ that rotates both the lattice and spin components around
the $x$-axis. The derivation of the group structure and irreps of
$\mathcal{G}_{D}$ is given in the [Supplemental Material]. The irreps of
$\mathcal{G}_{D}$ include all irreps of a regular $\mathrm{D}_{6h}\otimes
Z_{2}^{(T)}$ magnetic group, and an infinite number of two-dimensional irreps,
as presented in table 1.

Intuitively, one can expect that the continuous axial rotational symmetry
would separate the $z$-components of the spins from the planar components,
while leaving the latter degenerate. The decomposition of a $\mathbf{Q}=0$
magnetic structure consists of ten irreps, four of which involve only the
$z$-components of spins ($A_{2g}\oplus E_{2g}\oplus B_{1u}\oplus E_{1u}$), and
six two-dimensional planar irreps ($E_{g}^{(10)}\oplus E_{g}^{(12)}\oplus
E_{g}^{(14)}\oplus E_{u}^{(11)}\oplus E_{u}^{(13)}\oplus E_{u}^{(15)}$), as
shown in Fig. 2. Note that the order parameters corresponding to the two-
dimensional out-of-plane irreps ($E_{2g}$ and $E_{1u}$) do not have fixed
norms and therefore by themselves cannot be observed in a classical system
Essafi  _et al._ (2017). The labels $E_{a}^{(nm)}$ provide information about
the parity of the magnetic order parameters (label $a$), transformation
properties under spin rotations (subscript $n$), and the coupling between the
spin and spatial transformations (subscript $m$). In the case of the planar
irreps, the $m$ subscript can also be viewed as the “winding” of the spins
around the hexagon.

###  III.3 Intermediate coupling case

Finally, when the SOC is sufficiently strong, the SIA and bond-dependent
anisotropies become important, effectively reducing the symmetry of the spin
Hamiltonian to the magnetic group $\mathrm{D}_{6h}\otimes Z_{2}^{(T)}$. As
noted before, only those anisotropic interactions involving the planar spin
components explicitly break the global axial rotation symmetry, forcing the
spins to align with the local anisotropy axes $\mathbf{\hat{n}}_{i\alpha}$.
Therefore, these interactions are inherently Ising-like. This anisotropy
splits the $E_{g}^{(14)}$ and $E_{u}^{(11)}$ configurations each into two
singlets, while leaving the degeneracy of the remaining irreps unchanged (Fig.
2).

##  IV Self-duality transformations

In order to systematically describe the classical ground state properties of a
magnetic model, one must address the problem of the dimensionality of the
parameter space. In the general case considered in this work, there are more
than a dozen independent parameters, making the complete computational
analysis of the phase diagram forbiddingly expensive. As a result, it is
necessary to determine the means of reducing the parameter space. The most
common approach is to focus on a particular physical example (e.g. a family of
compounds) where the ranges of the coupling constants are approximately known
from either the experimental data or from ab initio calculations. This
approach allows one to bound the values of the parameters, and potentially
even ignore some of them. Although this method is of extreme utility for
explaining the properties of the known compounds, it may provide very limited
information for describing the functionality of novel compounds, since the
relevant parts of the parameter space may fall far beyond the explored
subspace. Moreover, even after the relevant parameter ranges are identified,
the dimensionality of the search space may be quite large and still require
extensive computations. For example, in the weak SOC limit, the magnetic model
considered in this work still contains 4-5 independent interactions.

Another approach, often neglected in the literature, is to determine hidden
relationships between the models with different coupling constants. It is
often true that the parameters in a given model are not completely
independent, and one can determine a set of transformations that map the
Hamiltonian onto itself, while changing the values of the coupling constants.
Such transformations are referred to as the self-duality transformations and
are the subject of this section. We note that while the two methods described
here are different in nature, one can and should use them together to achieve
a systematic yet physically relevant description of a spin model. In this
work, we constrain the values of the coupling constants, in particular those
originating from the SOC, by referring to the experimental and numerical
results for the $\mathrm{Mn}_{3}X$ compounds Soh  _et al._ (2020); Chen  _et
al._ (2020); Park  _et al._ (2018).

Self-duality transformations have played an important role in statistical
physics, an important example being a Kramers-Wannier duality that relates the
ordered and paramagnetic phases in the two-dimensional Ising model on a square
lattice Kramers and Wannier (1941); Onsager (1944); Wannier (1945); Wegner
(1971); Wu and Wang (1976); Savit (1980). Self-duality maps provide a natural
formulation of a renormalization flow and have therefore been used in studies
of critical phenomena. More recently, a different class of self-dual
transformations has been derived for Heisenberg-Kitaev models on honeycomb and
triangular lattices Chaloupka and Khaliullin (2015); Kimchi and Vishwanath
(2014); Maksimov  _et al._ (2019). These transformations have been referred to
in the literature as the Klein duality, since they form a group isomorphic to
the Klein group. The main interest in the self-duality maps has been the
search of accidental degeneracy points, where strongly anisotropic systems at
times display full rotational symmetry Chaloupka and Khaliullin (2015);
Maksimov  _et al._ (2019). Self-duality in two-dimensional kagome layers has
been considered before in Essafi  _et al._ (2016, 2017). There, it was used to
draw connections between different models that support spin liquid phases.

In this section, we determine the relevant self-duality transformations for
Hamiltonian (1). Since the number of possible self-dualities depends on the
types of magnetic interactions in the model, we separate the discussion into
three parts corresponding to the three types of spin Hamiltonians discussed in
Sec. III (see Fig. 2).

###  IV.1  Self-duality as a permutation of spin invariants

For the purposes of this work, we define a self-duality transformation as a
simultaneous transformation of the spin variables and model parameters that
leaves the Hamiltonian unchanged:

$\begin{split}\mu:\\{\mathbf{S}_{i}(\mathbf{r});\bm{\mathcal{A}}_{ij}(\bm{\delta})\\}\longrightarrow\\{\widetilde{\mathbf{S}}_{i};\bm{\widetilde{\mathcal{A}}}_{ij}(\bm{\delta})\\},\\\
\mathcal{H}\big{(}\\{\mathbf{S}_{i};\bm{\mathcal{A}}_{ij}(\bm{\delta})\\}\big{)}=\mathcal{H}\big{(}\\{\widetilde{\mathbf{S}}_{i};\bm{\widetilde{\mathcal{A}}}_{ij}(\bm{\delta})\\}\big{)}.\end{split}$
(17)

The transformation in Eq. (17) is implied to include all spins in the system,
as well as all symmetry-allowed coupling constants. We assume that the spin
transformations can be expressed as a site-dependent linear operation
$\widetilde{\mathbf{S}}_{i}(\mathbf{r})=\mathbf{M}_{i}^{T}(\mathbf{r})\mathbf{S}_{i}(\mathbf{r})$,
where $\mathbf{M}_{i}(\mathbf{r})$ is an orthogonal matrix. Then,

$\bm{\widetilde{\mathcal{A}}}_{ij}(\mathbf{r}-\mathbf{r}^{\prime})=\mathbf{M}_{i}^{T}(\mathbf{r})\bm{\mathcal{A}}_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\mathbf{M}_{j}(\mathbf{r}^{\prime}).$
(18)

We require that the matrices $\bm{\mathcal{A}}_{ij}(\bm{\delta})$ and
$\bm{\widetilde{\mathcal{A}}}_{ij}(\bm{\delta})$ only differ in the values of
the coupling parameters. Note that matrices $\mathbf{M}_{i}(\mathbf{r})$ are
not necessarily unique: we can define
$\widetilde{\mathbf{M}}_{i}(\mathbf{r})=\mathbf{R}\mathbf{M}_{i}(\mathbf{r})$
where $\mathbf{R}$ is a symmetry operation that leaves
$\bm{\mathcal{A}}_{ij}(\bm{\delta})$ invariant.

To demonstrate the relationship between the symmetry of the system and the
number of self-duality transformations, we first note that any quadratic spin
Hamiltonian of the form (3) can be written as a sum of bilinear spin
invariants:

$\mathcal{H}=\sum_{\lambda}\mathcal{A}^{(\lambda)}B^{(\lambda)},$ (19)

where $B^{(\lambda)}$ are the invariants, and $\mathcal{A}^{(\lambda)}$ are
the corresponding coupling constants. One can show [Supplemental Material]
that the bilinear spin invariants can be calculated by squaring the symmetry-
adapted order parameters corresponding to the irreps in the decomposition of
the spin structure (see Sec. III). The result can be written as

$B^{(\lambda)}=B_{kl}^{(\Gamma)}=\mathbf{S}_{k}^{(\Gamma)}\cdot\mathbf{S}_{l}^{(\Gamma)},$
(20)

where $\mathbf{S}_{k}^{(\Gamma)}$ is the symmetry-adapted order parameter
corresponding to irrep $\Gamma$, and $k,l$ label different order parameters
belonging to $\Gamma$. The number of components in this vector is equal to the
dimensionality of the $\Gamma$.

It can be shown that a transformation of the spins that results in a
permutation of the invariants,

$\mu\left(B^{(\lambda)}\right)=\sum_{\lambda}P_{\lambda\mu(\lambda)}B^{(\lambda)}=B^{(\mu(\lambda))},$
(21)

satisfies our definition of self-duality since

$\mu\left(\mathcal{H}\right)=\sum_{\lambda}\mathcal{A}^{(\lambda)}B^{(\mu(\lambda))}=\sum_{\lambda}\widetilde{\mathcal{A}}^{(\lambda)}B^{(\lambda)}.$
(22)

Here, $P_{\lambda\mu(\lambda)}$ permutes indices $\lambda$ and $\mu(\lambda)$,
and

$\widetilde{\mathcal{A}}^{(\lambda)}=\mu^{-1}\left(\mathcal{A}^{(\lambda)}\right).$
(23)

At the same time, a permutation of invariants occurs when we permute the order
parameters, and possibly change the sign of some of the components:

$\mu:S_{k\alpha}^{(\Gamma)}\longrightarrow\pm
S_{k^{\prime}\alpha}^{(\Gamma^{\prime})},$ (24)

where $\alpha$ labels the components of the order parameters, Therefore, we
look for the spin transformations $M_{i}(\mathbf{r})$ that correspond to
permutations between order parameters. Although this is by no means a rigorous
derivation, the self-dualities determined this way are sufficient to describe
most interesting properties observed in this paper 111We note that up to now
the discussion has been completely general and is applicable to any spin
system described by a set of symmetries. The derivation of the general
conditions for self-duality go beyond the scope of this work and are a subject
of an ongoing investigation.. In the remainder of this section we present the
relevant duality transformations in the context of the three SOC cases
discussed in the last section.

As a final note, we point out that the set of all self-duality maps resulting
from the permutations of the symmetry-adapted order parameters forms a group.
This can be deduced from the fact that we only allow orthogonal matrices as
the local transformations of spins. This fact significantly simplifies our
search, since the new transformations can be obtained by combining together
those already identified in the analysis.

###  IV.2 Decoupled case

In the decoupled limit, the coupling matrices are diagonal in the spin
components:

$\bm{\mathcal{A}}_{ij}(\bm{\delta})=J_{ij}(\bm{\delta})\bm{\mathbb{I}},$ (25)

where $\bm{\mathbb{I}}$ is an identity matrix. As a result, there is a single
self-duality transformation, which corresponds to transformations between
order parameters symmetric and antisymmetric under inversion operation:

$\gamma^{(-1)}:\begin{cases}T_{g}\longleftrightarrow T_{u},\\\
Q_{g}\longleftrightarrow Q_{u}.\end{cases}$ (26)

The corresponding group of spin transformations can be written as:

$\mathbf{M}_{i}^{(g)}=\begin{cases}g\bm{\mathbb{I}}&\text{if
}i\in\\{1,2,3\\},\\\ \phantom{-}\bm{\mathbb{I}}&\text{if
}i\in\\{4,5,6\\},\end{cases}$ (27)

where $g=\pm 1$. $\mathbf{M}_{i}^{(+1)}(\mathbf{r})$ is the identity, and
$\mathbf{M}_{i}^{(-1)}(\mathbf{r})$ flips all spins in layer A, while keeping
the layer B unchanged. While this transformation does not affect the in-plane
interactions, the out-of-plane couplings change sign

$\gamma^{(-1)}:\begin{cases}J_{1}\longrightarrow-J_{1},\\\
J_{2}\longrightarrow\phantom{-}J_{2},\\\ J_{4}\longrightarrow-
J_{4}.\end{cases}$ (28)

Thus, $\gamma^{(-1)}$ provides a map between models with ferro- and
antiferromagnetic out-of-plane interactions.

###  IV.3 Weak coupling case

Next, consider a model with weak SOC. The coupling matrix between two spins
(4) can then be simply written as

$\bm{\mathcal{A}}_{ij}(\bm{\delta})=\begin{bmatrix}J_{ij}(\bm{\delta})&D_{ij}(\bm{\delta})&0\\\
-D_{ij}(\bm{\delta})&J_{ij}(\bm{\delta})&0\\\
0&0&J_{ij}(\bm{\delta})+A_{ij}^{(z)}(\bm{\delta})\end{bmatrix}.$ (29)

We must formally include $A_{ij}^{(z)}(\bm{\delta})$ anisotropic interactions
in some cases in order to define proper self-duality transformations. However,
as will become clear in the next sections, this formality is not too
significant in practice, since the majority of the observed phases in this
work are restricted to the $xy$-plane. Therefore, many useful properties of
self-duality relating to the spin structures and phase diagrams remain exact
even if we ignore $A_{ij}^{(z)}(\bm{\delta})$, since the $z$-components of the
spins are decoupled from the planar components (Fig. 2), and so the duality
transformations are also decoupled.

We first focus on the permutations of the four out-of-plane order parameters.
There is a single non-trivial duality transformation corresponding to flipping
the $z$-components in the A layer. The group of transformations is then
written as

$\mathbf{M}_{i}^{(\eta)}=\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&\eta_{i}\end{bmatrix},$ (30)

where

$\eta_{i}=\begin{cases}\eta&\text{if }i\in\\{1,2,3\\},\\\ 1&\text{if
}i\in\\{4,5,6\\},\end{cases}$ (31)

and $\eta=\pm 1$. The parameters are mapped according to

$\mu^{(-1)}:\begin{cases}J_{i}&\longrightarrow\phantom{-}J_{i},\\\
A_{1}^{(z)}&\longrightarrow-A_{1}^{(z)}-2J_{1},\\\
A_{2}^{(z)}&\longrightarrow\phantom{-}A_{2}^{(z)},\\\
A_{4}^{(z)}&\longrightarrow-A_{4}^{(z)}-2J_{4}.\end{cases}$ (32)

Next, we identify the transformations that permute the six planar order
parameters. As discussed in the previous section, these are labeled by an
integer $m$ which represents the six-fold winding number. Therefore, the first
natural choice for a spin transformation is a local rotation written as

$\mathbf{M}^{(m)}_{i}=\begin{bmatrix}\cos\theta_{i}^{(m)}&-\sin\theta_{i}^{(m)}&0\\\
\sin\theta_{i}^{(m)}&\phantom{-}\cos\theta_{i}^{(m)}&0\\\ 0&0&1\end{bmatrix},$
(33)

where

$\theta_{i}^{(m)}=\frac{\pi ml_{i}}{3},$ (34)

$m\in\mathbb{Z}_{6}$ and $l_{i}$ label the positions of atoms on the hexagon
of the unit cell in the counter clock-wise direction:

$\\{l_{1},l_{2},l_{3},l_{4},l_{5},l_{6}\\}=\\{3,1,5,0,4,2\\}.$ (35)

However, there is another group of global transformations that leads to a
distinct duality:

$\mathbf{M}^{(\varepsilon)}=\begin{bmatrix}\varepsilon&0&0\\\ 0&1&0\\\
0&0&1\end{bmatrix},$ (36)

with $\varepsilon=\pm 1$. Note that the elements of $\mathbf{M}_{i}^{(m)}$ and
$\mathbf{M}^{(\varepsilon)}$ in general do not commute. We use elements
$\mathbf{M}^{(\varepsilon)}\mathbf{M}_{i}^{(m)}$ to define the self-duality
transformations $\mu_{m}^{(\varepsilon)}$, which transform the coupling
constants according to

$\mu_{m}^{(\varepsilon)}:\begin{dcases}J_{1}\longrightarrow
J_{1}\cos{\left(\frac{\pi m}{3}\right)}-D_{1}\varepsilon\sin{\left(\frac{\pi
m}{3}\right)},\\\ D_{1}\longrightarrow J_{1}\sin{\left(\frac{\pi
m}{3}\right)}+D_{1}\varepsilon\cos{\left(\frac{\pi m}{3}\right)},\\\
J_{2}\longrightarrow J_{2}\cos{\left(\frac{2\pi
m}{3}\right)}-D_{2}\varepsilon\sin{\left(\frac{2\pi m}{3}\right)},\\\
D_{2}\longrightarrow J_{2}\sin{\left(\frac{2\pi
m}{3}\right)}+D_{2}\varepsilon\cos{\left(\frac{2\pi m}{3}\right)},\\\
J_{4}\longrightarrow J_{4}(-1)^{m}.\end{dcases}$ (37)

In addition, if $\widetilde{J}_{i}$ is the new exchange coupling constant,
then

$\mu_{m}^{(\varepsilon)}:A_{i}^{(z)}\longrightarrow
A_{i}^{(z)}+J_{i}-\widetilde{J}_{i}.$ (38)

One can then prove the following relations:

$\displaystyle\mu_{0}^{(-1)}\mu_{m}^{(+1)}\mu_{0}^{(-1)}=\mu_{-m}^{(+1)},$
(39) $\displaystyle\mu_{m}^{(+1)}\mu_{n}^{(+1)}=\mu_{m+n}^{(+1)},$ (40)
$\displaystyle\mu_{m}^{(-1)}\mu_{m}^{(-1)}=\mu_{0}^{(+1)},$ (41)

where $m$ is implied to be a cyclic integer variable. These relations define
the structure of a dihedral group $\mathrm{D}_{6}$ with $\mu_{0}^{(+1)}$
serving as the identity element. The $\mu_{m}^{(+1)}$ transformations in this
case are equivalent to rotations and $\mu_{0}^{(-1)}$ is the reflection
operation. Since $\mu^{\eta}$ and $\mu_{m}^{(\varepsilon)}$ commute, the most
general duality transformation is written as
$\mu^{(\eta)}\mu_{m}^{(\varepsilon)}$ leading to group structure equivalent to
$\mathrm{D}_{6h}$. Thus, we obtain a very interesting situation where the
group of self-dualities is non-Abelian and is isomorphic to the point group of
the underlying lattice.

###  IV.4 Intermediate coupling case

When the SOC is sufficiently strong, the coupling matrix
$\bm{\mathcal{A}}_{ij}(\bm{\delta})$ takes the form of Eq. (4) and (5). The
symmetries of the Hamiltonian are reduced down to the paramagnetic group as
seen in Fig. 2. This has no affect on the out-of-plane order parameters, and
$\mu^{(\eta)}$ are still valid self-duality maps. However, the two-dimensional
$E_{g}^{(14)}$ and $E_{u}^{(11)}$ planar irreps each split into two one-
dimensional irreps ($\\{B_{1g},B_{2g}\\}$ and $\\{A_{1u},A_{2u}\\}$
respectively). As a result, the dualities obtained from permutations of these
irreps with other two-dimensional irreps are no longer exact. Out of the
twelve elements of $\mu_{m}^{(\varepsilon)}$, only four remain:
$\mu_{0}^{(+1)}$, $\mu_{1}^{(-1)}$, $\mu_{3}^{(+1)}$, and $\mu_{4}^{(-1)}$. It
is straight-forward to check that $\mu_{0}^{(+1)}$ and $\mu_{4}^{(-1)}$ leave
the anisotropic parameters unchanged,

$\mu_{0}^{(+1)},\mu_{4}^{(-1)}:\begin{cases}K_{\alpha}\phantom{-}\longrightarrow
K_{\alpha},\\\ A_{1}^{(xy)}\longrightarrow A_{1}^{(xy)},\\\
A_{2}^{(xy)}\longrightarrow A_{2}^{(xy)},\\\ A_{4}^{(xy)}\longrightarrow
A_{4}^{(xy)},\end{cases}$ (42)

while $\mu_{3}^{(+1)}$ and $\mu_{1}^{(-1)}$ change the sign of the out-of-
plane bond-dependent interactions,

$\mu_{3}^{(+1)},\mu_{1}^{(-1)}:\begin{cases}K_{\alpha}\phantom{-}\longrightarrow\phantom{-}K_{\alpha},\\\
A_{1}^{(xy)}\longrightarrow-A_{1}^{(xy)},\\\
A_{2}^{(xy)}\longrightarrow\phantom{-}A_{2}^{(xy)},\\\
A_{4}^{(xy)}\longrightarrow-A_{4}^{(xy)}.\end{cases}$ (43)

In addition to these transformations, global $C_{4}$ rotations around the
$z$-axis now lead to a valid self-duality map $\zeta$, since they interchange
the one-dimensional planar irreps, thus flipping the sign of the anisotropic
interactions:

$\zeta:\begin{cases}K_{\alpha}\phantom{-}\longrightarrow-K_{\alpha},\\\
A_{1}^{(xy)}\longrightarrow-A_{1}^{(xy)},\\\ A_{2}^{(xy)}\longrightarrow-
A_{2}^{(xy)},\\\ A_{4}^{(xy)}\longrightarrow-A_{4}^{(xy)}.\end{cases}$ (44)

The elements
$\\{\mu_{0}^{(+1)},\mu_{3}^{(+1)},\mu_{1}^{(-1)},\mu_{4}^{(-1)}\\}$ and
$\zeta$ commute (even though the corresponding spin transformation matrices
generally do not), and all together they form the Burnside group $B(3,2)$
$\mathbb{Z}_{2}\otimes\mathbb{Z}_{2}\otimes\mathbb{Z}_{2}$, which is Abelian.
Combining it with the $\mu^{(\eta)}$, the group of the self-duality
transformations in the intermediate SOC limit becomes the Burnside group
$B(4,2)$.

###  IV.5 Consequences of self-duality

The consequences of self-duality are far-reaching. In general, self-dual
transformations can be thought of as the symmetries of the parameter space
that do not change the energy of the system. As a result, the dual points in
the parameter space (self-dual images) will have a lot of the same physical
properties. These include the ground state energy and most thermodynamic
properties, since the self-duality applies also to the partition function.
These results are also not limited to classical systems either: the self-dual
images must have mostly the same quantum properties, within redefinition of
quantization axes. Therefore, given a description of a single phase, one can
immediately describe all phases that relate to it via self-duality
transformations. In the following, we discuss a few important properties of
self-duality relevant to this work.

From the properties of the duality relations it follows that in order to prove
that a certain spin structure is stable for some choice of physical
parameters, it is sufficient to show that one of its images is stable
somewhere in the parameter space. Given that even in the simple case of weak
SOC limit the parameter space is four-dimensional, the self-duality maps can
significantly reduce the time of computations since for every phase boundary
$f(\bm{\mathcal{A}})$, an image
$f\left(\mu\left(\bm{\mathcal{A}}\right)\right)$ must also correspond to a
phase boundary. At the phase boundaries, self-dualities become symmetry
operations, and the Hamiltonian remains invariant after the transformation.

$m$ | $\widetilde{J}_{1}$ | $\widetilde{D}_{1}$ | $\widetilde{A}^{(z)}_{1}$ | $\widetilde{J}_{2}$ | $\widetilde{D}_{2}$ | $\widetilde{A}^{(z)}_{2}$ | $\widetilde{J}_{4}$ | $\widetilde{A}^{(z)}_{4}$
---|---|---|---|---|---|---|---|---
0 | $J_{1}$ | 0 | 0 | $J_{2}$ | 0 | 0 | $J_{4}$ | 0
1 | $\frac{1}{2}J_{1}$ | $\frac{\sqrt{3}}{2}J_{1}$ | $-\frac{1}{2}J_{1}$ | $-\frac{1}{2}J_{2}$ | $\frac{\sqrt{3}}{2}J_{2}$ | $-\frac{3}{2}J_{2}$ | $-J_{4}$ | $-2J_{4}$
2 | $-\frac{1}{2}J_{1}$ | $\frac{\sqrt{3}}{2}J_{1}$ | $-\frac{3}{2}J_{1}$ | $-\frac{1}{2}J_{2}$ | $-\frac{\sqrt{3}}{2}J_{2}$ | $-\frac{3}{2}J_{2}$ | $J_{4}$ | 0
3 | $-J_{1}$ | 0 | $-2J_{1}$ | $J_{2}$ | 0 | 0 | $-J_{4}$ | $-2J_{2}$
4 | $-\frac{1}{2}J_{1}$ | $-\frac{\sqrt{3}}{2}J_{1}$ | $-\frac{3}{2}J_{1}$ | $-\frac{1}{2}J_{2}$ | $\frac{\sqrt{3}}{2}J_{2}$ | $-\frac{3}{2}J_{2}$ | $J_{4}$ | 0
5 | $\frac{1}{2}J_{1}$ | $-\frac{\sqrt{3}}{2}J_{1}$ | $-\frac{1}{2}J_{1}$ | $-\frac{1}{2}J_{2}$ | $-\frac{\sqrt{3}}{2}J_{2}$ | $-\frac{3}{2}J_{2}$ | $-J_{4}$ | $-2J_{4}$

Table 2: Points of accidental $\mathcal{G}_{J}$ symmetry obtained via self
duality maps $\mu_{m}^{(\varepsilon)}$. Note that the value of $\varepsilon$
does not change the number of points.

Self-dualities also allow one to quickly find the points of accidental
degeneracy. Consider, for example, Eq. (37) where the original DM parameters
are set to zero $D_{1}=D_{2}=0$. We get a set of new parameters, where
$\widetilde{D}_{1}$ and $\widetilde{D}_{2}$ are generally non-zero. However,
since the self-duality preserves the transformation properties and therefore
the symmetry group of the Hamiltonian, all self-dual images in this case must
be completely isotropic in spin space (i.e. have spin group
$\mathcal{G}_{J}$). Note that as stated before, in addition to the
$\widetilde{J}_{k}$ and $\widetilde{D}_{k}$ one must also specify
$\widetilde{A}^{(z)}_{k}=\widetilde{J}_{k}-J_{k}$ to ensure that the full
rotational symmetry is exact. With this, we get a set of points with
accidental $\mathcal{G}_{J}$ degeneracy, listed in table 2. At these points we
expect the spin-wave dispersion to obtain additional zero-energy pseudo-
Goldstone modes.

It is important to note, however, that the properties which depend directly on
the spin structure (e.g. average magnetization) are not always the same, since
the self-duality maps generally do not conserve the order parameter. An
important consequence of this is that the self-dual images will have the same
spin wave excitation spectra only when the corresponding spin transformations
produce proper local axes (i.e. when the local transformations are rotations).
This will be discussed in more detail in Sec. VII.

##  V Classical Phase Diagrams for Decoupled and Weakly Coupled Cases

###  V.1 Brief overview of the results

In this section we present the classical magnetic phases for models with
either negligible or weak SOC. The effects of intermediate SOC are discussed
in Sec. VIII. Before starting a detailed discussion of magnetic ground states,
we summarize the main results of this section. Based on our findings, we group
the magnetic structures into three categories: single-$\mathbf{Q}$
configurations, multi-$\mathbf{Q}$ structures, and incommensurate phases with
Ising-like ordering. We label these as $\bm{\Phi}_{m}^{(\mathbf{Q})}$,
$\bm{\Psi}_{m}^{(\mathbf{M})}$, and $\bm{\Lambda}_{m}$ respectively. The
integer label $m$ here refers to the different self-dual images produced by
the $\mu_{m}^{(\varepsilon)}$ elements in the weak SOC limit, as discussed in
the previous section. In addition to these states, in the decoupled limit,
some phase boundaries correspond to structures with degenerate wavevectors in
the Brillouin zone. The magnetic ordering wavevectors of all observed states
are shown in Fig. 4.

Figure 3: Ground state phase diagram for exchange-coupled spins on AB-stacked
kagome lattice. Here, the dashed and dash-dotted boundaries correspond to the
$\mathcal{M}_{\Sigma}$ and $\mathcal{M}_{[\mathbf{q}\mathbf{q}0]}$ manifold
states. Since the exchange interactions do not differentiate between the
different chirality values, the $m=1,5$ ($m=2,4$) states are written as linear
combinations. Figure 4: Magnetic ordering wavevectors corresponding to phases
studied in this paper shown in the first Brillouin zone of an AB-stacked
kagome lattice. $\mathcal{M}_{[\mathbf{q}\mathbf{q}0]}$ and
$\mathcal{M}_{\Sigma}$ represent degenerate ground state manifolds.
$\bm{\Phi}_{m}^{(\mathbf{Q})}$ and $\bm{\Psi}_{m}^{(\mathbf{M})}$ label the
single- and multi-$\mathbf{Q}$ states, and $\bm{\Lambda}_{m}$ labels the
incommensurate phases with Ising-type order, as described in the main text.

LT calculations were found to give correct spin configurations and energies
for the single-$\mathbf{Q}$ phases. For the remaining types of phases, the
analytical results gave non-normalized spin configurations and, as a result,
lower energies compared to the MC simulations. However, in most cases, LT
provided a decent estimation of the locations of the phase boundaries,
allowing us to optimize our numerical boundary search. Therefore, the exact
phase boundaries between $\bm{\Phi}_{m}^{(\mathbf{Q})}$-type states were
calculated using the LT method, and the locations of the remaining boundaries
were determined via numerical MC simulations.

We show that Heisenberg exchange interactions are sufficient to stabilize both
single- and multi-$\mathbf{Q}$ phases. The inclusion of DM interactions lifts
the degeneracy of the non-collinear phases by stabilizing configurations with
particular chirality, which further enriches the magnetic phases. For most
values of exchange constants, the $D_{1}-D_{2}$ phase diagrams display
multiple non-trivial ($\mathbf{Q}\neq 0$) structures, including unconventional
states with Ising-like ordering.

###  V.2 Decoupled case

Figure 5: Ground state phase diagrams on AB-stacked kagome lattice for
selected values of exchange and a range of DM values. In the middle diagram,
the brightened strips between the $\bm{\Lambda}_{2(4)}$ and
$\bm{\Phi}_{2(4)}^{(\bm{\Gamma})}$ phases are coexistence regions with many
stable multi-Q configurations.

First, we present phase diagrams for the model in Eq. (1) with
$\mathcal{H}=\mathcal{H}_{J}$ in Fig. 3. Both diagrams ($J_{2}>0$ and
$J_{2}<0$) display a clear “inversion” symmetry corresponding to flipping the
sign of both $J_{1}$ and $J_{4}$, as a consequence of the self-duality map
$\gamma^{(-1)}$, discussed in Sec. IV.2. The results yield eight distinct
phases, including single-, and multi-$\mathbf{Q}$ structures. When $J_{1}$ and
$J_{4}$ have the same sign and sufficiently large magnitudes, the ground state
spin configurations become collinear, with spins in the A and B layers either
parallel ($J_{1}<0$, $J_{4}<0$) or antiparallel ($J_{1}>0$, $J_{4}>0$) to each
other. On the other hand, opposing signs of the interactions introduce
frustration, which for large magnitudes of the couplings lead to uniform
120-degree states. Notably, the exchange interactions alone do not
differentiate between states with different chirality, leaving the structures
with $m=1,5$ ($m=2,4$) degenerate. Interestingly enough, for sufficiently
large magnitudes of $J_{1}$ there is a range of $J_{4}$ values that stabilizes
commensurate cycloidal configurations regardless of the sign of $J_{2}$. These
structures are characterized by a period three modulation ($\mathbf{Q}=K$)
with collinear spins in each unit cell (Fig. 7). The stability of these phases
extends indefinitely for large magnitudes of $J_{1}$. In the case of
antiferromagnetic nn in-plane interactions, one can also stabilize
multi-$\mathbf{Q}$ configurations for intermediate values of $J_{1}$ and
$J_{4}$. Similar to the uniform 120-degree phases, the
$\bm{\Psi}_{m}^{(\mathbf{M})}$ structures in Fig. 3 correspond to degenerate
configurations with opposite chiralities.

We also note that at certain phase boundaries structures with degenerate
wavevectors can be stabilized. In particular, the boundaries between
$\bm{\Phi}_{0}^{(\bm{\Gamma})}$ ($\bm{\Phi}_{3}^{(\bm{\Gamma})}$) and
$\bm{\Phi}_{0}^{(\mathbf{K})}$ ($\bm{\Phi}_{3}^{(\mathbf{K})}$) states
stabilize structures where all wavevectors with $Q_{z}=0$ are degenerate,
leading to a two-dimensional ground state manifold
$\mathcal{M}_{[\mathbf{q}\mathbf{q}0]}$ (dash-dotted lines in Fig. 3). Also,
the boundaries between $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ and
$\bm{\Psi}_{m}^{(\mathbf{M})}$ are degenerate along the $\Sigma$ lines in the
Brillouin zone and are labeled as $\mathcal{M}_{\Sigma}$ (dashed lines in Fig.
3). In this work we will focus more on the ordered states and leave the
description of the degenerate manifold states outside of the scope of this
work.

###  V.3 Weak coupling case

Next, we consider the effects of DM interactions. To reduce the number of
independent parameters, we focus on the case where $J_{4}=0$, and present the
diagrams that display the majority of the magnetic phases observed in this
study. The remaining phases can be obtained from the self-duality operations
$\mu_{m}^{\varepsilon}$ presented in Sec. IV.3 [Supplemental Material].

The phase diagrams are shown in Fig. 5. We see that DM interactions break the
degeneracy associated with the different chirality values of the non-collinear
spin structures. For large magnitudes of both $D_{1}$ and $D_{2}$, the ground
state eventually becomes one of the $\mathbf{Q}=\Gamma$ 120-degree structures.
For the intermediate values of these constants, the competition between the
exchange and DM interactions introduces a plethora of unconventional magnetic
phases. These include single-$\mathbf{Q}$ phases and multi-$\mathbf{Q}$
configurations already discussed, as well as incommensurate phases that
manifest themselves as Ising-like stripes. The latter are stable for a wide
range of the DM parameters and extend far beyond the ranges shown in Fig. 5.
These states were found to be hard to resolve in the simulations and required
a careful choice of parameters, as well as many MC steps in order to obtain
ordered stripe configurations. Furthermore, for the $J_{2}>0$
$J_{1}/|J_{2}|>\frac{3}{2}$ diagram, we identified thin coexistence regions
between the $\bm{\Lambda}_{2(4)}$ and $\bm{\Phi}_{2(4)}^{(\bm{\Gamma})}$
phases where various multi-$\mathbf{Q}$ configurations are stable.

##  VI Structure of the magnetic phases

###  VI.1 Single-$\mathbf{Q}$ phases

Figure 6: The $\bm{\Phi}^{(\bm{\Gamma})}_{m}$ spin configurations. Blue and
red arrows are used to differentiate the structures that are symmetric and
antisymmetric under the spatial inversion.

In this work, we encounter two types of single-$\mathbf{Q}$ phases: those with
$\mathbf{Q}=\Gamma$ and $\mathbf{Q}=K$, which we label as
$\bm{\Phi}^{(\bm{\Gamma})}_{m}$ and $\bm{\Phi}^{(\mathbf{K})}_{m}$
respectively. The spin configurations of these phases are shown in figures 6
and 7. Note that for exchange-only interactions, the spins generally also
possess out-of-plane components, since the order parameters in this case
belong to one of the four irreps in Sec. III.1.

When a weak SOC is turned on in the system, the in-plane and out-of-plane
components are no longer equivalent and correspond to different irreps (Sec.
III.2). Thus, in the decoupled limit, these structures can be parameterized as
$\mathbf{S}_{i}(\mathbf{r})=R(\phi_{z},\phi)R^{(\mathbf{Q})}(\mathbf{r})M_{i}^{(m)}\mathbf{\hat{S}}$,
where $\mathbf{\hat{S}}$ is an arbitrary in-plane unit vector, $M_{i}^{(m)}$
is defined in Eq. (33), and the remaining two rotation matrices are defined as

$\displaystyle R(\phi_{z},\phi)$
$\displaystyle=\begin{bmatrix}\cos(\phi)&-\sin(\phi)&0\\\
\cos(\phi_{z})\sin(\phi)&\cos(\phi_{z})\cos(\phi)&-\sin(\phi_{z})\\\
\sin(\phi_{z})\sin(\phi)&\sin(\phi_{z})\cos(\phi)&\cos(\phi_{z})\end{bmatrix},$
(45) $\displaystyle R^{(\mathbf{Q})}(\mathbf{r})$
$\displaystyle=\begin{bmatrix}\cos(\mathbf{Q}\cdot\mathbf{r})&-\sin(\mathbf{Q}\cdot\mathbf{r})&0\\\
\sin(\mathbf{Q}\cdot\mathbf{r})&\cos(\mathbf{Q}\cdot\mathbf{r})&0\\\
0&0&1\end{bmatrix},$ (46)

Matrix $R(\phi_{z},\phi)$ is used to define the globally phase of the spin
configurations. In the weakly-coupled limit, the parameterization is changed
to
$\mathbf{S}_{i}(\mathbf{r})=R(0,\phi)R^{(\mathbf{Q})}(\mathbf{r})M_{i}^{(m)}\mathbf{\hat{S}}$.
The corresponding classical energies of these spin configurations are
presented in table 3.

$m$ | $\mathcal{H}^{0}\left(\bm{\Phi}_{m}^{(\bm{\Gamma})}\right)$ | $\mathcal{H}^{0}\left(\bm{\Phi}_{m}^{(\mathbf{K})}\right)$
---|---|---
0 | $2J_{1}+2J_{2}+3J_{4}$ | $2J_{1}+\frac{J_{2}}{2}$
1 | $J_{1}-J_{2}-3J_{4}-\sqrt{3}(D_{1}+D_{2})$ | $J_{1}-\frac{J_{2}}{4}-\sqrt{3}(D_{1}+\frac{D_{2}}{4})$
2 | $-J_{1}-J_{2}+3J_{4}-\sqrt{3}(D_{1}-D_{2})$ | $-J_{1}-\frac{J_{2}}{4}-\sqrt{3}(D_{1}-\frac{D_{2}}{4})$
3 | $-2J_{1}+2J_{2}-3J_{4}$ | $-2J_{1}+\frac{J_{2}}{2}$
4 | $-J_{1}-J_{2}+3J_{4}+\sqrt{3}(D_{1}-D_{2})$ | $-J_{1}-\frac{J_{2}}{4}+\sqrt{3}(D_{1}-\frac{D_{2}}{4})$
5 | $J_{1}-J_{2}-3J_{4}+\sqrt{3}(D_{1}+D_{2})$ | $J_{1}-\frac{J_{2}}{4}+\sqrt{3}(D_{1}+\frac{D_{2}}{4})$
Table 3: Classical energies of the single-$\mathbf{Q}$ configurations. Figure
7: The $\bm{\Phi}^{(\mathbf{K})}_{m}$ spin configurations.

We note that while formally the $z$-components of spins are decoupled from the
planar ones in the weak SOC limit, the energy of the
$\bm{\Phi}_{0}^{(\mathbf{Q})}$ and $\bm{\Phi}_{3}^{(\mathbf{Q})}$ structures
is independent of the DM parameters leading to a full rotational degeneracy at
the classical level. This degeneracy originates from the collinear sublattice
structure in these phases. As will be discussed in Sec. VII, thermal
fluctuations break this classical degeneracy and orient the structures in the
out-of-plane direction.

###  VI.2 Multi-$\mathbf{Q}$ phases

The most prominent type of multi-$\mathbf{Q}$ states in the phase diagrams in
Fig. 3 and 5 is $\bm{\Psi}_{m}^{(\mathbf{M})}$. The structure factor of these
structures displays one peak at the $\Gamma$ point, as well as three peaks at
the $M$ points in the first Brillouin zone (forming the star of $M$), as shown
in Fig. 4. The spin configurations are shown in Fig. 8.

Figure 8: Magnetic unit cells of the $\bm{\Psi}_{m}^{(\mathbf{M})}$
configurations. Colors indicate spins that belong to the same orbit (see main
text).

The parameterization of these structures is not trivial. To make some further
progress, we must consider the symmetry of these configurations in order to
reduce the number of independent variables. The analysis below is performed
assuming the weak coupling case and then generalized to also describe the
decoupled case.

We first note that the peaks of the structure factor in Fig. 4 are completely
symmetric with respect to the $\mathrm{D}_{6h}$ point group operations.
Therefore, even though these spin structures break planar translations by a
single unit cell, the point group lattice transformations around the (0,0)
site are preserved, provided that we combine them with the necessary spin
rotations. Thus, we can still express the $\bm{\Psi}_{m}^{(\mathbf{M})}$
configurations in terms of the two-dimensional $E_{g}^{(1m)}$ and
$E_{u}^{(1m)}$ irreps of $\mathcal{G}_{D}$ in table 1. The general
representation of the spin structures in Fig. 8 decomposes as

$\bm{\Psi}^{(\mathbf{M})}=\bigoplus_{m}4E_{a_{m}}^{(1m)},$ (47)

where $a_{m}=g$ if $m$ is even and $a_{m}=u$ otherwise. Because all of the
symmetry operations transform spins globally (i.e. no site-dependent
transformations), it is impossible to have spin configuration described by
irreps with more than one value of the winding $m$, without inevitably
breaking the $\mathrm{D}_{6h}$ lattice symmetry. Therefore, a given spin
configuration must decompose into four irreps with the same winding number
$m$, $\bm{\Psi}_{m}^{(\mathbf{M})}=4E_{1}^{(1m)}$. As a result, the
parameterization of $\bm{\Psi}_{m}^{(\mathbf{M})}$ can be written in terms of
four planar vectors, which describe the spins belonging to the same orbit.
Here, we define an orbit as a set of all spin components related to each other
via lattice point group transformations (Fig. 8). The four orbits are given
below.

$\displaystyle\mathbf{O}^{(1)}_{i}$ $\displaystyle:\big{\\{}$
$\displaystyle\mathbf{S}_{1}(0,0),\mathbf{S}_{2}(0,0),\mathbf{S}_{3}(0,0),$
$\displaystyle\mathbf{S}_{4}(0,0),\mathbf{S}_{5}(0,0),\mathbf{S}_{6}(0,0)\big{\\}},$
(48) $\displaystyle\mathbf{O}^{(2)}_{i}$ $\displaystyle:\big{\\{}$
$\displaystyle\mathbf{S}_{1}(1,0),\mathbf{S}_{2}(1,1),\mathbf{S}_{3}(0,1),$
$\displaystyle\mathbf{S}_{4}(1,0),\mathbf{S}_{5}(1,1),\mathbf{S}_{6}(0,1)\big{\\}},$
(49) $\displaystyle\mathbf{O}^{(3)}_{i},\mathbf{O}^{(4)}_{i}$
$\displaystyle:\big{\\{}$
$\displaystyle\\{\mathbf{S}_{1}(1,1),\mathbf{S}_{2}(0,1),\mathbf{S}_{3}(1,0),$
$\displaystyle\mathbf{S}_{4}(1,1),\mathbf{S}_{5}(0,1),\mathbf{S}_{6}(1,0)\\},$
$\displaystyle\\{\mathbf{S}_{1}(0,1),\mathbf{S}_{2}(1,0),\mathbf{S}_{3}(1,1),$
$\displaystyle\mathbf{S}_{4}(0,1),\mathbf{S}_{5}(1,0),\mathbf{S}_{6}(0,1)\\}\big{\\}}.$
(50)

The positions of the spins inside of the magnetic unit cell are given by two
integers $(k,l)$, which correspond to displacement vector
$\mathbf{r}=k\mathbf{a}_{1}+l\mathbf{a}_{2}$. The last two orbits include
different components of spins on the same sites. Thus, the
$\bm{\Psi}_{m}^{(\mathbf{M})}$ order parameters can generally be written in
terms of the four vectors $\mathbf{O}_{i}^{(k)}$. However, due to the
constraint on the norm of the spin vectors, the orbits $\mathbf{O}^{(3)}_{i}$
and $\mathbf{O}^{(4)}_{i}$ are no longer independent and will share the
parameterization parameters. We can therefore write

$\bm{\Psi}_{m}^{(\mathbf{M})}=\min_{\phi_{k}}\sum_{k=1}^{4}R(0,\phi_{k})M_{i}^{(m)}\mathbf{O}_{i}^{(k)},$
(51)

where $M_{i}^{(m)}$ is given in Eq. (33), $R(\phi_{z},\phi)$ is defined in Eq.
(45), and $\phi_{k}$ are the angles that parameterize the spin structure. The
exact spin configurations are then obtained by minimizing the energy with
respect to the $\phi_{k}$. To simplify the minimization, we can remove the in-
plane rotations that do not change the energy by fixing the value of
$\phi_{1}$, and also use the equivalence of $\mathbf{O}_{i}^{(3)}$ and
$\mathbf{O}_{i}^{(4)}$ to deduce that $\phi_{3}=-\phi_{4}$. This leaves us
with two independent variables which can be determined using numerical
minimization.

Finally, we note that in the decoupled limit the above analysis still holds,
although the spin configurations may also possess out-of-plane components,
which can be parameterized by introducing a non-zero angle $\phi_{z}$.

###  VI.3  Ising-like phases

Figure 9: (a) A snapshot of the low-temperature $\bm{\Lambda}_{4}$ spin
configuration on a lattice with $18^{3}\times 6$ spins. This state was
obtained using $J_{1}=1$, $J_{2}=1$, $D_{1}=-1$, and $D_{2}=-0.5$. The order
parameters that define the Ising variables are shown in the top right corner.
(b) Structure factor calculated for a system with $30\times 30\times 2\times
6$ spins. The values were normalized with respect to the maxima. The
parameters used in the simulation are the same as in (a). (c) Fourier
transform of the effective Ising Hamiltonian in Eq. (60) obtained using the LT
approximation.

Lastly, we consider the $\bm{\Lambda}_{m}$ phases. The complete description of
these structures is complicated by the fact that there appears to be a large
number of local minima that are very close in energy to the ground state
configuration. Furthermore, the spin configurations depend strongly on the
size of the system, indicating an incommensurate nature of the magnetic order.
As a result, in order to resolve a single structure, it is necessary to
perform long ($\sim 10^{5}$ MC steps) simulations on large ($>3\cdot 10^{4}$
spins) systems. The resulting $\bm{\Lambda}_{m}$ spin configurations often
manifest in non-trivial patterns, where in each unit cell, the order parameter
is approximately equal to $\pm\bm{\Phi}_{m}^{(\bm{\Gamma})}$. The sign of this
effective order parameter alternates rapidly throughout the system, giving
rise to discrete domains. For this reason, we refer to these states as Ising-
like.

In the following, we will focus our attention on the description of this Ising
behavior by considering the $\bm{\Lambda}_{4}$ phase, although the discussion
equally applies to all $\bm{\Lambda}_{m}$ states by the virtue of the self-
duality transformations (Sec. IV). Since these phases do not appear in the
decoupled limit, the discussion below assumes weak SOC. Fig. 9 (a) shows a
typical low-temperature $\bm{\Lambda}_{4}$ spin configuration. It is important
to note that each kagome layer is identical, corresponding to $Q_{z}=0$, but
within the plane, the pattern can often appear seemingly random. The structure
factor (Fig. 9 (b)) displays a characteristic ring in the Brillouin zone with
a relatively well-defined radius. We note in this ring there are typically a
couple of peaks with larger intensity, corresponding to the dominant
orientations of the stripy domains. In our simulations, the spin deviations
from the $\bm{\Phi}_{4}^{(\bm{\Gamma})}$ state in each unit cell were found to
be restricted to the plane of the kagome (the out-of-plane deviations are
typically four orders of magnitude smaller). Therefore, for simplicity, we
express the spins as two-dimensional vectors

$\mathbf{S}_{i}(\mathbf{r})=\begin{bmatrix}\cos\left(\theta_{i}^{(4)}+\delta\tilde{\theta}_{i}(\mathbf{r})\right)\\\
\sin\left(\theta_{i}^{(4)}+\delta\tilde{\theta}_{i}(\mathbf{r})\right)\end{bmatrix},$
(52)

where $\theta_{i}^{(m)}$ are the angles in Eq. (34) that correspond to the
$\bm{\Phi}_{4}^{(\bm{\Gamma})}$ order parameter, and
$\delta\tilde{\theta}_{i}(\mathbf{r})$ are the deviations from the order
parameter. In order to account for the alternating sign of the
$\bm{\Phi}_{4}^{(\bm{\Gamma})}$ order parameter, we further rewrite these
deviations as

$\delta\tilde{\theta}_{i}(\mathbf{r})=-\frac{v(\mathbf{r})-1}{2}\pi+\delta\theta_{i}(\mathbf{r}),$
(53)

where $v(\mathbf{r})$ are Ising variables that take values $\pm 1$, and
$0\leq\delta\theta_{i}(\mathbf{r})<\pi$ are the residual angle deviations.
Substituting this parameterization into the Hamiltonian with exchange and DM
interactions is equivalent to performing a $\mu_{4}^{(+1)}$ gauge
transformation, from which we obtain

$\displaystyle\mathcal{H}=\sum_{\mathbf{r}\mathbf{r}^{\prime}}\Bigg{[}\sum_{\langle
ij\rangle}$
$\displaystyle\widetilde{J}_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\cos\Big{(}\delta\theta_{ij}(\mathbf{r},\mathbf{r}^{\prime})\Big{)}$
(54) $\displaystyle+$
$\displaystyle\widetilde{D}_{ij}(\mathbf{r}-\mathbf{r}^{\prime})\sin\Big{(}\delta\theta_{ij}(\mathbf{r},\mathbf{r}^{\prime})\Big{)}\Bigg{]}v(\mathbf{r})v(\mathbf{r}^{\prime}),$

where $\widetilde{J}_{ij}(\bm{\delta})$ and $\widetilde{D}_{ij}(\bm{\delta})$
are given in Eq. (37) ($m=4$ and $\varepsilon=+1$), and
$\delta\theta_{ij}(\mathbf{r},\mathbf{r}^{\prime})=\delta\theta_{i}(\mathbf{r})-\delta\theta_{j}(\mathbf{r}^{\prime})$.
Eq. (54) can be recast into an Ising model on a triangular lattice with non-
local nn interactions:

$\mathcal{H}=\sum_{\mathbf{r}\mathbf{r}^{\prime}}\mathcal{J}(\mathbf{r},\mathbf{r}^{\prime})v(\mathbf{r})v(\mathbf{r}^{\prime})+\mathcal{B}(\mathbf{r},\mathbf{r}^{\prime}).$
(55)

Here, both the Ising exchange coupling
$\mathcal{J}(\mathbf{r},\mathbf{r}^{\prime})$, and constant
$\mathcal{B}(\mathbf{r},\mathbf{r}^{\prime})$ depend on
$\delta\theta_{i}(\mathbf{r})$.

It is worth discussing the values of $\widetilde{J}_{ij}(\bm{\delta})$ and
$\widetilde{D}_{ij}(\bm{\delta})$ for which the $\bm{\Lambda}_{4}$ is stable.
In fig 9, the parameters are $J_{1}=J_{2}=1$, $D_{1}=-1$, and $D_{2}=-0.5$,
which corresponds to $\widetilde{J}_{1}\approx-1.37$,
$\widetilde{J}_{2}\approx-0.07$, $\widetilde{D}_{1}\approx-0.37$, and
$\widetilde{D}_{2}\approx 1.12$. Therefore, in the local coordinate frame, the
inter-plane interactions are dominated by the ferromagnetic exchange, whereas
the in-plane interactions almost exclusively come from the DM term. Notably,
$\bm{\Lambda}_{4}$ phase in this case (Fig. 5) appears to be stable for
arbitrarily large values of $D_{1}$, which leads to large negative values of
$\widetilde{J}_{1}$.

Using this information, we can expand Eq. (54) to second order in
$\delta\theta_{i}(\mathbf{r})$, while assuming that $\widetilde{J}_{2}$ is of
the same order of magnitude as the deviations. We also impose ferromagnetic
order along the $z$-direction by setting
$\delta\theta_{i}(\mathbf{r}\pm\mathbf{a}_{3})=\delta\theta_{i}(\mathbf{r})$
and $v(\mathbf{r}\pm\mathbf{a}_{3})=v(\mathbf{r})$. Under these assumptions,
Eq. (54) simplifies to the following quadratic equation:

$\mathcal{H}\approx\mathcal{H}^{(0)}+\sum_{\mathbf{r}}\sum_{i}g_{i}(\mathbf{r})\delta\theta_{i}(\mathbf{r})+\frac{1}{2}\sum_{\mathbf{r}}\sum_{ij}H_{ij}\delta\theta_{i}(\mathbf{r})\delta\theta_{j}(\mathbf{r}),$
(56)

where $\mathcal{H}^{(0)}$ is independent of $\delta\theta_{i}(\mathbf{r})$

$\displaystyle\mathcal{H}^{(0)}=\sum_{\mathbf{r}}6\widetilde{J}_{1}\Big{[}$
$\displaystyle 1+v(\mathbf{r})v(\mathbf{r}+\mathbf{a}_{3})\Big{]}$
$\displaystyle+2\widetilde{J}_{2}\Big{[}$ $\displaystyle
3+v(\mathbf{r})v(\mathbf{r}+\mathbf{a}_{1})+v(\mathbf{r})v(\mathbf{r}+\mathbf{a}_{2})+v(\mathbf{r})v(\mathbf{r}+\mathbf{a}_{1}+\mathbf{a}_{2})\Big{]},$
(57)

and $g_{i}(\mathbf{r})$ and $H_{ij}$ are the gradient vector and Hessian
matrix respectively:

$g_{i}(\mathbf{r})=-\widetilde{D}_{2}v(\mathbf{r})\begin{bmatrix}v(\mathbf{r}-\mathbf{a}_{2})-v(\mathbf{r}-\mathbf{a}_{1}-\mathbf{a}_{2})\\\
v(\mathbf{r}-\mathbf{a}_{1})-v(\mathbf{r}+\mathbf{a}_{2})\\\
v(\mathbf{r}+\mathbf{a}_{1}+\mathbf{a}_{2})-v(\mathbf{r}+\mathbf{a}_{1})\\\
v(\mathbf{r}+\mathbf{a}_{2})-v(\mathbf{r}+\mathbf{a}_{1}+\mathbf{a}_{2})\\\
v(\mathbf{r}+\mathbf{a}_{1})-v(\mathbf{r}-\mathbf{a}_{2})\\\
v(\mathbf{r}-\mathbf{a}_{1}-\mathbf{a}_{2})-v(\mathbf{r}-\mathbf{a}_{1})\end{bmatrix},H_{ij}=2\widetilde{J}_{1}\begin{bmatrix}-2&0&0&0&1&1\\\
0&-2&0&1&0&1\\\ 0&0&-2&1&1&0\\\ 0&1&1&-2&0&0\\\ 1&0&1&0&-2&0\\\
1&1&0&0&0&-2\end{bmatrix}.$ (58)

The Hessian matrix in Eq. (58) is positive semi-definite, with a single zero
in the eigenvalue spectrum corresponding to global rotations. After this mode
is integrated out by fixing the planar orientations of the spins, a straight-
forward minimization of the quadratic equation gives

$\delta\phi_{i}^{\text{min}}(\mathbf{r})=-\sum_{j}H^{-1}_{ij}g_{j}(\mathbf{r}).$
(59)

Thus, the minimized energy can be written as

$\mathcal{H}_{\text{min}}=\mathcal{H}^{(0)}+dE,$ (60)

where the energy contribution of the spin deviations is given by

$\displaystyle dE$ $\displaystyle=-\frac{1}{2}\sum_{\mathbf{r}}\sum_{ij}H$
${}^{-1}_{ij}\mu_{i}(\mathbf{r})\mu_{j}(\mathbf{r})$
$\displaystyle=\phantom{-}\mathcal{K}\sum_{\mathbf{r}}v(\mathbf{r})\Bigg{[}$
$\displaystyle\phantom{-}2v(\mathbf{r}+\mathbf{a}_{1})+2v(\mathbf{r}+\mathbf{a}_{2})+2v(\mathbf{r}+\mathbf{a}_{1}+\mathbf{a}_{2})$
$\displaystyle+v(\mathbf{r}+2\mathbf{a}_{1})+v(\mathbf{r}+2\mathbf{a}_{2})+v(\mathbf{r}+2\mathbf{a}_{1}+2\mathbf{a}_{2})$
$\displaystyle+2v(\mathbf{r}+\mathbf{a}_{1}-\mathbf{a}_{2})+2v(\mathbf{r}+2\mathbf{a}_{1}+\mathbf{a}_{2})+2v(\mathbf{r}+\mathbf{a}_{1}+2\mathbf{a}_{2})-15\Bigg{]},$
(61)

and the constant $\mathcal{K}$ is defined as

$\mathcal{K}=-\frac{\widetilde{D}_{2}^{2}}{12\widetilde{J}_{1}},$ (62)

Note that the positions of the Ising variables in Eq. (61) have been shifted
to better indicate the nature of the effective interactions. We see that the
spin deviations renormalize nearest-neighbor interactions between Ising
variables and also induce second- and third-neighbor in-plane
antiferromagnetic interactions (since $\widetilde{J}_{1}$ is negative). These
values of the coupling constants stabilize incommensurate phases in the
triangular Ising antiferromagnets Plumer  _et al._ (1993).

Qualitatively 222We note that a more systematic study of these phases is
currently underway., we can understand the nature of these unusual phases in
the following way. The $\widetilde{J}_{1}$ interactions establish
ferromagnetic ordering along the $z$-axis, ensuring $Q_{z}=0$.
$\widetilde{D}_{2}$ interactions couple small spin deviations to stabilize
incommensurate in-plane wavevectors with a fixed magnitude $Q_{I}$. LT
analysis of equation (60) gives a dispersion with a degenerate ring with
radius $Q_{c}$, in agreement with the MC results (Fig. 9 (c)). However, due to
the Ising nature of the ordering, we expect there to be a larger degeneracy in
the ground state than what is predicted by the LT method. Although it is not
evident from the analysis above, the phase diagrams in Fig. 5 indicate that
the larger values of $\widetilde{J}_{2}$ serve to tune $Q_{c}$, until it
either becomes zero ($\mathbf{Q}=\Gamma$) or reaches the zone boundary
($\mathbf{Q}=K,M$).

##  VII  Spin waves

###  VII.1 Dynamical matrix

In this section, we consider the elementary spin excitations in certain
magnetic phases presented in Sec. V. To do this, we solve the linearized
torque equation (11), as described in Sec. II.4. Once the appropriate local
coordinates are selected, the equations of motion can be generally written in
the matrix form as

$\frac{\text{d}}{\text{d}t}\begin{bmatrix}\widetilde{S}_{ix}(\mathbf{r},t)\\\
\widetilde{S}_{iy}(\mathbf{r},t)\end{bmatrix}=\sum_{\mathbf{r}^{\prime}}\sum_{j}\begin{bmatrix}\widetilde{\mathcal{A}}_{yxij}(\mathbf{r}-\mathbf{r}^{\prime})\widetilde{S}_{jx}(\mathbf{r}^{\prime},t)+\widetilde{\mathcal{A}}_{yyij}(\mathbf{r}-\mathbf{r}^{\prime})\widetilde{S}_{jy}(\mathbf{r}^{\prime},t)-\widetilde{\mathcal{A}}_{zzij}(\mathbf{r}-\mathbf{r}^{\prime})\widetilde{S}_{iy}(\mathbf{r},t)\\\
\widetilde{\mathcal{A}}_{zzij}(\mathbf{r}-\mathbf{r}^{\prime})S_{ix}(\mathbf{r})-\widetilde{\mathcal{A}}_{xyij}(\mathbf{r}-\mathbf{r}^{\prime})\widetilde{S}_{jy}(\mathbf{r}^{\prime},t)-\widetilde{\mathcal{A}}_{xxij}(\mathbf{r}-\mathbf{r}^{\prime})\widetilde{S}_{jx}(\mathbf{r}^{\prime},t)\end{bmatrix}$
(63)

where $\widetilde{\mathcal{A}}_{\alpha\beta
ij}(\mathbf{r}-\mathbf{r}^{\prime})$ are the elements of the coupling matrix
in the local coordinates, $\mathbf{r}$ and $\mathbf{r}^{\prime}$ determine the
positions of the magnetic unit cells, and $i$ and $j$ generally label the non-
equivalent magnetic sublattices. Eq. (63) is solved by first defining the
spatial Fourier transforms of the spin components as in Eq. (6), and then
Fourier transforming in the time domain
($\mathbf{S}_{i}(\mathbf{q},t)=\sum_{\omega}S_{i}(\mathbf{q},\omega)e^{-i\omega
t}$). The result is written as an eigenvalue equation:

$-i\omega(\mathbf{q})\begin{bmatrix}\mathbf{S}_{x}\\\
\mathbf{S}_{y}\end{bmatrix}=\begin{bmatrix}\mathbf{u}^{(xx)}(\mathbf{q})&\mathbf{u}^{(xy)}(\mathbf{q})\\\
\mathbf{u}^{(yx)}(\mathbf{q})&\mathbf{u}^{(yy)}(\mathbf{q})\end{bmatrix}\begin{bmatrix}\mathbf{S}_{x}\\\
\mathbf{S}_{y}\end{bmatrix}$ (64)

where

$\displaystyle u^{(xx)}_{ij}(\mathbf{q})$
$\displaystyle=\widetilde{\mathcal{A}}_{yxij}(\mathbf{q})$ (65) $\displaystyle
u^{(xy)}_{ij}(\mathbf{q})$
$\displaystyle=\widetilde{\mathcal{A}}_{yyij}(\mathbf{q})-\sum_{k}\widetilde{\mathcal{A}}_{yyjk}(\mathbf{q})\Delta_{ij}$
(66) $\displaystyle u^{(yx)}_{ij}(\mathbf{q})$
$\displaystyle=\sum_{k}\widetilde{\mathcal{A}}_{zzjk}(\mathbf{q})\Delta_{ij}-\widetilde{\mathcal{A}}_{xxij}(\mathbf{q})$
(67) $\displaystyle u^{(yy)}_{ij}(\mathbf{q})$
$\displaystyle=-\widetilde{\mathcal{A}}_{xyij}(\mathbf{q})$ (68)

Here, $\Delta_{ij}$ is a Kronecker delta function, and
$u^{(ab)}_{ij}(\mathbf{q})$ are the elements of the dynamical spin wave
matrix. For planar spin configurations, it is possible to choose local
coordinates such that $u^{(xx)}_{ij}(\mathbf{q})=u^{(yy)}_{ij}(\mathbf{q})=0$.
The frequencies $\omega(\mathbf{q})$ must generally be calculated numerically.
However, for some high symmetry points in the Brillouin zone, the dynamical
matrix can be diagonalized analytically.

###  VII.2  Symmetry properties in the decoupled and weak SOC limit

Before discussing the spin wave dispersions for the single-$\mathbf{Q}$
structures, it is useful to analyze the symmetry properties of the spin
configurations. We are interested in determining the groups of symmetries that
leave the magnetic structures unchanged, i.e. the stabilizer subgroups. The
stabilizers allow one to quickly determine the number and degeneracy of modes
at a given wavevector $\mathbf{q}$ in the Brillouin zone.

In the decoupled limit, the stabilizers are subgroups of $\mathcal{G}_{J}$,
which was derived in Sec. III.1. Out of all structures listed in Sec. VI.1,
$\bm{\Phi}_{0}^{(\bm{\Gamma})}$ (ferromagnetic) and
$\bm{\Phi}_{3}^{(\bm{\Gamma})}$ (collinear antiferromagnetic) configurations
deserve special attention, since any spin rotation around an axis collinear to
the spins leaves the spin configurations unchanged. Although these
configurations are very similar, their symmetry properties turn out to be
fundamentally different. The stabilizer group of the ferromagnetic state
contains all lattice symmetries and can therefore be written as a direct
product

$\mathcal{S}\left(\bm{\Phi}_{0}^{(\bm{\Gamma})}\right)=\mathrm{D}_{6h}^{(L)}\otimes\mathrm{SO}^{(S)}(2).$
(69)

On the other hand, the collinear antiferromagnetic state is not invariant
under lattice transformations that interchange A and B layers.

Figure 10: Spin wave dispersions for (a) collinear $\mathbf{Q}=\Gamma$, (b)
non-collinear $\mathbf{Q}=\Gamma$, and (c) $\mathbf{Q}=K$ single-$\mathbf{Q}$
structures stabilized in the decoupled limit. In (b), the labels
$\bm{\Phi}^{(\bm{\Gamma})}_{1,5}$ and $\bm{\Phi}^{(\bm{\Gamma})}_{2,4}$ are
used to indicate that the phases with $m=1,5$ (2,4) are degenerate and
therefore have the same spectra. The parameters used to calculate the
dispersions are $J_{2}<0,J_{1}/|J_{2}|=-1,J_{4}/|J_{2}|=-0.1$
($\bm{\Phi}^{(\bm{\Gamma})}_{0}$), $J_{2}<0,J_{1}/|J_{2}|=1,J_{4}/|J_{2}|=0.1$
($\bm{\Phi}^{(\bm{\Gamma})}_{3}$),
$J_{2}>0,J_{1}/|J_{2}|=-1,J_{4}/|J_{2}|=0.1$
($\bm{\Phi}^{(\bm{\Gamma})}_{1,5}$),
$J_{2}>0,J_{1}/|J_{2}|=1,J_{4}/|J_{2}|=-0.1$
($\bm{\Phi}^{(\bm{\Gamma})}_{2,4}$),
$J_{2}<0,J_{1}/|J_{2}|=-1.5,J_{4}/|J_{2}|=1$ ($\bm{\Phi}^{(\mathbf{K})}_{0}$),
and $J_{2}<0,J_{1}/|J_{2}|=1.5,J_{4}/|J_{2}|=-1$
($\bm{\Phi}^{(\mathbf{K})}_{3}$).

The remaining symmetry elements form a group $\mathrm{D}_{3h}^{(L)}$, and the
stabilizer can therefore be written as a semidirect product:

$\mathcal{S}\left(\bm{\Phi}_{3}^{(\bm{\Gamma})}\right)=\mathrm{D}_{3h}^{(L)}\otimes\left[\mathrm{SO}^{(S)}(2)\ltimes\mathrm{C}_{2}^{(SL)}\right],$
(70)

where $\mathrm{C}_{2}^{(SL)}$ contains a simultaneous $C_{2}$ rotation of the
lattice around $z$-axis and spins around an axis perpendicular to the spins.
This subtle difference in the structure of the stabilizers leads to
significant differences in the corresponding spin wave dispersion spectra.
Fig. 10 (a) demonstrates the differences between the dynamics of these two
spin configurations. The most notable difference is that in the case of
$\bm{\Phi}_{3}^{(\bm{\Gamma})}$, all branches are at least doubly degenerate,
leading to three distinct modes. However, as seen in the figure, in the case
of a ferromagnetic state, there are generally six non-degenerate modes. The
existence of the two-fold degeneracy in the antiferromagnetic case can be
attributed to the fact that at an arbitrary wavevector $\mathbf{q}$ the
dynamical matrix in Eq. (64) is invariant under a simultaneous lattice
inversion and reversal of the spin direction. This situation is equivalent to
the spin waves on a linear antiferromagnetic chain where two linearly-
polarized magnon modes have the same dispersion. As a result, the collinear
antiferromagnetic state has two linear Goldstone modes at the $\Gamma$ point,
corresponding to the in- and out-of-phase fluctuations of spins on A and B
sublattices. In contrast, the ferromagnetic state only has one (quadratic)
Goldstone mode.

In the case of the non-collinear single-$\mathbf{Q}$ configurations, the
symmetry analysis is simplified by the fact that the stabilizer subgroups are
now finite groups, since there are no spin rotations that leave these
structures invariant. One can see from Fig. 6 that for each non-collinear
$\bm{\Phi}_{m}^{(\bm{\Gamma})}$ configuration there is exactly one spin
rotation that can be combined with a given lattice symmetry transformation to
leave the spin configuration unchanged. Therefore, the stabilizer subgroups
for $\mathbf{Q}=\Gamma$ non-collinear phases are all isomorphic to
$\mathrm{D}_{6h}$. The situation is similar in the case of
$\bm{\Phi}_{m}^{(\mathbf{K})}$ structures. However, because certain lattice
symmetries are broken as a result of the cycloidal spin order, we are left
with the symmetry operations that belong to the group of the ordering
wavevector $\mathbf{Q}=K$. Thus, the stabilizers of
$\bm{\Phi}_{m}^{(\mathbf{K})}$ are all isomorphic to group $\mathrm{D}_{3h}$.
Figures 10 (b-c) show the dispersions for the non-collinear
single-$\mathbf{Q}$ spin configurations in phase diagrams in Fig. 3. Note that
the non-collinear nature of these phases implies that the Goldstone modes are
three-fold degenerate, since there are now three non-equivalent fluctuation
axes.

In the weakly-coupled limit, the stabilizers must be the subgroups of
$\mathcal{G}_{D}$. This does not change the symmetry of the non-collinear
single-$\mathbf{Q}$ states, although DM interactions lift the three-fold
degeneracy of the Goldstone modes, leaving only one zero-energy mode. Since
the out-of-plane spin rotations are no longer valid symmetry operations, the
stabilizers of the two collinear phases discussed above formally become
isomorphic to $\mathrm{D}_{6h}$, similar to the rest of the
$\bm{\Phi}_{m}^{(\bm{\Gamma})}$ configurations. We note that when one
considers the symmetry of the local spin components used for spin-wave
calculations, the stabilizer groups become the same (i.e. not just
isomorphic). Therefore, the spin wave eigenvectors in the local frame are
exactly the same for all six $\bm{\Phi}_{m}^{(\bm{\Gamma})}$
($\bm{\Phi}_{m}^{(\mathbf{K})}$) phases. This important property holds true
for all magnetic structures in the weak SOC limit and is one of the
consequence of the self-duality discussed in Sec. IV.

###  VII.3  Order by disorder and non-reciprocity

As discussed in Sec. VI.1, the classical energy of the collinear phases does
not depend on the DM constants, meaning that these spin configurations can be
rotated out-of-plane without any energy cost. Nevertheless, the energy of the
spin waves does, in fact, depend on the DM interactions. As a result, we
expect that the fluctuations would break the effective rotational symmetry by
selecting a particular orientation of the collinear structures via the order-
by-disorder mechanism Henley (1989); Chalker  _et al._ (1992). To prove
demonstrate predictions, we calculate the magnon free energy as function of
the out-of-plane angle $\theta$

$F(\theta)=\frac{1}{\beta}\sum_{\mathbf{q}}\ln\left(1-e^{-\beta\omega(\mathbf{q},\theta)}\right).$
(71)

The calculations of the frequencies $\omega(\mathbf{q},\theta)$ were performed
numerically using a $50\times 50\times 50$ discretized reciprocal lattice
grid. The result is shown in Fig. 11 (a).

Figure 11: (a) Magnon free energy as function of the out-of-plane angle
$\theta$ calculated numerically for the two collinear phases
($\bm{\Phi}_{0,3}^{(\bm{\Gamma})}$) with $J_{2}<0$, $J_{1}/|J_{2}|=0.5$,
$D_{1}=D_{2}=0.5$. The minima of the free energy occur at
$\theta=\frac{\pi}{2}$, indicating order by disorder. (b) Non-reciprocity in
the spin wave spectrum of the antiferromagnetic collinear state in the weak
SOC limit.

We can see that in both cases, thermal fluctuations rotate the spins
perpendicular to the kagome layers, thus preserving the continuous rotational
symmetry in these phases.

We note that in the case of the collinear antiferromagnet, the stabilization
of the out-of-plane order through thermal fluctuations introduces an
interesting new property in the spin wave spectrum, namely non-reciprocity
Rikken  _et al._ (2002); Zakeri  _et al._ (2010); Garcia-Sanchez  _et al._
(2014); Cheon  _et al._ (2018); dos Santos  _et al._ (2020); Borys  _et al._
(2021). Non-reciprocity implies that there is no unitary transformation that
could make the dynamical matrix (64) Hermitian. As a result, the spectrum is
no longer symmetric at $\omega(\pm\mathbf{q})$, as seen in Fig. 11 (b).

Finally, we note that the $\bm{\Phi}_{0,3}^{(\mathbf{K})}$ configurations also
possess the same accidental rotational degeneracy. This is evidenced from the
fact that the classical energy of these configurations is independent of the
DM constants. However, because the stabilizer groups of these phases are
significantly smaller than those of the collinear phases, the out-of-plane
rotations in this case generally lower the symmetry of the spin
configurations. Instead, the breaking of the continuous degeneracy occurs as a
result of small spin deviations that establish finite DM coupling.

###  VII.4 Excitations in the multi-$\mathbf{Q}$ phases

Next, we calculate the excitation spectra for the
$\bm{\Psi}_{m}^{(\mathbf{M})}$ phases. Unlike in the case of the
$\mathbf{Q}=K$ single-$\mathbf{Q}$ states, it is not possible to obtain all of
the spin wave modes by simply shifting the branches by $\pm K$, and one has to
construct a dynamical matrix where each of the four elements $u_{ij}^{(ab)}$
is a $24\times 24$ matrix (as the magnetic unit cell consists of 4
crystallographic unit cells). As a result, the excitation spectra for these
states consist of 24 modes.

Figure 12: Spin wave dispersion for the four out of six
$\bm{\Psi}_{m}^{(\mathbf{M})}$ configurations. The parameters used for the
calculation are $J_{2}>0$, $J_{1}/|J_{2}|=1$, $D_{1}/|J_{2}|=\mp 0.5$,
$D_{2}/|J_{2}|=\mp 0.5$ for $\bm{\Psi}_{2,4}^{(\mathbf{M})}$, and $J_{2}>0$,
$J_{1}/|J_{2}|=-1$, $D_{1}/|J_{2}|=\pm 0.5$, $D_{2}/|J_{2}|=\mp 0.5$ for
$\bm{\Psi}_{1,5}^{(\mathbf{M})}$.

Fig. 12 shows the dispersions for the $\bm{\Psi}_{m}^{(\mathbf{M})}$ phases
stabilized for $J_{2}>0$, $|J_{1}|=|J_{2}|$ in the weak SOC limit. As in the
case of the single-$\mathbf{Q}$ states, the DM interactions lead to a single
Goldstone mode, which in this case occurs at $\Gamma$ and M points, since the
period of the magnetic texture is double that of the lattice. Certain branches
in the spectra also exhibit other interesting features such as flat bands.
However further studies are needed for a full description of band topology in
these phases.

Figure 13: (a) A diagram illustrating the splitting of the classical energies
of the six planar irreps in the weak SOC limit. Here, the values of the
exchange and DM constants are chosen using the self-duality relations to give
the same energies in the zero anisotropy limit. Note that the energy scales
are not exact. (b) Distortion of the spin structures as a result of small SIA
(orange arrows) and bond-dependent anisotropy (grey arrows). The
$\bm{\Phi}_{1,4}^{(\bm{\Gamma})}$ remain unchanged since the spins are
collinear with the anisotropy axes. (c) Phase diagram in Fig. 5 (left) with
shifted phase boundaries (grey dashed lines) introduced by a small SIA
($K^{-}/|J_{2}|=0.05$).

##  VIII Effects of anisotropic interactions

###  VIII.1  Phase stability

Having described the important properties of the magnetic phases in the
decoupled and weak SOC limits, we now consider the case of the intermediate
SOC. We assume that in this limit, the SIA and bond-dependent anisotropy are
small but non-negligible.

As discussed in Sec. III, the anisotropic interactions break the axial
rotational symmetry, reducing the symmetry group of the Hamiltonian to the
paramagnetic group. As a result, the $\bm{\Phi}_{0}^{(\bm{\Gamma})}$
($\bm{\Phi}_{3}^{(\bm{\Gamma})}$) and $\bm{\Phi}_{2}^{(\bm{\Gamma})}$
($\bm{\Phi}_{5}^{(\bm{\Gamma})}$) states now belong to the same two-
dimensional irreps, while the remaining $\bm{\Phi}_{1}^{(\bm{\Gamma})}$ and
$\bm{\Phi}_{4}^{(\bm{\Gamma})}$ split into one-dimensional irreps (see Fig.
2). In general, the anisotropic interactions stabilize the spin configurations
that belong to the one-dimensional planar irreps, since in these structures
the spins are aligned collinear to the anisotropy axes. At the same time, the
configurations belonging to the same two-dimensional planar irreps remain
almost degenerate for small values of the anisotropy constants. This is
illustrated in Fig. 13 (a). We select the exchange and DM parameters using
self-duality relations in Sec. IV.3 such that in the absence of anisotropy the
energies of the six planar $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ configurations are
the same. The SIA and the in-plane bond-dependent interactions split the
energies by stabilizing $\bm{\Phi}_{1}^{(\bm{\Gamma})}$ and
$\bm{\Phi}_{4}^{(\bm{\Gamma})}$ (the sign of $K^{-}$ and $A_{2}^{(xy)}$
determines the orientations of spins). The out-of-plane bond anisotropy
further splits the energies based on the parity of the configurations (in the
figure, we choose $A_{1}^{(xy)}>0$, which stabilizes the structures that are
antisymmetric under space inversion).

Apart from the energy of the spin configurations, anisotropic interactions
have a significant impact on the spin structures themselves. Since some of the
planar phases now belong to the same two-dimensional irreps, the anisotropic
interactions couple these states, introducing small deviations from the irreps
in the weak SOC limit, as shown in Fig. 13 (b). This has been shown to play an
important effect on the physical properties of $\mathrm{Mn}_{3}X$ compounds,
where the anisotropic interactions couple the $\bm{\Phi}_{2}^{(\bm{\Gamma})}$
120-degree state and the ferromagnetic state, inducing a small magnetic moment
at zero field Soh  _et al._ (2020); Chen  _et al._ (2020); Zelenskiy  _et al._
(2021a). Our previous work in Ref. I has established that the relative
strengths of the SIA and the in-plane bond-anisotropy fix the overall
orientation of the spins as well as the direction of the spin deviations,
consequently determining both the magnitude and direction of the induced
magnetic moment in $\mathrm{Mn}_{3}X$ magnets.

Figure 14: (a) Diagram demonstrating the $\Gamma$-point splitting of the spin
wave modes in $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ phases in the three SOC limits.
Examples of broken symmetry in the excitations of (b)
$\bm{\Phi}_{2}^{(\bm{\Gamma})}$ and (c) $\bm{\Phi}_{5}^{(\bm{\Gamma})}$. The
solid blue and dashed red lines correspond to dispersions along the two paths
illustrated in the inset of figure (b). In both cases, we set
$K^{-}/|J_{2}|=0.1$, $J_{2}>0$, $J_{1}/|J_{2}|=1$, $D_{2}/|J_{2}|=-0.8$, and
set $D_{1}/|J_{2}|=0.8$ for $\bm{\Phi}_{2}^{(\bm{\Gamma})}$, and
$D_{1}/|J_{2}|=-0.8$ for $\bm{\Phi}_{5}^{(\bm{\Gamma})}$.

We note that the distortions of the spin structures introduced by the
anisotropic couplings may also provide a way for the experimental
quantification of these interactions. Techniques, such as elastic neutron
scattering, would allow one to determine the canting angles, which could in
turn be related to the values of the anisotropic parameters. Moreover, in the
case of the coupling between the ferromagnetic $\bm{\Phi}_{0}^{(\bm{\Gamma})}$
state and $\bm{\Phi}_{2}^{(\bm{\Gamma})}$, the induced magnetic moment would
serve as an excellent probe of the anisotropy in the system.

Although the discussion in this section so far has been in the context of the
$\mathbf{Q}=\Gamma$ phases, the general trends in the stabilization energy and
spin configurations mostly apply to the remaining types of phases, since every
spin structure discussed in this paper can be approximately constructed from
the rotated $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ states. As the deviations from the
$\bm{\Phi}_{m}^{(\bm{\Gamma})}$ local order parameters are small, we are still
able to make qualitative predictions about the stability of different phases
under small anisotropy. Therefore, we expect the anisotropic interactions to
stabilize those magnetic phases characterized by $m=1,4$, thus extending their
stability regions as compared to the weak SOC limit. This is clearly
demonstrated in Fig. 13 (c), where the phase boundaries are calculated using
$J_{2}>0$, $J_{1}/|J_{2}|=1$, $K^{-}/|J_{2}|=0.05$. The symmetry of the phase
diagram is broken, since $\mu_{0}^{(-1)}$ is no longer a valid self-duality
transformation. Notably, the in-plane bond-dependent anisotropy with
$A_{2}^{(xy)}=0.05$ produces a nearly identical shift of the phase boundaries.
The only exception to the established stability trend appears to be the
$\bm{\Lambda}_{4}$ phase which shrinks under the applied anisotropy. However,
this is not surprising given the result of Sec. VI.3 where the spin deviations
were shown to be the primary stabilizing factor of the Ising patterns.

###  VIII.2 Spin wave spectra

m | Stabilizer | irrep decomposition at $\Gamma$
---|---|---
0 | C${}_{2h}^{(SL)}$ | $A_{g}\oplus A_{u}\oplus 2B_{g}\oplus 2B_{u}$
1 | D${}_{6}^{(SL)}$ | $A_{2}\oplus B_{1}\oplus E_{1}\oplus E_{2}$
2 | C${}_{2h}^{(SL)}$ | $A_{g}\oplus A_{u}\oplus 2B_{g}\oplus 2B_{u}$
3 | D${}_{2}^{(SL)}$ | $A_{1}\oplus 2B_{1}\oplus B_{2}\oplus 2B_{3}$
4 | D${}_{3d}^{(SL)}$ | $2A_{2}^{\prime\prime}\oplus 2E^{\prime\prime}$
5 | D${}_{2}^{(SL)}$ | $A_{1}\oplus 2B_{1}\oplus B_{2}\oplus 2B_{3}$
Table 4: Stabilizer groups and irrep decomposition at the $\Gamma$ point for
the six $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ states in the intermediate anisotropy
limit.

The effects of anisotropic interactions are strongly manifested in the
excitations of the magnetic structures. Broken rotational symmetry implies
that the spin wave spectra are now completely gapped. Since the symmetry group
of the Hamiltonian is reduced down to the paramagnetic group, the stabilizer
subgroups of the spin configurations in the intermediate SOC limit are
typically small. In the case of the six $\mathbf{Q}=\Gamma$ phases, the
stabilizer groups and the corresponding irrep decompositions at the $\Gamma$
point are presented in table 4. We see that the stabilizers of states with
$m=0,2$ ($\mathrm{C}_{2h}$) and $m=3,5$ ($\mathrm{D}_{2}$) contain only one-
dimensional irreps, meaning that apart from the accidental mode crossings, the
spin wave spectra will generally only contain non-degenerate modes. However,
the remaining two configurations retain most of their symmetry from the weak
SOC limit and therefore have doubly-degenerate modes.

We now briefly summarize the effects of the SOC strength on the spin wave
spectra. Fig. 14 (a) demonstrates the evolution of the $\Gamma$ point
splitting of excitations in the six $\bm{\Phi}_{m}^{(\bm{\Gamma})}$ phases.
The transitions between the different SOC limits are accompanied by a
qualitative change in the spectra, associated with either a change in the
splitting or an opening of a gap. In addition, since the spin structures break
most of the symmetries of the Hamiltonian, the spectra in general will not be
symmetric in the Brillouin zone. This is demonstrated in Fig. 14 (b) and (c),
where two different paths in the Brillouin zone, related by a reflection yield
different values of $\omega(\mathbf{q})$.

These qualitative changes in the spin wave spectra provide a good probe for
the SOC strength in the AB-stacked kagome magnets. By comparing these results
to, for example, the inelastic neutron scattering or Raman scattering data,
one would be able to determine the appropriate SOC limit, and thus identify
the relevant spin interactions.

##  IX  Summary and Conclusions

In this work, we provide an extensive overview of the properties of the
general magnetic Hamiltonian (1) on hexagonal AB-stacked lattices. By studying
the symmetry of the model, we have determined the connection between the
strength of the SOC and the allowed spin symmetries in the system. We have
further identified three cases corresponding to decoupled, weak, and
intermediate SOC limits that yield different Hamiltonian symmetry groups. In
addition to the symmetries, we found a large number of self-duality
transformations that define the structure of the parameter space of the model.
Since these transformations directly depend on the symmetry of the physical
system, we identify three sets of duality transformations corresponding to the
three SOC limits.

The fundamental properties of the Hamiltonian allowed us to devise a strategy
for efficient exploration of the parameter space. We studied the ordered
phases in the decoupled and weak SOC limits, by constructing parameter-space
phase diagrams, using a combination of analytical LT and numerical MC methods.
We analyzed the structures of the resulting spin configurations, including
single- and multi-$\mathbf{Q}$ phases, and gave exact or nearly exact
parameterizations. Among the most interesting structures identified in this
work are the Ising-like patterns found in large pockets of the parameter space
in the weak and intermediate SOC limits. We determine that these states are
stabilized by small deviations from the idealized order parameters, and derive
the second order solution for the optimal spin canting, as well as the
stabilization energy.

Next, we calculated the spin-wave spectra of some of the single- and
multi-$\mathbf{Q}$ phases and analyzed the symmetry of the spin fluctuations.
We found that in the weak SOC limit the fluctuations drive the collinear
configurations out of the kagome plane, signifying order-by-disorder. As a
consequence of this rotation, the DM interactions make the excitations in the
collinear antiferromagnet phase non-reciprocal.

Finally, we study the effects of the intermediate SOC, manifested by the
additional SIA and bond-dependent interactions. We find that small amounts of
anisotropy can produce traceable changes in the structure and excitations of
the spin configurations, which opens the doors for potential experimental
quantification of these interactions.

As discussed throughout the paper, the results of this work are of direct
relevance to the $\mathrm{Mn}_{3}X$ family of compounds. For example, the
interplay of the anisotropic interactions has already been shown in our
previous work Ref. I to produce detectable changes in the static and dynamic
properties of these compounds. Another promising magnetic compound with the
AB-layered structure is the $\mathrm{Gd}_{3}\mathrm{Ru}_{4}\mathrm{Al}_{12}$
Chandragiri  _et al._ (2016); Matsumura  _et al._ (2019); Hirschberger  _et
al._ (2019). Although the magnetic parameters for these materials may not fall
in any of the non-trivial phases (such as multi-$\mathbf{Q}$ structures), one
may still be able to realize these states in systems where the SOC can be
tuned experimentally. For example, in thin-layer systems, the SOC close to the
surface can be enhanced by an addition of a heavy-metal capping layer Gradmann
(1986); Heinrich and Bland (1998). Importantly, this study has shown that the
unconventional phases often lie within the typical physical range of magnetic
parameters, including the DM and anisotropic couplings. In addition, the
existence of self-duality maps significantly increases the chances that a
compound with AB-layered kagome structure would exhibit non-trivial magnetic
phenomena. The last statement can be justified by recalling that in Sec. VI.3
the $\bm{\Lambda}$ phase was related via duality to a model with large
ferromagnetic out-of-plane exchange ($J_{1}\ll 0$), strong in-plane DM
interactions ($D_{2}\neq 0$) and negligible in-plane exchange ($J_{2}\approx
0$). Under normal circumstances, this situation would be extremely difficult
to realize, rendering the model as unphysical. However, the duality provides
us with a dual images that lie within a reasonable range of parameters.

In conclusion, we would like to point out that some of the assumptions made in
the beginning of this article could have an important impact on some of the
magnetic properties of $\mathrm{Mn}_{3}X$ systems. Since the interactions are
predominantly governed by the itinerant electrons, long-range interactions may
actually play an important role in stabilizing certain experimentally observed
phases. For example, further neighbor interactions along the
$\mathbf{\hat{z}}$-direction would be required in order to stabilize out-of-
plane spatial modulations through competition with the nn couplings Zelenskiy
_et al._ (2021b). This also means that the breathing anisotropy may also
affect certain properties of the magnetic phases. Furthermore, some recent
works have been dedicated to the importance of the coupling between the
elastic and magnetic degrees of freedom Theuss  _et al._ (2022); Reis  _et
al._ (2020); Sukhanov  _et al._ (2018, 2019). These are known to introduce
effective biquadratic spin interactions, which would further complicate the
theoretical analysis. As it stands however, this work provides an invitation
for further investigations of the rich magnetic phenomena offered by the AB-
stacked kagome compounds.

##  X  Acknowledgements

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC). The authors would like to thank J. S. R. McCoombs,
C Rudderham, D. Kalliecharan, and B. D. MacNeil for helpful discussions. We
acknowledge D. Kalliecharan’s suggestions for the visual presentation of this
paper. Finally, we thank S. H. Curnoe and J. P. F. LeBlanc for their useful
comments in regards to the symmetry and self-duality sections.

## References

  * Sadoc and Mosseri (1999) J.-F. Sadoc and R. Mosseri, Introduction to geometrical frustration, in  _Geometrical Frustration_, Collection Alea-Saclay: Monographs and Texts in Statistical Physics (Cambridge University Press, 1999) p. 1–13.
  * Diep (2005) H. T. Diep,  _Frustrated spin systems_ (World Scientific Publishing, 2005).
  * Moessner and Ramirez (2006) R. Moessner and A. P. Ramirez, Physics Today 59, 24 (2006), https://doi.org/10.1063/1.2186278 .
  * Lacroix  _et al._ (2011) C. Lacroix, P. Mendels, and F. Mila,  _Introduction to frustrated magnetism: materials, experiments, theory_ , Vol. 164 (Springer Science & Business Media, 2011).
  * Chalker  _et al._ (1992) J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Phys. Rev. Lett. 68, 855 (1992).
  * Harris  _et al._ (1992) A. B. Harris, C. Kallin, and A. J. Berlinsky, Phys. Rev. B 45, 2899 (1992).
  * Huse and Rutenberg (1992) D. A. Huse and A. D. Rutenberg, Phys. Rev. B 45, 7536 (1992).
  * Balents  _et al._ (2002) L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B 65, 224412 (2002).
  * Fu  _et al._ (2015) M. Fu, T. Imai, T.-H. Han, and Y. S. Lee, Science 350, 655 (2015).
  * Fujihala  _et al._ (2020) M. Fujihala, K. Morita, R. Mole, S. Mitsuda, T. Tohyama, S.-i. Yano, D. Yu, S. Sota, T. Kuwai, A. Koda, H. Okabe, H. Lee, S. Itoh, T. Hawai, T. Masuda, H. Sagayama, A. Matsuo, K. Kindo, S. Ohira-Kawamura, and K. Nakajima, Nature Communications 11, 3429 (2020).
  * Chen  _et al._ (2014) H. Chen, Q. Niu, and A. H. MacDonald, Phys. Rev. Lett. 112, 017205 (2014).
  * Nakatsuji  _et al._ (2015) S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212 (2015).
  * Kiyohara  _et al._ (2016) N. Kiyohara, T. Tomita, and S. Nakatsuji, Phys. Rev. Applied 5, 064009 (2016).
  * Nayak  _et al._ (2016) A. K. Nayak, J. E. Fischer, Y. Sun, B. Yan, J. Karel, A. C. Komarek, C. Shekhar, N. Kumar, W. Schnelle, J. Kübler, C. Felser, and S. S. P. Parkin, Science Advances 2, 10.1126/sciadv.1501870 (2016).
  * Ikhlas  _et al._ (2017) M. Ikhlas, T. Tomita, T. Koretsune, M.-T. Suzuki, D. Nishio-Hamane, R. Arita, Y. Otani, and S. Nakatsuji, Nature Physics 13, 1085 (2017).
  * Hong  _et al._ (2020) D. Hong, N. Anand, C. Liu, H. Liu, I. Arslan, J. E. Pearson, A. Bhattacharya, and J. S. Jiang, Phys. Rev. Materials 4, 094201 (2020).
  * Liu and Balents (2017) J. Liu and L. Balents, Phys. Rev. Lett. 119, 087202 (2017).
  * Park  _et al._ (2018) P. Park, J. Oh, K. Uhlířová, J. Jackson, A. Deák, L. Szunyogh, K. H. Lee, H. Cho, H.-L. Kim, H. C. Walker, D. Adroja, V. Sechovský, and J.-G. Park, npj Quantum Materials 3, 63 (2018).
  * Li  _et al._ (2019) X. Li, C. Collignon, L. Xu, H. Zuo, A. Cavanna, U. Gennser, D. Mailly, B. Fauqué, L. Balents, Z. Zhu, and K. Behnia, Nature Communications 10, 3021 (2019).
  * Reichlova  _et al._ (2019) H. Reichlova, T. Janda, J. Godinho, A. Markou, D. Kriegner, R. Schlitz, J. Zelezny, Z. Soban, M. Bejarano, H. Schultheiss, P. Nemec, T. Jungwirth, C. Felser, J. Wunderlich, and S. T. B. Goennenwein, Nature Communications 10, 5459 (2019).
  * Soh  _et al._ (2020) J.-R. Soh, F. de Juan, N. Qureshi, H. Jacobsen, H.-Y. Wang, Y.-F. Guo, and A. T. Boothroyd, Phys. Rev. B 101, 140411 (2020).
  * Chen  _et al._ (2020) Y. Chen, J. Gaudet, S. Dasgupta, G. G. Marcus, J. Lin, T. Chen, T. Tomita, M. Ikhlas, Y. Zhao, W. C. Chen, M. B. Stone, O. Tchernyshyov, S. Nakatsuji, and C. Broholm, Phys. Rev. B 102, 054403 (2020).
  * Zelenskiy  _et al._ (2021a) A. Zelenskiy, T. L. Monchesky, M. L. Plumer, and B. W. Southern, Phys. Rev. B 103, 144401 (2021a).
  * Kitaev (2006) A. Kitaev, Annals of Physics 321, 2 (2006).
  * Ran  _et al._ (2007) Y. Ran, M. Hermele, P. A. Lee, and X.-G. Wen, Phys. Rev. Lett. 98, 117205 (2007).
  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Yan  _et al._ (2011) S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
  * Messio  _et al._ (2012) L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. Lett. 108, 207204 (2012).
  * Bauer  _et al._ (2014) B. Bauer, L. Cincio, B. P. Keller, M. Dolfi, G. Vidal, S. Trebst, and A. W. W. Ludwig, Nature Communications 5, 5137 (2014).
  * Savary and Balents (2016) L. Savary and L. Balents, Reports on Progress in Physics 80, 016502 (2016).
  * Takagi  _et al._ (2019) H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Nature Reviews Physics 1, 264 (2019).
  * Bogdanov and Hubert (1994) A. Bogdanov and A. Hubert, Journal of Magnetism and Magnetic Materials 138, 255 (1994).
  * Rößler  _et al._ (2006) U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006).
  * Hayami and Motome (2021) S. Hayami and Y. Motome, Phys. Rev. B 103, 054422 (2021).
  * Leonov and Mostovoy (2015) A. O. Leonov and M. Mostovoy, Nature Communications 6, 8275 (2015).
  * Schulz  _et al._ (2012) T. Schulz, R. Ritz, A. Bauer, M. Halder, M. Wagner, C. Franz, C. Pfleiderer, K. Everschor, M. Garst, and A. Rosch, Nature Physics 8, 301 (2012).
  * Everschor-Sitte and Sitte (2014) K. Everschor-Sitte and M. Sitte, Journal of Applied Physics 115, 172602 (2014).
  * Neubauer  _et al._ (2009) A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Phys. Rev. Lett. 102, 186602 (2009).
  * Bruno  _et al._ (2004) P. Bruno, V. K. Dugaev, and M. Taillefumier, Phys. Rev. Lett. 93, 096806 (2004).
  * Kurumaji  _et al._ (2019) T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T. hisa Arima, and Y. Tokura, Science 365, 914 (2019).
  * He  _et al._ (2022) Y. He, S. Schneider, T. Helm, J. Gayles, D. Wolf, I. Soldatov, H. Borrmann, W. Schnelle, R. Schaefer, G. H. Fecher, B. Rellinghaus, and C. Felser, Acta Materialia 226, 117619 (2022).
  * Hirsch (1999) J. E. Hirsch, Phys. Rev. Lett. 83, 1834 (1999).
  * Chen and Byrnes (2019) T. Chen and T. Byrnes, Phys. Rev. B 99, 184427 (2019).
  * Zhang  _et al._ (2017) Y. Zhang, Y. Sun, H. Yang, J. Železný, S. P. P. Parkin, C. Felser, and B. Yan, Phys. Rev. B 95, 075128 (2017).
  * Hayami and Motome (2018) S. Hayami and Y. Motome, Phys. Rev. Lett. 121, 137202 (2018).
  * Akagi  _et al._ (2012) Y. Akagi, M. Udagawa, and Y. Motome, Phys. Rev. Lett. 108, 096401 (2012).
  * Ghosh  _et al._ (2016) S. Ghosh, P. O’Brien, C. L. Henley, and M. J. Lawler, Phys. Rev. B 93, 024401 (2016).
  * Moriya (1960) T. Moriya, Phys. Rev. 120, 91 (1960).
  * Luttinger and Tisza (1946) J. M. Luttinger and L. Tisza, Phys. Rev. 70, 954 (1946).
  * Zaliznyak and Zhitomirsky (1995) I. A. Zaliznyak and M. E. Zhitomirsky, J. Exp. Theor. Phys. 81, 579 (1995).
  * Maksimov  _et al._ (2019) P. A. Maksimov, Z. Zhu, S. R. White, and A. L. Chernyshev, Phys. Rev. X 9, 021017 (2019).
  * Glazer and Burns (2013) A. M. Glazer and G. Burns,  _Space groups for Solid State scientists_ (Elsevier, 2013).
  * Bradley and Cracknell (2010) C. J. Bradley and A. P. Cracknell,  _The mathematical theory of symmetry in solids: Representation theory for point groups and Space Groups_ (Clarendon Press, 2010).
  * Brinkman  _et al._ (1966) W. F. Brinkman, R. J. Elliott, and R. E. Peierls, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 294, 343 (1966).
  * Brinkman and Elliott (1966) W. Brinkman and R. J. Elliott, Journal of Applied Physics 37, 1457 (1966).
  * Brinkman (1967) W. Brinkman, Journal of Applied Physics 38, 939 (1967).
  * Litvin and Opechowski (1974) D. Litvin and W. Opechowski, Physica 76, 538 (1974).
  * Litvin (1977) D. B. Litvin, Acta Crystallographica Section A 33, 279 (1977).
  * Corticelli  _et al._ (2022) A. Corticelli, R. Moessner, and P. A. McClarty, Phys. Rev. B 105, 064430 (2022).
  * Liu  _et al._ (2022) P. Liu, J. Li, J. Han, X. Wan, and Q. Liu, Phys. Rev. X 12, 021016 (2022).
  * Essafi  _et al._ (2017) K. Essafi, O. Benton, and L. D. C. Jaubert, Phys. Rev. B 96, 205126 (2017).
  * Kramers and Wannier (1941) H. A. Kramers and G. H. Wannier, Phys. Rev. 60, 252 (1941).
  * Onsager (1944) L. Onsager, Phys. Rev. 65, 117 (1944).
  * Wannier (1945) G. H. Wannier, Rev. Mod. Phys. 17, 50 (1945).
  * Wegner (1971) F. J. Wegner, Journal of Mathematical Physics 12, 2259 (1971).
  * Wu and Wang (1976) F. Y. Wu and Y. K. Wang, Journal of Mathematical Physics 17, 439 (1976).
  * Savit (1980) R. Savit, Rev. Mod. Phys. 52, 453 (1980).
  * Chaloupka and Khaliullin (2015) J. c. v. Chaloupka and G. Khaliullin, Phys. Rev. B 92, 024413 (2015).
  * Kimchi and Vishwanath (2014) I. Kimchi and A. Vishwanath, Phys. Rev. B 89, 014414 (2014).
  * Essafi  _et al._ (2016) K. Essafi, O. Benton, and L. D. C. Jaubert, Nature Communications 7, 10297 (2016).
  * Note (1) We note that up to now the discussion has been completely general and is applicable to any spin system described by a set of symmetries. The derivation of the general conditions for self-duality go beyond the scope of this work and are a subject of an ongoing investigation.
  * Plumer  _et al._ (1993) M. L. Plumer, A. Mailhot, R. Ducharme, A. Caillé, and H. T. Diep, Phys. Rev. B 47, 14312 (1993).
  * Note (2) We note that a more systematic study of these phases is currently underway.
  * Henley (1989) C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Rikken  _et al._ (2002) G. L. J. A. Rikken, C. Strohm, and P. Wyder, Phys. Rev. Lett. 89, 133005 (2002).
  * Zakeri  _et al._ (2010) K. Zakeri, Y. Zhang, J. Prokop, T.-H. Chuang, N. Sakr, W. X. Tang, and J. Kirschner, Phys. Rev. Lett. 104, 137203 (2010).
  * Garcia-Sanchez  _et al._ (2014) F. Garcia-Sanchez, P. Borys, A. Vansteenkiste, J.-V. Kim, and R. L. Stamps, Phys. Rev. B 89, 224408 (2014).
  * Cheon  _et al._ (2018) S. Cheon, H.-W. Lee, and S.-W. Cheong, Phys. Rev. B 98, 184405 (2018).
  * dos Santos  _et al._ (2020) F. J. dos Santos, M. dos Santos Dias, and S. Lounis, Phys. Rev. B 102, 104401 (2020).
  * Borys  _et al._ (2021) P. Borys, O. Kolokoltsev, N. Qureshi, M. L. Plumer, and T. L. Monchesky, Phys. Rev. B 103, 144411 (2021).
  * Chandragiri  _et al._ (2016) V. Chandragiri, K. K. Iyer, and E. V. Sampathkumaran, Journal of Physics: Condensed Matter 28, 286002 (2016).
  * Matsumura  _et al._ (2019) T. Matsumura, Y. Ozono, S. Nakamura, N. Kabeya, and A. Ochiai, Journal of the Physical Society of Japan 88, 023704 (2019).
  * Hirschberger  _et al._ (2019) M. Hirschberger, T. Nakajima, S. Gao, L. Peng, A. Kikkawa, T. Kurumaji, M. Kriener, Y. Yamasaki, H. Sagayama, H. Nakao, K. Ohishi, K. Kakurai, Y. Taguchi, X. Yu, T.-h. Arima, and Y. Tokura, Nature Communications 10, 5831 (2019).
  * Gradmann (1986) U. Gradmann, Journal of Magnetism and Magnetic Materials 54-57, 733 (1986).
  * Heinrich and Bland (1998) B. Heinrich and J. A. C. Bland,  _Ultrathin Magnetic Structures_ , Vol. 1 & 2 (Springer Verlag, Berlin, 1998).
  * Zelenskiy  _et al._ (2021b) A. Zelenskiy, J. A. Quilliam, A. Y. Shapiro, and G. Quirion, Phys. Rev. B 103, 224422 (2021b).
  * Theuss  _et al._ (2022) F. Theuss, S. Ghosh, T. Chen, O. Tchernyshyov, S. Nakatsuji, and B. J. Ramshaw, Phys. Rev. B 105, 174430 (2022).
  * Reis  _et al._ (2020) R. D. d. Reis, M. Ghorbani Zavareh, M. O. Ajeesh, L. O. Kutelak, A. S. Sukhanov, S. Singh, J. Noky, Y. Sun, J. E. Fischer, K. Manna, C. Felser, and M. Nicklas, Phys. Rev. Materials 4, 051401 (2020).
  * Sukhanov  _et al._ (2018) A. S. Sukhanov, S. Singh, L. Caron, T. Hansen, A. Hoser, V. Kumar, H. Borrmann, A. Fitch, P. Devi, K. Manna, C. Felser, and D. S. Inosov, Phys. Rev. B 97, 214402 (2018).
  * Sukhanov  _et al._ (2019) A. S. Sukhanov, M. S. Pavlovskii, P. Bourges, H. C. Walker, K. Manna, C. Felser, and D. S. Inosov, Phys. Rev. B 99, 214445 (2019).