# Symmetry, Superposition and Fragmentation in Classical Spin Liquids: A
General Framework and Applications to Square Kagome Magnets

K.B. Yogendra, Suman Karmakar, and Tanmoy Das  Department of Physics, Indian
Institute of Science, Bangalore 560012, India

###### Abstract

Classical magnets offer glimpses of quantum-like features like spin liquids,
and fractionalization, promising an analogous construction of superposition
and projective symmetry in classical field theory. While models based on
system-specific spin-ice or soft-spin rules exist, a formal theory for general
classical magnets remains elusive. Here, we introduce a mutatis mutandis
symmetry group construction built from a vector field in a plaquette of
classical spins, demonstrating how classical spins superpose in irreducible
representations (irreps) of the symmetry group. The corresponding probability
amplitudes serve as order parameters and local spins as fragmented
excitations. The formalism offers a many-body vector field representation of
diverse ground states, including spin liquids and fragmented phases described
as degenerate ensembles of irreps. We apply the theory specifically to a
frustrated square Kagome lattice, where spin-ice or soft spin rules are inapt,
to describe spin liquids and fragmented phases, all validated through irreps
ensembles and unbiased Monte Carlo simulation. Our work sheds light on
previously unknown aspects of spin-liquid phases and fragmentation and
broadens their applications to other branches of field theory.

Classical spin models can potentially capture exotic phenomena like spin
liquid [1, 2, 3, 4, 5, 6], spin ice [7, 8, 9], and fragmentation [2, 10, 11,
12, 13], order by disorder [14, 15, 16, 17, 18], prethermal discrete time
crystals[19], and exciting progress lies in designing novel frameworks that
mimic quantum principles [6, 20, 21, 22, 23, 24, 25, 26]. One approach,
commonly known as the spin-ice rule, utilizes a class of Hamiltonians
expressed in a quadratic form of total spins within a unit cell, enabling the
depiction of a spin-zero degenerate manifold for spin ice/liquid ground states
[2, 6, 11, 22, 23]. On the other hand, within a soft-spin approximation,
analyzing the eigenenergy spectrum in the Fourier space of extended spin
states enables convenient exploration of global symmetry and topology of many-
body spins, but local spin conservation is sacrificed here[22, 23, 24, 25,
26]. Conversely, real-space studies of local spins offer complementary
advantages, accommodating local constraints, local symmetries, order
parameters, and monopole textures.[11, 4, 27, 28, 6] Both approaches
intertwine in capturing salient features of spin liquids: dispersionless
momentum-space energy mirrors extensive real-space spin degeneracy.
Furthermore, singular pinch-point features in the correlation function[11, 4,
29, 28, 6] corresponds to gapless points with singular wavefunction in
excitation energy dispersion.[24, 25] Another approach utilizes traditional
Landau’s coarse-grained magnetization fields, with or without enforcing local
constraints, which can be fragmented into components exhibiting distinct
correlation properties [30, 10, 11, 31].

Extensive research on these paradigms has explored frustrated lattices like
pyrochlore [32, 33, 34, 28, 6], triangular [35, 36, 37], and Kagome [38, 39,
29], with earlier work focussing on other lattices[15, 40, 18, 4]. Recent
excitement surrounds the square Kagome lattice (Fig. 1), spurred by
experimental hints of spin liquid phases [41, 42, 43] and supported by
theoretical investigations using variants of the Heisenberg model[44, 45, 46,
47, 48]. There are, however, indications of finite Dzyaloshinskii-Moriya (DM)
interaction in these materials [41, 42, 43] which cannot be adequately
captured within the spin-ice or soft-spin models. Moreover, a square Kagome
lattice, boasting several sublattices, presents a superior setting with an
enlarged degenerate manifold and increased fragmentation possibilities that
remained unexplored. We study a classical spin model with XXZ and DM
interactions on a two-dimensional square-Kagome lattice. Our approach
transcends a prior approach[28, 29], which initiates by defining the symmetry
group of the classical spin vector within a translationally invariant
plaquette. Consequently, irreps of the symmetry group form the local basis
states, enabling the plaquette field to superpose between them, and local
spins emerge as their fragmented entities. Expansion parameters, behaving as
Landau-like order parameters, transform, however, under ‘discrete’ spatial
rotations. Interestingly, the order parameters serve as spin’s ‘probability
amplitudes’ and ‘occupation densities’ to irreps state and energy levels.
Notably, the properties indeed echo quantum-like constructions. Through
unbiased classical Monte Carlo simulations, we observe that DM interactions
promote uniform or staggered ordering of irreps containing vortex or anti-
vortex, while CSL states emerge near their critical phase boundaries. In the
CSL phase, local spins remain either fully disordered if the ground state
consists of a randomly distributed irrep ensemble or are fragmented into
extended and point-like entities if the ground state scrambles order and
disorder irreps. Additionally, the spin-spin correlation function is analyzed
in each phase to distinguish between collective and fragmented excitations
with Bragg-like peaks in the order phase, coexisting with pinch-point
excitations in the liquid phases.

Figure 1: (a) A plaquette of a 2D square-Kagome lattice, belonging to the
$\mathsf{D}_{4}$ group, is shown with sublattices enumerated as $i=0-7$. (b)
Among five irreps with different multiplets, we show a few representative
irreps here, while others are shown in SMII.2. Each irrep consists of either
$S_{i}^{\perp}$ (horizontal arrow) or $S_{i}^{z}$ (open and filled dots for up
and down spins) components, with the sizes of the arrows or dots dictate their
magnitudes.

Mathematical foundation: Analogous to the product basis for the quantum case,
a many-body classical field can be expressed as a direct sum of local vector
spaces. We construct the local vector space from the irreps of a symmetry
group defined on a local network of spins within a plaquette $p$, invariant
under a lattice point group $\mathsf{G}$:

$\mathcal{S}_{p}=\bigoplus_{i\in p}S_{i}.$ (1)

$S_{i}=(S^{x}_{i}~{}S^{y}_{i}~{}S^{z}_{i})^{T}\in\mathsf{O}_{i}(3)$ at the
$i^{\rm th}$ site, and $\mathcal{S}_{p}\in\mathsf{O}_{p}(3n)$ where $n$ is the
number of sublattices in $p$. ($\mathsf{O}_{i}(n)$, $\mathsf{O}_{p}(n)$
distinguish the orthogonal symmetry of the site and plaquette fields,
respectively). We denote the irreps of $\mathsf{G}$ by
$m_{\alpha}\in\mathbb{R}$, and its vector representation by
$\mathcal{M}=\bigoplus_{\alpha}d_{\alpha}m_{\alpha}$, where $\alpha$ runs over
distinct irreps, and $d_{\alpha}\in\mathbb{Z}$ denotes their multiplets. The
transformation from the spin space to the irreps space involves an orthogonal
matrix, whose column vectors $\mathcal{V}_{\alpha}$ form the orthonormal basis
of the irreps representation. Expressing $\mathcal{S}_{p}$ in this irreps
space yields

$\mathcal{S}_{p}=\sum_{\alpha=1}^{3n}m_{\alpha}\mathcal{V}_{\alpha}.$ (2)

(The plaquette index is implicit in $m$, $\mathcal{V}$.) Interestingly,
$m_{\alpha}$ conforms to Landau’s order parameter as the coarse-grain average
of local fields, except, here it is invariant under discrete symmetry group in
a plaquette and is interpreted as the probability amplitude of vector field:
$m_{\alpha}=\mathcal{V}_{\alpha}^{T}\mathcal{S}_{p}$. The local spins are the
fragmented entities in the irreps space, defined by a rectangular projection
matrix $\mathcal{P}_{i\in p}$ as ${\bf S}_{i\in p}=\mathcal{P}_{i\in
p}\mathcal{S}_{p}=\sum_{\alpha}{m}_{\alpha}\mathcal{P}_{i\in
p}\mathcal{V}_{\alpha}$.

Reformulating the order parameters in terms of the irreps conveniently
decouples them in a symmetry invariant Hamiltonian, albeit the irreps’
multiples can mix among themselves. To account for the multiplets’ submanifold
and emergent symmetry, it is convenient to introduce an
$\mathsf{O}_{p}(d_{\alpha})$ ‘spinor’-like field
$\boldsymbol{m}_{\alpha}:=(m_{\alpha}^{(1)}~{}...~{}m_{\alpha}^{(d_{\alpha})})^{T}$
for the $\alpha$ irrep. Then, the eigenmodes are obtained by orthogonal
rotation
$\tilde{{\boldsymbol{m}}}_{\alpha}=e^{i\boldsymbol{\mathcal{L}}_{\alpha}\cdot{\boldsymbol{\phi}}_{\alpha}}\boldsymbol{{m}}_{\alpha}$,
where $\boldsymbol{\mathcal{L}}_{\alpha}$ are the corresponding generators for
the angle $\boldsymbol{\phi}_{\alpha}$. $\boldsymbol{\phi}_{\alpha}$ lives on
the Hamiltonian’s parameter space and assumes fixed values for the energy
eigenmodes. The orthonormal basis states ensure the constraint
$|\mathcal{S}_{p}|^{2}=\sum_{\alpha}d_{\alpha}|m_{\alpha}|^{2}=nS^{2}$,
$\forall p$, where $|S_{i}|=S$, $\forall i$ is an additional hardcore
constraint on the classical spins. Not all irreps necessarily adhere to the
local constraint, requiring them to collaborate with others for existence.
Such irreps ensembles may lead to non-analyticity and fragmentation into an
order-disorder mixed phase. Additionally, the collapse of the eigenmodes
$\tilde{{\boldsymbol{m}}}_{\alpha}$ into its constituent irrep
$\boldsymbol{{m}}_{\alpha}$ causes distinct fragmented excitation.

We have a $3nN$-dimensional vector space
$\mathcal{S}=\bigoplus_{p}\mathcal{S}_{p}$ for a generic N-unit cell lattice,
commencing a $3nN\times 3nN$-matrix valued quadratic-in-spin Hamiltonian.
However, thanks to nearest-neighbor interaction and discrete-translation-
invariance of the lattice, the Hamiltonian can be brought to a block-diagonal
form in terms of the plaquette Hamiltonian $H_{p}$:

$H_{p}=\frac{1}{2}\mathcal{S}_{p}^{T}\mathcal{H}_{p}\mathcal{S}_{p}.$ (3)

