# Non-Coplanar Magnetic Orders in Classical Square-Kagome Antiferromagnets

Martin Gembé  Institute for Theoretical Physics, University of Cologne, 50937
Cologne, Germany  Heinz-Jürgen Schmidt  Fachbereich Physik, Universität
Osnabrück, 49069 Osnabrück, Germany  Ciarán Hickey  Institute for Theoretical
Physics, University of Cologne, 50937 Cologne, Germany  School of Physics,
University College Dublin, Belfield, Dublin 4, Ireland  Centre for Quantum
Engineering, Science, and Technology, University College Dublin, Dublin 4,
Ireland  Johannes Richter  Institut für Physik, Otto-von-Guericke-Universität
Magdeburg, 39016 Magdeburg, Germany  Max-Planck-Institut für Physik Komplexer
Systeme, Nöthnitzer Straße 38, 01187 Dresden, Germany  Yasir Iqbal  Department
of Physics and Quantum Centers in Diamond and Emerging Materials (QuCenDiEM)
group, Indian Institute of Technology Madras, Chennai 600036, India  Simon
Trebst  Institute for Theoretical Physics, University of Cologne, 50937
Cologne, Germany

(August 28, 2024)

###### Abstract

Motivated by the recent synthesis of a number of Mott insulating square-kagome
materials, we explore the rich phenomenology of frustrated magnetism induced
by this lattice geometry, also referred to as the squagome or shuriken
lattice. On the classical level, square-kagome antiferromagnets are found to
exhibit extensive degeneracies, order-by-disorder, and non-coplanar ordering
tendencies, which we discuss for an elementary, classical Heisenberg model
with nearest-neighbor and cross-plaquette interactions. Having in mind that
upon introducing quantum fluctuations non-coplanar order can melt into chiral
quantum spin liquids, we provide detailed information on the multitude of non-
coplanar orders, including some which break rotational symmetry (possibly
leading to nematic quantum orders), as well as a number of (incommensurate)
spin spiral phases. Using extensive numerical simulations, we also discuss the
thermodynamic signatures of these phases, which often show multi-step thermal
ordering. Our comprehensive discussion of the classical square-kagome
Heisenberg model, often drawing comparisons to the conventional kagome
antiferromagnet, sets the stage for future explorations of quantum analogs of
the various phases, either conceptually such as in quantum spin-1/2
generalizations of our model or experimentally such as in the Cu-based
candidate materials.

##  I Introduction

Classical Heisenberg spin models with frustrated interactions are known to
host a rich variety of magnetic orders such as collinear, coplanar, or
helimagnetic states [1, 2, 3]. Of particular interest are non-Bravais lattices
(with more than one atom per unit cell) which offer the possibility of
stabilizing non-coplanar magnetic ordering [1, 4]. Such non-coplanar ground
states distinguish themselves from other magnetic orders by exhibiting a
scalar spin chirality. Notably, the spontaneous breaking of such a
$\mathbb{Z}_{2}$ (chiral) symmetry manifests itself in a finite-temperature
phase transition – even in two spatial dimensions [5], while other types of
magnetic order are, for the Heisenberg models of interest here, subject to the
Mermin-Wagner theorem [6]. The latter implies that, for the thermodynamic
limit of infinite system size, any fluctuation-driven phase transition in two
spatial dimensions, which only breaks the continuous spin symmetry, occurs at
zero temperature. Finite systems (often explored in numerical simulations)
will, however, exhibit a thermal crossover to magnetic order at some finite
temperature, with the accompanying entropy release leading to a peak in the
specific heat. This should be distinguished from cooperative paramagnetic
phases [7], which defy magnetic ordering tendencies due to the existence of
substantial residual entropies, even at temperatures orders of magnitude below
the coupling scales and for infinitely large systems. In classical Heisenberg
models, the thermodynamics of such cooperative paramagnetic phases, also
referred to as classical spin liquids, is typically signified by a plateau in
the specific heat [8].

The two themes, the formation of non-coplanar magnetic order and cooperative
paramagnetic phases, are conceptually tied when looking at their quantum
mechanical counterparts. By melting non-coplanar magnetic order via quantum
fluctuations, e.g. by going to small spins auch as $S=1/2$, one could possibly
restore spin rotational symmetry, i.e., realize a non-magnetic quantum ground
state. If the chiral symmetry breaking present in the parent classical
magnetic order would persist (at some finite temperature scale), one would
realize a much sought after chiral quantum spin liquid phase [9, 10, 11].
Similarly, the inclusion of quantum fluctuations on cooperative paramagnetic
ground states provides another promising route towards realizing
unconventional quantum phases such as quantum spin liquids, valence bond
crystals [12, 13, 14, 15], or spin and lattice nematics [16].

An ideal playground to explore this physics in experiment has come in the
arrival of materials based on the novel square-kagome lattice geometry [17],
whose potential to host intricately textured magnetic ground states or quantum
spin liquid phases is currently under much investigation [18]. Indeed, no sign
of long-range magnetic order down to 50 mK has been observed in the spin
$S=1/2$ Cu2+ based materials KCu6AlBiO4(SO4)5Cl [19] and
Na6Cu7BiO4(PO4)4[Cl,(OH)]3 [20] despite having large negative Curie-Weiss
temperatures of $-237\text{\,}\mathrm{K}$ and $-212\text{\,}\mathrm{K}$,
respectively. On the other hand, their sister compound KCu7(TeO4)(SO4)5Cl
presumably orders into a non-collinear antiferromagnetic structure [21]. In
general, the model Hamiltonians for these materials can host up to three
symmetry inequivalent couplings on the three sides of the elementary triangles
as well as potential longer-range Heisenberg couplings across the octagonal
plaquettes, see Fig. 1, whose presence could be resolved by ab-initio density
functional theory calculations. The latter are, in contrast to diagonal square
couplings, a key ingredient towards stabilizing non-coplanar magnetic orders
on the classical level and, potentially, chiral spin liquids in the quantum
realm. This is in a spirit similar to the diagonal couplings across hexagons
on the kagome lattice which are known to yield non-coplanar spin structures
dubbed cuboc orders [5, 22, 23, 1]. The details of these non-coplanar states
are, by their very nature, rather sensitive to the underlying lattice geometry
with unique features expected for the square-kagome lattice geometry at hand.

The purpose of this manuscript is to set the staging ground for future
explorations of square-kagome antiferromagnets by providing a comprehensive
discussion of their physics in the classical realm and identifying its unique
features. To this end, we investigate the ground state and thermodynamics of
the classical Heisenberg model on the square-kagome lattice in the presence of
nearest-neighbor $(J_{1},J_{2},J_{3})$ couplings as well as cross-plaquette
interactions inside the octagons, $J_{\times}$ and $J_{+}$, as indicated in
Fig. 1. Our analysis is based on extensive classical Monte Carlo simulations
and an analytical construction of ground states (beyond the Luttinger-Tisza
approach), which is shown to be rendered exact for some orders. As summarized
in the phase diagram of Fig. 5 below, we find a rich variety of non-coplanar
magnetic orders with cuboctohedral symmetry, including types of cuboc order
not found on the kagome lattice (or any other known lattice geometry). In
addition, we report a multitude of non-coplanar incommensurate spirals, in
addition to commensurate coplanar orders. Exploring the quantum analogs of
these phases in the future, either conceptually such as in quantum spin-1/2
generalizations of our model or experimentally such as in the Cu-based
candidate materials, might prove fruitful in identifying chiral quantum spin
liquids.

Figure 1: Square-kagome lattice and interactions. The square-kagome lattice,
also referred to as the squagome or shuriken lattice in the literature,
consists of two sets of topologically distinct sites – square sites and bow-
tie sites. Nearest-neighbor interactions are referred to as $J_{1}$ (square)
and $J_{2}$, $J_{3}$ (bow-tie), respectively. Additionally, we introduce
further neighbor cross-plaquette $J_{+}$-bonds ($J_{\times}$-bonds) as
indicated in blue (red). The unit cell contains six spins, namely the four
square sites and two bow-tie sites of a single shuriken star.

##  II Nearest-Neighbor Model

With just nearest-neighbor Heisenberg interactions the square-kagome model
shares much of the same physics as the nearest-neighbor kagome Heisenberg
model. Below, we briefly summarize some of the known results for the ground
states of the square-kagome model that can be inferred from the conventional
kagome antiferromagnet. We then move on to discuss its finite temperature
physics and explore critical fluctuations, going beyond what has been studied
for the conventional kagome scenario.

Figure 2: The nearest-neighbor model. (a) Zero-temperature phase diagram as a
function of nearest-neighbor couplings $J_{2}$ and $J_{3}$, reproduced from
Ref. [24]. (b) Specific heat traces for three representative points along the
$J_{2}=J_{3}$ diagonal in the phase diagram of panel (a). The one for the
isotropic point, $J_{1}=J_{2}=J_{3}=1$, closely resembles the well-studied
specific heat trace of the kagome AFM [25, 8]. (c) Specific heat scan across
the phase diagram of panel (a) at fixed, low temperature of $T=0.04$ (in units
of $J_{1}=1$).

###  II.1 Ground states

The Hamiltonian with nearest-neighbor couplings can be written as

$\displaystyle\mathcal{H}=\sum_{\left<i,j\right>\in
a}J_{a}\,\mathbf{S}_{i}\cdot\mathbf{S}_{j}\,,$ (1)

where $a=\\{1,2,3\\}$ runs over the three different types of bonds, as in Fig.
1. However, it can be more easily understood by rewriting it as

$\displaystyle\mathcal{H}=\sum_{i,j,k\in\triangle}$
$\displaystyle\left[\frac{1}{2}\left(\sqrt{\frac{J_{1}J_{3}}{J_{2}}}\mathbf{S}_{i}+\sqrt{\frac{J_{1}J_{2}}{J_{3}}}\mathbf{S}_{j}+\sqrt{\frac{J_{2}J_{3}}{J_{1}}}\mathbf{S}_{k}\right)^{2}\right.$
$\displaystyle\left.\,\,\,\,-\frac{1}{2}\left(\frac{J_{1}J_{3}}{J_{2}}+\frac{J_{1}J_{2}}{J_{3}}+\frac{J_{2}J_{3}}{J_{1}}\right)\right]\,,$
(2)

where the sum is now over all four types of elementary triangles of the
lattice and we have assumed all couplings to be antiferromagnetic, $J_{a}>0$
(such that all of the arguments of the square roots are positive). The
Hamiltonian of the kagome Heisenberg model can be written in the same form,
but with a key distinction being the number and nature of elementary triangles
that are summed over. The new form of the Hamiltonian allows us to easily
construct a special class of classical ground states. These are the states
that satisfy the constraint

$\left(\sqrt{\frac{J_{1}J_{3}}{J_{2}}}\mathbf{S}_{i}+\sqrt{\frac{J_{1}J_{2}}{J_{3}}}\mathbf{S}_{j}+\sqrt{\frac{J_{2}J_{3}}{J_{1}}}\mathbf{S}_{k}\right)=0\,.$
(3)

where $\mathbf{S}_{i},\mathbf{S}_{j}$ are the spins on the squares and
$\mathbf{S}_{k}$ is the spin on the bow-ties. To see that such states are
indeed ground states note that the Hamiltonian, for a given set of parameters,
is the sum of a squared term, which is always positive, and a constant term.
The minimal possible energy is thus obtained when the squared term is
precisely zero, i.e., the constraint above. However, there is the additional
constraint that the spins at each site are all properly normalized,
$\mathbf{S}_{i}=1\,\,\forall\,\,i$. Satisfying both of these constraints is
only possible within a restricted region of parameter space, namely

$-2\leq\frac{J_{2}J_{3}}{J_{1}^{2}}-\left(\frac{J_{2}}{J_{3}}+\frac{J_{3}}{J_{2}}\right)\leq+2\,.$
(4)

Thus, within this parameter region, all spin configurations satisfying the
constraint (3) on each and every triangle are guaranteed to be bonafide ground
state spin configurations. This includes both globally coplanar and non-
coplanar spin configurations as, though the three spins in each triangle are
constrained to lie within the same plane, it is not necessary that the planes
for different triangles are the same. The resulting unusual and highly
degenerate phase is the classical spin liquid indicated in the phase diagram
of Fig. 2(a) and previously reported in Ref. [24], which shares the same
qualitative physics as the classical spin liquid found in the distorted kagome
version of the model [26, 27, 28, 29].

At the isotropic point, $J_{1}=J_{2}=J_{3}$, the spins in each triangle lie
within the same plane at an angle of exactly $120^{\circ}$ away from one
another. On the other hand, away from the isotropic point, e.g. along the
diagonal line $J_{2}=J_{3}$, the coplanar spin configuration for a single
triangle obeying the constraint can be written as

