Non-Coplanar Magnetic Orders in Classical Square-Kagome Antiferromagnets

,


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 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 fluctuationdriven 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 such 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][16], or spin and lattice nematics [17].
An ideal playground to explore this physics in experiment has come in the arrival of materials based on the novel square-kagome lattice geometry [18], whose potential to host intricately textured magnetic ground states or quantum spin liquid phases is currently under much investigation [19].Indeed, no sign of long-range magnetic order down to 50 mK has been observed in the spin S = 1/2 Cu 2+ based materials KCu 6 AlBiO 4 (SO 4 ) 5 Cl [20] and Na 6 Cu 7 BiO 4 (PO 4 ) 4 [Cl,(OH)] 3 [21] despite having large negative Curie-Weiss temperatures of −237 K and −212 K, respectively.On the other hand, their sister compounds KCu 7 (TeO 4 )(SO 4 ) 5 Cl and NaCu 7 (TeO 4 )(SO 4 ) 5 Cl develop antiferromagnetic order [22,23], while related compounds Rb 7 (TeO 4 )(SO 4 ) 5 Cl and Cs 7 (TeO 4 )(SO 4 ) 5 Cl do not show sign of magnetic order down to 2K [23].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 [1,5,24,25].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 × 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.
The remainder of this manuscript is structured as follows.
To begin with, we discuss the nearest neighbor Heisenberg model on the square-kagome lattice in Sec.II, where we analyze its ground states, finite-temperature physics, and spinspin correlations.Afterwards, in Sec.III, we introduce additional further-neighbor cross-plaquette interactions and study the resulting rich phase diagram for fixed, antiferromagnetic nearest-neighbor interactions, looking at each phase separately in great detail.Finally, we briefly discuss results for both mixed and pure ferromagnetic nearest-neighbor interactions.

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 con-FIG.1. Square-kagome lattice and interactions.The squarekagome 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 J1 (square) and J2, J3 (bow-tie), respectively.Additionally, we introduce further neighbor cross-plaquette J+-bonds (J×-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.ventional kagome scenario.

A. Ground states
The Hamiltonian with nearest-neighbor couplings can be written as 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 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 where S i , S j are the spins on the squares and 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, |S i | = 1 ∀ i. Satisfying both of these constraints is only possible within a restricted region of parameter space, namely Thus, within this parameter region, all spin configurations satisfying the constraint (3) on each and every triangle are guaranteed to be bona fide 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. [26], which shares the same qualitative physics as the classical spin liquid found in the distorted kagome version of the model [28][29][30][31].
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 • 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 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 θ between S k (the spin on the bow-tie) and the other two spins (on the squares) increases from π/2 to π, passing through 2π/3 exactly at the isotropic point J 2 = J 1 [16].
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 [28].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, η a = ±1, which encodes whether the spins rotate clockwise or anti-clockwise as one goes from say i to j to k [32] (keep in mind that the angles between the spins in each triangle are fixed, e.g.all three angles are fixed to 2π/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 a η a = 0 for the four triangles surrounding a square, and a η a = 0, ±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π/3, but instead some angle incommensurate with respect to 2π.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 [28].
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-updown (UUD) state, sometimes also referred to as a Lieb ferrimagnet [33], 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).

B. 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 longrange 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 → 1 as T → 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) • 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 [27], 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 con-figurations, 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 → 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 foursite 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 → [6•(1/2)]/6 = 1/2 as T → 0 [34].(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 [35]).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 ≈ 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 entropicdriven order-by-disorder mechanism.This is accompanied by c V → [10 • (1/2) + 2 • (1/4)]/6 = 11/12 as T → 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 [27]).(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 → [8 • (1/2) + 4 • (1/4)]/6 = 10/12 as T → 0.
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 ≈ 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 [36].

C. 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 [37].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(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 S △ = 0 constraint of Eq. 3 and thus the onset of strong correlations be-tween 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 q = (4π, 0) and q = (2π, 2π) (and symmetry related points), respectively.Since one is still in the cooperative paramagnetic regime characterized by a dipolar ∼ 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 ≲ 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 • q = 0 type order [13,38], in contrast to the conventional kagome antiferromagnet which favors √ 3 × √ 3 correlations [39][40][41].While the kagome antiferromagnet develops long-range dipolar magnetic order of the √ 3 × √ 3 type in the limit T → 0 [41], it remains to be established on the square-kagome lattice whether true long-range q = 0 dipolar ordering of 120 • type asymptotically develops as T → 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 ≲ T ≲ 0.05), the pinch points are seen to be present, and with the principal spectral weight at the center of the lobes q = (2π, 2π) (and symmetry related points) which progressively increases, together with a relatively weaker signal at the centers of the squares q = (4π, 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 q = (4π, 0) type points, but in contrast to the isotropic point, these are not indicative of enhanced q = 0 correlations of the 120 • 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 ̸ = J 3 but still inside the degenerate manifold region, one observes that the S(q) are only rotationally invariant possibly reflective of the underlying symmetries of the incipient magnetic order in the limit T → 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(q) become more sharply peaked, as expected, and interestingly the S(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 ≥ 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π, 2π) (and symmetry related) points, and subdominant peaks at (2π, 0) (and symmetry related) points.

III. OCTAGON-PLAQUETTE INTERACTIONS
We now augment the nearest-neighbor model, discussed in the previous Section, with the cross octagon-plaquette interac- tions J + and J × (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 × , 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, obtained by Monte Carlo simulations for systems with a linear system size of L = 12 (N = 864), 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 • ordered phase (VI) in the lower left quadrant, i.e. for ferromagnetic J + and J × (indicated by dark gray in Fig. 5), as well as a distorted version of this 120 • order (VII) for slightly positive J × (indicated by light gray in the phase diagram).In the plain-vanilla 120 • ordered phase spins on elementary triangles form mutual angles of 2π/3, while this angle is increased beyond 2π/3 in the distorted phase as discussed in Sec.III B. 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 C. Somewhat more complex noncoplanar 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 [42][43][44][45][46][47][48], which we describe in Sec.III D.

A. Physics on the axes
But before we dive into the various magnetic orders that can be stabilized by the combined effects of cross octagonplaquette 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 × 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 S k in the constraint (3) is fixed in each and every triangle to be the same.Denoting this fixed direction by M, the remaining two spins (located on the squares) are thus subject to the local constraint on each triangle This results in a global spin configuration in which 1/3 of the spins point along M while the other 2/3 point along a ring at an angle of 2π/3 away from 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 • 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 → 11/12 as T → 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.For FM J × , the spins within each octagon coupled by J × 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 lowtemperature specific heat c V → [11 • (1/2) + 1 • (1/4)]/6 = 23/24 as T → 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(q) as shown in Fig. 7.The S(q) for (J + = −1.0,J × = 0.0) show maxima at positions where the Bragg peaks of the incipient q = 0 order of the 120 • type would show up as T → 0 [13].In contrast, for (J + = 0.0, J × = −1.0), the S(q) display maxima at the expected locations for √ 3 × √ 3 order [13].

B. 120 • order
Turning to the magnetically ordered states of our phase diagram in Fig. 5, we start with the 120 • order found in the lower left quadrant where both cross octagon-plaquette couplings are ferromagnetic.While at zero temperature the ground state is degenerate and coplanar states as well as noncoplanar states with the same energy exist, this degeneracy is lifted at small, finite temperatures and coplanar states are selected over non-coplanar states by a thermal order by disorder mechanism.Therefore, the following discussion concentrates on these coplanar states, whereas the non-coplanar states are discussed briefly in Appendix C.
For the square-kagome lattice geometry at hand, one can, in principle, distinguish three different types of 120 • coplanar 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π/3.As illustrated in Fig. 8, the magnetic unit cell coincides either with the geometric 6-site unit cell as in the regular q = 0 order shown in Fig. 8(a), or contains either 12 sites and is two times larger than the geometric unit cell for the type shown in Fig. 8(b), or 24 sites for the type in Fig. 8(c), with the difference between the two latter 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(b) forming a bow-tie ferromagnet, the bow-tie spins in the second type of 120 • order form stripy antiferromagnets with alternating rows/columns of up and down-pointing spins as shown in Fig. 8(c).The magnetization of the spins on the main diag-onal of the squares, m diag , is different for all three types of order with m diag = 1.0 for order Fig. 8(a), m diag = 0.5 for order Fig. 8(b), and m diag = 0.25 for order Fig. 8(c) (see also Fig. 9).A schematic of the static spin structure factors corresponding to these types of real-space 120 • order is shown for an extended Brillouin zone of the square-kagome lattice in Fig. 8(g), (h), and (i), respectively.All three 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 for varying strengths of the two couplings.
To elucidate the thermodynamics associated with these 120 • orders we show, in Fig. 9, specific heat traces for the point in the lower left corner, J + , J × = −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 FIG.9. Thermodynamics of 120 • order.(top) The specific heat of the 120 • 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 • order, whereas at the sharp feature a specific 120 • order is selected -in the example shown the structure factor at T = 0.1 (inset) coincides with the analytical structure factor in Fig. 8(h), whereas the intermediate structure factor can be obtained by averaging the real space correlations of all possible 120 • orders, i.e. those shown in Fig. 8(g)+(h)+(i) plus all possible rearrangements with the same order.The bottom panel shows the magnetization of the square-diagonal spins mdiag, i.e. the two spins on the main diagonal of the square of each unit cell (inset).The displayed square-diagonal 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 square-diagonal magnetization splits up into different branches that converge to the values of mdiag = 1.0, mdiag = 0.5, and mdiag = 0.25, corresponding to the 120 • orders shown in Fig. 8 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 • order -but without selecting one of the types of 120 • 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 possible 120 • order arrangements.Also the net magnetization of the square-diagonal 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 types 120 • orders, resulting in a sharp feature in the specific heat.Notably, all three types of orders appear with different probabilities, as indicated by histograms of m diag (see Appendix C).
A deformed version of the 120 • order, henceforth termed 120 • -d, exists in a small region touching the conventional 120 • order phase in the upper left quadrant of the phase diagram, indicated as phase VII in Fig. 5.As in the nondeformed case, there exist non-coplanar states (discussed in Appendix C) and coplanar states with the same energy, with the latter being selected by thermal order by disorder at small, but finite temperature.Therefore, the following discussion, again, concentrates on the coplanar deformed 120 • order.Its deformation -two of the three mutual angles between neighboring spins take a value of α > 2π/3, while the third one becomes smaller than 2π/3 -can be seen in the common origin plot of Fig. 5.It can be 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 α > 2π/3, while the third angle, between the red and the green sublattices, becomes 2π − 2α < 2π/3.The angle α is found to vary as cos α = −(2 + J × ) −1 and the ground-state energy per site takes the value 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 • -d = J× 3(2+J×) , which only depends on the values of J × (see also Fig. 25 in Appendix F).

C. Cuboc order
Let us now turn to the first instance of exclusively noncoplanar 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 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π/3, which corresponds to the cuboc1 state, discussed for the kagome lattice with crosshexagonal couplings [1,24].Note that there are eight possible FIG.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π/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 C4 symmetry.
ways to arrange cuboc1 order on the square-kagome lattice as illustrated in Fig. 22 of Appendix D, 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 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 D).
Turning to the thermodynamics of the cuboc phase, we show, in Fig. 11, the specific heat, the cuboc1 order parameter O, and its associated susceptibility χ O (both introduced in Appendix D) for different system sizes between L = 8 and L = 32.The specific heat displays a clearly visible doublepeak structure which, similar to our discussion of the coplanar 120 • 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 O (Eq.D1), which takes values of ±1 for different specific cuboc1 arrangements, builds up and the corresponding susceptibility χ O (Eq.D2) 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 × < 0 and J + ≥ 0, which can be derived by applying the generalized FIG.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 D).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. 22).Middle and bottom panels show the cuboc1 order parameter O (Eq.D1), and the corresponding susceptibility χO (Eq.D2).The lower inset shows how O ijklm is being calculated on a single skew bowtie.This quantity is then averaged over all skew bow-ties in order to calculate O.
Luttinger-Tisza method of Ref. [49], see also the Mathematica files in the Supplement [50].In this cuboc-d phase, the antipodal squares of the cuboctahedron are deformed into antipodal squares with modified z-values while the square in the equatorial plane remains undeformed.Its symmetry is thereby reduced to the symmetry of the dihedral symmetry group D s 4 where the superscript s denotes additional mirror symmetry on the xy-plane, i.e. s = diag(1, 1, −1).The ground state energy per site can be calculated to be The cuboc-d ground state still has zero magnetization, m = 0.