Here $\mathcal{H}_{p}$ is an orthogonal matrix-valued Hamiltonian, analogous
to the second quantized Hamiltonian, whose components consist of all possible
interactions between ${\bf S}_{i}$ and ${\bf S}_{j}$ for $\langle ij\rangle\in
p$. However, lattice symmetries restrict the allowed finite components in
$\mathcal{H}_{p}$, which we now consider for a square kagome lattice.

Realizations in a square-Kagome lattice: The square-Kagome lattice belongs to
the Dihedral $\mathsf{D}_{4}$ group with $n=8$ sublattice spins, giving a
$24$-dimensional vector representation. Denoting the group element
$\mathsf{g}\in\mathsf{D}_{4}$ in the $\mathcal{S}_{p}-$representation by the
matrix-valued operators $\mathcal{D}(\mathsf{g})$, we impose the symmetry
criterion that under a local symmetry transformation
$\mathcal{S}_{p}\rightarrow\mathcal{D}(\mathsf{g})\mathcal{S}_{p}$, the local
Hamiltonian $H_{p}$ is invariant if
$[\mathcal{D}(\mathsf{g}),\mathcal{H}_{p}]=0$, $\forall p,\mathsf{g}$. Since
local $\mathsf{O}_{i}(3)$ and sublattice symmetries are abandoned, the
plaquette symmetry allows us to have bond- and spin-dependent interactions
$J^{\mu\nu}_{ij}$ with six exchange and three DM interactions (see SM for the
details), leading to a bond-dependent XYZ-Heisenberg model with XY-DM
interaction. However, imposing bond-independent interactions, we consider an
XXZ model with DM interaction as more appropriate for realistic materials [41,
42, 43], $H=\sum_{\langle ij\rangle,\mu\nu}J^{\mu\nu}S_{i}^{\mu}S_{j}^{\nu}$.
This can, for future convenience, be expressed as:

$\displaystyle{}H=J\sum_{\langle ij\rangle,\tau=\pm}\left(D^{\tau}e^{{\rm
i}\tau(\Theta_{i}+\Theta_{j})}S^{\perp}_{i}S^{\perp}_{j}+\Delta
S_{i}^{z}S_{j}^{z}\right).$ (4)

Here $J^{\mu\nu}=J\delta_{\mu\nu}+JD\epsilon_{\mu\nu}$ for $\mu=x,y$, and
$J^{zz}=J\Delta$, $\delta_{\mu\nu}$ is the Kronecker delta and
$\epsilon_{\mu\nu}$ is the Levi-Civita tensor. $J$ is the exchange term,
$\Delta$ is the $z$-axis anisotropy ratio, and $JD$ is the XY DM interaction
strength. By diagonalizing the tensor $J^{\mu\nu}$, we define two ‘circularly
polarized’ fields: $S_{i}^{\tau}=|S^{\perp}_{i}|e^{{\rm
i}\tau\Theta_{i}}\in\mathsf{O}_{i}(2)\cong\mathsf{U}_{i}(1)$, where
$S^{\perp}_{i}=\sqrt{S^{2}-(S^{z}_{i})^{2}}$ is the coplanar spin magnitude
and $\Theta_{i}$ is the azimuthal angle in the spin space, which interact via
a complex (dimensionless) interaction $D^{\tau}=1+{\rm i}\tau D$.

Figure 2: Computed phase diagrams within the Monte Carlo simulation are shown
for (a) for AFM ($J=+1$) and (b) for the FM ($J=-1$) couplings. We highlight
spin textures in a randomly chosen four-plaquette setting for representative
phases. The displayed phases are (c) AFM CSL at $(J,\Delta,D)=(1,1,0)$, (d)
Anti-vortex order at $(1,0,-3)$, (e) Fragmented phase at $(1,4,-1)$ where
$\mathsf{S}_{i}^{z}$ values are random while $S_{i}^{\perp}$ are ordered in a
staggered AFM-anti-vortex texture, (f) Fragmented phase at $(-1,-2.5,0)$ where
$\mathsf{S}_{i}^{z}$ is disorder while $S_{i}^{\perp}$ exhibit collinear
ordering.

Irreps in square-Kagome lattice: There are five conjugacy classes in this non-
Abelian group, giving five irreps: $m_{\alpha}\equiv$
$\mathsf{A}^{(d_{\alpha})}_{1,2}$, $\mathsf{B}^{(d_{\alpha})}_{1,2}$, and a
two-dimensional $\mathsf{E}^{(d_{\alpha})}$, where the superscript denotes
their multiplicity $(d_{\alpha})=(2,4,3,3,6)$, respectively. Representative
irreps configurations are shown in Fig. 1(b).

We organize these irreps into an out-of-plane set
$\mathsf{m}_{z}:=\\{\mathsf{A}_{2}^{\rm(c,d)},\mathsf{B}_{1,2}^{\rm(c)},\mathsf{E}^{\rm(e,f)}\\}$,
and a coplanar set:$\mathsf{m}_{\perp}:=\mathsf{m}_{z}^{\complement}$.
Coplanar irreps $\mathsf{A}^{\rm{(a,b)}}_{1,2}$,
$\mathsf{B}_{1,2}^{\rm{(a,b)}}$ are even and odd under $\mathsf{C}_{4}$,
forcing $S_{i}^{\tau}$ to obey a homeomorphism $\Theta_{i\in
p}=Q_{p}\theta_{i}+\gamma_{p}$, where $\Theta_{i}$ and $\theta_{i}$ are the
(local) azimuthal angles in the spin and position manifolds, respectively,
$\gamma_{p}\in[0,\pi)$ is the (global in $p$) helicity angle, and
$Q_{p}\in\pi_{1}(\mathbb{S}^{1})\cong\mathbb{Z}$ is the topological charge. As
shown in Fig. 1(b), this leads to two concentric (anti-/) vortex substructures
in the outer and inner squares within each plaquette, which are unrelated by
symmetry and interact solely through $D^{\tau}$. $\mathsf{A}_{1}^{\rm{(a,b)}}$
consist of concentric vortices with the same/opposite helicities
($\gamma_{p}=\pm\pi/2$), while $\mathsf{A}_{2}^{\rm{(a,b)}}$, odd under
reflection, have $\gamma_{p}=\pm\pi$. $\mathsf{B}_{1,2}^{\rm{(a,b)}}$ irreps
(odd under $\mathcal{C}_{4}$) are similar, except they consist of anti-
vortices. The out-of-plane $\mathsf{A}_{2}^{\rm{(c,d)}}$ are colinear FM/AFM
irreps, while $\mathsf{B}_{1,2}^{\rm{(c)}}$ are colinear AFM irreps that do
not satisfy the local constraint. Finally, among the six-fold multiplets of
$\mathsf{E}$ irrep, $\mathsf{E}^{\rm(a-d)}$ are co-planer FM/ nematic/AFMorder
parameters, while $\mathsf{E}^{\rm(e,f)}$ are out-of-plane irreps that violate
the local constraints.

Eigen energies: The final task is to diagonalize the multiples of the irreps.
In our case, the irreps’ multiplets split as either
$\mathsf{O}_{p}(d_{\alpha})=\mathsf{O}_{p}(2)\oplus\mathsf{O}_{p}(2)\oplus...$,
or
$\mathsf{O}_{p}(d_{\alpha})=\mathsf{O}_{p}(2)\oplus\mathsf{Z}_{2}\oplus...$,
in which all $\mathsf{O}_{p}(2)$ operators have the same generator
$\mathcal{L}_{\alpha}=i\sigma_{y}$. $\phi_{\alpha}$ depends only on
$\rm{arg}(D^{\tau})$ in the eigenstates of $\mathcal{H}_{p}$. The resultant
diagonal Hamiltonian per plaquette is

$\displaystyle H_{p}=\sum_{\nu=1}^{3n}E_{\nu}|\tilde{{\bf m}}_{\nu}|^{2}.$ (5)

Here $|\tilde{{\bf m}}_{\nu}|^{2}$ serves as an ‘occupation density’ to the
${\nu}^{\rm th}$ order parameter’s energy level $E_{\nu}$. We, henceforth,
omit the tilde symbol for simplicity, and all irreps are taken to be
eigenmodes unless mentioned otherwise. Constrained by symmetry,
$E_{\nu\in\mathsf{m}_{\perp}}$ depends solely on $D^{\tau}$, while
$E_{\nu\in\mathsf{m}_{z}}$ are proportional to $\Delta$ 111Specifically, the
(anti-/) vortex irreps $\mathsf{B}_{1,2}^{\rm{(a,b)}}$,
$\mathsf{A}_{1,2}^{\rm{(a,b)}}$ are promoted by $\mp D$, while
$\mathsf{E}^{\rm(a,b)}$ do not depend on $D$.. One or more irrep (s) can form
a uniform (order) phase with a global energy minimum at $NE_{\nu}$ if they
satisfy the constraint and frustration; otherwise, they blends with other
irreps to form a degenerate ensemble, distributing randomly in the lattice.

Phase diagrams and correlation functions: Using classical Monte-Carlo
simulations of Eq.(̃4) under the imposed local constraint, we generate the
phase diagram presented in Fig. 2. Notably, across all phases (ordered,
disordered, and fragmented), the spin texture within each plaquette adheres to
the irreps, which permits us to construct a many-body ground state vector
field for all phases:

$\displaystyle\mathcal{S}_{\rm GS}$ $\displaystyle=$
$\displaystyle\bigoplus_{p}\sum_{\\{\nu_{p}\\}}m_{\nu_{p}}\mathcal{V}_{\nu_{p}}.$
(6)

The ordered phase harbors a summated state of a fixed irrep
$\bar{\nu}\in\\{\nu_{p}\\}$ (with $m_{\bar{\nu}}=\bar{m}$,
$m_{\nu\neq\bar{\nu}}=0$, $\forall p$); while the staggered phase features two
alternating but fixed irreps $\bar{\nu}_{p}$ and $\bar{\nu}_{q}$ in
neighboring plaquettes. The CSL state, on the other hand, combines a dynamic
ensemble of irreps $\\{\nu_{p}\\}$ within each plaquette $p$. Within this
ensemble, the probability amplitude $m_{\nu_{p}}$ may vary randomly, subject
to local constraints, for the same plaquette energy. The random distribution
of $m_{\nu_{p}}$ differs between plaquettes, resulting in an extensively
degenerate ground state.

In addition to the consistency between the unbiased Monte-Carlo simulation and
the irreps constructions, we also compare our results with a soft-spin
approximation in the Fourier space ([34, 50, 51, 4, 52, 6] and see SM). Given
that we have experimental access to the correlation function of local spins
${\bf S}_{i\in p}$, we report its correlation function. We project the
structure factor
$\chi(\textbf{k})=1/\mathcal{N}\sum_{i,j}\langle\textbf{S}_{i}\cdot\textbf{S}_{j}\rangle\exp{({{\rm
i}\textbf{k}\cdot(\textbf{r}_{i}-\textbf{r}_{j})})}$ to the irreps space as