$\displaystyle\mathbf{S}_{i}$
$\displaystyle=(-J_{2}/(2J_{1}),\,-\sqrt{4-(J_{2}/J_{1})^{2}}/2,\,0),$
$\displaystyle\mathbf{S}_{j}$
$\displaystyle=(-J_{2}/(2J_{1}),\,+\sqrt{4-(J_{2}/J_{1})^{2}}/2,\,0),$
$\displaystyle\mathbf{S}_{k}$ $\displaystyle=(1,\,0,\,0)\,,$ (5)

where we have fixed the spins to lie in the $xy$-plane for simplicity. As one
increases $J_{2}$ from $J_{2}=0$ to $J_{2}=2J_{1}$, the angle $\theta$ between
$\mathbf{S}_{k}$ (the spin on the bow-tie) and the other two spins (on the
squares) increases from $\pi/2$ to $\pi$, passing through $2\pi/3$ exactly at
the isotropic point $J_{2}=J_{1}$ [15].

Figure 3: The nearest-neighbor model II. (a) Specific heat traces for $0.0\leq
J_{2}=J_{3}\leq 1.0$ (in steps of 0.2). (b) Specific heat traces for $1.0\leq
J_{2}=J_{3}\leq 2.0$. Of special interest is the curve for $J_{2}=J_{3}=2.0$
which is the only one which converges to a value of $c_{V}=10/12$. (c) Finite
temperature specific heat along the diagonal $J_{2}=J_{3}$ illustrating the
$c_{V}=1$ plateau of the cooperative paramagnet (dark red area), the coplanar
phase with $c_{V}=11/12$ (light red area), the critical fan with $c_{V}\leq
10/12$ (grey area), and the paramagnetic phase (black region).

For globally coplanar spin configurations, which are the relevant
configurations at the lowest temperatures as we will see in the next section,
as well as the constraints already mentioned, there is one additional form of
constraint [26]. It is related to how the spins in the triangles around the
square and octagonal plaquettes of the lattice are arranged. For globally
coplanar ground state spin configurations, we can define on each triangle a
chirality variable, $\eta_{a}=\pm 1$, which encodes whether the spins rotate
clockwise or anti-clockwise as one goes from say $i$ to $j$ to $k$ [30] (keep
in mind that the angles between the spins in each triangle are fixed, e.g. all
three angles are fixed to $2\pi/3$ at the isotropic point). Now, starting from
an initial spin which points in some specific direction, if one travels in a
closed loop on the lattice then one must return back to that same initial spin
pointing in that same specific direction. However, as one travels along each
bond the chirality variables dictate in which direction the spins rotate, and
so in order to get back to the same initial spin there is a constraint on the
sum of the chirality variables along the closed loop. At the isotropic point,
we require $\sum_{a}\eta_{a}=0$ for the four triangles surrounding a square,
and $\sum_{a}\eta_{a}=0,\pm 6$ for the eight triangles surrounding an octagon.
Away from the isotropic point there will in general be more stringent
constraints as the angle between spins is no longer $2\pi/3$, but instead some
angle incommensurate with respect to $2\pi$. This results, for a given lattice
size, in a smaller number of allowed coplanar ground state spin
configurations, just as in the analogous kagome case [26].

Outside of the classical spin liquid region there are two distinct Néel states
in the phase diagram of Fig. 2(a), depending on whether $J_{2}$ or $J_{3}$
dominates. In each phase, spins are arranged antiferromagnetically along the
bonds with the dominant coupling, and ferromagnetically along the weaker
bonds. On the other hand, for dominant $J_{2}=J_{3}$, there is an up-up-down
(UUD) state, sometimes also referred to as a Lieb ferrimagnet [31], in which
the spins on the squares point along one direction while the spins on the bow-
ties point in the opposite direction (this can be seen from Eq. (5) with
$J_{2}=2J_{1}$, the critical point between the classical spin liquid and UUD
state).

###  II.2 Finite-temperature physics

The unusual nature of the classical spin liquid phase is naturally revealed by
examining how the specific heat behaves as a function of temperature. For two-
dimensional Heisenberg models with finite-range exchange interactions, true
long-range magnetic order cannot set in at any non-zero temperature, as laid
out in the Mermin-Wagner theorem [6]. However, it is possible for quasi-long-
range order to develop, signaled by a peak in the specific heat. Our
discussion, in the following, of such quasi-long-range “orders” is based on an
analysis of the symmetry of spin-spin correlations at distances shorter than
the correlation radius. In the case of non-coplanar, i.e., chiral orders,
however, a true thermal phase transition associated with the breaking of
discrete symmetries occurs.

At the lowest temperatures, one generically expects that, in the thermodynamic
limit, the specific heat per site $c_{V}\rightarrow 1$ as $T\rightarrow 0$.
This is because each spin is free to fluctuate about its ordered ground state
in two orthogonal directions on the unit sphere. These two quadratic modes
each contribute $(1/2)\cdot T$, as dictated by classical equipartition, to the
energy and thus $1/2$ to the specific heat (setting $k_{B}=1$). However, as
first discussed in the context of the kagome antiferromagnet [25], this simple
counting breaks down within the classical spin liquid region, as well as in an
extended finite temperature fan about the critical lines in the phase diagram.
The breakdown is due to the entropic selection of a subset of ground state
spin configurations, those which carry the largest entropy and thus the lowest
free energy at finite temperature. Classical fluctuations about this favored
subset include one or more zero modes at the harmonic level, which contribute
$1/4$ (i.e. they are quartic modes), rather than $1/2$, to the specific heat.
The deviation of the low-temperature specific heat from $c_{V}\rightarrow 1$
thus serves as a signature of this phenomenon of thermal order-by-disorder.

In Fig. 2(b), we show the specific heat as a function of temperature for three
special parameter points along the diagonal line $J_{2}=J_{3}$, with a number
of curves in between shown in Fig. 3(a) and (b). (i) First, starting with
$J_{2}=J_{3}=0$, we have the trivial limit of fully disconnected squares.
Alternatively, one can think of this limit as consisting of decoupled four-
site periodic Heisenberg chains. The spins order in a simple antiferromagnetic
arrangement within each square. There are 8 quadratic modes per square, minus
two due to the global rotational symmetry of the antiferromagnetic moment.
This leaves us with 6 independent quadratic modes, and a contribution to the
specific heat per site as $c_{V}\rightarrow[6\cdot(1/2)]/6=1/2$ as
$T\rightarrow 0$ [32]. (ii) At the isotropic point, there are three distinct
regimes, which share the same physics as the isotropic kagome model at finite
temperatures [8] (the similarity even extends to the spin-$1/2$ quantum case
[33]). There is the usual high-temperature paramagnetic region, followed by a
cooperative paramagnetic regime and finally a coplanar state at the lowest
temperatures. These three regimes can be observed throughout the classical
spin liquid phase. The cooperative paramagnet is clearly distinguished by a
plateau in the specific heat with $c_{V}\approx 1$. Within this temperature
window, the system fluctuates between the full (extensive) number of states
within the ground state manifold that satisfy the constraint in Eq. (3). At
lower temperatures, within the coplanar states, fluctuations select the subset
of globally coplanar states within the ground state manifold via the entropic-
driven order-by-disorder mechanism. This is accompanied by
$c_{V}\rightarrow[10\cdot(1/2)+2\cdot(1/4)]/6=11/12$ as $T\rightarrow 0$ due
to the presence of one zero mode per triangle (thus two zero modes per unit
cell) within the spectrum of classical harmonic fluctuations (very similar to
the conventional kagome case [25]). (iii) At $J_{2}=2J_{1}$ we have the
transition between the classical spin liquid and UUD state. Precisely at this
critical point an additional zero mode per triangle leads to
$c_{V}\rightarrow[8\cdot(1/2)+4\cdot(1/4)]/6=10/12$ as $T\rightarrow 0$.

Figure 4: Structure factors of the $J_{1}$-$J_{2}$-$J_{3}$ model. Structure
factors for a system with 864 Spins ($L=12$) for different temperatures
between T = 0.003 and T = 0.2 (left to right) and for different values of
$J_{2}$ and $J_{3}$ as indicated in the phase diagram in the left column (top
to bottom). The shown squares extent in reciprocal space from $-4\pi$ to
$4\pi$ in both dimensions. Rows (a) to (c) show structure factors on the
diagonal $J_{2}=J_{3}$ within the classical spin liquid regime, especially for
the isotropic point $J_{1}=J_{2}=J_{3}=1.0$ in row (b), which agrees well with
the result obtained via large-$N$ analysis in [13]. Rows (d) and (e) as well
show structure factors in the classical spin liquid regime, but off the
diagonal with different values for $J_{2}$ and $J_{3}$. All structure factors
within the classical spin liquid regime show sharp pinch point-like features
upon entering the cooperative paramagnetic regime at around $T=0.01$, which
broaden with increasing temperature. The statistical noise at lower
temperatures is due to freezing in Monte Carlo sampling. The lower three rows
(f) to (h) display structure factors outside the classical spin liquid regime,
namely directly on the transition line to the UUD phase, $J_{2}=J_{3}=2.0$, in
(f), within the UUD phase, $J_{2}=J_{3}=2.5$, in (g), and in the
$J_{1}$-$J_{3}$-Néel ordered phase, $J_{2}=0.25$ and $J_{3}=2.5$, in (h).

The evolution of the specific heat between these three special points, as one
increases $J_{2}=J_{3}$, is shown in Fig. 3. Starting from the trivial limit,
$J_{2}=J_{3}=0$, the width of the plateau that develops at $c_{V}\approx 1/2$
shrinks as the interaction scale that couples squares, $J_{2}=J_{3}$, grows.
At low temperatures, the characteristic two temperature regimes emerge, with
the cooperative paramagnet and coplanar regimes giving rise to plateaus at
$c_{V}=1$ and $11/12$, respectively. The $1/2$-feature disappears completely
at the isotropic point. As one moves past the isotropic point, the width of
the cooperative paramagnet plateau starts to decrease, until it eventually
disappears completely at the critical point $J_{2}=2J_{1}$. At the same time,
in the temperature regime above the cooperative paramagnet, there is a
substantial drop in the specific heat as a result of increasing fluctuations
on approaching the critical point (which damps the entropy loss in this
temperature regime and thus the specific heat). This can be more clearly seen
in Fig. 3(c), where the drop manifests itself as a finite temperature fan,
emanating from the zero-temperature critical point – reminiscent of the fan-
like structure at quantum critical points [34].

###  II.3 Spin-spin correlations

To explore the formation of quasi-long-range order and cooperative
paramagnetic phases in the absence of any true thermal phase transitions, we
turn to the static spin structure factor, i.e. the Fourier transform of the
equal-time real-space spin-spin correlations. Sharp peak-like features
indicate the formation of quasi-long-range order, while the tell-tale
signatures of cooperative paramagnets are pinch points in momentum space,
which map to algebraically decaying correlations in real space [35]. In Fig.
4, we provide a comprehensive overview of the static structure factor over a
wide range of temperatures and parameter points. In some cases, such as rows
(g) and (h), the ground state is quasi-long-range-ordered and sharp peaks are
clearly visible at the corresponding ordering wavevectors. In other cases,
i.e. within the classical spin liquid parameter region, we observe a complex
redistribution of weight as the system passes through the three distinct
finite temperature regimes.

For the isotropic point ($J_{2}=J_{3}=1.0$), the $S(\mathbf{q})$ at the
highest temperature shown in Fig. 4(b) only displays a broad diffuse profile
corresponding to a weakly correlated thermal paramagnet. Inside the
cooperative paramagnetic regime, at $T=0.1$ and $0.05$, pinch points (with a
finite-width set by the inverse correlation length) appear between the square
and lobed shaped regions of stronger relative intensity. Their presence
signifies the approximate fulfilment of the $\mathbf{S}_{\triangle}=0$
constraint of Eq. 3 and thus the onset of strong correlations between spins
within triangular plaquettes. Upon cooling further, at $T=0.03$ and $0.02$,
the intensity is redistributed to the centers of squares and lobes located at
$\mathbf{q}=(4\pi,0)$ and $\mathbf{q}=(2\pi,2\pi)$ (and symmetry related
points), respectively. Since one is still in the cooperative paramagnetic
regime characterized by a dipolar $\sim 1/r^{2}$ decay of spin-spin
correlations (at distances smaller compared to the correlation length), these
high intensity features do not correspond to Bragg peaks. Finally, once order-
by-disorder kicks in at $T\lesssim 0.01$ selecting coplanar states, we notice
the disappearance of spectral weight at the location of the pinch point as
well as the absence of narrow necks connecting the squares with the lobes. The
presence of well-defined maxima at the aforementioned points indicates
enhanced correlations of the $120^{\circ}$ $\mathbf{q}=\mathbf{0}$ type order
[36, 13], in contrast to the conventional kagome antiferromagnet which favors
$\sqrt{3}\times\sqrt{3}$ correlations [37, 38, 39]. While the kagome
antiferromagnet develops long-range dipolar magnetic order of the
$\sqrt{3}\times\sqrt{3}$ type in the limit $T\to 0$ [39], it remains to be
established on the square-kagome lattice whether true long-range
$\mathbf{q}=\mathbf{0}$ dipolar ordering of $120^{\circ}$ type asymptotically
develops as $T\to 0$.