D. 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 [42][43][44][45][46][47][48] 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 S i and S j point approximately in the same direction if S i • S j ≥ γ, where the exact value of γ slightly varies from case to case, but typically γ ≈ 0.995.
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, π) with R ∈ O(3) and π ∈ S M where S M is the permutation group of M elements, such that RS i = S π(i) for all spin vectors 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 I 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 × 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 F. In general, we find excellent agreement.
phase symmetry q vectors semi-analytical  I. Symmetry characterization of the ten ground-state phases of the phase diagram Fig. 5. Given are the ground state symmetry (second column), and the q vectors of each phase (if any, third column).For the six semi-analytically described phases (III-V, VIII-X), the compression of the parametrization of the spin spirals via clustering and symmetrization (Fig. 12) is given in the four columns on the right.Technically, this semi-analytical description is obtained by starting with a common origin plot with N = 864 spins sampled at T = 10 −4 for a linear system size L = 12, and is then given in terms of M , K, and the number of needed parameters (last column) to describe the phase.Let us illustrate this semi-analytical 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 × = −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 S i •S j ≥ 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 π/6 about the vertical symmetry axis as well as by a reflection at the equatorial plane, i.e.G = D s 12 .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 α, z 1 , z 2 in the following way With this compact representation at hand, the energy per site can now be explicitly calculated as Numerical minimization of this energy with J + = +1 and J × = −0.2 then leads to E III,semi-analytical = −1.34515which is in excellent agreement with and (as expected) slightly below the Monte Carlo result E III,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 E. It turns out that the six semi-analytically described phases can be divided into three regular ones that possess q vectors and three irregular ones, possessing no q vectors.A compact summary that characterizes all ten phases of the phase diagram Fig. 5 by the used method, the ground state symmetry, and corresponding q vectors, if any, is given in Table I.FIG. 13.Cuboctahedral order cuboc3.This variant of a noncoplanar cuboctahedral order is found when flipping the bow-tie interactions J2 and J3 of the nearest neighbor model to ferromagnetic (while keeping the square interactions J1 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π/3 on the squares and an angle of π/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 C4 symmetry.

E. 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 → −J 2 and J 3 → −J 3 .
Performing these local spin transformations on the orders of our original phase diagram yields qualitatively new types of orders: The original cuboctahedral order (Fig. 10) becomes a new cuboctahedral order cuboc3 (Fig. 13), where each neighboring pair of spins forms an angle of 2π/3 on the squares and an angle of π/3 on the triangles, which is different to cuboc1 order (where all neighboring pairs of spins form an angle of 2π/3) and to cuboc2 order (where all neighboring pairs of spins form an angle of π/3) [1].The 120 • orders (Fig. 8) are transformed into new coplanar orders in an analogous way, e.g.into a coplanar pentagonal order (Fig. 14), that results from flipping the bow-tie spins in the 120 • order Fig. 8(c).Note that, just as in the original 120 • case, at zero temperature there are both coplanar and non-coplanar states degenerate in energy but the coplanar states are selected at finite temperature via thermal order-by-disorder (see Appendix C for further details).

F. Octagonal and conical orders for FM interactions
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 × ≤ 0, which also extents to the other three quadrants in the phase diagram.Its ground state energy per site is given by It is bound by the hyperbola defined by (J × −1)(J + −1) = 1 2 .Bound by a second hyperbola, given by 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 (a) -(c)).Its ground state has zero magnetization, m = 0, and its energy per site is given by Between these two phases, which touch each other only at the point , there exist non-coplanar umbrellalike states that smoothly interpolate between octagonal and FM order (cf.Fig. 16(d)-(f)).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 π/2, whereas the second cone is formed by the four sublattices of bow-tie spins, again with mutual angles of π/2 but rotated π/4 with respect to the first cone.The two polar angles that describe these two cones depend on J + and J × , the ground state energy per site for those states can be calculated analytically and is given by where the plus sign applies to the left domain with 2 , and the minus sign to the right domain with 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 × .Analytically, one finds Here, the plus sign applies to the right domain with 2 , and the minus sign to the left domain with J + < 1 − 1 √ 2 , and w = 4J 2 × J 2 + − 3J × 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 ̸ = 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 orderby-disorder mechanism.We show that upon introduction of ferromagnetic cross-plaquette interactions of just one type, i.e., either J + or J × , 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 → 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 [51], which hosts bulk semion excitations and a chiral gapless edge mode [52].Furthermore, in our phase diagrams, the noncoplanar orders break lattice symmetries in such a way that simply restoring spin rotational symmetry would not fully restore all lattice symmetries [53], potentially leading to descendent nematic chiral liquids for small values of spin [54].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 [55][56][57].For the non-coplanar umbrella-like double cone states with m ̸ = 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 [58][59][60].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 [61] 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 manybody approaches constitute an important direction of future research.
In the context of material realizations, it is worth noting that in Na 6 Cu 7 BiO 4 (PO 4 ) 4 [Cl,(OH)] 3 [21], the J × 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 × interaction akin to the scenario realized in the kagome based materials kapellasite and haydeeite [62].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 [63].In KCu 6 AlBiO 4 (SO 4 ) 5 Cl [20], it is the sulfate SO 2− 4 that occupies the octagon centers and one may consider the possibility of substituting it with a selenate group SeO 2− 4 to enhance both the cross-plaquette couplings.These would constitute interesting future explorations on the material synthesis front.FIG.19.Three families of non-coplanar orders arising from coplanar 120 • order.In the 120 • phase, there are three one-parameter families of non-coplanar states, each described by a single angle α, β, or γ, whose limiting cases yield the coplanar orders presented in Fig. 8, as sketched on the left.The first family, here depicted in (a)+(d)+(g) by its real space arrangement, a schematic common origin plot, and its structure factor, respectively, for α = π/2, interpolates between the coplanar q = 0 state in Fig. 8(a)+(d)+(g) (limit α = 0) and the coplanar state Fig. 8(b)+(e)+(h) (limit α = π).The second family, here depicted in (b)+(e)+(h) by its real space arrangement, a schematic common origin plot, and its structure factor, respectively, for β = π/4, interpolates between the coplanar state in Fig. 8(b)+(e)+(h) (limit β = 0) and the coplanar state Fig. 8(c)+(f)+(i) (limit β = π/2).The ratio of the weight of the subdominant peaks λ to the weight of the dominant peaks Λ in the structure factor is λ/Λ = 40%.The third family, here depicted in (c)+(f)+(i) by its real space arrangement, a schematic common origin plot, and its structure factor, respectively, for γ = π/4, has the coplanar state in Fig. 8(c)+(f)+(i) for both limits γ = 0 and γ = π/2 and non-coplanar states in between.The ratio of the weight of the subdominant peaks λ to the weight of the dominant peaks Λ in the structure factor is λ/Λ ≈ 39%. a ground state.However, in addition the above three coplanar ground states, these rules are also satisfied by certain 1parameter families of, in general, non-coplanar spin configurations.These can be further characterized by the fact that the sequence of spin vectors in the squares of unit cells is either of the form (a, b, a, b) (1st family) or (a, b, a, c) (2nd family) or (a, b, c, d) (3rd family).It turns out that the 1st family connects the coplanar ground states (a) and (b) of Fig. 8, the 2nd one the coplanar ground states (b) and (c), and the third family loops (c) with itself, see Figure 19.Moreover, in the region of the deformed 120 • phase, there is another one-parameter family of non-coplanar ground states.For more details see the Mathematica files in the Supplement [50].