$\displaystyle\langle\textbf{S}_{i}\cdot\textbf{S}_{j}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\nu_{p}\nu_{q}}{m}_{\nu_{p}}{m}_{\nu_{q}}\langle\mathcal{{V}}_{\nu_{p}}^{T}\mathcal{P}^{T}_{i}\mathcal{P}_{j}\mathcal{{V}}_{\nu_{q}}\rangle,$
(7)

with ri is the $i^{\rm th}$ spin’s position in $p$ and $j\in q$ plaquette.

Figure 3: Simulated static structure factors are plotted in the momentum space
for the four phases discussed in Fig. 2. (a) AFM CSL at
$(J,\Delta,D)=(1,1,0)$, where red dots are plotted separately to signify
additional strong Bragg-like peaks that overwhelm the spectral density of the
disordered pattern. (b-c) Fragmented phase at $(1,4,-1)$ where the plots for
the ordered $S_{i}^{\perp}$ and disordered $S_{i}^{z}$ components are
separated in (b) and (c), respectively. (d) Anti-vortex order at $(1,0,-3)$
showing Bragg peaks similar to $S_{i}^{\perp}$ components in (b). (e-f)
Fragmented phase at $(-1,-2.5,0)$ with FM ordered $S_{i}^{\perp}$ and disorder
$S_{i}^{z}$ are separated in (e) and (f). Panels (a) and (f) host pinch-points
around $(\pi,3\pi)$ and its equivalent points.

The phase diagram reveals a predominance of (uniform or staggered) order
phases in both $J<0$ (frustration inactive) and $J>0$ (frustration active)
regions, with a CSL phase emerging at the critical line of $D\rightarrow 0$,
otherwise, it turns into distinct mixed phases for $2|D|/\Delta<1$. For
$D\rightarrow 0$, $J>0$, three distinct CSL phases emerge with varying
$\Delta$. At $D=0$, the coplanar irreps $\mathsf{A}_{1,2}^{\rm(a,b)}$,
$\mathsf{B}_{1,2}^{\rm(a,b)}$ become degenerate at $-2J$, while
$\mathsf{E^{\rm(c,d)}}$, satisfying the constraint but not frustration, have
the lowest energy at $-4J$. As $\Delta\rightarrow 0$, the Hamiltonian (first
term in Eq. 4) is constrained by a local $\mathsf{O}_{i}(2)$ symmetry of the
$S_{i}^{\tau}$ fields, and the structure factor $\chi({\bf k})$ receives only
finite contribution from $S_{i}^{\tau}$ and no allocation to $S_{i}^{z}$.
Moreover, $\chi({\bf k})$ displays a characteristic disorder pattern without
any Bragg-like peak but with a prominent pinch-point around ${\bf
k}=(\pm\pi,\pm 3\pi)$. The pinch-point characterizes an algebraic correlation
between the topological charge of the $\mathsf{O}_{p}(2)$ multiplets. At
$\Delta=1$, the Hamiltonian is subject to a full $\mathsf{O}_{i}(3)$ symmetry
constraint per site, resulting in symmetry-allowed access to the entire
ensemble $\\{m_{\nu_{p}}\\}\subseteq\mathsf{m}_{\perp}\cup\mathsf{m}_{z}$. For
example,
$\\{m_{\nu}\\}\in\\{\mathsf{A}_{1,2}^{\rm(a,b,c,d)},\mathsf{B}_{1,2}^{\rm(a,b)}\\}$
are degenerate at $E_{\nu}=-2J$ and
$\\{m_{\nu^{\prime}}\\}\in\\{\mathsf{B}_{1,2}^{\rm(c)},\mathsf{E}^{\rm(c,d)}\\}$
at $E_{\nu^{\prime}}=-4J$, making a larger CSL ensemble degenerate at energy
$E_{p}=m_{\nu}^{2}E_{\nu}+m_{\nu^{\prime}}^{2}E_{\nu^{\prime}}=-4J$ for
$m_{\nu}=\sqrt{2}m_{\nu^{\prime}}$. Consequently, $\chi({\bf k})$ displays
pinch-point correlations among both $S_{i}^{\tau}$ and $S_{i}^{z}$. Finally,
as $\Delta\rightarrow\infty$, the Hamiltonian (last term in Eq. 4) retains a
residual local $\mathsf{Z}_{2}$ symmetry constraint, and the disorder ground
state solely stems from the $\\{m_{\nu_{p}}\\}\subseteq\mathsf{m}_{z}$
ensemble. $\chi({\bf k})$ is contributed solely by $S_{i}^{z}$ with pinch-
points at $\bf k=(\pm\pi,\pm 3\pi)$. Based on their distinct local
constraints, it is convenient to refer to these phases as $\mathsf{O}(2)$,
$\mathsf{O}(3)$, and $\mathsf{Z}_{2}$ CSLs, respectively, without implying a
Landau-type phase boundary between them.

Any finite $D$ steers the CSL phase into either order or fragmented (mixed)
phases. For weak out-of-plane anisotropy $\Delta<2|D|$, $D^{\tau}$ $(\tau={\rm
sgn}(D))$ interaction uniformly selects a coplanar irrep
$m_{\mathsf{B}_{1}/\mathsf{A}_{1}}$ containing concentric (anti-/) vortices,
respectively (red/magenta regions in Fig. 2). Here, the (anti-/) vortices of
the same topological charge are staggered between the neighboring plaquettes
with a $\gamma_{p}=\pi$ phase shift, with Bragg-like peaks at ${\bf
k}=(\pi,\pi)$ in $\chi({\bf k})$. However, for strong $\Delta>2|D|$ (and
$J>0$), the homogenous coplanar order becomes scrambled with disordered out-
of-plane irreps: $\\{m_{\nu_{p}}\\}_{\rm mix}\subseteq
m_{\mathsf{A}_{1}/\mathsf{B}_{1}}\cup\mathsf{m}_{z}$, in Fig. 2(e). Their
interplay yields an interesting fragmentation feature in which the outer
vortex maintains coplanarity, while the inner vortex mixes with the
$\mathsf{B}_{1}^{{\rm(c)}}\in\mathsf{m}_{z}$ irrep in each plaquette. The
combination $m_{\mathsf{A}_{1}}=-m_{\mathsf{B}_{1}}$ produces a novel AFM-
vortex/AFM-anti-vortex texture within the inner square where neighboring spins
possess opposite easy axes 222This AFM-vortex topology is homotopically
distinct from the known AFM skyrmion[AFMskyrmion], and has not been predicted
previously. Consequently, $\mathsf{O}_{i}(3)$ field fragments into its
$S_{i}^{z}$ components becoming non-interacting and fail to order or exhibit
any significant correlation, while the $S_{i}^{\tau}$ fields exhibit long-
range order with Bragg-like peaks in the structure factor, see Fig. 3(b). In
essence, this is a unique spin liquid-crystal-like phase arising from a
coordinated spatial distribution of the probability amplitude ($m_{\alpha}$)
of the contributing irreps.

The interplay between the FM interaction, $J=-1$, and strong AFM out-of-plane
anisotropy $\Delta>2|D|$ generates fragmented phases of distinct
characteristics. Here, the in-plane FM irrep $\mathsf{E}^{\rm{(a,b)}}$ pairs
with the out-of-plane AFM $\mathsf{E}^{\rm{(f)}}$ counterpart. The latter
violates the local constraint, leading to an intriguing fragmented structure
in $\chi({\bf k})$, showing a FM ordering in $S_{i}^{\tau}$, but a pinch-point
disorder in $S^{z}_{i}$. This irreps ensemble satisfies
$m_{\mathsf{E}^{\rm(a)}}=m_{\mathsf{E}^{\rm(b)}}=m_{\mathsf{E}^{\rm(f)}}/\sqrt{2}$,
but an extensive degeneracy arises from the possibilities of the four-fold
orientations of the $\mathsf{E}^{\rm(f)}$ irrep, see Fig. 1(b). The DM
interaction disfavors this mixed phase for $D>2\Delta$, leading to a
transition to similar in-plane orders of (anti-/) vortices observed in the
$J=1$ phase diagram. The remaining two phases are readily identifiable: a
uniform coplanar FM order with $\bar{m}_{\mathsf{E}^{\rm(a_{y})}}$ irrep at
$\Delta\rightarrow 0$, and an out-of-plane FM order with
$\bar{m}_{\mathsf{A}_{2}^{\rm(c)}}$) for $J\Delta\rightarrow\infty$.

Conclusions and outlook. We followed an analysis that draws parallels between
quantum and classical field theories in the context of spin liquids and
fractionalization. While the distinction lies in the quantum statistics
manifesting as a direct product basis versus a direct summated field, the
concept of superposition and symmetry group representation remain central to
both. This shared concept underpins the emergence of fragmentation and spin
liquid ground states. Discussions on their excitations and phase transition
are merited. Among the ordered phases, the (anti-/) vortex order phases (red
and magenta) exhibit novel collective excitations. As their descriptions are
equivalent, we discuss the vortex case here. Gapless collective excitations
emerge from the long-wavelength fluctuation of the helicity angle $\gamma_{p}$
across the lattice, protected by the topology of the irreps space through the
charge $Q_{p}\in\mathbb{Z}$. These modes, termed helicity phase modes or
phasons, possess novel characteristics. The two concentric vortices per
plaquette are coupled by interaction but not symmetry. Frustration affects
only the outer vortex, resulting in the fragmentation of the excitation
spectrum into a collective mode for the ordered fields and local excitations
for the disordered components. The Mermin-Wagner theorem dictates the
instability of ordered states to gapless magnons or phason modes, while
disorder phases tend to order via thermal fluctuations according to the order-
by-disorder paradigm [14, 16, 17, 18]. Notably, the two ordered phases of
vortex and anti-vortices for $\pm D$ consists of different irreps, i.e.,
distinct conjugacy classes that do not couple in the Hamiltonian. Hence their
phase boundary at $D=0$ signifies a topological phase transition, associated
with a spin liquid phase at the critical point, reminiscence of the deconfined
critical point [54]. The CSL critical point can be extended by applying a
magnetic field in the $z$-direction (see SM). Moreover, transitions between
ordered and fragmented phases, or within fragmented phases, offer intriguing
avenues for studying non-Landau-type phase transitions.