Away from the isotropic point, but still on the line $J_{2}=J_{3}$, we first
note that the cooperative paramagnetic and coplanar temperature regimes are
pushed down to smaller $T$ and shrink in extent [see Fig. 3]. At
$J_{2}=J_{3}=0.5$, within the cooperative paramagnetic regime ($0.008\lesssim
T\lesssim 0.05$), the pinch points are seen to be present, and with the
principal spectral weight at the center of the lobes $\mathbf{q}=(2\pi,2\pi)$
(and symmetry related points) which progressively increases, together with a
relatively weaker signal at the centers of the squares $\mathbf{q}=(4\pi,0)$
(and symmetry related points) [see Fig. 4(a)]. However, upon entering the
coplanar regime, at $T=0.003$, an equally strong maximum develops at the
$\mathbf{q}=(4\pi,0)$ type points, but in contrast to the isotropic point,
these are not indicative of enhanced $\mathbf{q}=\mathbf{0}$ correlations of
the $120^{\circ}$ type. Similar observations hold true at $J_{2}=J_{3}=1.5$
[see Fig. 4(c)] with the noticeable difference being the presence of a finite
spectral weight at the Brillouin zone centre, being more pronounced in the
intermediate temperature, i.e., cooperative paramagnetic regime. Moving away
from the symmetric line, i.e., $J_{2}\neq J_{3}$ but still inside the
degenerate manifold region, one observes that the $S(\mathbf{q})$ are only
rotationally invariant possibly reflective of the underlying symmetries of the
incipient magnetic order in the limit $T\to 0$. Upon cooling, a progressive
redistribution of spectral weight occurs leading to the appearance of new soft
maxima in the coplanar regime at $(J_{2},J_{3})=(1.0,2.5)$ [see Fig. 4(d)]
while at $(J_{2},J_{3})=(1.7,0.8)$ [see Fig. 4(e)] interestingly a similar
intensity distribution prevails across all temperatures. Finally, inside the
magnetically ordered regions, the $S(\mathbf{q})$ become more sharply peaked,
as expected, and interestingly the $S(\mathbf{q})$ at
$(J_{2},J_{3})=(2.0,2.0)$ [see Fig. 4(f)] and $(J_{2},J_{3})=(2.5,2.5)$ [see
Fig. 4(g)] resemble those inside the disordered regime at
$(J_{2},J_{3})=(1.5,1.5)$ seen in Fig. 4(c) along the $J_{2}=J_{3}$ axis as if
pre-empting the UUD order which onsets for $J_{2}=J_{3}\geq 2.0$. Inside the
Néel phase at $(J_{2},J_{3})=(0.25,2.5)$, in Fig. 4(h) there is no spectral
weight at the centre of the Brillouin zone as expected, and instead we find
dominant peaks at the $(2\pi,2\pi)$ (and symmetry related) points, and
subdominant peaks at $(2\pi,0)$ (and symmetry related) points.

Figure 5: Phase diagram in the presence of cross-octagonal couplings. At the
center we show the phase diagram with ten different phases (labeled I-X) as a
function of $J_{+}$ and $J_{\times}$ for a nearest neighbor model with
isotropic, antiferromagnetic interactions, i.e. $J_{1}=J_{2}=J_{3}=1$. (A
companion phase diagram for ferromagnetic nearest neighbor couplings is shown
in Fig. 15 below.) Besides indicating the phase boundaries (solid, dashed and
dotted lines), the phases are described by symmetry (with $D_{n}$ referring to
the dihedral group of order $n$ and $O_{h}$ to full octahedral symmetry),
coplanarity, and magnetization of their ground states. As a quick
visualization of their distinct nature, we show common origin plots and spin
structure factors for each phase left and right around the phase diagram. The
only coplanar orders are found in form of $120^{\circ}$ order (VI) in the
lower left quadrant, i.e. ferromagnetic $J_{+}$ and $J_{\times}$ (darker gray
phase), as well as a distorted version of this (VII) for sligthly positive
$J_{\times}$ (lighter gray phase). All other phases are non-coplanar, which
besides a rigid phase with cuboctahedral order (I) in the upper right quadrant
and a distorted version of this (II) includes a variety of different spirals
(III-V and VIII-X). Phase boundaries of special interest are found upon
exiting the coplanar $120^{\circ}$ order, indicated by the dashed and dotted
lines, along which the low temperature specific heat $c_{V}(T\rightarrow 0)$
takes values of $11/12$ and $23/24$, respectively (cf. Fig. 6).

##  III Octagon-plaquette Interactions

We now augment the nearest-neighbor model, discussed in the previous Section,
with the cross octagon-plaquette interactions $J_{+}$ and $J_{\times}$ (cf.
Fig. 1), which microscopically arise upon the inclusion of longer-range
Heisenberg couplings. Conceptually, they are interesting as they are expected
to stabilize non-coplanar magnetic orders, akin to the cross-hexagonal
interactions in the conventional kagome case, and in distinction to the
square-diagonal couplings.

Indeed, we find that the cross octagon-plaquette interactions induce a
plethora of non-coplanar orders as summarized in the global phase diagram of
Fig. 5. One can, in fact, distinguish ten different phases, indicated by the
different colors in the phase diagram, as one varies the relative coupling
strength of the two couplings, $J_{+}$ and $J_{\times}$, starting from an
isotropic nearest-neighbor model (i.e. for fixed $J_{1}=J_{2}=J_{3}=1$). The
distinct nature of these phases can be easily visualized by representative
common origin plots for each phase, i.e. the collapse of an extended ground-
state real-space configuration of spins to a single unit sphere by placing all
spins at a joint origin. These are shown around the phase diagram. In
addition, we show their respective spin structure factors. Let us briefly go
through these phases here, before providing a much more detailed description
in the remainder of this Section. The only phases with coplanar order come in
the form of a $120^{\circ}$ ordered phase (VI) in the lower left quadrant,
i.e. for ferromagnetic $J_{+}$ and $J_{\times}$ (indicated by dark gray in
Fig. 5), as well as a distorted version of this $120^{\circ}$ order (VII) for
slightly positive $J_{\times}$ (indicated by light gray in the phase diagram).
In the plain-vanilla $120^{\circ}$ ordered phase spins on elementary triangles
form mutual angles of $2\pi/3$, while this angle is increased beyond $2\pi/3$
in the distorted phase as discussed in Sec. III.2. All other phases of the
phase diagram exhibit non-coplanar order, which come in commensurate and
incommensurate forms. The simpler, commensurate variant of such non-coplanar
order is the extended phase with cuboctahedral order (I) in the upper right
quadrant and a distorted version of this (II), which we discuss in depth in
Sec. III.3. Somewhat more complex non-coplanar orders come in incommensurate
spin spiral order, which we find for the remaining six phases (III-V and
VIII-X). Remarkably, however, these can still be described by a semi-
analytical ground-state construction, which we describe in Sec. III.4.

###  III.1 Physics on the axes

But before we dulge into the various magnetic orders that can be stabilized by
the combined effects of cross octagon-plaquette interactions, we first
consider their exclusive effect, i.e. we consider the horizontal and vertical
axes in the middle of the phase diagram of Fig. 5. Adding only one of the two
cross octagon-plaquette interactions ($J_{\times}$ or $J_{+}$ only), it turns
out that it is still possible to locally satisfy the constraint of Eq. (3)
provided that the added interaction is ferromagnetic in nature.

For FM $J_{+}$, the spins located on the bow-ties become locked together and
are all ferromagnetically aligned, meaning that $\mathbf{S}_{k}$ in the
constraint (3) is fixed in each and every triangle to be the same. Denoting
this fixed direction by $\mathbf{M}$, the remaining two spins (located on the
squares) are thus subject to the local constraint on each triangle

$\mathbf{S}_{i}+\mathbf{S}_{j}=-\mathbf{M}.$ (6)

This results in a global spin configuration in which $1/3$ of the spins point
along $\mathbf{M}$ while the other $2/3$ point along a ring at an angle of
$2\pi/3$ away from $\mathbf{M}$. This removes the possibility of having an
extensive number of both globally coplanar and globally non-coplanar ground
states, and thus there are only two distinct finite temperature regimes. In
other words, entropic-driven selection of globally coplanar configurations
would result in a regular $120^{\circ}$ ordered state. At high temperatures
there is the usual paramagnetic regime, and at low temperatures a crossover
into the ground state manifold with an accompanying specific heat
$c_{V}\rightarrow 11/12$ as $T\rightarrow 0$, due to the continued existence
of one zero mode per triangle. This is precisely what we find in finite-
temperature Monte Carlo simulations as shown in Fig. 6.

Figure 6: Specific heat and ground states of the extended model on the axes.
For either only $J_{+}$ or $J_{\times}$ turned on and being ferromagnetic, one
again finds special values of $c_{V}(T\rightarrow 0)$ that differ from 1. For
For $J_{\times}=0.0,J_{+}<0.0$, the value is 11/12 (cf. purple curve with
$J_{+}=-1.0$), and for $J_{+}=0.0,J_{\times}<0.0$, we find a value of 23/24
(cf. mint curve with $J_{\times}=-1.0$). The upper insets show common origin
plots of the corresponding ground states. Figure 7: Structure factors of the
extended model on the axes. Structure factors for different temperatures
between T = 0.003 and T = 0.1 (left to right) on the axes
($J_{+}=-1.0,J_{\times}=0.0$ on the top and $J_{+}=0.0,J_{\times}=-1.0$ on the
bottom). Note the sharp maxima at $(-4\pi,0)$ and symmetry related momenta in
the top row (indicated by red circle). The shown squares extent in reciprocal
space from $-4\pi$ to $4\pi$ in both dimensions. For $T=0.1$, the first
Brillouin zone and the extended Brillouin zone are indicated.

For FM $J_{\times}$, the spins within each octagon coupled by $J_{\times}$
become locally ferromagnetically aligned, but the spins from one octagon to
the next are not. This locks neighboring triangles together, resulting in one
zero mode per unit cell, as opposed to one zero mode per triangle, and thus a
low-temperature specific heat
$c_{V}\rightarrow[11\cdot(1/2)+1\cdot(1/4)]/6=23/24$ as $T\rightarrow 0$.
This, again, is in perfect agreement with finite-temperature Monte Carlo
simulations as shown in Fig. 6.

In comparison to the spin structure factors, summarized for the nearest-
neighbor model in Fig. 4, the spin-spin correlations discussed above lead to a
deformation of $S({\bf q})$ as shown in Fig. 7. The $S(\mathbf{q})$ for
($J_{+}=-1.0,J_{\times}=0.0$) show maxima at positions where the Bragg peaks
of the incipient $\mathbf{q}=\mathbf{0}$ order of the $120^{\circ}$ type would
show up as $T\to 0$ [13]. In contrast, for ($J_{+}=0.0,J_{\times}=-1.0$), the
$S(\mathbf{q})$ display maxima at the expected locations for
$\sqrt{3}\times\sqrt{3}$ order [13].

###  III.2 Coplanar $\bf 120^{\circ}$ order

Turning to the magnetically ordered states of our phase diagram in Fig. 5, we
start with the coplanar $120^{\circ}$ order found in the lower left quadrant
where both cross octagaon-plaquette couplings are ferromagnetic. For the
square-kagome lattice geometry at hand, one can, in principle, distinguish two
different types of $120^{\circ}$ orders, each with three sublattices where the
spins on each sublattice point to a different corner of an equilateral
triangle such that each neighboring pair of spins forms an angle of $2\pi/3$.
As illustrated in Fig. 8, the magnetic unit cell contains either 12 sites and
is two times larger than the geometric unit cell for the type shown in Fig.
8(a), or 24 sites for the type in Fig. 8(d), with the difference between the
two types coming down to the bow-tie spins. While these bow-tie spins are all
pointing in the same direction in the order of Fig. 8(a) forming a bow-tie
ferromagnet, the bow-tie spins in the second type of $120^{\circ}$ order form
stripy antiferromagnets with alternating rows/columns of up and down-pointing
spins as shown in Fig. 8(d). A schematic of the static spin structure factors
corresponding to these two types of real-space $120^{\circ}$ order is shown
for an extended Brillouin zone of the square-kagome lattice in Fig. 8(c) and
(f), respectively. Both spin configurations exhibit the symmetry of the
dihedral group $D_{3}$. Their net magnetization vanishes, $m=0$, and the
energy per site can be calculated as

$E_{\text{120}^{\circ}}=-1+\frac{1}{3}(J_{+}+J_{\times})\,,$ (7)

for varying strengths of the two couplings.