Order-by-disorder selection of coplanar order
The energy of the three families of non-coplanar ground states, as well as the three coplanar ground states, are all identical.However, as these spin configurations are not related by a symmetry of the Hamiltonian, their free energy, the relevant quantity that needs to be minimized at finite temperature, will not be identical.The selection of a subset of ground states via such a thermal order-by-disorder mechanism is a well-studied phenomenon within the context of frustrated magnets.By expanding to quadratic order in classical fluctuations about a particular ordered ground state one can compute the lowest-order correction to the entropy.With the free energy F = E − T S, the states with the highest entropy will be selected at finite temperatures which, as shown in Fig. 20, are precisely the three coplanar states.Interestingly, these all possess the exact same correction at this lowest quadratic order, meaning that the ultimate selection of the ground state must be via higher-order non-linear corrections.The calculation of such corrections are well beyond the scope of the current manuscript so we instead apply a numerical approach to tease out the preferred ground state.
The three different coplanar 120 • orders, shown in Fig. 8, each have a distinct value of m diag .In Fig. 21, we show histograms from Monte Carlo timeseries showing the equilibrium distribution of m diag within the 120 • phase (J + = J × = −1.0)for different system sizes at different temperatures.Below the phase transition at T = 0.27 (cf.Fig. 9), the distribution of m diag splits up into three distinct peaks, corresponding to the three different coplanar 120 • orders.At the lowest temperature, T = 0.05, and for L = 24, the probability of q = 0 order (Fig. 8(a)) is 2%, whereas the orders Fig. 8(b) and (c) have probabilities of 88% and 10%, respectively.We thus suspect that the order in (b) is the likely spin configuration favored by thermal fluctuations at finite temperatures, on the largest system sizes.Only for smaller system sizes, q = 0 order occurs with significant probability.Overall, the probabilities only slightly change between T = 0.05 and T = 0.15.