Acknowledgements: TD acknowledges research funding from the S.E.R.B.
Department of Science and Technology, India, under Core Research Grant (CRG)
Grant No. CRG/2022/00341 and acknowledges the computational facility at
S.E.R.C. Param Pravega under NSM grant No. DST/NSM/R&D HPC
Applications/2021/39.

## References

  * Balents [2010] L. Balents, Nature 464, 199 (2010).
  * Henley [2010] C. L. Henley, Annual Review of Condensed Matter Physics 1, 179 (2010).
  * Moessner and Chalker [1998] R. Moessner and J. T. Chalker, Phys. Rev. Lett. 80, 2929 (1998).
  * Rehn  _et al._ [2016] J. Rehn, A. Sen, K. Damle, and R. Moessner, Phys. Rev. Lett. 117, 167201 (2016).
  * Benton  _et al._ [2016] O. Benton, L. Jaubert, H. Yan, and N. Shannon, Nature Communications 7, 10.1038/ncomms11572 (2016).
  * Benton and Moessner [2021] O. Benton and R. Moessner, Phys. Rev. Lett. 127, 107202 (2021).
  * Castelnovo  _et al._ [2012] C. Castelnovo, R. Moessner, and S. Sondhi, Annual Review of Condensed Matter Physics 3, 35–55 (2012).
  * Nisoli  _et al._ [2013] C. Nisoli, R. Moessner, and P. Schiffer, Rev. Mod. Phys. 85, 1473 (2013).
  * Ross  _et al._ [2011] K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Phys. Rev. X 1, 021002 (2011).
  * Brooks-Bartlett  _et al._ [2014] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X 4, 011007 (2014).
  * Powell [2015] S. Powell, Phys. Rev. B 91, 094431 (2015).
  * Petit  _et al._ [2016] S. Petit, E. Lhotel, B. Canals, M. Ciomaga-Hatnean, J. Ollivier, H. Mutka, E. Ressouche, A. R. Wildes, M. R. Lees, and G. Balakrishnan, Observation of magnetic fragmentation in spin ice (2016), arXiv:1603.05008 [cond-mat.str-el] .
  * Canals  _et al._ [2016] B. Canals, I.-A. Chioar, V.-D. Nguyen, M. Hehn, D. Lacour, F. Montaigne, A. Locatelli, T. O. Menteş, B. S. Burgos, and N. Rougemaille, Nature Comm. 7, 10.1038/ncomms11446 (2016).
  * Villain, J.  _et al._ [1980] Villain, J., Bidaux, R., Carton, J.-P., and Conte, R., J. Phys. France 41, 1263 (1980).
  * Henley [1989] C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Chandra  _et al._ [1990] P. Chandra, P. Coleman, and A. I. Larkin, Phys. Rev. Lett. 64, 88 (1990).
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Phys. Rev. Lett. 68, 855 (1992).
  * Shender and Holdsworth [1996] E. F. Shender and P. C. W. Holdsworth, Order by disorder and topology in frustrated magnetic systems, in  _Fluctuations and Order: The New Synthesis_, edited by M. Millonas (Springer US, New York, NY, 1996) pp. 259–279.
  * Pizzi  _et al._ [2021] A. Pizzi, A. Nunnenkamp, and J. Knolle, Phys. Rev. Lett. 127, 140602 (2021).
  * Placke  _et al._ [2023] B. Placke, O. Benton, and R. Moessner, Ising fracton spin liquid on the honeycomb lattice (2023), arXiv:2306.13151 [cond-mat.str-el] .
  * Dahlbom  _et al._ [2023] D. Dahlbom, H. Zhang, Z. Laraib, D. M. Pajerowski, K. Barros, and C. Batista, Renormalized classical theory of quantum magnets (2023), arXiv:2304.03874 [cond-mat.str-el] .
  * Nutakki  _et al._ [2023] R. P. Nutakki, L. D. C. Jaubert, and L. Pollet, SciPost Phys. 15, 040 (2023).
  * Davier  _et al._ [2023] N. Davier, F. A. Gómez Albarracín, H. D. Rosales, and P. Pujol, Phys. Rev. B 108, 054408 (2023).
  * Yan  _et al._ [2023a] H. Yan, O. Benton, R. Moessner, and A. H. Nevidomskyy, Classification of classical spin liquids: Typology and resulting landscape (2023a), arXiv:2305.00155 [cond-mat.str-el] .
  * Yan  _et al._ [2023b] H. Yan, O. Benton, A. H. Nevidomskyy, and R. Moessner, Classification of classical spin liquids: Detailed formalism and suite of examples (2023b), arXiv:2305.19189 [cond-mat.str-el] .
  * Fang  _et al._ [2023] Y. Fang, J. Cano, A. H. Nevidomskyy, and H. Yan, Classification of classical spin liquids: Topological quantum chemistry and crystalline symmetry (2023), arXiv:2309.12652 [cond-mat.str-el] .
  * Essafi  _et al._ [2016] K. Essafi, O. Benton, and L. Jaubert, Nature Communications 7, 10.1038/ncomms10297 (2016).
  * Yan  _et al._ [2017] H. Yan, O. Benton, L. Jaubert, and N. Shannon, Phys. Rev. B 95, 094422 (2017).
  * Essafi  _et al._ [2017] K. Essafi, O. Benton, and L. D. C. Jaubert, Phys. Rev. B 96, 205126 (2017).
  * Bramwell [2012] S. T. Bramwell, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 5738 (2012).
  * Savary and Balents [2012] L. Savary and L. Balents, Phys. Rev. Lett. 108, 037202 (2012).
  * Siddharthan  _et al._ [1999] R. Siddharthan, B. S. Shastry, A. P. Ramirez, A. Hayashi, R. J. Cava, and S. Rosenkranz, Phys. Rev. Lett. 83, 1854 (1999).
  * Bramwell and Gingras [2001] S. T. Bramwell and M. J. P. Gingras, Science 294, 1495 (2001).
  * Isakov  _et al._ [2004] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 93, 167204 (2004).
  * Mól  _et al._ [2012] L. A. S. Mól, A. R. Pereira, and W. A. Moura-Melo, Phys. Rev. B 85, 184410 (2012).
  * Li  _et al._ [2020] Y. Li, P. Gegenwart, and A. A. Tsirlin, Journal of Physics: Condensed Matter 32, 224004 (2020).
  * Liu  _et al._ [2016] D. T. Liu, F. J. Burnell, L. D. C. Jaubert, and J. T. Chalker, Phys. Rev. B 94, 224413 (2016).
  * Iqbal  _et al._ [2015] Y. Iqbal, H. O. Jeschke, J. Reuther, R. Valentí, I. I. Mazin, M. Greiter, and R. Thomale, Phys. Rev. B 92, 220404 (2015).
  * Mizoguchi  _et al._ [2017] T. Mizoguchi, L. D. C. Jaubert, and M. Udagawa, Phys. Rev. Lett. 119, 077207 (2017).
  * Yildirim  _et al._ [1998] T. Yildirim, A. B. Harris, and E. F. Shender, Phys. Rev. B 58, 3144 (1998).
  * 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, _et al._ , Nat. Commun. 11, 3429 (2020).
  * Yakubovich  _et al._ [2021] O. V. Yakubovich, L. V. Shvanskaya, G. V. Kiriukhina, A. S. Volkov, O. V. Dimitrova, and A. N. Vasiliev, Inorganic Chemistry 60, 11450 (2021), pMID: 34264636.
  * Markina  _et al._ [2023] M. M. Markina, P. S. Berdonosov, T. M. Vasilchikova, K. V. Zakharov, A. F. Murtazoev, V. A. Dolgikh, A. V. Moskin, V. N. Glazkov, A. I. Smirnov, and A. N. Vasiliev, Static and resonant properties of decorated square kagome lattice compound $\ce{KCu_{7}(TeO_{4})(SO_{4})_{5}Cl}$ (2023), arXiv:2212.11623 [cond-mat.str-el] .
  * 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, Phys. Rev. B 104, L220408 (2021).
  * Richter  _et al._ [2022] J. Richter, O. Derzhko, and J. Schnack, Phys. Rev. B 105, 144427 (2022).
  * Niggemann  _et al._ [2023] N. Niggemann, N. Astrakhantsev, A. Ralko, F. Ferrari, A. Maity, T. Müller, J. Richter, R. Thomale, T. Neupert, J. Reuther, Y. Iqbal, and H. O. Jeschke, Phys. Rev. B 108, L241117 (2023).
  * Richter and Schnack [2023] J. Richter and J. Schnack, Phys. Rev. B 107, 245115 (2023).
  * Gembé  _et al._ [2023] M. Gembé, H.-J. Schmidt, C. Hickey, J. Richter, Y. Iqbal, and S. Trebst, Phys. Rev. Res. 5, 043204 (2023).
  * Note [1] Specifically, the (anti-/) vortex irreps $\mathsf{B}_{1,2}^{\rm{(a,b)}}$, $\mathsf{A}_{1,2}^{\rm{(a,b)}}$ are promoted by $\mp D$, while $\mathsf{E}^{\rm(a,b)}$ do not depend on $D$.
  * Conlon and Chalker [2009] P. H. Conlon and J. T. Chalker, Phys. Rev. Lett. 102, 237206 (2009).
  * Conlon and Chalker [2010] P. H. Conlon and J. T. Chalker, Phys. Rev. B 81, 224413 (2010).
  * Rehn  _et al._ [2017] J. Rehn, A. Sen, and R. Moessner, Phys. Rev. Lett. 118, 047201 (2017).
  * Note [2] This AFM-vortex topology is homotopically distinct from the known AFM skyrmion[AFMskyrmion], and has not been predicted previously.
  * Senthil  _et al._ [2004] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004).

##  I Supplementary Material

##  II Detailed derivation of the Symmetry properties

Here, we provide further details of the relevant mathematical constructions
that are used in the main text. We start with a system of $\mathcal{N}$ spins.
Much like how one starts in the quantum case with a direct product state basis
to construct exotic entangled states, here we can also start with a many-body
$3\mathcal{N}$-dimensional vector field as a direct sum basis:
$\mathcal{S}=\oplus_{i}^{\mathcal{N}}{\bf S}_{i}$, where ${\bf
S}_{i=1}\in\mathsf{O}(3)$. Then, the most general two-spin interaction
Hamiltonian is written as $H=\mathcal{S}^{T}\mathcal{H}\mathcal{S}$, where
$\mathcal{H}$ is the $3\mathcal{N}\times 3\mathcal{N}$ matrix-valued
Hamiltonian. Short-range interaction and (discrete) translational symmetry
drastically simplifies this Hamiltonian, giving a block-diagonal one.