Figure 8: $\bf 120^{\circ}$ orders. (a)+(d) Real space arrangement of spins in
two different coplanar $120^{\circ}$ ordered states. There are each three
sublattices (each one corresponding to one color), with a 12-site magnetic
unit cell (large gray square), which is two times larger than the geometric
unit cell (small square), and a 24-site magnetic unit cell, respectively. The
magnetization of the bow-tie spins, $m_{\text{bow-tie}}$, is $1.0$ for
arrangement (a) and $0.5$ for arrangement (d). (b)+(e) Each sublattice of
spins points towards a different corner of an equilateral triangle such that
each neighboring pair of spins forms an angle of $2\pi/3$. (c)+(f) First and
extended Brillouin zones of the square-kagome lattice showing the positions of
the corresponding dominant and subdominant Bragg peaks. In (f), the ratio of
the weight of the subdominant peaks $\lambda$ to the weight of the dominant
peaks $\Lambda$ is $\lambda/\Lambda\approx 85\%$.

To elucidate the thermodynamics associated with these $120^{\circ}$ orders we
show, in Fig. 9, specific heat traces for the point in the lower left corner,
$J_{+},J_{\times}=-1$ (i.e. deep in the phase), for different system sizes
between $L=8$ and $L=32$. Next to a sharp peak at $T=0.27(3)$, there is a
second, more subtle feature slightly at $T=0.33(2)$. This smaller bump can, in
fact, be associated with the build-up of quasi-long-range $120^{\circ}$ order
– but without selecting one of the two types of $120^{\circ}$ order.
Accordingly, the structure factor for the regime between the two features in
the specific heat $c_{V}$, e.g. for $T=0.3$, can be obtained by averaging the
real space correlations of the two possible $120^{\circ}$ order arrangements.
Also the net magnetization of the bow-tie spins, shown in the lower panel of
Fig. 9, builds up at this feature at higher temperatures. The lower
temperature peak at $T=0.27(3)$ then corresponds to the spontaneous selection
of one of the two types $120^{\circ}$ orders, resulting in a sharp feature in
the specific heat. Notably, both type of orders appear with the same
probability, which is signified in the bow-tie magnetization. For small enough
temperatures and large enough system system sizes a bimodal distribution
around $m_{\text{bow-tie}}=1$ and $m_{\text{bow-tie}}=0.5$ is sampled, whereas
for higher temperatures and smaller system sizes the Monte Carlo sampling is
able to overcome the free energy barrier between these two states and samples
them in an ergodic fashion, resulting in a net bow-tie magnetization of
$m_{\text{bow-tie}}=0.75$.

Figure 9: Thermodynamics of $\bf 120^{\circ}$ order. (top) The specific heat
of the $120^{\circ}$ phase shows a sharp feature at $T=0.27(3)$ and a subtle
bump slightly at $T=0.33(2)$. The subtle bump can be associated with the
build-up of $120^{\circ}$ order, whereas at the sharp feature a specific
$120^{\circ}$ order is selected – in the example shown the structure factor at
$T=0.1$ (inset) conincides with the analytical structure factor in Fig. 8 (f),
whereas the intermediate structure factor can be obtained by averaging the
real space correlations of all possible $120^{\circ}$ orders, i.e. both shown
in Fig. 8 plus all possible rearrangements with the same order. The bottom
panel shows the magnetization of the bow-tie spins. The displayed bow-tie
magnetization starts to build up right at the bump-like feature slightly above
$T=0.3$. Below $T=0.3$, at the sharp peak in the specific heat, the bow-tie
magnetization splits up into different branches that converge to the values of
$m_{\text{bow-tie}}=0.5$, and $m_{\text{bow-tie}}=1.0$, respectively.

A deformed version of the $120^{\circ}$ order, henceforth termed
$120^{\circ}$-d, exists in a small region touching the conventional
$120^{\circ}$-order phase in the upper left quadrant of the phase diagram,
indicated as phase VII in Fig. 5. Its deformation – two of the three mutual
angles between neighboring spins take a value of $\alpha>2\pi/3$, while the
third one becomes smaller than $2\pi/3$ – can be seen in the common origin
plot of Fig. 5. It can derived by elementary geometric considerations: For
instance, let the angle between the blue and the red sublattices and between
the blue and the green sublattices shown in Fig. 8 increase to
$\alpha>2\pi/3$, while the third angle, between the red and the green
sublattices, becomes $2\pi-2\alpha<2\pi/3$. The angle $\alpha$ is found to
vary as $\cos\alpha=-(2+J_{\times})^{-1}$ and the ground-state energy per site
takes the value

$E_{\text{120}^{\circ}\text{-d}}=-\frac{6+J_{\times}^{2}-2J_{+}-J_{\times}(J_{+}-4)}{3(J_{\times}+2)}\,.$
(8)

The ground state symmetry is reduced to the symmetry of the dihedral symmetry
group $D_{1}$ and has a non-zero, yet small magnetization
$m_{120^{\circ}\text{-d}}=\tfrac{J_{\times}}{3(2+J_{\times})}$, which only
depends on the values of $J_{\times}$ (see also Fig. 20 in Appendix C).

###  III.3 Cuboc order

Let us now turn to the first instance of non-coplanar order and consider the
cuboctahedral (cuboc) order found in the upper right quadrant of the phase
diagram, where both cross-octagon plaquette couplings are antiferromagnetic.
In real-space order this order is described by 12 sublattices where the spins
on each sublattice point towards a different corner of a cuboctahedron, as
illustrated in Fig. 10(b). The underlying magnetic unit cell contains 24 sites
and is thus four times larger than the geometric unit cell, see Fig. 10(a).
All neighboring spins form an angle of $2\pi/3$, which corresponds to the
cuboc1 state, discussed for the kagome lattice with cross-hexagonal couplings
[22, 1]. Note that there are eight possible ways to arrange cuboc1 order on
the square-kagome lattice as illustrated in Fig. 17 of Appendix A, each of
which breaks $C_{4}$ rotation symmetry. The remaining symmetry of this ordered
state is the full octahedral symmetry group $O_{h}$. It has zero
magnetization, $m=0$, and its energy per site is given by

$E_{\text{cuboc}}=-1-\frac{1}{3}(J_{+}+J_{\times})\,.$ (9)

Figure 10: Cuboctahedral order cuboc1. (a) Real space arrangement of spins in
a cuboctahedral (cuboc) ordered state. There are 12 sublattices (each one
corresponding to one color) with a 24-site magnetic unit cell (large gray
square), which is four times larger than the geometric unit cell (small
square). (b) Each sublattice of spins points towards a different corner of a
cuboctahedron such that each neighboring pair of spins forms an angle of
$2\pi/3$. This order corresponds to the cuboc1 state in [1]. (c) First and
extended Brillouin zones of the square-kagome lattice showing the positions
and the fractions of total spectral weight of the corresponding Bragg peaks.
The order breaks $C_{4}$ symmetry.

The real-space correlations of the cuboc1 order give rise to a spin structure
factor as schematically visualized for one of the eight possible real-space
configurations in Fig. 10 (c), which shows the positions of the corresponding
Bragg peaks within the extended Brillouin zone of the square-kagome lattice
together with the associated fraction of total spectral weight. Averaging over
the real-space correlations of all eight possible cuboc1 arrangements,
restores $C_{4}$ symmetry in the structure factor (cf. Appendix A).

Figure 11: Thermodynamics of cuboc1 phase. The specific heat (top panel) of
the cuboc1 phase displays a double-peak structure. The high-temperature peak
at $T=0.371$ can be associated with the initial build-up of coexisting cuboc1
orders. At the low-temperature peak (at $T=0.256$), one specific realization
of cuboc1 order is selected (cf. Appendix A). This can be seen from the
structure factors (inset). While the structure factor at $T=0.1$ corresponds
to one specific realization of cuboc1 order, the intermediate structure factor
at $T=0.3$ conincides with the analytical structure factor of the
superposition of all possible cuboc1 order realizations (see also Fig. 17).
Middle and bottom panels show the cuboc1 order parameter $\mathcal{O}$ (Eq.
17), and the corresponding susceptibility $\chi_{\mathcal{O}}$ (Eq. 18). The
lower inset shows how $\mathcal{O}_{ijkl}$ is being calculated on a single
skew bow-tie. This quantity is then averaged over all skew bow-ties in order
to calculate $\mathcal{O}$.

Turning to the thermodynamics of the cuboc phase, we show, in Fig. 11, the
specific heat, the cuboc1 order parameter $\mathcal{O}$, and its associated
susceptibility $\chi_{\mathcal{O}}$ (both introduced in Appendix A) for
different system sizes between $L=8$ and $L=32$. The specific heat displays a
clearly visible double-peak structure which, similar to our discussion of the
coplanar $120^{\circ}$ order, can be rationalized by the coexistence of
multiple possible cuboc1 arrangements on the square-kagome lattice: At the
high-temperature peak in $c_{V}$ (at $T=0.371$), the system builds up cuboc1
order, but does not select a specific arrangement out of the eight possible
realization, as can be seen from the $C_{4}$ symmetric structure factor in the
inset. At the low-temperature peak (at $T=0.256$) a specific cuboc1 order is
then spontaneously selected. At this temperature, the order parameter
$\mathcal{O}$ (Eq. 17), which takes values of $\pm 1$ for different specific
cuboc1 arrangements, builds up and the corresponding susceptibility
$\chi_{\mathcal{O}}$ (Eq. 18) diverges as shown in the lowest panel of Fig.
11.

A deformed version of cuboc1 order (denoted as phase II in our phase diagram),
termed cuboc-d, extends to a part of the lower right quadrant in the phase
diagram Fig. 5 with $J_{\times}<0$ and $J_{+}\geq 0$, which can be derived by
applying the generalized Luttinger-Tisza method of Ref. [40], see also the
Mathematica files in the Supplement [41]. In this cuboc-d phase, the antipodal
squares of the cuboctahedron are deformed into non-coplanar closed polygon
chains with equal edges, while the square in the equatorial plane is deformed
into a coplanar rectangle. Its symmetry is thereby reduced to the symmetry of
the dihedral symmetry group $D_{4}^{s}$ where the superscript $s$ denotes
additional mirror symmetry on the $xy$-plane, i.e. $s=\text{diag}(1,1,-1)$.
The ground state energy per site can be calculated to be

$E_{\text{cuboc-d}}=\frac{3+J_{\times}^{2}+J_{+}-J_{\times}(J_{+}+3)}{3(J_{\times}-1)}.$
(10)

The cuboc-d ground state still has zero magnetization, $m=0$.

###  III.4 Spiral orders

The upper left and lower right quadrant of our phase diagram are occupied by
spin spiral phases, coming in the form of six different variants (labeled
III-V and VIII-X, respectively). The complexity of these incommensurate, non-
coplanar orders becomes immediately clear when looking at their common origin
plots, whose intricate patterns point to magnetic unit cells of hundreds of
spins. This renders any direct analytical description of these phases rather
elusive, but it turns out that one can, in fact, deduce a semi-analytical
description of these phases from low-temperature numerical simulation data. As
we will discuss below, this approach provides us with a symmetry-optimized
description of these spin spirals including explicit expressions of their
ground-state energy as function of the coupling parameters. The latter then
allows us to establish sharp phase boundaries between these complex spin
spiral phases as depicted in the phase diagram of Fig. 5.

#### Semi-analytical approach

The starting point of our semi-analytical approach is numerical data in the
form of a common origin plot of the $N$ spin vectors of a ground-state spin
configuration sampled in Monte Carlo simulations at ultra-low temperatures
$T=10^{-4}$, typically explored in conjunction with a parallel tempering
scheme. Example input for the spin spiral phase III in the lower right
quadrant is shown on the left in the schematic illustration of Fig. 12.

In a second step, we then identify a smaller number of $M<N$ unique spin
vectors by grouping spins in the initial common origin plot that point
approximately in the same direction, see the middle panel of Fig. 12. In
practice, we say that two spins $\mathbf{S}_{i}$ and $\mathbf{S}_{j}$ point
approximately in the same direction if
$\mathbf{S}_{i}\cdot\mathbf{S}_{j}\geq\gamma$, where the exact value of
$\gamma$ slightly varies from case to case, but typically $\gamma\approx
0.995$.