Appendix D: 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. 22, 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 where Appendix E: 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 D. For all scenarios the nearestneighbor 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 .
a. Spiral phase IV with D s,t 6 symmetry We consider the ground state for J + = J × = −0.2.On the left of Fig. 23  Carlo result E IV,MC = −1.1051(4)for the same parameters.For more details see the explicit description in the Mathematica files of the Supplement [50].
Specific heat traces for this spiral phase for different system sizes are shown in the top panel of Fig. 24.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. Spiral phase V with {id, s} symmetry For J + = 0.2 and J × = −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. 23

(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 {id,s}.The energy per spin then is a function of 48 parameters and yields, after minimization, E V,semi-analytical = −1.36719 in accordance with the Monte Carlo result E V,MC = −1.367(1)for the same parameters.For more details see the explicit description in the Mathematica files of the Supplement [50].Specific heat traces for this spiral phase for different system sizes are shown in the bottom panel of Fig. 24.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.For J + = −0.2 and J × = 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π/3 about the symmetry axis, i.e.G = D σ 3 , where the superscript σ implies additional point reflection symmetry, i.e. σ = diag(−1, −1, −1).The energy per spin is then given by a function of 26 variables and minimization yields E IX,semi-analytical = −1.3272as compared to the Monte Carlo result E IX,MC = −1.327(1)for the same parameters.For more details see the explicit description in the Mathematica files of the Supplement [50].
e. Spiral phase X with D s 6 symmetry For J + = −0.05 and J × = 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 π/3 about the symmetry axis.The energy per site can be expressed as a function of three variables α 1 , α 2 , z as  Fig. 26.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.

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

FIG. 4 .
FIG. 4. Structure factors of the J1-J2-J3 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 J2 and J3 as indicated in the phase diagram in the left column (top to bottom).The shown squares extent in reciprocal space from −4π to 4π in both dimensions.Rows (a) to (c) show structure factors on the diagonal J2 = J3 within the classical spin liquid regime, especially for the isotropic point J1 = J2 = J3 = 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 J2 and J3.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, J2 = J3 = 2.0, in (f), within the UUD phase, J2 = J3 = 2.5, in (g), and in the J1-J3-Néel ordered phase, J2 = 0.25 and J3 = 2.5, in (h).

FIG. 5 .
FIG. 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× for a nearest neighbor model with isotropic, antiferromagnetic interactions, i.e.J1 = J2 = J3 = 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 Dn 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, obtained from Monte Carlo simulations for systems with linear system size L = 12 (N = 864).The only coplanar orders are found in form of 120 • order (VI) in the lower left quadrant, i.e. ferromagnetic J+ and J× (darker gray phase), as well as a distorted version of this (VII) for sligthly positive J× (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 • order, indicated by the dashed and dotted lines, along which the low temperature specific heat cV (T → 0) takes values of 11/12 and 23/24, respectively (cf.Fig. 6).

FIG. 6 .FIG. 7 .
FIG.6.Specific heat and ground states of the extended model on the axes.For either only J+ or J× turned on and being ferromagnetic, one again finds special values of cV (T → 0) that differ from 1.For For J× = 0.0, J+ < 0.0, the value is 11/12 (cf.purple curve with J+ = −1.0),and for J+ = 0.0, J× < 0.0, we find a value of 23/24 (cf.mint curve with J× = −1.0).The upper insets show common origin plots of the corresponding ground states.

FIG. 8 .
FIG. 8.Coplanar 120 • orders.(a), (b), and (c) show real space arrangements of spins in the three different coplanar 120 • ordered states, with (a) corresponding to regular q = (0, 0) order, (b) corresponding to q = (0, 0) and q = (π, π) and (c) having no qvectors at all.The states have either a 6-site magnetic unit cell that coincides with the geometric unit cell (a), a 12-site magnetic unit cell (large gray rectangle), which is two times larger than the geometric unit cell (small square) (b), or a 24-site magnetic unit cell (c), respectively.The magnetization of the square diagonal spins, mdiag, is 1.0 for state (a), 0.5 for state (b), and 0.25 for state (c).(d)+(e)+(f) There are each three sublattices (each one corresponding to one color).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π/3.(g)+(h)+(i) First and extended Brillouin zones of the square-kagome lattice showing the positions of the corresponding dominant and subdominant Bragg peaks.In (i), the ratio of the weight of the subdominant peaks λ to the weight of the dominant peaks Λ is λ/Λ ≈ 85%.Arising from these three coplanar orders, there are three one-parameter families of non-coplanar orders, whose details are shown in Fig.19.
FIG.9.Thermodynamics of 120 • order.(top) The specific heat of the 120 • 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 • order, whereas at the sharp feature a specific 120 • order is selected -in the example shown the structure factor at T = 0.1 (inset) coincides with the analytical structure factor in Fig.8(h), whereas the intermediate structure factor can be obtained by averaging the real space correlations of all possible 120 • orders, i.e. those shown in Fig.8(g)+(h)+(i) plus all possible rearrangements with the same order.The bottom panel shows the magnetization of the square-diagonal spins mdiag, i.e. the two spins on the main diagonal of the square of each unit cell (inset).The displayed square-diagonal 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 square-diagonal magnetization splits up into different branches that converge to the values of mdiag = 1.0, mdiag = 0.5, and mdiag = 0.25, corresponding to the 120 • orders shown in Fig.8(a), (b), and (c), respectively.

FIG. 12 .
FIG.12.Semi-analytical scheme for spin spiral phases.Starting point is a common origin plot of the N spin vectors of a groundstate spin configuration sampled in Monte Carlo simulations at ultralow 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× = −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 π/6 about the vertical symmetry axis as well as by a reflection at the equatorial plane, i.e.G = D s 12 , and therefore K = 4.

FIG. 14 .
FIG. 14. Pentagonal order.Coplanar pentagonal order arises when flipping the bow-tie interactions J2 and J3 of the nearest neighbor model to ferromagnetic (while keeping the square interactions J1 antiferromagnetic).It is obtained from the original 120 • order (Fig. 8(c)) 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 24site 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 and subdominant Bragg peaks.The ratio of the weight of the subdominant peaks λ to the weight of the dominant peaks Λ is λ/Λ ≈ 85%.

FIG. 15 . 1 √ 2 .
FIG. 15.FM phase diagram.For fixed J1 = J2 = J3 = −1, the model shows a large ferromagnetic phase (light gray), a coplanar phase with octagonal order (dark gray) and non-coplanar umbrellalike double cone states (turquoise) that are an interpolation between octagonal and ferromagnetic states.Ferromagnetic and coplanar phases touch at the single point J+ = J× = 1 − 1 √ 2 .On the right hand side are corresponding common origin plots and structure factors from Monte Carlo simulations.

FIG. 21 .
FIG. 21.Square diagonal magnetization distribution in the 120 • phase.Histograms from Monte Carlo timeseries showing the equilibrium distribution of mdiag in phase VI (J+ = J× = −1.0)for different system sizes at different temperatures.
D2) 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, O yields +1 (−1).O and its absolute value, |O|, from Monte Carlo simulations are shown in Fig. 11 as well as the associated susceptibility FIG. 23.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.

FIG. 26 .
FIG. 26.Energy cuts for AFM nearest neighbor interactions.Energy cuts from Monte Carlo simulations for fixed J× = ±1.0(top) and J+ = ±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.

FIG. 27 .
FIG.27.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.15in 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.