We assume that there exists a unit cell with sublattices that are invariant
under a point group symmetry $\mathsf{G}$. The spins sitting at the cell
coordinates interact with the spins from the neighboring cells. This
interaction term is translated back to a periodically equivalent interaction
between the spins within the cell. This allows us to define a plaquette
containing $n$ sublattices (counting the sites fully that are shared with the
neighboring cells, and hence, the number of sublattices in a plaquette is
larger than that in a periodic cell). In this prescription, the Hamiltonian
$\mathcal{H}$ becomes block diagonal into a $3n\times 3n$ plaquette
Hamiltonian $\mathcal{H}_{p}$, and the many-body spin vector field splits as
$\mathcal{S}=\oplus_{p=1}^{N=\mathcal{N}/n}\mathcal{S}_{p}$, where
$\mathcal{S}_{p}$ the vector field in the plaquette.

Here, we focus on the square Kagome lattice, which has $n=8$ sites in a
plaquette, giving a 24-dimensional reducible representation $\mathcal{S}_{p}$,
as shown in Fig. X. Our first job is to find the irreducible representation of
the Dihedral group $\mathsf{D}_{4}$ group in this vector field representation.
The group elements are denoted by
$\mathsf{D}_{4}=\\{e,\mathsf{C}_{4},\mathsf{C}_{4}^{2},\mathsf{C}_{4}^{3},\sigma_{v}^{x},\sigma_{v}^{y}=\mathsf{C}_{4}^{-1}\sigma_{v}^{x}\mathsf{C}_{4},\sigma_{v}^{xy},\sigma_{v}^{yx}=\mathsf{C}_{4}^{-1}\sigma_{v}^{xy}\mathsf{C}_{4}\\}$,
where $\mathsf{C}_{4}$ is the four-fold rotation, $\sigma_{v}$ are the
reflection with respect to the verticle plane passing through the $x,y-$ axis,
or diagonal ($xy/yx$), as shown in Fig. X. In this
$\mathcal{S}_{p}$-representation, we can split each of the $\mathsf{D}_{4}$
group elements as successive transformations on how the onsite spin ${\bf
S}_{i}\in\mathsf{O}(3)$ undergoes an internal spin rotation, followed by how
each component $S^{\mu}_{i=1-8}$ of the 8 sublattices reorders in the
plaquette vector $\mathcal{S}_{p}$. Noticeably further, the inner and outer
squares of the square kagome lattice are decoupled from each other in terms of
the $\mathsf{D}_{4}$ symmetries and give a trivial transformation between the
two concentric squares of four sublattices. In what follows, if we denote the
$\mathcal{S}_{p}$-representation of the group elements
$\mathsf{g}\in\mathsf{D}_{4}$ as $\mathcal{D}(\mathsf{g})$, then it can be
decomposed into a direct product of three symmetries:
$\mathcal{D}(\mathsf{g})=\mathcal{R}_{I}(\mathsf{g})\otimes\mathcal{R}_{L}(\mathsf{g})\otimes\mathcal{R}_{S}(\mathsf{g})$,
where $\mathcal{R}_{S}(\mathsf{g})$ are the $3\times 3$ rotational matrices of
the local $\mathsf{O}_{i}(3)$ spin, $\mathcal{R}_{L}(\mathsf{g})$ are the
$4\times 4$ rotational matrices of the four sublattices, and
$\mathcal{R}_{I}(\mathsf{g})$ is the $2\times 2$ transformation between the
inner and outer squares.

$\displaystyle\mathcal{D}(\mathsf{C}_{4})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(1)}(\mathsf{C}_{4})+\tau_{x}\otimes\mathcal{R}_{L}^{(2)}(\mathsf{C}_{4})\right]\otimes\mathcal{R}_{S}(\mathsf{C}_{4}),$
$\displaystyle\mathcal{D}(\mathsf{C}^{2}_{4})$ $\displaystyle=$
$\displaystyle\tau_{x}\otimes\mathbb{I}_{4\times
4}\otimes\mathcal{R}_{S}(\mathsf{C}_{4}^{2}),$
$\displaystyle\mathcal{D}(\mathsf{C}^{3}_{4})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(2)}(\mathsf{C}_{4})+\tau_{x}\otimes\mathcal{R}_{L}^{(1)}(\mathsf{C}_{4})\right]\otimes\mathcal{R}_{S}(\mathsf{C}_{4}^{3}),$
$\displaystyle\mathcal{D}(\sigma_{v}^{x})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(1)}(\sigma_{v}^{x})+\tau_{x}\otimes\mathcal{R}_{L}^{(2)}(\sigma_{v}^{x})\right]\otimes\mathcal{R}_{S}(\sigma_{v}^{x}),$
$\displaystyle\mathcal{D}(\sigma_{v}^{y})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(2)}(\sigma_{v}^{x})+\tau_{x}\otimes\mathcal{R}_{L}^{(1)}(\sigma_{v}^{x})\right]\otimes\mathcal{R}_{S}(\sigma_{v}^{y}),$
$\displaystyle\mathcal{D}(\sigma_{v}^{xy})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(1)}(\sigma_{v}^{xy})+\tau_{x}\otimes\mathcal{R}_{L}^{(2)}(\sigma_{v}^{xy})\right]\otimes\mathcal{R}_{S}(\sigma_{v}^{xy}),$
$\displaystyle\mathcal{D}(\sigma_{v}^{yx})$ $\displaystyle=$
$\displaystyle\left[\tau_{0}\otimes\mathcal{R}_{L}^{(2)}(\sigma_{v}^{xy})+\tau_{x}\otimes\mathcal{R}_{L}^{(1)}(\sigma_{v}^{xy})\right]\otimes\mathcal{R}_{S}(\sigma_{v}^{yx}).$
(8)

Here $\tau_{0},\tau_{x}$ are Pauli matrices defining the internal symmetry
$\mathcal{D}_{4}(\mathsf{g})$, and

$\displaystyle\mathcal{R}_{L}^{(1)}(\mathsf{C}_{4})=\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ 1&0&0&0\\\ 0&1&0&0\\\
\end{pmatrix},\mathcal{R}_{L}^{(2)}(\mathsf{C}_{4})=\begin{pmatrix}0&0&1&0\\\
0&0&0&1\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix},\mathcal{R}_{L}^{(1)}(\sigma_{v}^{x})=\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&1\\\
\end{pmatrix},\mathcal{R}_{L}^{(2)}(\sigma_{v}^{x})=\begin{pmatrix}0&0&1&0\\\
0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&0\\\ \end{pmatrix},$
$\displaystyle\mathcal{R}_{L}^{(1)}(\sigma_{v}^{xy})=\begin{pmatrix}1&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix},\mathcal{R}_{L}^{(2)}(\sigma_{v}^{xy})=\begin{pmatrix}0&0&0&0\\\
0&0&0&1\\\ 0&0&1&0\\\ 0&1&0&0\\\ \end{pmatrix}.$

Under $\mathsf{C}_{4}$, the continuous $\mathsf{O}_{i}(3)$ symmetry simply
becomes a discrete angle of rotation by $2\pi/4$ with $L_{z}$ being the
angular momentum, while under the mirror, spin is rotated as an axial vector.
This gives

$\displaystyle\mathcal{R}_{S}(\mathsf{C}_{4})=\begin{pmatrix}0&-1&0\\\
1&0&0\\\
0&0&1\end{pmatrix},~{}~{}~{}\mathcal{R}_{S}(\sigma_{v}^{x})=\begin{pmatrix}1&0&0\\\
0&-1&0\\\
0&0&-1\end{pmatrix},~{}~{}~{}\mathcal{R}_{S}(\sigma_{v}^{xy})=\begin{pmatrix}0&-1&0\\\
-1&0&0\\\ 0&0&-1\end{pmatrix},$ (9)

and
$\mathcal{R}_{S}(\mathsf{C}_{4}^{2})=(\mathcal{R}_{S}(\mathsf{C}_{4}))^{2}$,
$\mathcal{R}_{S}(\mathsf{C}_{4}^{3})=(\mathcal{R}_{S}(\mathsf{C}_{4}))^{3}$,
$\mathcal{R}_{S}(\sigma_{v}^{y})=\mathcal{R}_{S}(\mathsf{C}_{4})^{-1}\mathcal{R}_{S}(\sigma_{v}^{x})\mathcal{R}_{S}(\mathsf{C}_{4})$,
and
$\mathcal{R}_{S}(\sigma_{v}^{yx})=\mathcal{R}_{S}(\mathsf{C}_{4})^{-1}\mathcal{R}_{S}(\sigma_{v}^{xy})\mathcal{R}_{S}(\mathsf{C}_{4})$.

###  II.1 Symmetry of the Hamiltonian

The generic plaquette Hamiltonian is expressed in the main text as
$H_{p}=\frac{1}{2}\mathcal{S}_{p}^{T}\mathcal{H}_{p}\mathcal{S}_{p}$, where
$\mathcal{H}_{p}$ is the $24\times 24$ symmetric matrix containing all
possible nearest neighbor interactions. The symmetry constraints make many
terms vanish or be identical to other terms. Under a symmetry, the vector
field transforms to
$\mathcal{S}_{p}^{\prime}=\mathcal{D}(\mathsf{g})\mathcal{S}$,
$\forall\mathsf{g}\in\mathsf{D}_{4}$, and if the Hamiltonian to $H_{p}$ is
invariant, then the Hamiltonian matrix transforms as
$\mathcal{D}^{T}(\mathsf{g})\mathcal{H}_{p}\mathcal{D}(\mathsf{g})=\mathcal{H}_{p}$,
$\forall p$.

Under these conditions, we find that the interaction terms among the four
triangles are related to each other by symmetry, while those within a triangle
are independent of each other; see Fig. 1 (a). Consider the one independent
triangle at sites $i=\\{0,1,2\\}$ in Fig. 1 (a), and we obtain three distinct
$3\times 3$ matrices between sites $i$ and $j$ :

$\displaystyle(\mathcal{H}_{p})_{01}=\begin{pmatrix}J^{xx}&D^{xy}&0\\\
D^{yx}&J^{yy}&0\\\
0&0&J^{zz}\end{pmatrix},(\mathcal{H}_{p})_{12}=\begin{pmatrix}J^{xx}&-D^{yx}&0\\\
-D^{xy}&J^{yy}&0\\\ 0&0&J^{zz}\end{pmatrix},\text{ and
}(\mathcal{H}_{p})_{20}=\begin{pmatrix}J^{\prime xx}&D^{\prime xy}&0\\\
-D^{\prime xy}&J^{\prime yy}&0\\\ 0&0&J^{\prime zz}\end{pmatrix}.$ (10)