Figure 12: Semi-analytical scheme for spin spiral phases. Starting point is a common origin plot of the $N$ spin vectors of a ground-state spin configuration sampled in Monte Carlo simulations at ultra-low temperatures $T=10^{-4}$, as shown on the left. In a second step, we identify $M<N$ unique spin directions by grouping spins that point approximately in the same direction (middle). From these unique spin vectors, we then identify symmetries that further reduce the number of unique spin vectors to $K<M$ and allow us to describe the spiral phase analytically (right). The data shown is for spin spiral phase III of the lower right quadrant with couplings $J_{+}=+1$ and $J_{\times}=-0.2$. The initial common origin plot has $N=864$ spins corresponding to a system size of $L=12$. The initial reduction leads to $M=72$ points on five circles (as indicated in the middle panel). The symmetry group $G$ of the example is generated by rotations of $\pi/6$ about the vertical symmetry axis as well as by a reflection at the equatorial plane, i.e. $G=D_{12}^{s}$, and therefore $K=4$. spiral phase | symmetry | $N$ | $M$ | $K$ | # parameters
---|---|---|---|---|---
III | $D_{12}^{s}$ | 864 | 72 | 4 | 3
IV | $D_{6}^{s,t}$ | 864 | 228 | 13 | 19
V | $\\{{\rm id,s}\\}$ | 864 | 48 | 24 | 47
VIII | $D_{12}$ | 864 | 37 | 3 | 3
IX | $D_{3}^{\sigma}$ | 864 | 108 | 18 | 26
X | $D_{6}^{s}$ | 864 | 32 | 4 | 3
Table 1: Semi-analytical approach. Compression of the parametrization of the
spin spirals via clustering and symmetrization (Fig. 12), starting with a
common origin plot with $N=864$ spins sampled at $T=10^{4}$ for a linear
system size $L=12$.

For these $M$ spin directions we then identify all possible symmetries, which
allows us to further reduce the number of unique spin vectors to $K<M$. A
symmetry in the aforementioned sense is a tuple $(R,\pi)$ with $R\in O(3)$ and
$\pi\in S_{M}$ where $S_{M}$ is the permutation group of $M$ elements, such
that $R\mathbf{S}_{i}=\mathbf{S}_{\pi(i)}$ for all spin vectors
$\mathbf{S}_{i}$. From the remaining $K$ spin vectors, all spin vectors can be
generated by applying these symmetries.

In total, this approach allows us to describe the ground state by at most
$2K-1$ parameters – maximally two parameters per spin minus one parameter due
to a global rotation around the symmetry axis, but less if some polar or
azimuthal angles of the ground state assume fixed values. This compact
representation is summarized in Table 1 for all six spin spiral phases of our
phase diagram. Having such an analytical representation at hand, we can then
explicitly calculate various observables such as the magnetization or ground-
state energy for arbitrary couplings $J_{+}$ and $J_{\times}$ for all phases,
which in turn allows us to determine the phase boundaries shown in the phase
diagram Fig. 5. Several cross-checks can be used to validate this approach,
including a comparison of the analytical ground-state energy and the Monte
Carlo result as well as the determination of the phase boundaries, which we
compare to scans of derivatives of the Monte Carlo energy, as shown in
Appendix C. In general, we find excellent agreement.

#### Example: Spiral phase III with $D_{12}^{s}$ symmetry

Let us illustrate this semi-classical approach and its validity for an
explicit example, picking the spiral phase III in the lower right quadrant of
the phase diagram Fig. 5 for $J_{+}=+1$ and $J_{\times}=-0.2$. Our schematic
illustration of the semi-classical approach in Fig. 12 also uses this example.
The common origin plot on the left consists $N=864$ spin vectors of the
numerical ground state (at $T=10^{-4}$) of a system of linear length $L=12$.
By grouping spins that point to the same direction, using the criterion
$\mathbf{S}_{i}\cdot\mathbf{S}_{j}\geq 0.999$, we find that there are only
$M=72$ unique spin directions which are shown in the middle panel. Performing
a symmetry analysis, one finds that the symmetry group $G$ of this state is
generated by rotations of $\pi/6$ about the vertical symmetry axis as well as
by a reflection at the equatorial plane, i.e. $G=D_{12}^{s}$. This leaves us
with just $K=4$ representative spins that describe the entire spin spiral
configuration – a significant reduction compared to the $N=864$ spins in the
original real-space configuration. The $K=4$ representative spins can be
written as functions of three parameters $\alpha,z_{1},z_{2}$ in the following
way

$\displaystyle{\mathbf{S}}_{1}=$
$\displaystyle\left(0,\sqrt{1-z_{1}^{2}},z_{1}\right)\;,$
$\displaystyle{\mathbf{S}}_{2}=$
$\displaystyle\left(\frac{\sqrt{1-z_{2}^{2}}}{\sqrt{2}},\frac{\sqrt{1-z_{2}^{2}}}{\sqrt{2}},z_{2}\right)\;,$
$\displaystyle{\mathbf{S}}_{3,4}=$
$\displaystyle\left(\cos{\left(\frac{\pi}{4}\pm\alpha\right)},\sin{\left(\frac{\pi}{4}\pm\alpha\right)},0\right)\;.$
(11)

With this compact representation at hand, the energy per site can now be
explicitly calculated as

$\displaystyle E_{\text{III}}=$
$\displaystyle\frac{1}{12}\bigg{[}J_{\times}\cos
2\alpha+\sqrt{3}J_{\times}\sin 2\alpha$
$\displaystyle-2\Big{(}\sqrt{3}J_{+}-(\sqrt{3}-2)J_{+}z_{1}^{2}+4z_{1}z_{2}$
$\displaystyle+(\sqrt{6}-\sqrt{2})\sqrt{(1-z_{1}^{2})(1-z_{2}^{2})}+J_{\times}(2z_{2}^{2}-1)$
$\displaystyle+2(\sqrt{3-3z_{2}^{2}}-\sqrt{2-2z_{1}^{2}})\sin\alpha$
$\displaystyle+2(\sqrt{2-2z_{1}^{2}}+\sqrt{1-z_{2}^{2}})\cos\alpha\Big{)}\bigg{]}\,.$
(12)

Numerical minimization of this energy with $J_{+}=+1$ and $J_{\times}=-0.2$
then leads to $E_{\text{III},\text{semi-analytical}}=-1.34515$ which is in
excellent agreement with and (as expected) slightly below the Monte Carlo
result $E_{\text{III},\text{MC}}=-1.345$, which is shifted upwards by finite-
temperature fluctuations commensurate with a temperature of $T=10^{-4}$.

We present explicit results of this semi-classical approach for all other spin
spiral phases in Appendix B.

Figure 13: Cuboctahedral order cuboc3. This variant of a non-coplanar
cuboctahedral order is found when flipping the bow-tie interactions $J_{2}$
and $J_{3}$ of the nearest neighbor model to ferromagnetic (while keeping the
square interactions $J_{1}$ antiferromagnetic). It is obtained from the
original cuboc order (Fig. 10) via a local spin transformation (see main
text). (a) Real space arrangement of spins in the cuboc3 ordered state. There
are 12 sublattices (each one corresponding to one color) with a 24-site
magnetic unit cell (large gray square), four times larger than the geometric
unit cell (small square). (b) Each sublattice of spins points towards a
different corner of a cuboctahedron such that each neighboring pair of spins
forms an angle of $2\pi/3$ on the squares and an angle of $\pi/3$ on the
triangles, making it different from the cuboc1 and cuboc2 orders discussed in
the literature [1]. (c) First and extended Brillouin zones of the square-
kagome lattice showing the positions of the corresponding Bragg peaks. The
order breaks $C_{4}$ symmetry. Figure 14: Pentagonal order. Coplanar
pentagonal order arises when flipping the bow-tie interactions $J_{2}$ and
$J_{3}$ of the nearest neighbor model to ferromagnetic (while keeping the
square interactions $J_{1}$ antiferromagnetic). It is obtained from the
original $120^{\circ}$ order (Fig. 8) via a local spin transformation (see
main text). (a) Real space arrangement of spins in the coplanar pentagonal
state. There are five sublattices (each one corresponding to one color) with a
24-site magnetic unit cell (large gray square), four times larger than the
geometric unit cell (small square). (b) The spins on the five sublattices
point to five out of the six corners of a hexagon. (c) First and extended
Brillouin zones of the square-kagome lattice showing the positions of the
corresponding dominant a nd subdominant Bragg peaks. The ratio of the weight
of the subdominant peaks $\lambda$ to the weight of the dominant peaks
$\Lambda$ is $\lambda/\Lambda\approx 85\%$.

###  III.5 Cuboc3 and pentagonal order for mixed interactions

Finally, we note that one could also consider variations of the model at hand
where one changes the sign of the interactions in the original nearest-
neighbor $(J_{1},J_{2},J_{3})$ Heisenberg model (see Fig. 1 for a reminder of
the coupling geometries). Flipping the sign of the bow-tie interactions to
ferromagnetic couplings, i.e. $J_{2}=J_{3}=-1$, while keeping the square
interactions antiferromagnetic, i.e. $J_{1}=1$, yields exactly the same phase
diagram as shown in Fig. 5, up to local spin transformations. Specifically,
since the spins on the bow-ties are coupled via $J_{2}$ and $J_{3}$ to their
nearest neighbors and via $J_{+}$ to other bow-tie spins, the total energy
remains unchanged if all bow-tie spins are inverted simultaneously while the
sign of the triangular couplings is changed, $J_{2}\rightarrow-J_{2}$ and
$J_{3}\rightarrow-J_{3}$.

Performing these local spin transformations on the orders of our original
phase diagram yields two qualitatively new types of orders: The original
cuboctahedral and $120^{\circ}$ orders (Figs. 10 and 8) turn into a non-
coplanar cuboc3 order (Fig. 13) and a coplanar pentagonal order (Fig. 14),
respectively. In the new cuboctahedral order cuboc3, each neighboring pair of
spins forms an angle of $2\pi/3$ on the squares and an angle of $\pi/3$ on the
triangles, which is different to cuboc1 order (where all neighboring pairs of
spins form an angle of $2\pi/3$) and to cuboc2 order (where all neighboring
pairs of spins form an angle of $\pi/3$) [1].

###  III.6 Octagonal and conical orders for FM interactions

Figure 15: FM phase diagram. For fixed $J_{1}=J_{2}=J_{3}=-1$, the model shows
a large ferromagnetic phase (light gray), a coplanar phase with octagonal
order (dark gray) and non-coplanar umbrella-like double cone states
(turquoise) that are an interpolation between octagonal and ferromagnetic
states. Ferromagnetic and coplanar phases touch at the single point
$J_{+}=J_{\times}=1-\tfrac{1}{\sqrt{2}}$. On the right hand side are
corresponding common origin plots and structure factors from Monte Carlo
simulations. Figure 16: Octagonal and double-conical orders. (a) Real space
arrangement of spins in the coplanar octagonal ordered state stabilized in the
purely ferromagnetic nearest-neighbor model. There are eight sublattices with
a 24-site magnetic unit cell (large gray square), four times larger than the
geometric unit cell (small square). (b) Each sublattice of spins points
towards a different corner of an octagon such that each pair of neighboring
spins forms an angle of $\pi/4$. (c) First and extended Brillouin zones
showing the positions and the fractions of total spectral weight of the
corresponding Bragg peaks. (d) Real space arrangement of spins in the non-
coplanar double cone state. The eight sublattices coincide with those of the
octagonal state shown in (b). (e) The directions in which the spins on the
sublattices point form two cones. The first cone is formed by the four
sublattices of spins located on the squares with mutual angles of $\pi/2$,
whereas the second cone is formed by the four sublattices of bow-tie spins,
again with mutual angles of $\pi/2$ but rotated $\pi/4$ with respect to the
first cone. (f) First and extended Brillouin zones showing the positions of
the corresponding Bragg peaks. The ratio of the weight of the subdominant
peaks $\lambda$ to the weight of the dominant peaks $\Lambda$ depends on the
two polar angles of the cones and therefore on $J_{+}$ and $J_{\times}$.

Flipping the sign of all three couplings to ferromagnetic in the underlying
nearest neighbor Heisenberg model, i.e. $J_{1}=J_{2}=J_{3}=-1$, the phase
diagram for varying cross octagon-plaquette interactions changes its topology
substantially. As depicted in Fig. 15 there are only three distinct phases.
Trivially, there is a large ferromagnetic phase when $J_{+},J_{\times}\leq 0$,
which also extents to the other three quadrants in the phase diagram. Its
ground state energy per site is given by

$E_{\text{FM}}=-2+\frac{J_{\times}+J_{+}}{3}\,.$ (13)

It is bound by the hyperbola defined by
$(J_{\times}-1)(J_{+}-1)=\tfrac{1}{2}$. Bound by a second hyperbola, given by
$(J_{\times}+\tfrac{1}{\sqrt{2}})(J_{+}+\tfrac{1}{\sqrt{2}})=1$, there is a
coplanar ordered phase with eight sublattices where the spins on each
sublattice point to the corners of an octagon. The magnetic unit cell contains
24 sites and is four times larger than the geometric unit cell (cf. Fig. 16
(d) - (f)). Its ground state has zero magnetization, $m=0$, and its energy per
site is given by

$E_{\text{octagonal}}=-\frac{1}{3}(2+2\sqrt{2}+J_{+}+J_{\times})\,.$ (14)

Between these two phases, which touch each other only at the point
$J_{+}=J_{\times}=1-\tfrac{1}{\sqrt{2}}$, there exist non-coplanar umbrella-
like states that smoothly interpolate between octagonal and FM order. The
magnetic unit cell of this phase coincides with the one of the octagonal
phase, but the spins are no longer coplanar, but rather the directions in
which the spins on the sublattices point form two cones. The first cone is
formed by the four sublattices of spins located on the squares with mutual
angles of $\pi/2$, whereas the second cone is formed by the four sublattices
of bow-tie spins, again with mutual angles of $\pi/2$ but rotated $\pi/4$ with
respect to the first cone (cf. Fig. 16 (h)). The two polar angles that
describe these two cones depend on $J_{+}$ and $J_{\times}$, the ground state
energy per site for those states can be calculated analytically and is given
by

$\displaystyle E_{\text{conical}}=$
$\displaystyle-\frac{1}{6J_{+}J_{\times}}\bigg{(}J_{+}+J_{\times}+4J_{+}J_{\times}$
$\displaystyle\pm(J_{+}-J_{\times})\sqrt{1-12J_{+}J_{\times}+4J_{+}^{2}J_{\times}^{2}}\bigg{)}\,,$
(15)

where the plus sign applies to the left domain with
$J_{+}<1-\tfrac{1}{\sqrt{2}}$, and the minus sign to the right domain with
$J_{+}>1-\tfrac{1}{\sqrt{2}}$. While the ferromagnetic phase has uniform
magnetization and the octagonal phase has zero magnetization, the
magnetization in the conical phase depends on $J_{+}$ and $J_{\times}$.
Analytically, one finds

$\displaystyle m_{\text{conical}}$ $\displaystyle=[(2(J_{\times}+4)J_{+}\pm
w)+1][24J_{+}|J_{\times}|]^{-1}$
$\displaystyle\times\sqrt{\frac{\left(2J_{\times}^{2}+1\right)w\pm\left(-4J_{\times}^{3}J_{+}+6J_{\times}^{2}+6J_{\times}J_{+}-1\right)}{w}}\,.$
(16)

Here, the plus sign applies to the right domain with
$J_{+}>1-\tfrac{1}{\sqrt{2}}$, and the minus sign to the left domain with
$J_{+}<1-\tfrac{1}{\sqrt{2}}$, and
$w=\sqrt{4J_{\times}^{2}J_{+}^{2}-3J_{\times}J_{+}+1}$.

##  IV Discussion and Outlook

We have unveiled a rich variety of magnetic textures that span the phase
diagram of an extended classical Heisenberg model on the square-kagome
lattice. Motivated by the possibility of stabilizing non-coplanar magnetic
orders, we show that a minimal set of interactions needed to realize these
involves introducing cross-plaquette interactions on top of the nearest-
neighbor Heisenberg model, with either antiferromagnetic or ferromagnetic
nearest-neighbor interactions, or a combination of both (ferromagnetic on
bowtie bonds and antiferromagnetic on square bonds). A thorough classical
Monte Carlo analysis reveals a plethora of non-coplanar states including a new
type of cuboc order (dubbed cuboc3, Fig. 13), as well as highly intricate
incommensurate non-coplanar spirals. The underlying magnetic unit cells of
these spin spiral states feature a highly complex structure and large sizes
but, remarkably, we are able to provide a semi-analytical construction of
these phases based on a symmetry-optimized parameterization. This renders
possible obtaining explicit expressions (depending only on a small number of
parameters) for their ground-state energy as a function of the coupling
strengths, which in turn enables us to establish, with high precision, the
phase boundaries between these complex spiral phases.

Besides the ground state, we also study the thermodynamics of non-coplanar
states employing classical Monte Carlo simulations. By virtue of being chiral,
non-coplanar states are expected to feature a symmetry breaking phase
transition at $T\neq 0$. In particular, for the cuboc order we present the
temperature evolution of the specific heat, the chiral order parameter, and
its susceptibility, for different system sizes which manifestly exhibits
signatures of a chiral phase transition. For the elementary model with only
three symmetry inequivalent nearest-neighbor couplings, we show that, besides
the isotropic point, within the entire region occupied by an extensively
degenerate manifold of ground states one can always identify three distinct
temperature regimes, namely, a high-temperature thermal paramagnet, an
intermediate-temperature cooperative paramagnet, and a low-temperature
coplanar regime selected via an order-by-disorder mechanism. We show that upon
introduction of ferromagnetic cross-plaquette interactions of just one type,
i.e., either $J_{+}$ or $J_{\times}$, the extensive degeneracy is reduced, but
there still persists one zero mode per per triangle or per unit cell,
respectively. This results in a $T\to 0$ limiting value of specific heat which
is less than one, however, the intermediate-temperature cooperative
paramagnetic temperature regime disappears.

Since classical non-coplanar magnetic orders are characterized by a finite
scalar spin chirality, it is plausible that, if quantum fluctuations are
turned on and are successful in restoring spin rotational symmetry, e.g., in
the extreme quantum limit of small spin $S=1/2$, the chiral symmetry breaking
still persists and carries over in the resulting non-magnetic ground state [9,
10, 11]. Such quantum melting would give rise to novel chiral paramagnetic
phases, characterized by a spontaneous breaking of time reversal and lattice
symmetries, while preserving their product. One such example is the chiral
quantum spin liquid, first elucidated by Kalmeyer and Laughlin [42], which
hosts bulk semion excitations and a chiral gapless edge mode [43].
Furthermore, in our phase diagrams, the non-coplanar orders break lattice
symmetries in such a way that simply restoring spin rotational symmetry would
not fully restore all lattice symmetries 111In other words, to use the
language of [1], the non-coplanar orders found here are not regular magnetic
orders on the square-kagome lattice. A regular magnetic order would be a
classical spin state which respects all lattice symmetries modulo global
$O(3)$ spin transformations. Thus, restoring spin rotational symmetry would
result in a fully symmetric state., potentially leading to descendent nematic
chiral liquids for small values of spin [45]. We thus provide a detailed
symmetry analysis for the myriad of non-coplanar orders, paving the way for
the systematic classification of descendant chiral spin liquid states. It
should be noted here that, unlike the kagome lattice, the square-kagome has an
even number of sites per unit cell. This leaves open the possibility of fully
trivial paramagnetic phases appearing, as the Lieb-Schultz-Mattis-Hastings-
Oshikawa theorem does not preclude a fully symmetric, yet topologically
trivial, gapped phase [46, 47, 48]. For the non-coplanar umbrella-like double
cone states with $m\neq 0$, interpolating between the ferromagnetic and
octagonal orders, the proximity to ferromagnetism opens the exciting
possibility of realizing in the corresponding quantum model the elusive spin
nematic orders [49, 50, 51]. The spin rotational symmetry breaking in these
quadrupolar ordered states is described by a time-reversal invariant order
parameter, given by a symmetric traceless rank-2 tensor in contrast to a
conventional dipolar order parameter which breaks time-reversal symmetry. The
exploration of the aforementioned exotic phases in the corresponding quantum
models employing state-of-the-art numerical quantum many-body approaches
constitute an important direction of future research.

In the context of material realizations, it is worth noting that in
Na6Cu7BiO4(PO4)4[Cl,(OH)]3 [20], the $J_{\times}$ interactions are mediated by
Cu-O-Na-O-Cu superexchange pathways, while the $J_{+}$ bonds pass through a
chloride group but distances that prevent any Cu-Cl hybridization, and thus
not triggering a superexchange. The presence of the nonmagnetic Na ions in the
center of the octagons is likely to trigger a finite $J_{\times}$ interaction
akin to the scenario realized in the kagome based materials kapellasite and
haydeeite [52]. One could think of possible chemical substitutions which are
likely to enhance this interactions, e.g., by replacing Na with Cs which has a
larger ionic radius. Another route towards strengthening the cross-plaquette
couplings would involve preparing the corresponding sulfide version instead of
oxides leading to Cu-S-Na-S-Cu type superexchange pathways, similar to the
breathing chromium spinels [53]. In KCu6AlBiO4(SO4)5Cl [19], it is the sulfate
SO${}_{4}^{2-}$ that occupies the octagon centers and one may consider the
possibility of substituting it with a selenate group SeO${}_{4}^{2-}$ to
enhance both the cross-plaquette couplings. These would constitute interesting
future explorations on the material synthesis front.

###### Acknowledgements.

We thank Harald O. Jeschke for helpful discussions. The Cologne group
acknowledges partial funding from the DFG within Project-ID 277146847, SFB
1238 (projects C02, C03). S.T. thanks the Center for Computational Quantum
Physics at the Flatiron Institute, New York, for hospitality during the
initial stages of this project. M.G. thanks the Bonn-Cologne Graduate School
of Physics and Astronomy (BCGS) for support. The numerical simulations were
performed on the JUWELS cluster at the Forschungszentrum Juelich and the
CHEOPS cluster at RRZK Cologne. Y. I. acknowledges support from the Department
of Science and Technology (DST), India through the MATRICS Grant No.
MTR/2019/001042, CEFIPRA Project No. 64T3-1, the ICTP through the Associates
Programme and from the Simons Foundation through grant number 284558FY19. This
research was supported in part by the National Science Foundation under Grant
No. NSF PHY-1748958, IIT Madras through the Institute of Eminence (IoE)
program for establishing the QuCenDiEM group (Project No.
SB20210813PHMHRD002720), the International Centre for Theoretical Sciences
(ICTS), Bengaluru, India during a visit for participating in the program
“Frustrated Metals and Insulators” (Code: ICTS/frumi2022/9). Y. I.
acknowledges the use of the computing resources at HPCE, IIT Madras.