Therefore, we have nine independent parameters: three exchange interactions
$J^{\mu\mu}$, $J^{\prime\mu\mu}$. and three DM interactions $D^{xy}$,
$D^{yx}$, and $D^{\prime xy}$. Due to in-plane inversion symmetry, no in-plane
DM interaction is allowed. We take a simpler XXZ + DM interaction model in
which $J^{\mu\mu}=J^{\prime\mu\mu}$, $J^{xx}=J^{yy}=J^{zz}/\Delta=J$, and
$D^{xy}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-}D^{yx}=D^{\prime
xy}=JD$. This gives us three independent parameters, among which the global
energy scaling by $J$ is removed, except its sign $\pm$ is considered in the
main text.

###  II.2 Irreducible spin configurations

Finally, we find the irreducible representation of the $\mathcal{S}_{p}$
vector. There are five classes in the group $\mathsf{D}_{4}$ denoted by
$E=\\{\mathsf{e}\\},C_{4}=\\{\mathsf{C}_{4},\mathsf{C}_{4}^{3}\\},C_{2}=\\{\mathsf{C}_{4}^{2}\\},C_{2}^{\prime}=\\{\sigma_{v}^{xy},\sigma_{v}^{yx}\\},C_{2}^{\prime\prime}=\\{\sigma_{v}^{x},\sigma_{v}^{y}\\}$.
The character table for this symmetry group is given in Table S1.

$\mathsf{D}_{4}$ | $d_{\alpha}$ | $E$ | $2C_{4}$ | $2C_{2}^{\prime\prime}$ | $C_{2}$ | $2C_{2}^{\prime}$
---|---|---|---|---|---|---
$\mathsf{A}_{1}$ | 2 | 1 | 1 | 1 | 1 | 1
$\mathsf{A}_{2}$ | 4 | 1 | 1 | -1 | 1 | -1
$\mathsf{B}_{1}$ | 3 | 1 | -1 | 1 | 1 | -1
$\mathsf{B}_{2}$ | 3 | 1 | -1 | -1 | 1 | 1
$\mathsf{E}$ | 6 | 2 | 0 | 0 | -2 | 0
$\mathsf{S}_{p}$ |  | 24 | 0 | -2 | 0 | -2
Table S1: Character table of the group $\mathsf{D}_{4}$. The last row
corresponds to the characters of the reducible representation
$\mathcal{S}_{p}$ for each class. N${}_{\textit{k}}$C${}_{\textit{k}}$ notion
is used in the first row. Nk is the number of elements in each conjugacy
class, Ck.

We have five irreps, which we denote by $m_{\alpha}$ for $\alpha=1-5$. Then
the vector representation of the irreps is a direct sum of the irreps
$\mathcal{M}=\bigoplus_{\alpha}d_{\alpha}m_{\alpha}$ with $d_{\alpha}$ giving
the number of times the $\alpha$-th irrep appears in the sum. $d_{\alpha}$ is
calculated from orthogonality relation with the characters:
$\chi_{m_{\alpha}}(C_{\textit{k}}),\chi_{\mathcal{M}}(C_{\textit{k}})$ of the
24-dimensional representations $m_{\alpha}(C_{k})$, $\mathcal{M}(C_{k})$
respectively, for each conjugacy class $C_{k}$, where $k$ runs over the five
conjugacy classes:

$d_{\alpha}=\frac{1}{h}\sum_{k}N_{k}\chi_{m_{\alpha}}(C_{k})\bar{\chi}_{\mathcal{M}}(C_{\textit{k}})$
(11)

where $h=8$ is the order of the group $\textbf{D}_{4}$, and $N_{\textit{k}}$
is the number of elements in $C_{k}$ conjugacy class. The values of
$d_{\alpha}$ are given in the second column in Table S1.

The final task in this section is to find the basis functions
$\mathcal{V}_{\alpha}$ of each irrep. We denote the basis vectors as
$|\mathcal{V}_{\alpha}^{\mu}\rangle$, where $\alpha=1$ for one-dimensional
irreps, and $\mu=1,2$ (which are relabelled as $x,y$ in Fig. S1) for the two-
dimensional irrep $\mathsf{E}$. The basis vectors follow a relation :
$\mathcal{D}(\mathsf{g})|\mathcal{V}_{\alpha}^{\mu}\rangle=\sum_{\mu^{\prime}}(U_{\alpha}(\mathsf{g}))_{\mu\mu^{\prime}}|\mathcal{V}_{\alpha}^{\mu^{\prime}}\rangle$,
$\forall\mathsf{g}$. $(U_{\alpha}(\mathsf{g}))_{\mu\mu^{\prime}}$ are the
$\mu\times\mu$-matrix for the $\mu$-dimensional irrep $\alpha$ defined for the
group element $\mathsf{g}$. For the one-dimensional irreps $\mathsf{A}_{1,2}$
and $\mathsf{B}_{1,2}$, $U_{\alpha}(\mathsf{g})$ simply gives the character of
the group, and then $|\mathcal{V}_{\alpha}^{\mu}\rangle$ are the simultaneous
eigenvectors of the group elements with the character being the eigenvalue.
They can be solved easily and the corresponding basis functions for the one-
dimensional irreps are shown in Fig. S1(a-d). For the two-dimensional
$\mathsf{E}$ irrep, the orthogonal condition of the basis vector simplifies
the above equation to
$(U_{\alpha}(\mathsf{g}))_{\mu\mu^{\prime}}=\langle\mathcal{V}_{\alpha}^{\mu}|\mathcal{D}(\mathsf{g})|\mathcal{V}_{\alpha}^{\mu^{\prime}}\rangle$.
We solve this matrix for the $\mathsf{E}$ irrep for each group elements, which
comes out to be $U_{\mathsf{E}}(\mathsf{e})=\mathbb{I}_{2\times 2}$,
$U_{\mathsf{E}}(\mathsf{C}_{4})=-i\tau_{y}$,
$U_{\mathsf{E}}(\mathsf{C}_{4}^{2})=-\mathbb{I}_{2\times 2}$,
$U_{\mathsf{E}}(\mathsf{C}_{4}^{3})=i\tau_{y}$,
$U_{\mathsf{E}}(\mathsf{\sigma}_{v}^{x})=\tau_{z}$,
$U_{\mathsf{E}}(\mathsf{\sigma}_{v}^{y})=-tau_{z}$,
$U_{\mathsf{E}}(\mathsf{\sigma}_{v}^{xy})=\tau_{x}$,
$U_{\mathsf{E}}(\mathsf{\sigma}_{v}^{yx})=-\tau_{x}$. $\tau_{\mu}$ are the
$2\times 2$ Pauli matrices.

We have the multiplets as $d_{\alpha}=2,4,3,3$ for the four one-dimensional
irreps $\mathsf{A}_{1}$, $\mathsf{A}_{2}$, $\mathsf{B}_{1}$, $\mathsf{B}_{2}$,
giving 12 basis vectors, while the two-dimensional irrep with multiplicity
$d_{\mathsf{E}}=6$ gives another 12 basis vectors, as shown in Fig. S1(e).
Among them, sixteen are in-plane, defined in the set $\mathsf{m}_{\perp}$, and
eight are out-of-plane, defined in the set $\mathsf{m}_{z}$ in the main text.
Among them, six out-of-plane irreps do not satisfy the local constraint of
$S=1$ per site.

Figure S1: We plot all the irreps’ basis functions. The verticle dashed line
demarcates the out-of-plane irreps $\mathsf{m}_{z}$ on the right-hand side,
among which only the top row satisfies the local constraint while the others
do not. The horizontal arrows dive the spin direction for $S_{i}^{\perp}$,
while the filled and open dots correspond to $S_{i}^{z}$. The size of the dots
corresponds to $|S_{i}^{z}|$. For $\mathsf{A}_{2}^{\rm{(c,d)}}$, the size of
the dots is adjusted for $|S_{i}^{z}|=1$, while for
$\mathsf{B}_{1,2}^{\rm{(c)}}$, sites with symbols give $|S_{i}^{z}|=\sqrt{2}$,
while sites without symbols have $|S_{i}|=1$. Similar consideration is used
for the $\mathsf{E}$ irreps that do not meet the local constraint.

###  II.3 XXZ and DM interactions

In the plaquette Hamiltonian, after substituting
$\mathcal{S}_{p}=\sum_{\alpha=1}^{3n}m_{\alpha}\mathcal{V}_{\alpha}$, we
obtain a Hamiltonian that is block diagonal between the irreps but contains
cross-terms along the multiplicity within an irrep. So we define a
$d_{\alpha}$-dimensional spinor field for each irrep as
$\boldsymbol{{m}}_{\alpha}:=({m}_{\alpha}^{(1)}~{}...~{}{m}_{\alpha}^{(d_{\alpha})})^{T}\in\mathsf{O}_{p}(d_{\alpha})$,
in which the plaquette Hamiltonian splits as

$\displaystyle
H_{p}=\sum_{\alpha=1}^{5}\boldsymbol{{m}}_{\alpha}^{T}\mathcal{H}_{\alpha}\boldsymbol{{m}}_{\alpha},$
(12)

where we have suppressed the plaquette index on the right-hand side.
$\mathcal{H}_{\alpha}$ is a $d_{\alpha}\times d_{\alpha}$ matrix. The
$\mathsf{O}_{p}(d_{\alpha})$ symmetry of each irrep breaks into
$\mathsf{O}_{p}(2)$ and $\mathsf{Z}_{2}$ symmetry as follows.

For $\alpha=1$, the $\mathsf{A}_{1}$ irrep with $d_{1}=2$ multiplets follows
an $\mathsf{O}_{p}(2)$ symmetry.

For $\alpha=2$, the $\mathsf{A}_{2}$ irrep with $d_{2}=4$, we have an emergent
$\mathsf{O}_{p}(2)\times\mathsf{O}_{p}(2)$ symmetry among the multiplets,
giving
$\mathcal{H}_{\mathsf{A}_{2}}=\mathcal{H}_{\mathsf{A}_{2}^{\rm{(a,b)}}}\oplus\mathcal{H}_{\mathsf{A}_{2}^{\rm{(c,d)}}}$.
This is obvious because $\mathsf{A}_{2}^{\rm{(a,b)}}$ consists of coplanar
spins while $\mathsf{A}_{2}^{\rm{(c,d)}}$ are the two out-of-plane spins.