## References

  * Messio  _et al._ [2011] L. Messio, C. Lhuillier, and G. Misguich, Lattice symmetries and regular magnetic orders in classical frustrated antiferromagnets, Phys. Rev. B 83, 184401 (2011).
  * Ghosh  _et al._ [2019a] P. Ghosh, T. Müller, F. P. Toldin, J. Richter, R. Narayanan, R. Thomale, J. Reuther, and Y. Iqbal, Quantum paramagnetism and helimagnetic orders in the Heisenberg model on the body centered cubic lattice, Phys. Rev. B 100, 014420 (2019a).
  * Balla  _et al._ [2020] P. Balla, Y. Iqbal, and K. Penc, Degenerate manifolds, helimagnets, and multi-$Q$ chiral phases in the classical Heisenberg antiferromagnet on the face-centered-cubic lattice, Phys. Rev. Res. 2, 043278 (2020).
  * Iqbal  _et al._ [2019] Y. Iqbal, T. Müller, P. Ghosh, M. J. P. Gingras, H. O. Jeschke, S. Rachel, J. Reuther, and R. Thomale, Quantum and Classical Phases of the Pyrochlore Heisenberg Model with Competing Interactions, Phys. Rev. X 9, 011005 (2019).
  * Domenge  _et al._ [2005] J.-C. Domenge, P. Sindzingre, C. Lhuillier, and L. Pierre, Twelve sublattice ordered phase in the ${J}_{1}-{J}_{2}$ model on the kagomé lattice, Phys. Rev. B 72, 024433 (2005).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
  * Villain [1979] J. Villain, Insulating Spin Glasses, Z. Phys. B 33, 31 (1979).
  * Zhitomirsky [2008] M. E. Zhitomirsky, Octupolar ordering of classical kagome antiferromagnets in two and three dimensions, Phys. Rev. B 78, 094423 (2008).
  * Bieri  _et al._ [2016] S. Bieri, C. Lhuillier, and L. Messio, Projective symmetry group classification of chiral spin liquids, Phys. Rev. B 93, 094437 (2016).
  * Hickey  _et al._ [2016] C. Hickey, L. Cincio, Z. Papić, and A. Paramekanti, Haldane-Hubbard Mott Insulator: From Tetrahedral Spin Crystal to Chiral Spin Liquid, Phys. Rev. Lett. 116, 137202 (2016).
  * Hickey  _et al._ [2017] C. Hickey, L. Cincio, Z. Papić, and A. Paramekanti, Emergence of chiral spin liquids via quantum melting of noncoplanar magnetic orders, Phys. Rev. B 96, 115115 (2017).
  * Schmoll  _et al._ [2023] P. Schmoll, A. Kshetrimayum, J. Naumann, J. Eisert, and Y. Iqbal, Tensor network study of the spin-$\frac{1}{2}$ Heisenberg antiferromagnet on the shuriken lattice, Phys. Rev. B 107, 064406 (2023).
  * Astrakhantsev  _et al._ [2021] N. Astrakhantsev, F. Ferrari, N. Niggemann, T. Müller, A. Chauhan, A. Kshetrimayum, P. Ghosh, N. Regnault, R. Thomale, J. Reuther, T. Neupert, and Y. Iqbal, Pinwheel valence bond crystal ground state of the spin-$\frac{1}{2}$ Heisenberg antiferromagnet on the shuriken lattice, Phys. Rev. B 104, L220408 (2021).
  * Ralko and Rousochatzakis [2015] A. Ralko and I. Rousochatzakis, Resonating-valence-bond physics is not always governed by the shortest tunneling loops, Phys. Rev. Lett. 115, 167202 (2015).
  * Rousochatzakis  _et al._ [2013] I. Rousochatzakis, R. Moessner, and J. van den Brink, Frustrated magnetism and resonating valence bond physics in two-dimensional kagome-like magnets, Phys. Rev. B 88, 195109 (2013).
  * Lugan  _et al._ [2019] T. Lugan, L. D. C. Jaubert, and A. Ralko, Topological nematic spin liquid on the square Kagome lattice, Phys. Rev. Research 1, 033147 (2019).
  * Siddharthan and Georges [2001] R. Siddharthan and A. Georges, Square kagome quantum antiferromagnet and the eight-vertex model, Phys. Rev. B 65, 014417 (2001).
  * Yakubovich  _et al._ [2021] O. V. Yakubovich, L. V. Shvanskaya, G. V. Kiriukhina, A. S. Volkov, O. V. Dimitrova, and A. N. Vasiliev, Hydrothermal synthesis and a composite crystal structure of Na6Cu7BiO4(PO4)4[Cl(OH)]3 as a candidate for quantum spin liquid, Inorg. Chem. 60, 11450 (2021).
  * 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, Gapless spin liquid in a square-kagome lattice antiferromagnet, Nat. Commun. 11, 3429 (2020).
  * Liu  _et al._ [2022] B. Liu, Z. Zeng, A. Xu, Y. Sun, O. Yakubovich, L. Shvanskaya, S. Li, and A. Vasiliev, Low-temperature specific-heat studies on two square-Kagome antiferromagnets, Phys. Rev. B 105, 155153 (2022).
  * [21] M.M. Markina, P.S. Berdonosov, T.M. Vasilchikova, K.V. Zakharov, A.F. Murtazoev, V.A. Dolgikh, A.V. Moskvin, V.N. Glazkov, A.I. Smirnov, A.N. Vasiliev, Static and resonant properties of decorated square kagome lattice compound KCu7(TeO4)(SO4)5Cl, arXiv:2212.11623 .
  * Janson  _et al._ [2008] O. Janson, J. Richter, and H. Rosner, Modified Kagome Physics in the Natural Spin-$1/2$ Kagome Lattice Systems: Kapellasite ${\mathrm{Cu}}_{3}\mathrm{Zn}(\mathrm{OH}{)}_{6}{\mathrm{Cl}}_{2}$ and Haydeeite ${\mathrm{Cu}}_{3}\mathrm{Mg}(\mathrm{OH}{)}_{6}{\mathrm{Cl}}_{2}$, Phys. Rev. Lett. 101, 106403 (2008).
  * Janson  _et al._ [2009] O. Janson, J. Richter, and H. Rosner, Intrinsic peculiarities of real material realizations of a spin-1/2 kagomé lattice, J. Phys. Conf. Ser. 145, 012008 (2009).
  * Morita and Tohyama [2018] K. Morita and T. Tohyama, Magnetic phase diagrams and magnetization plateaus of the spin-1/2 antiferromagnetic Heisenberg model on a square-Kagome lattice with three nonequivalent exchange interactions, J. Phys. Soc. Jpn. 87, 043704 (2018).
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Hidden order in a frustrated system: Properties of the Heisenberg Kagomé antiferromagnet, Phys. Rev. Lett. 68, 855 (1992).
  * Wang  _et al._ [2007] F. Wang, A. Vishwanath, and Y. B. Kim, Quantum and classical spins on the spatially distorted kagomé lattice: Applications to volborthite ${\mathrm{Cu}}_{3}{\mathrm{V}}_{2}{\mathrm{O}}_{7}{(\mathrm{O}\mathrm{H})}_{2}\cdot 2{\mathrm{H}}_{2}\mathrm{O}$, Phys. Rev. B 76, 094421 (2007).
  * Yavors’kii  _et al._ [2007] T. Yavors’kii, W. Apel, and H.-U. Everts, Heisenberg antiferromagnet with anisotropic exchange on the kagomé lattice: Description of the magnetic properties of volborthite, Phys. Rev. B 76, 064430 (2007).
  * Schnyder  _et al._ [2008] A. P. Schnyder, O. A. Starykh, and L. Balents, Spatially anisotropic Heisenberg kagome antiferromagnet, Phys. Rev. B 78, 174420 (2008).
  * Hering  _et al._ [2022] M. Hering, F. Ferrari, A. Razpopov, I. I. Mazin, R. Valentí, H. O. Jeschke, and J. Reuther, Phase diagram of a distorted kagome antiferromagnet and application to Y-kapellasite, npj Comput. Mater. 8, 10 (2022).
  * P. Chandra  _et al._ [1993] P. Chandra, P. Coleman, and I. Ritchey, The anisotropic kagome antiferromagnet: a topological spin glass?, J. Phys. I France 3, 591 (1993).
  * Lieb and Mattis [1962] E. Lieb and D. Mattis, Ordering Energy Levels of Interacting Spin Systems, Journal of Mathematical Physics 3, 749 (1962).
  * [32] H.-J. Schmidt and C. Schröder, Thermodynamics of the spin square, arXiv:2210.13957 .
  * Richter  _et al._ [2022] J. Richter, O. Derzhko, and J. Schnack, Thermodynamics of the spin-half square kagome lattice antiferromagnet, Phys. Rev. B 105, 144427 (2022).
  * Vojta [2003] M. Vojta, Quantum phase transitions, Rep. Prog. Phys. 66, 2069 (2003).
  * Henley [2010] C. L. Henley, The “Coulomb Phase” in Frustrated Systems, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).
  * Richter  _et al._ [2009] J. Richter, J. Schulenburg, P. Tomczak, and D. Schmalfuß, The Heisenberg antiferromagnet on the square-kagome lattice, Cond. Matter Phys. 12, 507 (2009).
  * Huse and Rutenberg [1992] D. A. Huse and A. D. Rutenberg, Classical antiferromagnets on the Kagomé lattice, Phys. Rev. B 45, 7536 (1992).
  * Harris  _et al._ [1992] A. B. Harris, C. Kallin, and A. J. Berlinsky, Possible Néel orderings of the Kagomé antiferromagnet, Phys. Rev. B 45, 2899 (1992).
  * Chern and Moessner [2013] G.-W. Chern and R. Moessner, Dipolar Order by Disorder in the Classical Heisenberg Antiferromagnet on the Kagome Lattice, Phys. Rev. Lett. 110, 077201 (2013).
  * Schmidt and Richter [2022] H.-J. Schmidt and J. Richter, Classical ground states of spin lattices, J. Phys. A: Math. Theor. 55, 465005 (2022).
  * [41] See the ancillary files of this arXiv submission for two Mathematica notebooks that provide a detailed description of (i) the deformed cuboctahedral ground states (phase II) applying the generalized Luttinger-Tisza theory of Ref. 40, and (ii) the semi-analytical approach for all spin-spiral phases.
  * Kalmeyer and Laughlin [1987] V. Kalmeyer and R. B. Laughlin, Equivalence of the resonating-valence-bond and fractional quantum Hall states, Phys. Rev. Lett. 59, 2095 (1987).
  * Wen [1990] X. G. Wen, Topological orders in rigid states, International Journal of Modern Physics B 04, 239 (1990).
  * Note [1] In other words, to use the language of [1], the non-coplanar orders found here are not regular magnetic orders on the square-kagome lattice. A regular magnetic order would be a classical spin state which respects all lattice symmetries modulo global $O(3)$ spin transformations. Thus, restoring spin rotational symmetry would result in a fully symmetric state.
  * Lu [2016] Y.-M. Lu, Symmetric ${Z}_{2}$ spin liquids and their neighboring phases on triangular lattice, Phys. Rev. B 93, 165113 (2016).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * Oshikawa [2000] M. Oshikawa, Commensurability, Excitation Gap, and Topology in Quantum Many-Particle Systems on a Periodic Lattice, Phys. Rev. Lett. 84, 1535 (2000).
  * Hastings [2004] M. B. Hastings, Lieb-Schultz-Mattis in higher dimensions, Phys. Rev. B 69, 104431 (2004).
  * Shannon  _et al._ [2006] N. Shannon, T. Momoi, and P. Sindzingre, Nematic Order in Square Lattice Frustrated Ferromagnets, Phys. Rev. Lett. 96, 027213 (2006).
  * Iqbal  _et al._ [2016] Y. Iqbal, P. Ghosh, R. Narayanan, B. Kumar, J. Reuther, and R. Thomale, Intertwined nematic orders in a frustrated ferromagnet, Phys. Rev. B 94, 224403 (2016).
  * [51] S. Jiang, J. Romhányi, S. R. White, M. E. Zhitomirsky, and A. L. Chernyshev, The Elusive Spin Nematic, arXiv:2209.00010 .
  * Iqbal  _et al._ [2015] Y. Iqbal, H. O. Jeschke, J. Reuther, R. Valentí, I. I. Mazin, M. Greiter, and R. Thomale, Paramagnetism in the kagome compounds $(\mathrm{Zn},\mathrm{Mg},\mathrm{Cd}){\mathrm{Cu}}_{3}{(\mathrm{OH})}_{6}{\mathrm{Cl}}_{2}$, Phys. Rev. B 92, 220404 (2015).
  * Ghosh  _et al._ [2019b] P. Ghosh, Y. Iqbal, T. Müller, R. T. Ponnaganti, R. Thomale, R. Narayanan, J. Reuther, M. J. P. Gingras, and H. O. Jeschke, Breathing chromium spinels: a showcase for a variety of pyrochlore Heisenberg Hamiltonians, npj Quantum Materials 4, 63 (2019b).

Figure 17: Variants of cuboc1 order on the square-kagome lattice. (a)-(d) show
four different possible arrangements of cuboc1 order on the square-kagome
lattice, each with its corresponding structure factor. The arrangement in (a)
corresponds to the arrangement shown in Fig. 10 in the main text and has a
clockwise placement of the building blocks highlighted in blue, green, yellow,
and red, in the magnetic unit cell. The arrangement in (c) also has this
clockwise placement, but in a different order, and has the same structure
factor as (a). The other two arrangements, (b) and (d), on the other hand,
have an anti-clockwise placement of blue, green, yellow, and red, resulting in
a structure factor that is rotated by $\pi/2$ (two additional arrangements are
possible in each case, but not explicitly shown). (e) Averaging the real space
correlations of a clockwise and an anti-clockwise arrangement results in a
structure factor with twice as much subdominant peaks. (f) The averaged real
space correlations of all possible cuboc1 orders yields a structure factor
with only main peaks on the edges of the extended Brillouin zone without
further subdominant peaks.

##  Appendix A Different cuboc1 order realizations

There are several different possibilities to set up cuboc1 order on the
square-kagome lattice. These can be best described by dividing the 24-site
magnetic unit cell into four building blocks, each encompassing one geometric
unit cell, with some specific fixed arrangement of spins – indicated as blue,
green, yellow, and red squares in Fig. 17, where four possible such
arrangements are shown. These building blocks can be arranged clockwise or
anti-clockwise, each in four different ways. The order parameter to
distinguish between clockwise and anti-clockwise arrangements of blue-green-
yellow-red is

$\mathcal{O}=\frac{1}{M}\sum_{\scalerel*{\includegraphics{skew_bow_tie.png}}{\sum}}\mathcal{O}_{ijkl}$
(17)

where

$\mathcal{O}_{ijkl}=\frac{1}{\sqrt{2}}\vec{S}_{i}\cdot\left[(\vec{S}_{j}-\vec{S}_{k})\times(\vec{S}_{l}-\vec{S}_{m})\right]$
(18)

is summed over all skew bow-ties as shown in the inset of the lower panel of
Fig. 11 and $M$ is the number of skew bow-ties. For cuboc1 order with a
clockwise (anti-clockwise) arrangement of blue-green-yellow-red, $\mathcal{O}$
yields $+1$ $(-1)$. $\mathcal{O}$ and its absolute value, $|\mathcal{O}|$,
from Monte Carlo simulations are shown in Fig. 11 as well as the associated
susceptibility

$\chi_{\mathcal{O}}=\frac{M}{T}(\langle\mathcal{O}^{2}\rangle-\langle|\mathcal{O}|\rangle^{2})\,.$
(19)

##  Appendix B Spiral phases for AFM interactions

In this section we will briefly characterize the remaining spiral phases which
are not described in the main text. In all cases, the procedure is based on
the semi-analytical method described in Section III.4. For all scenarios the
nearest-neighbor couplings are set to $J_{1}=J_{2}=J_{3}=1.0$ and the
numerical Monte Carlo ground states are taken at a temperature of $T=10^{-4}$.

####  B.0.1 Spiral phase IV with $D_{6}^{s,t}$ symmetry

We consider the ground state for $J_{+}=J_{\times}=-0.2$. On the left of Fig.
18 (a), a common origin plot of the $N=864\ (L=12)$ spin vectors of the
numerical ground state (at $T=10^{-4}$) is shown. After grouping those spins
that point to the same direction, we are left with $M=228$ unique spin
directions which are shown on the right of Fig. 18 (a). The symmetry group $G$
of this state is generated by two reflections and rotations of $\pi/3$ about
the symmetry axis i.e. $G=D_{6}^{s,t}$, where the superscripts imply
additional mirror symmetries on the $xy$-plane and the $yz$-plane, i.e.
$s=\text{diag}(1,1,-1)$ and $t=\text{diag}(-1,1,1)$. The energy per site can
eventually be described as a function of 19 parameters and numerical
minimization of this energy with $J_{+}=-0.2$ then leads to
$E_{\text{IV},\text{semi-analytical}}=-1.10525$ in accordance with the Monte
Carlo result $E_{\text{IV},\text{MC}}=-1.1051(4)$ for the same parameters. For
more details see the explicit description in the Mathematica files of the
Supplement [41].

Specific heat traces for this spiral phase for different system sizes are
shown in the top panel of Fig. 19. However, due to the incommensurability of
these spin spirals there are significant finite-size effects, making it hard
to narrow down the location of the finite-temperature transition.

Figure 18: Spin spiral phases. Shown are common origin plots (left column) of
the numerical Monte Carlo ground states of spiral phases IV (a), V (b), VIII
(c), IX (d), and X (e), at $T=10^{-4}$. Each common origin plot depicts the
spin direction of all $N=864$ spins on the unit sphere. The right column shows
the result of reducing the number of spin directions by grouping spins that
point to the same direction which eventually allows for the symmetry analysis
and analytical description outlined in the text.

####  B.0.2 Spiral phase V with $\\{\text{id},s\\}$ symmetry

For $J_{+}=0.2$ and $J_{\times}=-1.0$, there are $M=47$ unique spin vectors
left after grouping the $N=864$ ground state spin vectors according to unique
directions (cf. Fig. 18 (b)). There is no exact rotation symmetry but a
reflection symmetry with respect to the equatorial plane, thus the symmetry
group of this phase is $\\{\text{id,s}\\}$. The energy per spin then is a
function of 48 parameters and yields, after minimization,
$E_{\text{V},\text{semi-analytical}}=-1.36719$ in accordance with the Monte
Carlo result $E_{\text{V},\text{MC}}=-1.367(1)$ for the same parameters. For
more details see the explicit description in the Mathematica files of the
Supplement [41].

Figure 19: Specific heat of spiral phases IV and VIII. To exemplify the
complicated behaviour of the specific heat of the spiral phases, we show
$c_{V}$ traces for the spiral phase IV (top) and VIII (bottom). In both cases,
the specific heat qualitatively changes drastically with the system size which
indicates that the spiral ground state depends very sensitively on the system
size.

####  B.0.3 Spiral phase VIII with $D_{12}$ symmetry

The ground state for $J_{+}=-1.0$ and $J_{\times}=1.0$ (Fig. 18 (c)) can be
constructed depending on two polar angles $\alpha$ and $\beta$ and one
azimuthal angle $\gamma$ as follows. Let $\vec{S}_{i}\ (i=1,\dots,6)$ be the
six spin vectors of the primitive unit cell with $\vec{S}_{1}$ and
$\vec{S}_{2}$ being the bow-tie spins. One can define

$\displaystyle\vec{S}_{1}$ $\displaystyle=(0,0,1),$ $\displaystyle\vec{S}_{2}$
$\displaystyle=(\sin\alpha,0,\cos\alpha),$ $\displaystyle\vec{S}_{3}$
$\displaystyle=(\sin\beta\cos\gamma,\sin\beta\sin\gamma,\cos\beta),$
$\displaystyle\vec{S}_{4}$
$\displaystyle=(\sin\beta\cos\gamma,-\sin\beta\sin\gamma,\cos\beta),$
$\displaystyle\vec{S}_{5}$
$\displaystyle=(-\sin\beta\cos\gamma,\sin\beta\sin\gamma,\cos\beta),$
$\displaystyle\vec{S}_{6}$
$\displaystyle=(-\sin\beta\cos\gamma,-\sin\beta\sin\gamma,\cos\beta).$ (20)

With $\mu=(m,n,i)$ denoting the general index of a spin site on the lattice,
any spin of the ground state can thus be written as
$\vec{S}_{\mu}=\sigma^{n}\rho^{m}\vec{S}_{i}$ with

$\sigma=\begin{pmatrix}-1&0&0\\\ 0&-1&0\\\
0&0&1\end{pmatrix},\quad\rho=\begin{pmatrix}\cos\tfrac{\pi}{6}&\sin\tfrac{\pi}{6}&0\\\
-\sin\tfrac{\pi}{6}&\cos\tfrac{\pi}{6}&0\\\ 0&0&1\end{pmatrix}.$ (21)

The ground state energy per spin of this state is given as

$\displaystyle E_{\text{VIII}}=\frac{1}{12}$
$\displaystyle\bigg{[}2J_{+}(1+\cos^{2}\alpha)+8\cos\beta(1+\cos\alpha)$
$\displaystyle+4(2+J_{\times})\cos^{2}\beta+\sqrt{3}J_{+}\sin^{2}\alpha$
$\displaystyle+8\cos\gamma\sin\alpha\sin\beta+2\sqrt{3}J_{\times}\cos^{2}\gamma\sin^{2}\beta$
$\displaystyle+4J_{\times}\cos\gamma\sin^{2}\beta\sin\gamma-8\sin^{2}\beta\sin^{2}\gamma$
$\displaystyle-2\sqrt{3}J_{\times}\sin^{2}\beta\sin^{2}\gamma\bigg{]}\,.$ (22)

Specific heat traces for this spiral phase for different system sizes are
shown in the bottom panel of Fig. 19. However, due to the incommensurability
of these spin spirals there are significant finite-size effects, making it
hard to narrow down the location of the finite-temperature transition.

####  B.0.4 Spiral phase IX with $D_{3}^{\sigma}$ symmetry

For $J_{+}=-0.2$ and $J_{\times}=1.0$, grouping the $N=864$ ground state
vectors according to their direction leaves us with $M=108$ unique spin
vectors. These can be reduced further by identifying the symmetry group $G$ of
this state, which is generated by a point reflection and rotations of $2\pi/3$
about the symmetry axis, i.e. $G=D_{3}^{\sigma}$, where the superscript
$\sigma$ implies additional point reflection symmetry, i.e.
$\sigma=\text{diag}(-1,-1,-1)$. The energy per spin is then given by a
function of 26 variables and minimization yields $E_{\text{IX},\text{semi-
analytical}}=-1.3272$ as compared to the Monte Carlo result
$E_{\text{IX},\text{MC}}=-1.327(1)$ for the same parameters. For more details
see the explicit description in the Mathematica files of the Supplement [41].

####  B.0.5 Spiral phase X with $D_{6}^{s}$ symmetry

For $J_{+}=-0.05$ and $J_{\times}=0.4$, there are $M=32$ unique spin
directions left after identifying the numerical ground state spin vectors that
are close to each other. The symmetry group of this state is generated by a
point reflection and by rotations of $\pi/3$ about the symmetry axis. The
energy per site can be expressed as a function of three variables
$\alpha_{1},\alpha_{2},z$ as

$\displaystyle E_{\text{X}}=\frac{1}{36}$
$\displaystyle\bigg{[}8J_{\times}\left(z^{2}-1\right)\sin\left(\pi\left(-\alpha_{1}+\alpha_{2}+\frac{1}{6}\right)\right)$
$\displaystyle+2J_{\times}\left(z^{2}-1\right)\cos\left(\pi\left(\alpha_{1}-\alpha_{2}-\frac{5}{3}\right)\right)$
$\displaystyle+2J_{\times}\left(z^{2}-1\right)\cos\left(\pi\left(\alpha_{1}-\alpha_{2}+\frac{1}{3}\right)\right)$
$\displaystyle-12J_{\times}z^{2}-3J_{+}-2\sqrt{1-z^{2}}\cos\left(\pi\left(\alpha_{1}+\frac{2}{3}\right)\right)$
$\displaystyle+2\sqrt{1-z^{2}}\cos\left(\pi\left(\alpha_{1}+\frac{5}{3}\right)\right)$
$\displaystyle+2\left(z^{2}-1\right)\cos\left(\frac{1}{3}\pi(3\alpha_{1}-3\alpha_{2}-4)\right)$
$\displaystyle+6\left(z^{2}-1\right)\cos\left(\frac{1}{3}\pi(3\alpha_{1}-3\alpha_{2}+2)\right)$
$\displaystyle+2\left(z^{2}-1\right)\cos\left(\pi\left(\alpha_{1}-\alpha_{2}-\frac{4}{3}\right)\right)$
$\displaystyle+2\left(z^{2}-1\right)\cos\left(\pi\left(\alpha_{1}-\alpha_{2}+\frac{2}{3}\right)\right)$
$\displaystyle+4\sqrt{1-z^{2}}\left(\sqrt{3}\sin(\pi\alpha_{1})+\cos(\pi\alpha_{1})\right)$
$\displaystyle+12\sqrt{1-z^{2}}\cos(\pi\alpha_{2})+12z^{2}-24z-12\bigg{]}.$
(23)

Minimization of this expression gives $E_{\text{X},\text{semi-
analytical}}=-1.12813$ in agreement with the Monte Carlo result
$E_{\text{X},\text{MC}}=-1.1280(3)$ for the same parameters.

##  Appendix C Energy and magnetization cuts

We close by showing Monte Carlo data of energy and magnetization for selected
horizontal and vertical cuts through the phase diagrams of Figs. 5 (AFM
nearest-neighbor model) and 15 (FM nearest-neighbor model) in the main text.

###  C.1 AFM nearest-neighbor model

Fig. 20 shows Monte Carlo magnetization cuts for fixed $J_{+}=-1.0$, and
$J_{\times}=1.0$, respectively. The only two phases with non-zero
magnetization are the $120^{\circ}$-d phase (VII) and the spiral phase VIII.
In the deformed $120^{\circ}$ phase, the magnetization only depends on
$J_{\times}$ and is given by
$m_{120^{\circ}\text{-d}}=J_{\times}/(3(2+J_{\times}))$.

In order to demonstrate the good agreement between the (semi-)analytically
derived phase boundaries and Monte Carlo results, we show cuts of $E/N$ for
fixed values of $J_{+}=\pm 1.0$, and $J_{\times}=\pm 1.0$, respectively, from
Monte Carlo simulations along with the corresponding second derivatives in
Fig. 21. The second derivatives of $E/N$ show sharp features exactly at all
(semi-)analytically determined phase boundaries.

Figure 20: Magnetization cuts for AFM nearest neighbor interactions.
Magnetization cuts from Monte Carlo simulations for fixed $J_{+}=-1.0$ (left)
and $J_{\times}=1.0$ (right) with the colorcode of the phases shown in Fig. 5
in the background. Only the deformed $120^{\circ}$ phase (light gray) and the
spiral phase with $D_{12}$ symmetry (yellow) have non-zero magnetization. For
the deformed $120^{\circ}$ phase the magnetization as a function of
$J_{\times}$ can be expressed analytically as $m=J_{\times}/(3(2+J_{\times}))$
(red curve). Figure 21: Energy cuts for AFM nearest neighbor interactions.
Energy cuts from Monte Carlo simulations for fixed $J_{\times}=\pm 1.0$ (top)
and $J_{+}=\pm 1.0$ (bottom). The inset shows the second derivatives of the
energy per spin with the phases shown in the phase diagram Fig. 5 in the
background with the phase boundaries (dashed lines) obtained
(semi-)analytically as described in the main text. The second derivatives of
$E/N$ show sharp features exactly at all phase boundaries.

###  C.2 FM nearest-neighbor model

Similarly good agreement between analytically determined phase boundaries and
Monte Carlo results we find in the case of FM interactions. Fig. 22 shows cuts
of $E/N$ for fixed values of $J_{+}=0.0$, and $J_{+}=0.5$, respectively, from
Monte Carlo simulations along with the corresponding derivatives in Fig. 21.
Since the phase transition is one order higher than usual we consider the
third derivatives of $E/N$. These show sharp features exactly at all
(semi-)analytically determined phase boundaries.

Figure 22: Energy cuts for FM nearest neighbor interactions. Energy cuts from
Monte Carlo simulations for fixed $J_{+}=0.0$ (blue) and $J_{+}=0.5$ (red).
The inset shows the third derivatives of the energy per spin with the phases
shown in Fig. 15 in the background with phase boundaries (dashed lines)
obtained analytically as described in the main text. The third derivatives of
$E/N$ show sharp features exactly at all phase boundaries.