For both $\alpha=3,d$, the $\mathsf{B}_{1,2}$ irreps with $d_{3,4}=3$, we have
an emergent $\mathsf{O}_{p}(2)\times\mathsf{Z}_{2}$ symmetry with
$\mathcal{H}_{\mathsf{B}_{1,2}}=\mathcal{H}_{\mathsf{B}_{1,2}^{\rm{(a,b)}}}\oplus\mathcal{H}_{\mathsf{B}_{1,2}^{\rm{(c)}}}$.
Here, the $\mathsf{B}^{\rm{(a,b)}}_{1,2}$ multiplets are coplanar spins
forming $\mathsf{O}(2)$ symmetry, while $\mathsf{B}^{\rm{(c)}}_{1,2}$ consists
of out-of-plane spins that do not obey local constraints.

For $\alpha=5$, the two-dimensional $\mathsf{E}$ irrep with $d_{5}=6$, each
component of each multiplicity gives emergent $\mathsf{O}_{p}(2)$ rotation as
$\mathcal{H}_{\mathsf{E}}=\mathcal{H}_{\mathsf{E}^{\rm{(a,b)}}}\oplus\mathcal{H}_{\mathsf{E}^{\rm{(c,d)}}}\oplus\mathcal{H}_{\mathsf{E}^{\rm{(e,f)}}}$.

All the $\mathsf{O}_{p}(2)$ invariant $2\times 2$ Hamiltonian matrices for all
irreps have this general form

$\displaystyle(\mathcal{H}_{\alpha})_{k,k^{\prime}}=\epsilon_{\alpha}^{(k+)}\sigma_{0}+\epsilon_{\alpha}^{(k-)}\sigma_{z}+\lambda_{\alpha}^{(kk^{\prime})}\sigma_{x},$
(13)

where $k,k^{\prime}=1,2\in$ (a,b) or (c,d) or (e,f), and
$\epsilon^{k\pm}_{\alpha}=[\epsilon^{(k)}_{\alpha}\pm\epsilon^{(k^{\prime})}_{\alpha}]/2$
and $\epsilon^{(k)}_{\alpha}$ is the onsite energy for the $k^{\rm th}$
multiplet of the $\alpha$-irrep, and $\lambda_{\alpha}^{(kk^{\prime})}$ is the
‘hopping energy’ between the $k$ and $k^{\prime}$ multiples. The onsite
energies of the two vortices with different helicities are
$\epsilon_{\mathsf{A}_{1}^{\rm(a)}}=\epsilon_{\mathsf{A}_{2}^{\rm(a)}}=2\sqrt{2}+2(\sqrt{2}-1)D$,
$\epsilon_{\mathsf{A}_{1}^{\rm(b)}}=\epsilon_{\mathsf{A}_{2}^{\rm(b)}}=-2\sqrt{2}-2(\sqrt{2}+1)D$,
while the energy cost to change the helicity angle is
$\lambda_{\mathsf{A}_{1}^{\rm(a,b)}}=\lambda_{\mathsf{A}_{2}^{\rm(a,b)}}=-4D$.
The same for the two anti-vortices are:
$\epsilon_{\mathsf{B}_{1}^{\rm(a)}}=\epsilon_{\mathsf{B}_{2}^{\rm(a)}}=-2\sqrt{2}+2(\sqrt{2}+1)D$,
$\epsilon_{\mathsf{B}_{1}^{\rm(b)}}=\epsilon_{\mathsf{B}_{2}^{\rm(b)}}=2\sqrt{2}-2(\sqrt{2}-1)D$,
$\lambda_{\mathsf{B}_{1}^{\rm(a,b)}}=\lambda_{\mathsf{B}_{2}^{\rm(a,b)}}=-4D$.
The out-of-plane irreps with parallel and anti-parallel spins and spin-flip
energies between them as $\epsilon_{\mathsf{A}_{2}^{\rm(c)}}=6\Delta$,
$\epsilon_{\mathsf{A}_{2}^{\rm(d)}}=-2\Delta$,
$\lambda_{\mathsf{A}_{2}^{\rm(c,d)}}=4\Delta$. The two irreps with only inner
and out-square out-of-plane spins have the onsite energy:
$\epsilon_{\mathsf{B}_{1}^{\rm(c)}}=\epsilon_{\mathsf{B}_{2}^{\rm(c)}}=-4\Delta$.
Each two-dimensional irreps is degenerate. The in-plane FM $\mathsf{E}$ irreps
have the energies $\epsilon_{\mathsf{E}^{\rm(a)}}=6$,
$\epsilon_{\mathsf{E}^{\rm(b)}}=-2$, and their hopping energy
$\epsilon_{\mathsf{E}^{\rm(a,b)}}=4$. The in-plane AFM $\mathsf{E}$ irreps
have the energies $\epsilon_{\mathsf{E}^{\rm(c)}}=4D-2$,
$\epsilon_{\mathsf{E}^{\rm(d)}}=-4D-2$, and
$\epsilon_{\mathsf{E}^{\rm(c,d)}}=-4$. The two out-of-plane $\mathsf{E}$
irreps that do not mix have the energies
$\epsilon_{\mathsf{E}^{\rm(e)}}=2\sqrt{2}$,
$\epsilon_{\mathsf{E}^{\rm(b)}}=-2\sqrt{2}$. All energies are multiplied with
$J$.

The explicit form of Hamiltonian in terms of the matrix elements in the basis
of the irrep order parameter is

${}\begin{split}H_{p}=\sum_{\alpha=\mathsf{A}_{1,2},\mathsf{B}_{1,2}}\sum_{k,k^{\prime}}(\mathcal{H}_{\alpha})_{k,k^{\prime}}m_{\alpha}^{(k)}m_{\alpha}^{(k^{\prime})}+\sum_{k,k^{\prime}}(\mathcal{H}_{\mathsf{E}})_{k,k^{\prime}}{\bf
m}_{\mathsf{E}}^{(k)}\cdot{\bf
m}_{\mathsf{E}}^{(k^{\prime})}+\sum_{\alpha=\mathsf{B}_{1,2},k={\rm
c}}(\mathcal{H}_{\mathsf{E}})_{k,k}(m_{\alpha}^{k})^{2}.\end{split}$ (14)

where $k,k^{\prime}=$a,b for all irreps, and in addition, we have
$k,k^{\prime}=$c,d for $\mathsf{A}_{2}$ and $k,k^{\prime}=$c,d, and
$k,k^{\prime}=$e,f for $\mathsf{E}$.

Then, for all $\mathsf{O}_{p}(2)$ order parameters, diagonalize the
corresponding $2\times 2$ $\mathcal{H}_{\alpha}$ matrices by the orthogonal
transformation:

$\begin{pmatrix}\tilde{m}_{\alpha}^{(k)}\\\
\tilde{m}_{\alpha}^{(k^{\prime})}\end{pmatrix}=\left[\sigma_{0}\cos\phi_{\alpha}^{(k,k^{\prime})}-i\sigma_{y}\sin\phi_{\alpha}^{(k,k^{\prime})}\right]\begin{pmatrix}m_{\alpha}^{(k)}\\\
m_{\alpha}^{(k^{\prime})}\end{pmatrix}$

where $\phi_{\alpha}^{(k,k^{\prime})}$ is a fixed angle of rotation that
diagonalizes the corresponding irrep multiplets. Eventually, we obtain a fully
diagonal Hamiltonian as

${}H_{p}=\sum_{\nu=(\alpha,k=1,d_{\alpha})}E_{\nu}|\tilde{m}_{\nu}|^{2}.$ (15)

We have abandoned the $\alpha$ and $k$ symbols for the irreps and multiplicity
and combined them into a single symbol $\nu$ which runs from 1 to $3n$ in the
eigenmodes, for simplicity. Here
$E_{\nu}=\epsilon_{\alpha}^{+}\pm\sqrt{(\epsilon_{\alpha}^{-})^{2}+\lambda_{\alpha}^{2}}$
for each $\mathsf{O}_{p}(2)$ multipltes of $\alpha$-irreps. Their explicit
forms are

$\displaystyle E_{\nu=1,2}$ $\displaystyle=$ $\displaystyle-2D\pm
2\sqrt{D^{2}+(1+D)^{2}},\qquad\rm{for}~{}\alpha=\mathsf{A}_{1}^{\rm(a,b)},$
$\displaystyle E_{\nu=3,4}$ $\displaystyle=$ $\displaystyle
E_{\nu=1,2},\qquad\text{\hskip
82.51282pt}\rm{for}~{}\alpha=\mathsf{A}_{2}^{\rm(a,b)},$ $\displaystyle
E_{\nu=5,6}$ $\displaystyle=$ $\displaystyle
2\Delta(1\pm\sqrt{5}),\qquad\text{\hskip
59.75095pt}\rm{for}~{}\alpha=\mathsf{A}_{2}^{\rm(c,d)},$ $\displaystyle
E_{\nu=7,8}$ $\displaystyle=$ $\displaystyle 2D\pm
2\sqrt{D^{2}+2(1-D)^{2}},\text{\hskip
22.76228pt}\rm{for}~{}\alpha=\mathsf{B}_{1}^{\rm(a,b)},$ $\displaystyle
E_{\nu=9}$ $\displaystyle=$ $\displaystyle-4\Delta,\text{\hskip
110.96556pt}\rm{for}~{}\alpha=\mathsf{B}_{1}^{\rm(c)},$ $\displaystyle
E_{\nu=10-12}$ $\displaystyle=$ $\displaystyle E_{\nu=7-9}\text{\hskip
102.43008pt}\rm{for}~{}\alpha=\mathsf{B}_{2}^{\rm(a,b,c)},$ $\displaystyle
E_{\nu=13,14}$ $\displaystyle=$ $\displaystyle 2\pm 2\sqrt{5},\text{\hskip
93.89418pt}\rm{for}~{}\alpha=\mathsf{E}^{\rm(a,b)},$ $\displaystyle
E_{\nu=15,16}$ $\displaystyle=$ $\displaystyle-2\pm
2\sqrt{1+4D^{2}},\text{\hskip
54.06006pt}\rm{for}~{}\alpha=\mathsf{E}^{\rm(c,d)},$ $\displaystyle
E_{\nu=17,18}$ $\displaystyle=$ $\displaystyle\pm 2\sqrt{2}\Delta\text{\hskip
99.58464pt}\rm{for}~{}\alpha=\mathsf{E}^{\rm(e,f)}.$ (16)

All the energies are defined with respect to $J$. The values of the angle
$\phi$ are:

$\displaystyle\phi_{\mathsf{A}_{1}^{\rm(a,b)}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\tan^{-1}\left(\frac{D}{\sqrt{2}(1+D)}\right),\text{\hskip
14.51074pt}\phi_{\mathsf{A}_{2}^{\rm(a,b)}}=\phi_{\mathsf{A}_{1}^{\rm(a,b)}},\text{\hskip
14.51074pt}\phi_{\mathsf{A}_{1}^{\rm(c,d)}}=-\frac{1}{2}\tan^{-1}\left(\frac{1}{2}\right),$
$\displaystyle\phi_{\mathsf{B}_{1}^{\rm(a,b)}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\tan^{-1}\left(\frac{D}{\sqrt{2}(1-D)}\right),\text{\hskip
14.51074pt}\phi_{\mathsf{B}_{2}^{\rm(a,b)}}=\frac{1}{2}\tan^{-1}\left(\frac{-D}{\sqrt{2}(1-D)}\right),$
$\displaystyle\phi_{\mathsf{E}^{\rm(a,b)}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\tan^{-1}\left(\frac{1}{2}\right),\text{\hskip
14.51074pt}\phi_{\mathsf{E}^{\rm(c,d)}}=\frac{1}{2}\tan^{-1}\left(\frac{1}{2D}\right).$
(17)

##  III Details of Classical Monte Carlo

In the classical Monte Carlo calculation, the final temperature is achieved by
annealing from the high temperature at each step with $8\times 10^{5}$ Monte
Carlo steps. The expectation values of the observables are calculated by
taking the average over the last $7\times 10^{5}$ configurations of a total
$8\times 10^{5}$ Monte Carlo steps with system size $N=6L^{2}$, with $L$
number of unit cells. All the static structure factor averages are performed
over system size, $L=20$ at temperature $10^{-3}$. The position vectors of
each sublattice (denoted with indices 0,1, … in Fig. 1(a) of main text) are
taken as considering the origin at the center of the square,

$\displaystyle\delta_{0}=\left(\frac{-1}{4},\frac{-1}{4}\right),\text{\hskip
7.14183pt}\delta_{1}=\left(\frac{1}{4},\frac{-1}{4}\right),\text{\hskip
7.14183pt}\delta_{2}=\left(\frac{1}{4},\frac{1}{4}\right),\text{\hskip
7.14183pt}\delta_{3}=\left(\frac{-1}{4},\frac{1}{4}\right),\text{\hskip
7.14183pt}\delta_{4}=\left(0,\frac{-1}{2}\right),\text{\hskip
7.14183pt}\text{\hskip 7.14183pt}\delta_{5}=\left(\frac{1}{2},0\right)$ (18)

##  IV Structure Factor Plots

In this section, we list the real space spin configurations of all the phases
and their respective structure factors. As defined in the main text, the
different structure factors are

$\displaystyle\chi(\textbf{k})$ $\displaystyle=$ $\displaystyle
1/\mathcal{N}\sum_{i,j}\langle\textbf{S}_{i}\cdot\textbf{S}_{j}\rangle\exp{({{\rm
i}\textbf{k}\cdot(\textbf{r}_{i}-\textbf{r}_{j})})}$
$\displaystyle\chi^{\perp}(\textbf{k})$ $\displaystyle=$ $\displaystyle
1/\mathcal{N}\sum_{i,j}\langle\textbf{S}^{\perp}_{i}\textbf{S}^{\perp}_{j}\rangle\exp{({{\rm
i}\textbf{k}\cdot(\textbf{r}_{i}-\textbf{r}_{j})})}$
$\displaystyle\chi^{z}(\textbf{k})$ $\displaystyle=$ $\displaystyle
1/\mathcal{N}\sum_{i,j}\langle\textbf{S}^{z}_{i}\textbf{S}^{z}_{j}\rangle\exp{({{\rm
i}\textbf{k}\cdot(\textbf{r}_{i}-\textbf{r}_{j})})}$ (19)

Figure S2: The real spin configurations (left panel) and the corresponding
structure factor (right panel) are plotted for various phases for the AFM
coupling $J=+1$. (a) Order phase (red region in the phase diagram) with
staggered anti-vortices between the neighboring sites, showing Bragg-like
peaks at a finite but preferential wavevector. (b) Mixed or fragmented phase
where the inner anti-vortices turn into an AFM-anti-vortex with opposite
$S_{i}^{z}$ components, while $S_{i}^{z}=0$ for the outer anti-vortex. The
$S_{i}^{z}$ values, however, take random values and show disorder features in
the corresponding structure factor without any pinch-point correlation. This
is expected as the inner vortices become decoupled from each other, lacking
any significant correlation between them. (c) A CSL phase (close to the
$\mathsf{Z}_{2}$ CSL phase) showing larger spectral weight the $S_{i}^{z}$
correlation function with pinch-points. (d) The mixed or fragmented phase for
$D<0$ which is similar to the mixed phase for $D>0$ except here vortices
replace the anti-vortices. (e) Orderd phase for D<0, similar to the $D>0$ case
in (a), with vortices replacing anti-vortices. (f) A collinear out-of-plane FM
phase arising in the limit of strong our-of-plane anisotropy term
$\Delta\rightarrow-\infty$.  Figure S3: Similar to Fig. S2, but for the FM
interaction $J=-1$. All three phases shown here are the fragmented phases at
different values of $D$, and $\Delta$, showing pinch-point in the $S_{i}^{z}$
correlation function, but FM ordering in the in-plane component.

##  V Soft-spin Approximation

In this section, we analyze the Hamiltonian in Eq. 4 with ’soft spin’
approximation i.e. spin length constraint (|Si|${}^{2}=1$) is softened from
exact value of 1 per site to the global value of
$\sum_{i}^{N}|\textbf{S}_{i}|^{2}=NS$. Because of the global constraint, we
have a uniform (fixed) chemical potential (Lagrangian multiplier) in the
theory. Then, following Ref. 25, we have diagonalized the Hamiltonian in the
Fourier space of the spin. There, a spin vector is defined per unit cell, not
in the plaquette, which means we have six sublattices as
$\mathcal{S}_{i}$=(S${}_{0}^{x}$, S${}_{0}^{y}$, S${}_{0}^{z}$, S${}_{1}^{x}$,
…, S${}_{5}^{z}$). We Fourier transform the spin vector as
$\mathcal{S}(\textbf{q})=\frac{1}{\sqrt{N}}\sum_{i}\mathcal{S}_{i}e^{-i\textbf{q}\cdot\textbf{r}_{i}}$,
where $\textbf{r}=a\textbf{n}_{1}+b\textbf{n}_{2}$ with integers a, b and unit
vectors $\textbf{n}_{1}=(1,0),\textbf{n}_{2}=(0,1)$. The Hamiltonian is then
diagonal in the momentum space as

$H=\sum_{\bf
q}\mathcal{S}^{T}(\textbf{q})^{T}\mathcal{H}(\textbf{q})\mathcal{S}(\textbf{q}),$
(20)

where $\mathcal{H}(\textbf{q})$ is a $18\times 18$ matrix. We can now
diagonalize the $\mathcal{H}(\textbf{q})$ matrix, which gives the energy
eigenvalues $E_{\nu}({\bf q})$. The lowest energy state is the ground state,
and then we plot a few low-energy excited states in Fig. S4.

We note that the analysis on the Fourier basis leads to a violation of the
local constraint and hence, inconsistency is expected between the real-space
model and the Fourier space one, especially in the spin liquid phase. In the
CSL phase, we find an extremely flat band as the lowest energy state,
suggesting extensive degeneracy as expected here. We see the flat band in all
the mixed phases as well. In addition, the spectrum is gapless in both phases,
with gapless points present at $(\pm\pi,\pm\pi),(\pm\pi,\pm 3\pi)$, and $(\pm
3\pi,)\pm 3\pi)$, as shown in Fig. S4. The band degeneracy, denoted with $d$
in the spectrum at each region is different: $d$=4(2) for $\Delta<1(>1)$,
$d$=6 at $\Delta=1$ in the CSL phase where $D$=0; and $d$=2 for mixed phases
both for $J$ = +1 and -1. Hence, there is no simple positive sum of the
constrainer rule here; the direct matching of singular/non-singular bands to
emergent gauge fields/fragility is not possible.

Figure S4: Energy dispersion of the Hamiltonian $\mathcal{H}({\bf q})$ at four
with re)spective degeneracy of flat bands d, for (a) $\Delta=1.0,D=0.0$ (CSL),
(b) $\Delta=4.0,D=1.0$ (Mixed phase) for $J=+1$ and (c) $\Delta=-2.5,D=0.0$
(d) $\Delta=4.0,D=1.0$ (Mixed phases) for J = -1

As discussed rigorously in the main text, the spin liquids (cyon(/black)
colored phase for J=+1(/-1)) phase has pinch points belonging to the algebraic
class of CSLs with ’emergent’ low-energy gauge field excitations. The mixed
(black-colored phase for $J$=+1) phase has no pinch points and, hence, belongs
to the fragile class of CSLs. All the other ordered phase regions have
dispersive bands.

###  V.1 Finite Magnetic field

Figure S5: Phase diagram at D=0, as a function of h and $\Delta$. For
$\Delta<1$, the phase is a mixed phase and spin liquid for another case. The
mixed phase here is unstable for any finite value of D; the phase becomes
ordered in and out-of-plane for non-zero D value.

The external magnetic field is applied along the $z$-axis to the Hamiltonian,
now written as

$H_{mag}=H_{\rm XXZ-DM}-h\sum_{i}S_{i}^{z}.$ (21)

The phase diagram as a function of $h$ and $\Delta$ is presented in Fig. S5
for $D=0$. A mixed phase of disordered in-plane spins with ordered out-of-
plane components is observed at $D=0$ for $\Delta<1$ with increasing $h$. The
in-plane disordered spins exhibit a coexisting Bragg-like leak at $(0,4\pi)$,
and pinch points at $(\pm\pi,\pm 3\pi)$. The ordering along the $z$ components
is FM type. This phase is unstable for any finite value of $D$. A finite value
of $D$ gives an ordered phase depending on the sign of the $D$ value, where
the in-plane spins form an ordered supercell structure and the out-of-plane
spins are ferromagnetically ordered. As $\Delta>1$, the spins become
disordered both in in-plane and out-of-plane components. This phase also has
pinch-points in the correlation function, indicating power-law correlations.
This phase survives at finite values of $D$. Therefore, we conclude that, by
applying the external magnetic field, the spin liquid phase can be stabilized
in these materials.