# Classical spin models of the windmill lattice and their relevance for
PbCuTe2O6

Anna Fancelli  Helmholtz-Zentrum Berlin für Materialien und Energie, Hahn-
Meitner-Platz 1, 14109 Berlin, Germany  Dahlem Center for Complex Quantum
Systems and Fachbereich Physik, Freie Universität Berlin, Arnimallee 14, 14195
Berlin, Germany  Johannes Reuther  Helmholtz-Zentrum Berlin für Materialien
und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Dahlem Center for
Complex Quantum Systems and Fachbereich Physik, Freie Universität Berlin,
Arnimallee 14, 14195 Berlin, Germany  Department of Physics and Quantum Center
for Diamond and Emergent Materials (QuCenDiEM), Indian Institute of Technology
Madras, Chennai 600036, India  Bella Lake  Helmholtz-Zentrum Berlin für
Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Institut
für Festkörperforschung, Technische Universität Berlin, 10623 Berlin, Germany

###### Abstract

We investigate classical Heisenberg models on the distorted windmill lattice
and discuss their applicability to the spin-$1/2$ spin liquid candidate
PbCuTe2O6. We first consider a general Heisenberg model on this lattice with
antiferromagnetic interactions $J_{n}$ ($n=1,2,3,4$) up to fourth neighbors.
Setting $J_{1}=J_{2}$ (as approximately realized in PbCuTe2O6) we map out the
classical ground state phase diagram in the remaining parameter space and
identify a competition between $J_{3}$ and $J_{4}$ that opens up interesting
magnetic scenarios. Particularly, these couplings tune the ground states from
coplanar commensurate or non-coplanar incommensurate magnetically ordered
states to highly degenerate ground state manifolds with subextensive or
extensive degeneracies. In the latter case, we uncover an unusual classical
spin liquid defined on a lattice of corner sharing octahedra. We then focus on
the particular set of interaction parameters $J_{n}$ that has previously been
proposed for PbCuTe2O6 and investigate the system’s incommensurate magnetic
ground state order and finite temperature multistage ordering mechanism. We
perform extensive finite temperature simulations of the system’s dynamical
spin structure factor and compare it with published neutron scattering data
for PbCuTe2O6 at low temperatures. Our results demonstrate that thermal
fluctuations in the classical model can largely explain the signal
distribution in the measured spin structure factor but we also identify
distinct differences. Our investigations make use of a variety of different
analytical and numerical approaches for classical spin systems, such as
Luttinger-Tisza, classical Monte Carlo, iterative minimization, and molecular
dynamics simulations.

##  I Introduction

Magnetic frustration plays a central role in determining the collective
behavior of interacting spins at low temperatures. In particular, it can
prevent conventional magnetic long-range order and give rise to exotic phases,
with spin liquids perhaps representing the most notable ones [1]. From a
classical point of view, frustrated spin systems are often characterized by a
large ground state degeneracy at zero temperature [2, 3]. A typical situation
where this is the case arises when the lattice consists of clusters of sites,
in such a way that neighboring clusters only have one site in common [4, 5, 6,
7, 8, 9]. This is, for example, realized in the two-dimensional kagome [4] and
in the three-dimensional pyrochlore [5] lattices, where the corner-sharing
units are given by triangles and tetrahedra, respectively. If, additionally,
antiferromagnetic Heisenberg interactions couple all spins within a cluster,
the classical ground states are determined by the condition that the spins
within each cluster have to sum up to zero. Depending on details of the spin
degrees of freedom (number of spin components) and the precise lattice
geometry (number of sites in a cluster and number of clusters per unit cell)
these constraints typically leave sufficient freedom for the relative spin
orientations within a cluster, such that the constraints can be satisfied by
an infinite number of spin configurations [10, 11]. This creates a situation
where, on the one hand, the classical spins are free to fluctuate, but on the
other hand are constrained in their collective behaviors. A system with such a
ground state was called by Villain a cooperative paramagnet [12] and is
commonly identified as a classical spin liquid.

Adding quantum fluctuations enables tunneling between the classical ground
states, such that the new (and now possibly unique) ground state is a
macroscopic superposition of the formerly degenerate classical states, a
situation that is particularly promising for producing a quantum spin liquid
[13, 14, 15]. Because of the rich variety of phenomena arising from this
construction, the lattices hosting corner-sharing units have been of central
interest in the field of frustrated magnetism. The most celebrated examples
are the Heisenberg models on the kagome [16, 17, 18] and pyrochlore [10, 5,
11, 19] lattices which, classically, both have an extensive ground state
degeneracy, that means the number of ground states scales exponentially with
the number of lattice sites. When turning on quantum fluctuations the kagome
Heisenberg model is widely believed to realize a quantum spin liquid [20, 21,
22]. Quantum spin systems on the pyrochlore lattice are, generally, also good
candidates for quantum spin liquids, however, this is well established only in
the case when the Ising model is perturbed by small transverse interactions
[14, 13, 23, 24], while latest results for the pyrochlore Heisenberg
antiferromagnet rather indicate a symmetry broken ground state [25, 26, 27,
28].

Besides the well-known kagome and pyrochlore networks, the distorted windmill
lattice [29, 30, 31] is an alternative and less explored corner-sharing
lattice geometry that, likewise, appears very promising in the context of
possible spin liquid behavior. In fact, the distorted windmill lattice
represents a family of lattices that can be constructed from the hyperkagome
lattice. The latter is a three-dimensional arrangement of corner-sharing
triangles, where each site participates in two triangles and which can be
thought of as a three-dimensional generalization of the kagome lattice.
Interestingly, there exists a particular way of deforming the hyperkagome
lattice – described by one real parameter – such that its point group remains
unchanged. By tuning this parameter, one can realize different corner-sharing
geometries [31, 7]. Moreover, several material realizations are known [32, 33,
34, 35, 36, 37, 38], providing concrete possibilities for obtaining quantum
spin liquids. Among them, $\mathrm{PbCuTe_{2}O_{6}}$ with spin-1/2 Cu2+ ions
has recently attracted the greatest attention. This material realizes a
distorted windmill lattice in which the nearest neighbors form isolated
triangles, while the second neighbors form a hyperkagome network. Since the
theoretical predictions for the Heisenberg interactions indicate that first
and second-neighbor couplings are antiferromagnetic and almost equally strong
[36], the dominant couplings form a network of corner-sharing triangles, where
each spin is shared by three triangles. The classical Heisenberg model for
this lattice (which has recently been dubbed the hyper-hyperkagome model [36])
admits an infinitely large ground state degeneracy, which, however, is only
subextensive, i.e., the number of ground states scales exponentially only in
the linear system size [30, 31]. Couplings beyond second neighbors which are
also present in the Heisenberg Hamiltonian for $\mathrm{PbCuTe_{2}O_{6}}$ lift
this degeneracy [36].

The experimental findings are consistent with a spin liquid phase at low
temperatures. Particularly, thermodynamic probes do not find any signature of
symmetry breaking via a magnetic phase transition [37] and muon spin
relaxation experiments show no signs of static magnetism [38]. Moreover,
inelastic neutron scattering on single crystals shows a broad dispersionless
continuum of magnetic excitations [36] that can be interpreted as resulting
from fractional spinon quasiparticles which are characteristic for quantum
spin liquids. The diffuse spin structure factor measured by neutron scattering
is compatible with either a $U(1)$ gapless or a ${\mathbb{Z}}_{2}$ gapped
quantum spin liquid according to a fermionic parton mean field theory [39].

Inspired by the experimental results on $\mathrm{PbCuTe_{2}O_{6}}$, this work
adds a different perspective on this material by performing theoretical
investigations of the classical version of the model Hamiltonian in Ref. [36].
A first focus in Sec. II is the reexamination of the lattice structure and the
unusual change of connectivity of lattice bonds upon varying the tuning
parameter for the site positions. Section III is dedicated to understanding
the precise role of the interactions $J_{n}$ ($n=1,2,3,4$) up to fourth
neighbors in determining the magnetic ground state. The classical ground state
phase diagram, which has so far only been investigated up to third neighbor
interactions [40], is mapped out including all four interactions and setting
$J_{1}=J_{2}$ (as is approximately realized in $\mathrm{PbCuTe_{2}O_{6}}$).
Surprisingly, the seemingly inconspicuous fourth neighbor coupling is found to
have significant impact on the system’s magnetic properties and can tune the
network of interacting spins towards an interesting and previously unexplored
lattice of corner-sharing octahedra. We identify an extensive ground state
degeneracy in this system giving rise to an unusual type of classical spin
liquid. Moreover, the phase diagram contains regions with subextensive ground-
state degeneracies as well as commensurate and incommensurate ground-state
magnetic orders.

In Sec. IV, we focus on the particular set of couplings that have previously
been proposed to be realized in $\mathrm{PbCuTe_{2}O_{6}}$ [36]. We
investigate in detail the magnetic properties of this system such as the
nature of its incommensurate ground state and the sequence of two finite
temperature phase transitions at which this order builds up. A second focus in
Sec. IV lies on investigating the classical model Hamiltonian for
$\mathrm{PbCuTe_{2}O_{6}}$ in the paramagnetic regime, i.e., above the
ordering transitions. Calculating the system’s dynamical spin structure factor
and comparing it with measured neutron scattering data, we investigate whether
thermal fluctuations can mimic the effects of quantum fluctuations in
$\mathrm{PbCuTe_{2}O_{6}}$. We indeed find that the overall shape of the
simulated spin structure factor agrees well with the measured data, except for
a feature of strong intensity in our simulations that can be associated with
the magnetic long-range order below the critical temperature. These results
shed light on the nature of the observed spin fluctuations in
$\mathrm{PbCuTe_{2}O_{6}}$ and reveal a partial quantum-to-classical
correspondence.

##  II Lattice and Hamiltonian

Table 1: Positions of the twelve atoms in the cubic unit cell (with the lattice constant set to unity) of the distorted windmill lattice, parameterized by the real parameter $y\in\mathds{R}$. For $\mathrm{PbCuTe_{2}O_{6}}$ this value is given by $y=-0.2258$ [36]. All other sites of the lattice are obtained by adding integer multiples of the lattice vectors $\bm{\hat{x}}=(1,0,0)$, $\bm{\hat{y}}=(0,1,0)$, $\bm{\hat{z}}=(0,0,1)$. Sublattice | Position
---|---
1 | (3/4 + $y$, 3/8, 1 - $y$)
2 | (1/2 + $y$, 1/4 - $y$, 7/8)
3 | (5/8, 1/2 - $y$, 3/4 - $y$)
4 | (1/2 - $y$, 3/4 - $y$, 5/8)
5 | (3/4 - $y$, 5/8, 1/2 - $y$)
6 | (7/8, 1/2 + $y$, 1/4 - $y$)
7 | (1 - $y$, 3/4 + $y$, 3/8)
8 | (1/4 - $y$, 7/8, 1/2 + $y$)
9 | (3/8, 1 - $y$, 3/4 + $y$)
10 | (1/4 + $y$, 1/8, $y$)
11 | (1/8, $y$, 1/4 + $y$)
12 | ($y$, 1/4 + $y$, 1/8)

The magnetic behavior of $\mathrm{PbCuTe_{2}O_{6}}$ is determined by
$\mathrm{Cu^{2+}}$ ions with spin $S=1/2$, and these ions form a distorted
windmill lattice [37]. This lattice has twelve atoms in its cubic unit cell,
located at the $12d$ Wyckoff positions of the $P4_{1}32$ space group [29].
These positions depend on a real parameter $y\in\mathds{R}$, as reported in
Table 1. By comparing the crystallography data on $\mathrm{PbCuTe_{2}O_{6}}$
with the $12d$ Wyckoff positions of the $P4_{1}32$ space group, one obtains
$y=-0.2258$ for this compound [36, 39].

The Hamiltonian for $\mathrm{PbCuTe_{2}O_{6}}$ is given by

$H=\sum_{n=1}^{4}J_{n}\sum_{{\langle
i<j\rangle}_{n}}\bm{S}_{i}\cdot\bm{S}_{j},$ (1)

where $\langle\ldots\rangle_{n}$ indicates the sum over the $n^{th}$-nearest
neighbours [36]. The estimates for the interaction strengths $J_{n}$ are
reported in Table 2. Since we focus on the classical model, the spins are
treated as three-component vectors
${\bm{S}}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i})$ with unitary norm
$\lvert{\bm{S}}_{i}\rvert=1$.

Table 2:  Exchange parameters of the Heisenberg Hamiltonian (1) for $\mathrm{PbCuTe_{2}O_{6}}$ from DFT calculations [36].  $J_{1}$ | $J_{2}$ | $J_{3}$ | $J_{4}$
---|---|---|---
1.13 meV | 1.07 meV | 0.59 meV | 0.12 meV

In the distorted windmill lattice, bonds of the same length often form
connected or disconnected arrangements of equilateral triangles. The special
feature of the distorted windmill lattice is that, by tuning the value of $y$,
the number and connectivity of these triangles change. To visualize this, in
Fig. 1(a) we plot the side lengths of the smallest equilateral triangles
created by connecting $n^{th}$-nearest neighbors, as a function of the lattice
parameter $y$. Other bond types, which do not form equilateral triangles are
omitted in Fig. 1(a). Since the change $y\rightarrow y+1$ leaves the lattice
invariant, it is sufficient to consider $y$ in the interval $y\in[0,1]$ only.
The triangles with the smallest side length for a fixed value of $y$ also
correspond to the first nearest neighbours among all bonds. For
$y^{*}<y<1-y^{*}$, $y\neq 1/2$, with $y^{*}=(9-\sqrt{33})/16\simeq 0.2035$,
the first nearest neighbours create isolated triangles, see Fig. 1(b). For
$0\leq y<y^{*}$, $1-y^{*}<y<1$, and $y=1/2$ the first nearest neighbours
create a network of corner-sharing triangles, where each spin is shared by two
triangles, see Fig. 1(c). This lattice has been referred to as the hyperkagome
lattice [7]. At the special points $y=y^{*}$ and $y=1-y^{*}$ the corner-
sharing and isolated triangles have equal size and each spin belongs to three
triangles [30]. At $y=3/8$ and $y=5/8$ the isolated triangles formed by the
first nearest neighbors collapse into one point and the positions in the unit
cell become three times degenerate. The resulting lattice for $y=3/8$ and
$y=5/8$ with four sites per unit cell corresponds to the so-called trillium
lattice [41].

Various discontinuities in the side lengths as a function of $y$ can be
observed in Fig. 1(a). For example, at $y=0.5$ a hyperkagome network is
created with much smaller bond lengths than at slightly smaller or larger
values of $y$. This is explained by the observation that at this point, the
two smallest isolated triangles have equal size and, when combining them,
become a hyperkagome network. The smallest triangles that are truly isolated
then occur with larger bond lengths outside the plotted region. Other
discontinuities in Fig. 1(a) have a similar origin.

Figure 1: (a) Side lengths of the smallest triangles in the distorted windmill
lattice that are created by connecting bonds of the same length, as a function
of the position parameter $y$. Only the two smallest types of isolated
triangles (labeled ‘isolated 1’ and ‘isolated 2’) are shown. In these isolated
triangles, each spin only takes part in one triangle [as shown in (b)], such
that the coordination number is two. The hyperkagome triangles correspond to a
network of corner-sharing triangles [as shown in (c)], with the coordination
number four. At the special points $y=y^{*}$ and $y=1-y^{*}$ with
$y^{*}=(9-\sqrt{33})/16\simeq 0.2035$ the hyperkagome and isolated triangles
have equal size and the system has six nearest neighbors. At $y=3/8$ and
$y=5/8$ the smallest isolated triangles collapse into one point and the
positions in the unit cell become three times degenerate, resulting in a
lattice with four sites per unit cell. The classical nearest neighbor
Heisenberg models with $y=0.125$ and $y=y^{*}$ (dotted gray lines) have
already been studied in Ref. [7] and Refs. [31, 30], respectively. (b)-(e)
First to fourth nearest neighbor bonds of the distorted windmill lattice with
$y=-0.2258$ as realized in $\mathrm{PbCuTe_{2}O_{6}}$. The black and grey dots
correspond to sites inside and outside the cubic unit cell, respectively. (b)
First nearest neighbors form isolated triangles. (c) Second nearest neighbors
create a network of corner-sharing (hyperkagome) triangles. (d) Third nearest
neighbours form chains parallel to the $\bm{\hat{x}}$, $\bm{\hat{y}}$,
$\bm{\hat{z}}$ Cartesian directions. (e) Fourth nearest neighbors form chains
along the body diagonals.

For $y=-0.2258$ realized in $\mathrm{PbCuTe_{2}O_{6}}$, the first and second
nearest neighbours create, respectively, a network of isolated and corner-
sharing triangles [Fig. 1(b), (c)]. This value [red line in Fig. 1(a)] lies
close to the special point $1-y^{*}$. In fact, it was found that $J_{1}\simeq
J_{2}$ (Table 2), reflecting the similar geometrical distance between the
first and second neighbor sites. As a result, each spin interacts almost
equally with six other spins, forming a network of corner-sharing triangles
where each site contributes to three triangles. This network has the same
connectivity as the one at $y=y^{*}$ and $y^{*}-1$ formed by the first nearest
neighbors. Because of the higher connectivity with respect to the hyperkagome
lattice, the present case has been referred to as a hyper-hyperkagome lattice
[36].

The magnetic Hamiltonian for $\mathrm{PbCuTe_{2}O_{6}}$ also contains
antiferromagnetic interactions between the third and fourth nearest neighbors
(Table 2). These form chains along the $\bm{\hat{x}}$, $\bm{\hat{y}}$,
$\bm{\hat{z}}$ Cartesian directions [Fig. 1(d)] and the body diagonals [Fig.
1(e)], respectively. In the next section, we will focus on the role of all
these interactions in determining the magnetic ground state.

##  III Hyper-hyperkagome model and magnetic phase diagram

In this section, we present the $T=0$ classical phase diagram for the
$J_{1}$-$J_{2}$-$J_{3}$-$J_{4}$ model at $J_{1}=J_{2}$, obtained by combining
the Luttinger-Tisza (LT) [42] and the iterative minimization (IM) [43]
techniques. Particularly, we also discuss the role of the $J_{4}$ coupling
which has not been considered in previous works [40]. Details of the LT and IM
methods and their comparison are given, respectively, in the Appendices A and
B. The LT method is an analytic approach where the classical energy is
minimized under the approximation of a partially relaxed spin length
constraint where only the total spin but not necessarily the individual spins
are normalized. The IM approach, on the other hand, consists of minimizing the
classical energy numerically, by aligning the spins to the effective magnetic
field created by the surrounding spins, starting from a random configuration.
As an iterative technique, it is prone to detecting local energy minima
instead of global ones. Our procedure to identify the ground state magnetic
order consists of comparing results from IM and LT. A first insight into the
ground state spin configuration is obtained by the magnetic ordering wave
vector $\bm{Q}$ that minimizes the classical energy within LT. This wave
vector can also be calculated with IM, as it corresponds to the maxima of the
equal-time spin structure factor for spins belonging to the same sublattice,
summed over the sublattices [40]

$\mathcal{S}_{\text{sub}}(\bm{q})=\frac{12}{N}\sum_{I,J}^{L\times L\times
L}\sum_{\alpha=1}^{12}\bm{S}_{\alpha,I}\bm{S}_{\alpha,J}e^{i\bm{q}(\bm{R}_{I}-\bm{R}_{J})},$
(2)

where $I,J$ index the unit cell, $\alpha$ the sublattice and $\bm{R}_{I}$
denotes the real space position of the unit cell $I$. Furthermore, $N$ is the
total number of spins in the simulated system, which is a cube of $L\times
L\times L$ crystallographic unit cells ($N=12L^{3}$). Since
$\mathcal{S}_{\text{sub}}(\bm{q})$ only contains spatial Fourier transforms of
spin correlations on the same sublattice, it describes how spins rotate
between different unit cells but neglects information on how spins are
correlated within a unit cell. Using $\mathcal{S}_{\text{sub}}(\bm{q})$ for
characterizing magnetic order can be useful when the lattice has a large unit
cell. Besides that, it can be directly compared to the LT results, by checking
if $\bm{Q}$ is one of the wave vectors that minimizes the energy within LT.
The quantity $\mathcal{S}_{\text{sub}}(\bm{q})$ should be distinguished from
the equal-time spin structure factor

$\mathcal{S}(\bm{q})=\frac{1}{N}\sum_{i,j}\bm{S}_{i}\bm{S}_{j}e^{i\bm{q}(\bm{r}_{i}-\bm{r}_{j})},$
(3)

which uses the actual site positions $\bm{r}_{i}$ and considers correlations
between all sublattices and which will also be used below.

We consider the set of interactions described in Fig. 1, fixing
$J_{1}=J_{2}=J$ (hyper-hyperkagome case) and letting $J_{3}$ and $J_{4}$ vary
in the range $0\leq J_{3},J_{4}\leq J$. Setting $J_{1}=J_{2}$ is motivated by
the physical situation in $\mathrm{PbCuTe_{2}O_{6}}$ where $y\simeq y^{*}$,
i.e., the isolated nearest neighbor triangles and the hyperkagome triangles
are approximately of the same size, and the same is true for the corresponding
interactions.

Figure 2: Classical zero temperature phase diagram for the Heisenberg model on
the distorted windmill lattice with fixed $J_{1}=J_{2}=J$, as a function of
$J_{3}$ and $J_{4}$. Four different cases can be distinguished. (1) In the
green region where $J_{4}<J_{3}$ the ground state is magnetically ordered with
an incommensurate wave vector $\bm{Q}$. (2) In the red region where
$J_{3}<J_{4}$ the ground state has coplanar magnetic order with
$\bm{Q}_{0}=0$. (3) Along the $J_{3}=J_{4}$ line with the point
$J_{3}=J_{4}=0$ included and the point $J_{3}=J_{4}=J$ excluded, the
degenerate ground states form a sub-extensive (1D) manifold
($\mathcal{M}^{1}$), equal to the one found for the $J_{3}=J_{4}=0$ case [30].
(4) At the point $J_{3}=J_{4}=J$ the degenerate ground states form an
extensive (3D) manifold ($\mathcal{M}^{3}$), as further explained in the main
text. The point associated with PbCuTe2O6 is drawn by approximating $J_{2}$ by
$J_{1}$. It corresponds to the ratios $J_{3}/J_{1}$ and $J_{4}/J_{1}$, with
the values of $J_{i}$ reported in Table 2. Figure 3: (a), (b) Left: Spin
configurations in one unit cell corresponding to (a) the $\bm{Q}_{0}=(0,0,0)$
coplanar order and (b) the $\bm{Q}_{1/3}=2\pi(1/3,1/3,1/3)$ coplanar order.
The coloring of the bonds with interactions $J_{1}$ (blue), $J_{2}$ (orange),
$J_{3}$ (green) and $J_{4}$ (red) matches Fig. 1(b)-(e). Spins with the same
orientation have the same color. Spins with different colors enclose an angle
of 120$\degree$ around each triangle. Right: Corresponding construction of
spin states over the full lattice. The $\bm{Q}_{0}=(0,0,0)$ unit cell
configuration illustrated by red boxes in (a) is repeated in each unit cell.
The spins in the $\bm{Q}_{1/3}=2\pi(1/3,1/3,1/3)$ unit cell undergo a rotation
about the axis perpendicular to the plane they span, by an angle
$\bm{Q}_{1/3}\cdot\bm{R}$ when proceeding to the other unit cells, where
$\bm{R}$ denotes the position of the unit cell. Consequently, there are three
distinct spin configurations in one unit cell, indicated by different colors,
that are arranged over the full lattice as depicted in (b). (c) Possible
ground state configuration for the $J_{3}=J_{4}$ case obtained by mixing the
$\bm{Q}_{0}=(0,0,0)$ (red) and the $\bm{Q}_{1/3}=2\pi(1/3,1/3,1/3)$ unit cells
(from blue to green) with periodic boundary conditions. The obtained
configuration consists of an alternating stacking of $\bm{Q}_{1/3}$ and
$\bm{Q}_{0}$ unit cells along the $(1,1,1)$ direction.

As shown in Fig. 2, the phase diagram contains two different extended phases
and different types of ground state degeneracies at the boundaries, described
in more detail below.

(1) $J_{4}<J_{3}<J$: Incommensurate $\bm{Q}$ order
In the region where $J_{4}<J_{3}<J$ the system orders magnetically in a non-
coplanar spin state with incommensurate wave vector $\bm{Q}$, which varies as
$J_{3}$ and $J_{4}$ vary. The same result was found in Ref. [40] for the case
$J_{4}=0$. This means that, besides a continuous change of $\bm{Q}$, turning
on $J_{4}$ does not change the type of magnetic order, as long as
$J_{4}<J_{3}$ holds. The procedure used to identify the incommensurate order
is explained in the next section, where we discuss the system’s ground state
and finite temperature properties for the set of interactions of
$\mathrm{PbCuTe_{2}O_{6}}$ (black dot in the phase diagram in Fig. 2).

(2) $J_{3}<J_{4}<J$: Coplanar $\bm{Q}=0$ order
In the region where $J_{3}<J_{4}<J$ the system orders magnetically in a
coplanar state with wave vector $\bm{Q}_{0}=0$. This means that the spin
arrangement, as depicted in Fig. 3(a), is the same in each unit cell. Pairs of
spins coupled by $J_{1}$, $J_{2}$, and $J_{4}$ form an angle of 120$\degree$
with each other. On the other hand, the chains formed by $J_{3}$ have a
ferromagnetic spin alignment. As a consequence, increasing the value of
$J_{3}$ increases the level of frustration, as this bond gives a positive
contribution to the energy, see Table 3 for $\bm{Q}_{0}$. This spin
configuration is unique up to global spin rotations. A special subgroup of
such rotations consists of permutations of the three spin directions which
form this state.

(3) $J_{3}=J_{4}\neq J$: One-dimensional degenerate ground state manifold
Along the $J_{3}=J_{4}\neq J$ line, the system is characterized by a sub-
extensive ground state manifold, equal to the one found for the
$J_{3}=J_{4}=0$ case [30], where the number of degenerate ground states scales
exponentially in the linear system size $L$. We will first briefly review the
$J_{3}=J_{4}=0$ case and then show how the results can be generalized to the
$J_{3}=J_{4}\neq J$ case. When $J_{3}=J_{4}=0$, there are two different types
of coplanar ground states, characterized by $\bm{Q}_{0}=0$ and
$\bm{Q}_{1/3}=2\pi(\pm 1/3,\pm 1/3,\pm 1/3)$, where the signs $\pm$ in the
three components of $\bm{Q}_{1/3}$ can be chosen independently. The
$\bm{Q}_{0}=0$ ground state corresponds to the one found in the region (2). An
example for a $\bm{Q}_{1/3}$ state with the choice of signs given by
$\bm{Q}_{1/3}=2\pi(1/3,1/3,1/3)$ is depicted in Fig. 3(b) for one unit cell.
(Note that here and in the following, the term ‘unit cell’ refers to the
crystallographic unit cell, not to the magnetic one.) Other combinations of
signs correspond to other spin configurations not shown. The configurations in
the neighboring unit cells are found by rotating the spins about the axis
perpendicular to the plane in which the spins lie, by an angle given by
$\bm{Q}_{1/3}\cdot\bm{R}$, where $\bm{R}$ denotes the position of the unit
cell.

Table 3:  Bond contributions $c_{n}=\frac{1}{N}\sum_{{\langle i<j\rangle}_{n}}\bm{S}_{i}\cdot\bm{S}_{j}$ for the coplanar $\bm{Q}_{0}$, $\bm{Q}_{1/3}$ and mixed states, defined such that $H=N\sum_{n=1}^{4}J_{n}c_{n}$. The mixed state is obtained by mixing $a$ unit cells with $\bm{Q}_{0}$ order and $b$ unit cells with $\bm{Q}_{1/3}$ order in such a way that $J_{1}$ and $J_{2}$ triangles remain in energetically optimal $120\degree$ spin configurations, see main text for details. As indicated in the last row, when $J_{3}=J_{4}$, the energy contribution from $J_{3}$ and $J_{4}$ bonds given by $\sim J_{3}c_{3}+J_{4}c_{4}=J_{3}(c_{3}+c_{4})$ is independent of the parameters $a$ and $b$ and equal to the $\bm{Q}_{0}$ and $\bm{Q}_{1/3}$ states. |  $\bm{Q}_{0}$ state |  $\bm{Q}_{1/3}$ state | mixed state
---|---|---|---
$c_{1}$ | -0.5 | -0.5 | -0.5
$c_{2}$ | -1 | -1 | -1
$c_{3}$ | 1 | -0.125 | $\frac{a-0.125b}{a+b}$
$c_{4}$ | -0.5 | 0.625 | $\frac{-0.5a+0.625b}{a+b}$
$c_{3}+c_{4}$ | 0.5 | 0.5 | 0.5

These rotations always involve an angle of 120$\degree$. Therefore, for each
combination of signs, there are three possible arrangements of spins in the
unit cell denoted $A$, $B$, $C$ that are cyclically or anti-cyclically
repeated along each Cartesian direction. Particularly, for the signs as given
in $\bm{Q}_{1/3}=2\pi(1/3,1/3,1/3)$, a cyclic progression $A\rightarrow
B\rightarrow C\rightarrow A$ occurs along all three Cartesian directions which
results in stacked layers of different unit cells along the $(111)$ direction
[see Fig. 3(b)], while other sign choices lead to stacks along other body
diagonals.

Both the $\bm{Q}_{0}$ and the $\bm{Q}_{1/3}$ coplanar ground states feature
spins in the isolated and hyperkagome triangles that have angles of
120$\degree$ between each other, see Fig. 2(b) and Fig. 2(c), respectively.
Having the same mutual arrangement of spins in the triangles allows the
construction of an infinite set of ground states, by mixing stacked layers of
$A$, $B$, $C$ unit cells (corresponding to a $\bm{Q}_{1/3}$ order) with
stacked layers of $\bm{Q}_{0}$ order along the respective body diagonal. Such
a composition of unit cells along the body diagonal has to be done under a
precise set of rules, explained in Ref. [30], to ensure that the triangles
shared by neighboring unit cells maintain the 120$\degree$ angles. An example
of a possible ground state configuration is depicted in Fig. 3(c). The
existence of mixed ground states gives rise to an exponentially large ground
state degeneracy $\sim 2^{L}/6$, where the exponent scales linearly with the
system size, see Ref. [30].

Remarkably, all the degenerate ground states in the case $J_{3}=J_{4}=0$
remain ground states when $J_{3}=J_{4}\neq 0$. This can be understood by
examining the energy contributions from the individual $J_{1}$, $J_{2}$,
$J_{3}$ and $J_{4}$ bonds in the $\bm{Q}_{0}$, $\bm{Q}_{1/3}$ and mixed
states, respectively, reported in Table 3. Particularly, when $J_{3}=J_{4}$
the energy contributions from the chain interactions are equal in all three
states and also independent of the precise mixing of $\bm{Q}_{0}$ and
$\bm{Q}_{1/3}$ orders. Furthermore, the corresponding energies agree with the
ones found within IM and LT, indicating that these degenerate states are also
ground states. Thus, the ground state manifold of the $J_{3}=J_{4}<J$ case is
identical to the one in the $J_{3}=J_{4}=0$ case.

A finite-temperature order-by-disorder transition towards the $\bm{Q}_{0}$
state was previously observed in the $J_{3}=J_{4}=0$ case [30]. To investigate
whether an analogous behavior also occurs for $J_{3}=J_{4}\neq 0$ we perform
finite-temperature classical Monte Carlo (MC) simulations of $N=12\times
L\times L\times L$ spins with periodic boundary conditions, interacting
according to the Hamiltonian in Eq. (1). Details of the simulations are given
in Appendix C. We fix $J_{1}=J_{2}$ and change the value of $J_{3}=J_{4}$. The
specific heat as a function of temperature shows a peak for each value of
$J_{3}=J_{4}$, indicating the presence of an order-by-disorder transition
[Fig. 4(a)]. Below the transition the system is always found to be in the
$\bm{Q}_{0}$ state, in analogy with the $J_{3}=J_{4}=0$ case. The temperature
of the order-by-disorder transition decreases as $J_{3}=J_{4}$ increases. This
can be understood by investigating the energy bands from LT which are obtained
by diagonalizing the coupling matrix $J_{ij}$. As shown in Fig. 4(b)
increasing the value of $J_{3}=J_{4}$ causes a flattening of the lowest band
around the minimum (which forms a line in momentum space), making the system
more prone to explore different configurations driven by thermal fluctuations.
As a consequence, the order-by-disorder transition shifts towards lower
temperatures.

Figure 4: (a) Specific heat as a function of temperature obtained by MC
simulations with $L=8$, for $J_{1}=J_{2}=J$ and four different values of
$J_{3}=J_{4}\in\\{0,0.25,0.5,0.75\\}$. For every value of $J_{3}=J_{4}$, the
system undergoes an order-by-disorder transition, as indicated by the peak. By
increasing the value of $J_{3}=J_{4}$, the temperature of the transition
decreases. (b) Energy ${\bm{e}}_{0}$ of the lowest band from the LT method,
along the momentum cut $[h,0,0]$ that contains the energy minimum
${\bm{e}}_{min}$ at $\bm{Q}_{0}=(0,0,0)$. For better comparison of the
behaviors at low energies, the bands are shifted by the respective minimum
$\bm{e}_{min}$ for each value of $J_{3}=J_{4}$. The flattening of the band
with increasing $J_{3}=J_{4}$ causes a downward shift of the order-by-disorder
transition temperature.

(4) $J_{3}=J_{4}=J$: Three-dimensional degenerate ground state manifold
At the point $J_{3}=J_{4}=J$ the system features an extensive ground state
manifold where the number of degenerate ground states scales exponentially in
the number of sites $N$. This is due to the fact that when
$J_{1}=J_{2}=J_{3}=J_{4}=J$, the Hamiltonian can be rewritten as a sum over
spin clusters that have the shape of irregular octahedra,

$H=\frac{J}{2}\sum_{\text{octa}}(\bm{S}_{\text{octa}})^{2}+\text{const.},$ (4)

where $\bm{S}_{\text{octa}}$ is the sum of the spins in each irregular
octahedral cluster. More precisely, in this specific case, each spin is
coupled to ten other spins with the same interaction strength. These ten spins
can be divided into two groups of five spins, where each group together with
the reference spin forms an irregular octahedra. This property holds for each
spin in the lattice, creating thus a network of corner-sharing clusters of six
spins. Within each cluster, all pairs of spins are coupled with the same
interaction $J$. We call these clusters irregular octahedra because they have
six vertices and eight faces, however, the triangles composing the faces are
not equilateral (as they would be in a regular octahedron) even though all
interaction strengths are the same [see Fig. 5(a)]. The underlying lattice
created by the center points of the octahedral clusters is a trillium lattice
with four sites in the cubic unit cell where each site is shared by three
triangles [41], see Fig. 5(b). Due to this property, we will refer to this
network as the dual-trillium lattice.

The present case is similar to the well-known classical Heisenberg model on
the pyrochlore lattice, in which the spins create a network of corner-sharing
tetrahedra. As for the pyrochlore lattice, the origin of the extensive
degeneracy on the dual trillium lattice can be understood from a Maxwellian
counting argument [10, 11]. The number of degrees of freedom can be estimated
as $F=\frac{N_{\text{c}}}{2}(n-1)q$, where $n$ is the number of spin
components, $q$ is the number of spins in a cluster (for the dual trillium
lattice $q=6$), and $N_{\text{c}}$ is the number of clusters. The factor $n-1$
comes from the fact that the spin’s length constraint fixes one of its
components, while the factor $1/2$ is due to the fact that each spin belongs
to two clusters. The number of constraints $\bm{S}_{\text{octa}}=0$ that
determine the ground state manifold is given by $K=nN_{\text{c}}$. If the
constraints are independent the dimension of the ground state manifold is
given by $D=F-K=\frac{N_{\text{c}}}{2}[n(q-2)-q]$, which for $q=6$ is
extensive for Heisenberg spins with $n=3$ (where $D=3N_{\text{c}}$) and for XY
spins with $n=2$ (where $D=N_{\text{c}}$). Interestingly, these estimates for
the dimension of the ground state manifolds are even larger than for the
pyrochlore lattice where $q=4$ and, consequently, $D=N_{\text{c}}$ for
Heisenberg spins and $D=0$ for XY spins.

It is also instructive to compare the magnetic properties of the dual trillium
Heisenberg model and the hyperkagome Heisenberg model, since the latter,
likewise, features an extensive ground state degeneracy [7]. As mentioned
earlier the hyperkagome network is a special type of distorted windmill
lattice where only $J_{2}$ bonds but no other $J_{1}$, $J_{3}$, or $J_{4}$
bonds contribute. In the hyperkagome lattice, the spins create a network of
corner-sharing triangles and the dual lattice that is formed by connecting the
centers of these triangles corresponds to a lattice with two interconnected
trillium lattices, that we call the double-trillium lattice, see Fig. 5(c) and
(d) and Ref. [44]. This lattice has a doubled unit cell with respect to the
trillium lattice, and each site is connected to three sites of the respective
other trillium sublattice (consequently the double-trillium lattice is
bipartite). The classical hyperkagome Heisenberg model is characterized by
dipolar spin correlations, which manifest as pinch points in the equal time
spin structure factor [7]. On the other hand, however, the system shows
partial order at low temperatures $T/J\lesssim 0.001$ [7], where an order-by-
disorder transition drives the system into a phase with coplanar spin
configurations.

The classical dual trillium Heisenberg model obtained at
$J_{1}=J_{2}=J_{3}=J_{4}=J$ shows significant differences compared to the
hyperkagome Heisenberg model. First, the equal time spin structure factor [see
Eq. (3)] does not show any pinch points, neither within the large-$N$
approximation [17] nor within classical MC, see Fig. 6(a) for a comparison of
the equal time spin structure factors in the $[h,h,l]$ plane obtained by both
methods. As explained in Ref. [45], the absence of pinch points can be traced
back to the fact that the parent trillium lattice (that is created by
connecting the centers of the irregular octahedra) is non-bipartite. In
contrast, the double trillium lattice (that is created by connecting the
centers of the hyperkagome triangles) is bipartite. A non-bipartite nature of
a lattice of corner-sharing clusters gives rise to fascinating physical
phenomena, most strikingly, these systems have been proposed to host
“classical $\mathbb{Z}_{2}$ spin liquids” [45], which can be understood as the
classical counterparts of quantum $\mathbb{Z}_{2}$ spin liquids (see also the
recent works in Refs. [46, 47, 48] for general classification schemes for
classical spin liquids). In the classical case, these phases are characterized
by exponentially decaying spin correlations and an absence of pinch points in
the spin structure factor.

Figure 5: (a) Distorted windmill lattice with $y=-0.2258$. The colored bonds
represent the $J_{1}$, $J_{2}$, $J_{3}$, and $J_{4}$ couplings which are all
of equal strength $J$. We refer to this network as the dual-trillium lattice.
All bonds within an irregular octahedron are illustrated with the same color.
(b) Centers of the irregular octahedra plotted in (a). The centers of the
octahedra that share a site are connected by a black line, creating a nearest
neighbor trillium network. The trillium lattice has four sites per unit cell
which are represented by four different colors of the sites. These colors
match the ones used for the octahedral bonds in (a). (c) Hyperkagome lattice
where the colored bonds correspond to the hyperkagome couplings $J_{2}$ shown
in Fig. 1(c). The same color is used for the hyperkagome triangles
corresponding to the same trillium sublattice of the double-trillium lattice
formed by the centers of the triangles plotted in (d). The centers of the
hyperkagome triangles that share a site are connected by a black line in (d),
creating a double-trillium network. The double-trillium lattice is bipartite
and composed of two trillium sublattices highlighted here by the two colors
that match the colors of the corresponding hyperkagome triangles in (c).

Furthermore, the specific heat from MC simulations of the classical dual
trillium Heisenberg model does not show any sign of an order-by-disorder
transition [see Fig. 6(b)].

Figure 6: (a) Equal time spin structure factor [Eq. (3)] of the classical dual
trillium Heisenberg model with $J_{1}=J_{2}=J_{3}=J_{4}=J$ in the $[h,h,l]$
plane at $T/J=0.01$ obtained by means of MC simulations (right), compared with
the one obtained by means of the large-$N$ approach [17] (left). (b) Specific
heat per spin for the same model as a function of temperature obtained by MC
simulations. No peak is observed down to $T/J=0.01$.

This is in agreement with the argument given in Ref. [19], according to which
a system is expected to show an order-by-disorder transition if
$n<(q+2)/(q-2)$. In the present case with six-site clusters ($q=6$) and
Heisenberg spins ($n=3$) this condition is not satisfied. Our MC calculations
predict a specific heat per spin that approaches the value $c_{v}=k_{B}/2$ at
low temperatures. This value can again be explained within the aforementioned
Maxwellian counting argument and is also consistent with the absence of an
order-by-disorder transition in the following way. The $D$ modes in the ground
state manifold (which correspond to flat bands in LT) do not contribute to the
specific heat. The other $F-D=K=nN_{c}$ modes contribute with $k_{B}/2$ to the
specific heat, provided that these are all harmonic modes. Hence the specific
heat per spin is expected to be $c_{v}=nN_{c}k_{B}/(2N)$. Since the number of
sites $N$ and the number of octahedral clusters $N_{c}$ are related by
$N=3N_{c}$ and assuming Heisenberg spins ($n=3$) one indeed reproduces the
value $c_{v}=k_{B}/2$ from MC calculations. For this agreement, it is
essential that all $K$ non-flat modes outside the ground-state manifold are
harmonic. If that was not the case, particularly, if a subset of the $K$ modes
were quartic [4, 11] (which would then contribute $k_{B}/4$ to the specific
heat), these modes would be selected within an order-by-disorder mechanism, in
contradiction with our observations.

The dual trillium lattice constitutes an interesting platform for
investigating novel types of classical spin liquids. We defer a more detailed
investigation of this system to future work and now focus on the relevance of
the distorted windmill lattice for $\mathrm{PbCuTe_{2}O_{6}}$.

##  IV Classical study of $\mathbf{PbCuTe_{2}O_{6}}$

In this section, we focus on the classical version of the Heisenberg
Hamiltonian for $\mathrm{PbCuTe_{2}O_{6}}$ in Eq. (1) with the couplings given
in Table 2 and the site positions from Table 1 with $y=-0.2258$. Besides an
investigation of the magnetic properties of this model, we discuss to which
extent previous experimental results on $\mathrm{PbCuTe_{2}O_{6}}$ can already
be explained on a classical level, despite the $S=1/2$ quantum nature of this
compound. We simulate $N=12\times L\times L\times L$ classical spins with
periodic boundary conditions, using different analytical and numerical
techniques. To characterize the magnetic order at $T=0$, we proceed similarly
to the previous section by applying IM and LT. We then study the finite
temperature behavior by means of MC simulations. Finally, we investigate the
magnetic excitations, i.e., the dynamical spin structure factor in the
paramagnetic regime using molecular dynamics (MD) simulations.

###  IV.1 Static magnetic order at zero and finite temperatures

As indicated in the phase diagram in Fig. 2, the classical Hamiltonian for
$\mathrm{PbCuTe_{2}O_{6}}$ has a magnetically ordered ground state with an
incommensurate wave vector $\bm{Q}$. To identify the incommensurate order by
means of IM we proceeded similarly to Ref. [40], by investigating different
linear system sizes $L$. As reported in Fig. 7(a), the ground state energy and
the corresponding ordering wave vector $\bm{Q}$ depend on the simulated linear
size $L$. The lowest energy is found for $L=11$ and the ordering wave vector
for this system size is given by $\bm{Q}=2\pi(\pm 3/11,\pm 3/11,\pm 3/11)$,
where all combinations of signs occur and yield identical energies. For $L=7$
and $L=14$ we find only slightly larger ground state energies which differ
from the one at $L=11$ only after the fourth decimal digit. The corresponding
ordering wave vector for $L=7$ and $L=14$ is given by $\bm{Q}=2\pi(\pm 2/7,\pm
2/7,\pm 2/7)$. Similar energies for these two wave vectors are expected, since
$2/7$ and $3/11$ only differ by approximately $5\%$. These results are also
consistent with LT, where the ordering wave vector $\bm{Q}$ is found to be
$\bm{Q}=2\pi(\pm\xi,\pm\xi,\pm\xi)$ where $\xi=0.285756(1)$, i.e., very close
to $2/7\simeq 0.285714$. The dependence of these results on the linear size
and the presence of two almost degenerate estimates of the ground state energy
indicate that the exact ground state is not reached by finite-size
simulations. This is the case when $\bm{Q}$ approaches an irrational number,
which corresponds to a ground state characterized by incommensurate magnetic
order.

Next, we study the finite temperature behavior by means of MC simulations. The
energy per spin as a function of temperature exhibits a drop that becomes
steeper with increasing linear size $L$ [Fig. 7(b)], which is a clear
signature of a phase transition. The position of the energy drop, i.e., the
critical temperature $T_{c}$ of the transition again strongly depends on the
linear system size $L$. This size dependence of $T_{c}$ is expected since
different linear sizes realize different magnetic orders. By comparing the
critical temperatures of the linear sizes characterized by the same ordering
wave vector $\bm{Q}$, one recognizes another type of dependence on the linear
size: The critical temperature $T_{c}$ decreases with increasing system sizes,
see, for example, the graphs for $L=4,8,12$ in Fig. 7(b). Such a behavior
could be due to thermal hysteresis effects [30]. To verify this possibility,
we perform simulations starting from a low-temperature configuration and
gradually heating up the system. The curves obtained, respectively, by heating
and cooling the system do not match around the critical temperature,
confirming the presence of a hysteresis effect [Fig. 7(b) inset]. The energy
jump at the critical temperature and the thermal hysteresis effect are
indications that the phase transition is of first order [49].

Because of the strong dependence of the results on the linear system size, it
is difficult to give a precise estimate of the transition temperature.
Suitable system sizes $L$ for quantifying the transition temperature are
$L=11$ and $L=14$ since they yield the lowest ground state energies in IM and,
additionally, $L=14$ is the largest simulated system size that realizes the
ordering wave vector $\bm{Q}=2\pi(\pm 2/7,\pm 2/7,\pm 2/7)$. Thus, we locate
the critical temperature around the transitions observed upon cooling at
$T_{c}=3.3K$ for $L=11$ and $T_{c}=3.25K$ for $L=14$.

Figure 7: (a) Ground state energy per spin $e_{\text{GS}}$ (left) and ordering
wave vector $\bm{Q}=2\pi(\pm m/L,\pm m/L,\pm m/L)$ (right) of the classical
Heisenberg Hamiltonian for $\mathrm{PbCuTe_{2}O_{6}}$ obtained within IM as a
function of the linear system size $L$. The horizontal blue line marks the
optimal wave vector $\bm{Q}=2\pi(\pm\xi,\pm\xi,\pm\xi)$ with $\xi=0.285756(1)$
found by means of LT. (b) Average energy per spin $e$ for different linear
system sizes $L$ as a function of the temperature $T$, obtained by means of MC
simulations. The energy drop indicates a phase transition towards a
magnetically ordered state. The curve for $L=13$ is almost completely covered
by the curve for $L=12$. Inset: Hysteresis loop for $L=11$. The arrow pointing
down (up) indicates that the curve is obtained by decreasing (increasing) the
temperature during a MC simulation. Figure 8: (a) Specific heat for $L=11$ as
a function of temperature for the classical Heisenberg Hamiltonian for
$\mathrm{PbCuTe_{2}O_{6}}$, obtained by means of MC simulations. The peak at
$T_{c}=3.3K$ corresponds to the onset of the magnetic long-range order. The
peak at $T_{c_{2}}=1.05K$ is associated with additional discrete lattice
symmetry breaking, as discussed in the main text. (b)
$\mathcal{S}_{\text{sub}}(\bm{q})$ at $\bm{q}=\bm{Q}_{1}$ and
$\bm{q}=\bm{Q}_{2}$ as a function of the temperature, obtained by means of MC
simulations at finite temperature and with IM at $T=0$. The vectors
$\bm{Q}_{1}$ and $\bm{Q}_{2}$ are two distinct combinations of signs in the
ordering wave vector $\bm{Q}=2\pi(\pm 3/11,\pm 3/11,\pm 3/11)$. The gray line
marks the total intensity
$\sum_{\bm{q}\in\text{HBZ}}\mathcal{S}_{\text{sub}}(\bm{q})=N/2$ where the
summation is only done over half of the first Brillouin zone (HBZ) since the
other half is identical due to
$\mathcal{S}_{\text{sub}}(\bm{q})=\mathcal{S}_{\text{sub}}(-\bm{q})$. This
value $N/2$ results from a sum rule and is independent of the temperature. (c)
$\mathcal{S}_{\text{sub}}(\bm{q})$ in the $[h,k,k]$ plane, obtained from one
MC run at $T=3K$ (left) and $T=0.5K$ (right). At $T_{c_{2}}<T=3K<T_{c}$ a peak
emerges in $\mathcal{S}_{\text{sub}}(\bm{q})$ at $\bm{Q}_{1}=\pm
2\pi(-3/11,3/11,3/11)$. At $T=0.5K<T_{c_{2}}$ a second less intense peak
emerges in $\mathcal{S}_{\text{sub}}(\bm{q})$ at $\bm{Q}_{2}=\pm
2\pi(3/11,3/11,3/11)$.

To study the magnetically ordered low-temperature regime in more detail, we
focus on $L=11$ for the rest of this subsection. According to our previous
results, this value provides the lowest ground state energy among all
simulated system sizes. The specific heat as a function of temperature
displays two peaks, one at $T_{c}=3.3K$ and another, less pronounced peak at
$T_{c_{2}}=1.05K$ [Fig. 8(a)]. The peak at $T_{c}$ corresponds to the position
of the energy drop in Fig. 7(b) and marks the onset of magnetic long-range
order. To investigate the origin of the second peak in the specific heat at
$T_{c_{2}}$ we consider the Fourier transform of the spin-spin correlations
$\mathcal{S}_{\text{sub}}(\bm{q})$ as defined in Eq. (2). (Since
$\mathcal{S}_{\text{sub}}(\bm{q})$ only takes into account the spin-spin
correlations within the same sublattice, it is periodic in momentum space with
respect to the first Brillouin zone.) Particularly insightful is
$\mathcal{S}_{\text{sub}}(\bm{q}=\bm{Q})$ as a function of temperature at the
ordering wave vector $\bm{Q}$, which is $\bm{Q}=2\pi(\pm 3/11,\pm 3/11,\pm
3/11)$ for $L=11$, where we now explicitly distinguish between different
combinations of signs in this vector [see Fig. 8(b)]. The results at $T>0$
($T=0$) in Fig. 8(b) are obtained by means of MC calculations (IM
calculations). We find that the onset of magnetic order at $T_{c}$ is
associated with the selection of a particular sign combination in
$\bm{Q}=2\pi(\pm 3/11,\pm 3/11,\pm 3/11)$ (due to the random initialization
and spin updates within MC, this selection is also random in each MC run).
Denoting this chosen wave vector as $\bm{Q}_{1}$, below $T_{c}$ we observe the
onset of a finite weight in $\mathcal{S}_{\text{sub}}(\bm{q}=\bm{Q}_{1})$.
Note that, since Eq. (2) is invariant with respect to momentum inversion
$\bm{q}\rightarrow-\bm{q}$ the same signal is also found at $-\bm{Q}_{1}$,
i.e.,
$\mathcal{S}_{\text{sub}}(\bm{q}=\bm{Q}_{1})=\mathcal{S}_{\text{sub}}(\bm{q}=-\bm{Q}_{1})$.
Below $T_{c_{2}}$, we find that, additionally,
$\mathcal{S}_{\text{sub}}(\bm{q}=\bm{Q}_{2})$ acquires a finite weight, where
$\bm{Q}_{2}$ corresponds to another selected combination of signs in
$\bm{Q}=2\pi(\pm 3/11,\pm 3/11,\pm 3/11)$. The sum
$\mathcal{S}_{\text{sub}}(\bm{Q}_{1})+\mathcal{S}_{\text{sub}}(\bm{Q}_{2})$
exhibits a small kink at $T_{c_{2}}$ and accounts for almost the total
intensity at $T=0$. This indicates that the changes which the system undergoes
at $T_{c_{2}}$ do not just correspond to a redistribution of signal in
$\mathcal{S}_{\text{sub}}(\bm{q})$ but rather the generation of additional
signal, pointing towards another magnetic ordering transition associated with
the formation of enhanced static magnetic moments. To further illustrate the
presence of the additional magnetic Bragg peak, in Fig. 8(c) we show
$\mathcal{S}_{\text{sub}}(\bm{q})$ in a plane that contains both $\bm{Q}_{1}$
and $\bm{Q}_{2}$, at two different temperatures $T_{c_{2}}<T<T_{c}$ and
$T<T_{c_{2}}$, from a single MC run. The fact that we observe this transition
for all simulated system sizes, indicates that this is an intrinsic feature of
the system, rather than a finite size effect.

The physical picture that emerges is the following: The spins undergo a two-
stage ordering mechanism, in which at $T_{c}$ they first develop magnetic
order characterized by a single ordering wave vector $\bm{Q}_{1}$. At the
lower temperature $T_{c_{2}}$ additional long-range correlations described by
another (symmetry related) wave vector $\bm{Q}_{2}$ set in. This corresponds
to the loss of certain types of fluctuations that are still allowed for
$T_{c_{2}}<T<T_{c}$. Since the presence of two peaks of unequal height at
$\bm{Q}_{1}$ and $\bm{Q}_{2}$ represents a weight distribution in
$\mathcal{S}_{\text{sub}}(\bm{q})$ with lower momentum space symmetries
compared to the presence of just one peak at $\bm{Q}_{1}$, we expect that the
magnetic order that develops at $T_{c_{2}}$ is associated with additional
lattice symmetry breaking. A possible explanation for the occurrence of two
consecutive phase transitions is provided by the LT approach. Within this
method a single-$\bm{Q}$ magnetic ordering in the ground state is allowed, but
with the caveat that the length constraint of individual spins is not
fulfilled. At sufficiently large temperatures (i.e., in the temperature range
$T_{c_{2}}<T<T_{c}$) a single-$\bm{Q}$ state can nevertheless be realized
since thermal fluctuations reduce the static magnetic moment such that the
spin’s length constraint is effectively relaxed. With decreasing temperature
and the loss of thermal fluctuations, however, the fulfillment of the
individual length constraints eventually requires an additional ordering wave
vector $\bm{Q}_{2}$ which leads to the observed second phase transition at
$T_{c_{2}}$. Such a multi-$\bm{Q}$ ordering in the ground state is not common
but can occur when the interactions are particularly complex and generate a
high level of frustration [50, 51].

In summary, the classical Heisenberg model corresponding to
$\mathrm{PbCuTe_{2}O_{6}}$ orders magnetically in a state with incommensurate
$\bm{Q}$. This behavior is very different from the one obtained experimentally
for $\mathrm{PbCuTe_{2}O_{6}}$. In fact, both thermodynamic probes and
inelastic neutron scattering (INS) experiments show no sign of magnetic order
down to $T=0.01K$ [37, 36]. This indicates that quantum effects, completely
neglected in the classical model, play an important role in determining the
behavior of the real material at low temperatures.

###  IV.2 Dynamics

As pointed out in the previous section, the low-temperature magnetic behavior
of the classical Heisenberg model in Eq. (1) differs strongly from the one of
the real material $\mathrm{PbCuTe_{2}O_{6}}$: While the model Hamiltonian
exhibits magnetic order below $T_{c}\simeq 3.3K$, the material shows no sign
of a magnetic phase transition down to $T=0.01K$. These differences are also
reflected in the low-temperature properties of the dynamical spin structure
factor. A magnetically ordered classical spin system gives rise to well-
defined dispersive spin wave features in the spin structure factor, in
contrast to the dispersionless and broad features found in INS experiments on
$\mathrm{PbCuTe_{2}O_{6}}$ [36]. Thus, classical simulations performed at the
same temperature as INS experiments (i.e., below the critical temperature
$T_{c}$ of the classical model), do not reproduce the dynamical spin structure
factors from experiments. On the other hand, in the paramagnetic regime at
$T>T_{c}$ the classical model also shows broad and dispersionless features,
due to the absence of magnetic order caused by thermal fluctuations.
Consequently, in this regime, it is possible to attempt a comparison between
the simulated classical model at suitably chosen temperatures $T>T_{c}$ and
the experimental data. In the following, we will investigate whether thermal
fluctuations at $T>T_{c}$ can indeed mimic the strong quantum fluctuations
observed in $\mathrm{PbCuTe_{2}O_{6}}$ at low temperatures.

We compute the dynamical spin structure factor at $T>T_{c}$ by means of
molecular dynamics (MD) simulations. Due to the absence of magnetic order in
this temperature regime, possible finite size effects resulting from the
incompatibility of the finite system and incommensurate ordering wave vectors
are not present. As a consequence, we can fix the linear system size to a
moderate (i.e., not too large) value $L=8$. For this system size, the critical
temperature is given by $T_{c}=3.25K$. MD simulations consist of numerically
solving the classical equations of motion for the spins, where the starting
spin configuration is taken from a snapshot of a MC run at a given
temperature, when the system has reached thermal equilibrium. Further
numerical details are given in Appendix D. Within MD simulations the dynamical
spin structure factor $\mathcal{S}(\bm{q},\omega)$ is accessible via the
Fourier transform of the time-dependent spin-spin correlation function

$\mathcal{S}(\bm{q},\omega)=\frac{\beta\omega\lvert{\mathcal{F}(\bm{q})}\rvert^{2}}{N\sqrt{N_{t}}}\sum_{n_{t}=0}^{N_{t}}e^{-i\omega
n_{t}\delta
t}\sum_{i,j}^{N}e^{i\bm{q}(\bm{r}_{i}-\bm{r}_{j})}\bm{S}_{i}(t)\bm{S}_{i}(0),$
(5)

where $\bm{r}_{i}$ is the position of site $i$. Furthermore, $N_{t}$ is the
number of time steps and $\delta t$ is the time step size. The factors
$\mathcal{F}(\bm{q})$ and $\beta\omega$ correspond respectively to the form
factor for the $\mathrm{Cu^{2+}}$ ions [52] and a commonly used factor to take
into account the differences between quantum and classical statistics [53].

In the specific case of an isotropic Heisenberg model as the one under exam,
$\mathcal{S}(\bm{q},\omega)$ in Eq.(5) corresponds to the perpendicular
structure factor. The latter is used to compare theoretical calculations with
INS data [54, 53], because it takes into account the impossibility of neutrons
to probe excitations with momentum parallel to the scattering direction. Also
note that, since $\mathcal{S}(\bm{q},\omega)$ contains the actual site
positions $\bm{r}_{i}$ which are non-integer multiples of the lattice vectors
in direct space, this quantity is not periodic in the first Brillouin zone. To
take into account the spin length $S=1/2$ of the $\mathrm{Cu^{2+}}$ ions in
our comparisons between theory and experiment, we rescale the frequency
$\omega$ by a factor of $2$. This can be justified by linear spin wave theory
where the spin wave energies are proportional to $S$ and match the classical
ones from MD simulations (where spins are normalized as $|\bm{S}_{i}|=1$) when
$S=1$. Consequently, to adjust the energy scales of a classical MD simulation
with the energy scales of a spin-wave approach at $S=1/2$, we need to rescale
our MD results according to $\omega\rightarrow\omega/2$.

We compute the dynamical spin structure factor by means of MD simulations for
four different temperatures $T=4K,5K,6K,7K$ in the paramagnetic regime and
compare it with INS data from Ref. [36].

Figure 9: (a) Dynamical spin structure factor from MD in the $[h,k,0]$ plane,
integrated in the energy range $0.415\text{meV}\leq\omega\leq
0.604\text{meV}$, at $T=4K$ and $T=7K$ (top) compared with the one at $T=3K$
(bottom). (b) Dynamical structure factor at $T=0.1K$ from INS experiments in
the $[h,k,0]$ plane, integrated in the energy range
$0.4\text{meV}\leq\omega\leq 0.6\text{meV}$ [36]. (c) Dynamical spin structure
factor from INS and MD as a function of $h$ along the line cut $[h,1.69,0]$
where the signal is integrated in the perpendicular direction $[0,k,0]$
respectively over $1.125\leq k\leq 2$ (MD) and $1.1\lesssim k\lesssim 2$
(INS), as indicated by the orange shaded region in (a) and (b). Subfigures (d)
and (e) show the same as (a) and (b) but now in the $[h,h,l]$ plane. The
colored arrows in the simulated data at $T=3K$ mark the Bragg peaks discussed
in the main text. The experimental data are additionally integrated along the
perpendicular direction $[h,-h,0]$ in the range $-0.1\leq h\leq 0.1$. (f)
Dynamical spin structure factor from INS and MD as a function of $h$ along the
line cut $[h,h,1.69]$ where the MD signal is integrated in the perpendicular
direction $[0,0,l]$ respectively over $1.25\leq l\leq 2.125$ (MD) and
$1.2\lesssim k\lesssim 2.1$ (INS), as indicated by the orange shaded region in
(d) and (e). (g) From left to right: dynamical spin structure factor along the
$[0.3,0.3,l]$ direction and as a function of frequency $\omega$, computed by
means of MD simulations at $T=4K,5K,6K,7K$. The white dashed lines mark the
positions of the Bragg peaks below $T_{c}$ at $(\pm 0.25,\pm 0.25,\pm 1.75)$
[green arrow in (d)] and $(\pm 0.25,\pm 0.25,\pm 3.25)$ [yellow arrow in (d)].
The simulated data are integrated along the perpendicular direction $[h,h,0]$
in the range $0.25\leq h\leq 0.375$ [white shaded area in (d)]. The finite
resolution in energy of the INS data ($\delta E=0.18meV$) is added to the
simulated data by means of a Gaussian broadening. (h) Dynamical spin structure
factor along the $[0.3,0.3,l]$ direction from INS experiments at $T<0.1K$
[36]. The intensity at low frequencies is due to non-magnetic incoherent
background and is not reproduced in MD simulations. In all plots, the
simulated intensity is rescaled such that the maximum matches the one from the
INS experiments.

In Fig. 9(a) and (b) the comparison is shown for the $[h,k,0]$ plane while
Fig. 9(d) and (e) present the comparison in the $[h,h,l]$ plane. The MD
results in these plots are obtained at $T=4K$ and $T=7K$, while the INS data
has been measured at $T<0.1K$. The simulated inelastic spin structure factor
is also shown at $T=3K<T_{c}$ to highlight the changes which the system
undergoes when entering the magnetically ordered phase.

In the $[h,k,0]$ plane [see Fig. 9(a), (b)], the MD simulations are able to
reproduce the ring-shape features found in experiments. However, the regions
of highest intensities on this ring do not occur around $(h,k,0)=(\pm 2,0,0)$
and $(0,\pm 2,0)$, as seen in the experimental data, but are shifted towards
$(\pm 1,\pm 1,0)$. At $T=3K$, i.e., below the critical temperature, the
inelastic structure factor continues to display similar broad and ring-shaped
features. This is because the $[h,k,0]$ plane does not contain magnetic Bragg
peaks and, consequently, the signal at $T=3K$ stems from the fluctuating part
of the spins.

In the $[h,h,l]$ plane [see Fig. 9(d), (e)], the simulated inelastic spin
structure factors at $T=4K$ and $T=7K$ again reproduce the overall ring-like
distribution of signal observed in INS experiments, except for the region
around $(\pm 1,\pm 1,0)$ which appears overemphasized in our simulations, as
already mentioned above. Interestingly, even weaker features outside the ring
are contained in our MD results in the $[h,h,l]$ plane such as two radial
streaks at $l>2$ and a feature of enhanced intensity at $h\simeq 2$. Below
$T_{c}$ the simulated magnetic structure factor shows multiple magnetic Bragg
peaks in the $[h,h,l]$ plane, the one with the highest intensity is located at
$(h,h,l)=(\pm 1.25,\pm 1.25,\pm 0.25)$ [pink arrow in Fig. 9(d)]. This peak
position is consistent with our findings for the system’s ground state order
from IM in Fig. 7(a) where an ordering wave vector $\bm{Q}=2\pi(\pm m/L,\pm
m/L,\pm m/L)$ with $m/L=1/4$ was found for $L=8$. More precisely, the dominant
magnetic Bragg peak at $(\pm 1.25,\pm 1.25,\pm 0.25)$ can be written as
$(\pm(1+m/L),\pm(1+m/L),\pm m/L)$ where the integer parts correspond to
reciprocal lattices vectors. This dominant peak also explains the strong
simulated signal at $T>T_{c}$ in the nearby region around $(h,h,l)=(\pm 1,\pm
1,0)$ which can be interpreted as a molten remnant of the magnetic long-range
order. In contrast, the INS data does not show a particularly strong signal
around $(h,h,l)=(\pm 1.25,\pm 1.25,\pm 0.25)$. This indicates that despite the
obvious similarities between the calculated and the measured data, the spin
structure factor from INS cannot be simply interpreted as displaying smeared
remnants of classical magnetic Bragg peaks.

The plots in Fig. 9(c) and (f) show the same comparison between MD results and
INS data but now along line cuts in the $[h,k,0]$ and $[h,h,l]$ planes as
indicated by the shaded regions in Fig. 9(a), (b) and Fig. 9(d), (e),
respectively. Furthermore, to match the processing of the INS data, the
calculated spin structure factor has been integrated in a direction
perpendicular to the line cut, as indicated by the width of the shaded
regions. These plots demonstrate that in the temperature range $4K\leq T\leq
7K$ the width of simulated features in the spin structure factor approximately
agrees with the typical broadening from quantum fluctuations in the INS data
(overall, however, the changes of our results within this temperature range
are rather small). They also confirm our observation that the best agreement
between MD and INS data occurs in the $[h,h,l]$ plane away from the $(\pm
1,\pm 1,0)$ region, where the simulated results along the investigated line
cut reproduce the experimental data even within error bars [Fig. 9(f)]. On the
other hand, for the plotted line cut in the $[h,k,0]$ plane which contains
contributions from the strong signal at $(\pm 1,\pm 1,0)$, larger deviations
between theory and experiment are observed [Fig. 9(c)].

Finally, in Fig. 9(g) we report the simulated dynamical spin structure factor
in a plane spanned by the frequency $\omega$ (vertical axis) and the
$[0.3,0.3,l]$ momentum direction (horizontal axis). At all simulated
temperatures $T=4K,5K,6K,7K$, the MD results exhibit a rather featureless
region of strong signal. With decreasing temperature, the features in the
simulated data become slightly more distinct and the intensity shifts towards
lower frequencies, reminiscent of spin waves. The wave vectors where the
intensities are largest at low frequencies correspond to the magnetic Bragg
peaks at $(\pm m/L,\pm m/L,\pm(2-m/L))=(\pm 0.25,\pm 0.25,\pm 1.75)$ and $(\pm
m/L,\pm m/L,\pm(3+m/L))=(\pm 0.25,\pm 0.25,\pm 3.25)$ (where $m/L=1/4$ for the
simulated system size $L=8$), see, respectively, green and yellow arrows in
Fig. 9(d). The second peak is clearly less intense than the first peak, also
due to the form factor in $\mathcal{S}(\bm{q},\omega)$ that reduces the
intensity at large wave vectors. Comparing these observations with
experimental results, INS data [Fig. 9(h)] show a similar dispersionless and
smeared feature, whose intensity reaches up to maximal energies
$1.5\text{meV}\lesssim\omega\lesssim 2\text{meV}$, in agreement with our
calculations. The MD simulations are also able to reproduce this maximum
intensity region in the correct momentum space position. More precisely, the
INS data show large intensities at $(\pm 0.3,\pm 0.3,-1.67)$ which is
compatible with the Bragg peak position at $(\pm 0.25,\pm 0.25,\pm 1.75)$
found in the classical simulations below $T_{c}$.

Overall, our classical simulations at finite temperatures correctly reproduce
the main features of the dynamical spin structure factor of
$\mathrm{PbCuTe_{2}O_{6}}$ measured at much lower temperatures and subject to
strong quantum fluctuations. An exception is the region around $(h,k,l)=(\pm
1,\pm 1,0)$ where our calculated signal is comparably large. This region in
momentum space is close to the dominant Bragg peak position observed in the
magnetically ordered low-temperature regime of the classical model. This
indicates that despite the demonstrated similar effects of thermal
fluctuations in the classical model and quantum fluctuations in
$\mathrm{PbCuTe_{2}O_{6}}$, the spin structure factor of
$\mathrm{PbCuTe_{2}O_{6}}$ cannot be purely explained by molten remnants of
classical magnetic long-range order. On the other hand, the signal
distribution in the measured spin structure factor of
$\mathrm{PbCuTe_{2}O_{6}}$ can be consistently interpreted in terms of spin
liquid behavior, as already demonstrated in previous works [36, 39].

##  V Conclusions

In the first part of this work, we investigated the classical $T=0$ phase
diagram of the hyper-hyperkagome $J_{1}=J_{2}$ Heisenberg model with varying
chain interactions $J_{3}$ and $J_{4}$. We identified a rich variety of ground
states in this phase diagram, with a subextensive degenerate manifold along
the $J_{3}=J_{4}$ line, and an extensive manifold for the point
$J_{1}=J_{2}=J_{3}=J_{4}$. We expect that, when extending the phase diagram
for $J_{3},J_{4}>J_{1}=J_{2}$, further subextensively degenerate lines can be
found, which connect to the $J_{1}=J_{2}=J_{3}=J_{4}$ point. In fact, highly
degenerate points are often conjunctions of subextensively degenerate lines in
parameter space [43, 55]. Our preliminary LT calculations, indeed, show that
for the line cuts $J_{1}=J_{2}=J_{3}$, $J_{4}>J_{1}$ and $J_{1}=J_{2}=J_{4}$,
$J_{3}>J_{1}$ the minima of the lowest energy band form a surface in momentum
space.

It remains an open question how the presented phase diagram changes when
quantum fluctuations are taken into account. The latter usually lift classical
ground state degeneracies and stabilize new phases that do not have a
classical counterpart [56, 57, 58]. The experimental results on
$\mathrm{PbCuTe_{2}O_{6}}$ can be considered as a first indication of how
quantum fluctuations affect the classical system. In fact, the results of the
classical study presented in this work strongly differ from the quantum
picture that has emerged from experiments. First, the classical model has a
phase transition at around $T_{c}=3.3K$ towards a magnetically ordered state,
that is absent in the real material. It is worth mentioning that further
experimental works on $\mathrm{PbCuTe_{2}O_{6}}$ identified a ferroelectric
transition at $T\simeq 1K$ [59, 60]. For a rough estimate of whether this
ferromagnetic transition could be related to the magnetic transition we found
at $T_{c}=3.3K$, one needs to adjust the temperature (energy) scales of our
classical simulation with spins normalized to one ($|\bm{S}|=1$) and the
$S=1/2$ degrees of freedom of the $\mathrm{Cu^{2+}}$ magnetic ions. Since
critical temperatures of magnetic transitions are expected to scale as $\sim
S^{2}$, this means that the ferroelectric transition temperature has to be
divided by $S^{2}=1/4$ to compare it with our simulated transition temperature
(assuming that the ferroelectric transition was accompanied by magnetic
order). This would give a critical temperature of $T=4K$, that is higher than
$T_{c}=3.3K$ found for the classical model. If the observed ferroelectric
transition was accompanied by magnetic order occurring at the same transition
temperature, we would expect the classical model to predict a higher critical
temperature than observed in experiments. This is because quantum fluctuations
are expected to lower magnetic transition temperatures. Indeed, there is no
evidence for magnetic long-range order at any temperature in
$\mathrm{PbCuTe_{2}O_{6}}$ and quantum fluctuations have been experimentally
demonstrated to play an important role in suppressing order.

Finally, our simulated dynamical spin structure factors in the paramagnetic
regime show an overall good agreement with experimental data. This indicates
that thermal fluctuations are able to reproduce the broad features in the
experimental dynamical spin structure factor caused by quantum fluctuations.
There are already existing examples for materials that show a good agreement
between the dynamical structure factors from classical simulations and INS
experiments at different temperatures [53, 61]. However, in these cases, the
behavior of the system does not change drastically when quantum fluctuations
are taken into account. More specifically, no phase transition occurs between
the temperature used in the simulations and the one at which the INS
experiments are carried out. The temperature for the classical simulations is
tuned to get the best match between simulated and measured dynamical structure
factors. This usually results in a higher temperature with respect to the
experimental one, in agreement with the fact that the materials are also
subject to quantum fluctuations. In the present case, at the temperature of
the INS experiments, the classical model is magnetically ordered, while the
real material does not show any magnetic phase transition. Nevertheless, the
dynamical spin structure factor in the paramagnetic regime of the classical
model is able to reproduce the main features found by means of INS
experiments. Although it is not possible to draw specific conclusions on the
nature of quantum fluctuations based on this agreement, identifying and
explaining such quantum-to-classical correspondences [62] will remain an
interesting research direction that may give new insights into the quantum
ground states of strongly frustrated spin systems.

###### Acknowledgements.

We thank S. Chillal, A. R. N. Hanna, and R. Moessner for their useful
insights. The simulations were performed on the Tron cluster at the Physics
Department of Freie Universität Berlin. B.L. acknowledges the support of
Deutsche Forschungsgemeinschaft (DFG) through project B06 of SFB 1143:
Correlated Magnetism: From Frustration To Topology (ID 247310070).

##  Appendix A Luttinger-Tsiza (LT)

To apply the LT method we proceed as follows. First we rewrite the Hamiltonian
in Eq. (1) as

$H=\frac{1}{2}\sum_{\alpha,\beta}\sum_{i\in\alpha,j\in\beta}J_{i,j}(\Delta\bm{r}_{i,j})\bm{S}_{i}\bm{S}_{j},$
(6)

where $i,j$ index the sites and $\alpha,\beta=1\ldots 12$ index the
sublattices. Furthermore, $\bm{r}_{i}$ is the position of the site $i$ and
$\Delta\bm{r}_{i,j}=\bm{r}_{i}-\bm{r}_{j}$.

Then we define the Fourier transform of the spins on the individual
sublattices

$\tilde{\bm{S}}_{\alpha}(\bm{q})=\frac{1}{\sqrt{N/12}}\sum_{i\in\alpha}e^{-i\bm{q}\bm{r}_{i}}\bm{S}_{i},$
(7)

and its inverse

$\bm{S}_{i\in\alpha}=\frac{1}{\sqrt{N/12}}\sum_{\bm{q}}e^{i\bm{q}\bm{r}_{i}}\tilde{\bm{S}}(\bm{q}).$
(8)

Inserting Eq. (8) into Eq. (6) we obtain

$H=\sum_{\bm{q}}\sum_{\alpha,\beta}\tilde{J}_{\alpha,\beta}(\bm{q})\tilde{\bm{S}}_{\alpha}(\bm{q})\tilde{\bm{S}}_{\beta}(\bm{q}),$
(9)

where

$\tilde{J}_{\alpha,\beta}(\bm{q})=\frac{1}{2}\sum_{i\in\alpha,j\in\beta}J_{i,j}(\Delta\bm{r}_{i,j})e^{i\bm{q}(\bm{r}_{i}-\bm{r}_{j})},$
(10)

is the Fourier transform of the interaction matrix.

The ground state energy and the magnetic ordering vector correspond
respectively to the minimum in $\bm{q}$ space of the lowest eigenvalue of the
matrix $\tilde{J}_{\alpha,\beta}(\bm{q})$ in Eq. (10) and the wave vector
$\bm{q}$ where this minimum is located. In the case of a non-Bravais lattice
(such as our distorted windmill lattice), the LT method only provides a lower
bound to the ground state energy and the predicted $\bm{Q}$ ordering may not
correspond to the exact ground state order [63]. In fact, this procedure
ensures that the spin length constraint $\lvert{\bm{S}}_{i}\rvert=1$ (the so-
called strong constraint) is satisfied only globally,
$\sum_{i}|\bm{S}_{i}|^{2}=N$ (which is referred to as the weak constraint). As
a consequence, the obtained ground state could be ruled out when the strong
constraint is also imposed.

In this work, we performed the diagonalization of
$\tilde{J}_{\alpha,\beta}(\bm{q})$ numerically: each component of
$\bm{q}=(q_{x},q_{y},q_{z})$ in Eq. (10) is discretised in finite steps
between $-\pi$ and $+\pi$, and for each value of $\bm{q}$ the diagonalization
of the $12\times 12$ matrix in Eq. (10) is carried out.

##  Appendix B Iterative Minimisation (IM)

The IM method consists of finding the classical ground state, based on the
property that in a ground state, all spins are aligned to the effective
magnetic field created by the surrounding spins with which they interact. More
precisely, we start by rewriting the Hamiltonian in Eq. (1) as

$H=-\sum_{i}\bm{B}_{i}^{\text{eff}}\bm{S}_{i},$ (11)

where

$\bm{B}_{i}^{\text{eff}}=-\frac{1}{2}\sum_{j}J_{i,j}\bm{S}_{j}$ (12)

is the effective local magnetic field to which the spin $\bm{S}_{i}$ is
subject, due to the Heisenberg interactions with other spins. The procedure is
simple, we first initialize a finite number of spins, in our case $N=12\times
L\times L\times L$ spins, in a random configuration. We then select a random
site and align it to its effective magnetic field

$\bm{S}_{i}\rightarrow\frac{\bm{B}_{i}^{\text{eff}}}{\lvert\bm{B}_{i}^{\text{eff}}\rvert},$
(13)

where the magnetic field $\bm{B}_{i}^{\text{eff}}$ on the right-hand side is
divided by its norm to keep the spin length equal to one. We iterate this
procedure until the energy per spin does not change up to the $10^{th}$
decimal digit. The $\bm{Q}$ ordering wave vector is obtained by the wave
vector $\bm{q}$ at which $\mathcal{S}_{\text{sub}}(\bm{q})$ in Eq. (2)
calculated from the resulting spin configuration has its maximum.

Coplanar configurations are characterized by the nematic order parameter [18]

$Q^{\nu,\mu}=\frac{1}{N}\sum_{i}\left(\bm{S}^{\nu}_{i}\bm{S}^{\mu}_{i}-\frac{1}{3}\delta^{\nu,\mu}\right),$
(14)

where $\nu,\mu={x,y,z}$. In particular, the trace of the second moment of (14)

$\sum_{\mu,\nu}Q^{\mu,\nu}Q^{\nu,\mu}=\frac{1}{N^{2}}\sum_{i,j}\left((\bm{S}_{i}\bm{S}_{j})^{2}-\frac{1}{3}\right),$
(15)

is equal to $1/6$ in a coplanar configuration [18]. This quantity has been
used to determine whether a spin configuration is coplanar or not.

We remark that IM does not have any control over whether it reaches the real
ground state or a metastable state. Thus, one has to perform different
attempts with the same linear system size $L$, and then check for convergence
to the same result when varying the linear size $L$.

##  Appendix C Classical Monte Carlo (MC)

We simulated the Heisenberg model given in Eq. (1) on the distorted windmill
lattice with periodic boundary conditions using different linear sizes $L$,
which correspond to a total number of spins $N=12\times L\times L\times L$. A
single MC step includes $N$ local heat-bath moves [64], followed by $3$
overrelaxation moves [65]. Each run is initialized with a random spin
configuration and gradually cooled down from $T\simeq 2J$ to $T\simeq 0.02J$,
where $J$ corresponds to the largest coupling. Physical quantities are
averaged over $2\cdot 10^{5}$ MC steps, after $3\cdot 10^{5}$ MC steps for
thermalization.

##  Appendix D Molecular Dynamics (MD)

We performed molecular dynamics simulations to compute the Fourier transform
of the time-dependent spin-spin correlation function
$\mathcal{S}(\bm{q},\omega)$ defined in Eq. (5) at a fixed temperature. The
starting configuration is taken from the MC simulations. The spins are then
evolved in time according to the equation of motion

$\frac{d\bm{S}_{i}}{dt}=\bm{S}_{i}\times\bm{B}_{i}^{\text{eff}},$ (16)

where $\bm{B}_{i}^{\text{eff}}$ is defined in Eq. (12). This equation
describes the precession of the spins around their local effective field. We
numerically integrate Eq. (16) with a fourth-order Runge-Kutta method. This
procedure is repeated for different starting configurations at a fixed
temperature and then averaged over 50 independent runs. To compute
$\mathcal{S}(\bm{q},\omega)$, we employed the approach described in Ref. [53].
This consists of accumulating the time Fourier transform of each spin
component during the simulation, rather than computing the time-dependent
correlations and Fourier transform afterward.

## References

  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Ramirez  _et al._ [1999] A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. Shastry, Zero-point entropy in ‘spin ice’, Nature 399, 333 (1999).
  * Chalker [2011] J. T. Chalker, Geometrically frustrated antiferromagnets: Statistical mechanics and dynamics, in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 3–22.
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Hidden order in a frustrated system: Properties of the Heisenberg Kagomé antiferromagnet, Phys. Rev. Lett. 68, 855 (1992).
  * Reimers [1992] J. N. Reimers, Absence of long-range order in a three-dimensional geometrically frustrated antiferromagnet, Phys. Rev. B 45, 7287 (1992).
  * Petrenko and McK. Paul [2000] O. A. Petrenko and D. McK. Paul, Classical Heisenberg antiferromagnet on a garnet lattice: A Monte Carlo simulation, Phys. Rev. B 63, 024409 (2000).
  * Hopkinson  _et al._ [2007] J. M. Hopkinson, S. V. Isakov, H.-Y. Kee, and Y. B. Kim, Classical Antiferromagnet on a Hyperkagome Lattice, Phys. Rev. Lett. 99, 037201 (2007).
  * Henley [2010] C. L. Henley, The “Coulomb Phase” in Frustrated Systems, Annual Review of Condensed Matter Physics 1, 179 (2010).
  * Rehn  _et al._ [2016] J. Rehn, A. Sen, K. Damle, and R. Moessner, Classical Spin Liquid on the Maximally Frustrated Honeycomb Lattice, Phys. Rev. Lett. 117, 167201 (2016).
  * Reimers  _et al._ [1991] J. N. Reimers, A. J. Berlinsky, and A.-C. Shi, Mean-field approach to magnetic ordering in highly frustrated pyrochlores, Phys. Rev. B 43, 865 (1991).
  * Moessner and Chalker [1998a] R. Moessner and J. T. Chalker, Properties of a Classical Spin Liquid: The Heisenberg Pyrochlore Antiferromagnet, Phys. Rev. Lett. 80, 2929 (1998a).
  * Villain [1979] J. Villain, Insulating spin glasses, Zeitschrift für Physik B Condensed Matter 33, 31 (1979).
  * Hermele  _et al._ [2004] M. Hermele, M. P. A. Fisher, and L. Balents, Pyrochlore photons: The $U(1)$ spin liquid in a $S=\frac{1}{2}$ three-dimensional frustrated magnet, Phys. Rev. B 69, 064404 (2004).
  * Benton  _et al._ [2012] O. Benton, O. Sikora, and N. Shannon, Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice, Phys. Rev. B 86, 075154 (2012).
  * Shannon  _et al._ [2012] N. Shannon, O. Sikora, F. Pollmann, K. Penc, and P. Fulde, Quantum Ice: A Quantum Monte Carlo Study, Phys. Rev. Lett. 108, 067204 (2012).
  * Huse and Rutenberg [1992] D. A. Huse and A. D. Rutenberg, Classical antiferromagnets on the Kagomé lattice, Phys. Rev. B 45, 7536 (1992).
  * Garanin and Canals [1999] D. A. Garanin and B. Canals, Classical spin liquid: Exact solution for the infinite-component antiferromagnetic model on the kagomé lattice, Phys. Rev. B 59, 443 (1999).
  * Zhitomirsky [2008] M. E. Zhitomirsky, Octupolar ordering of classical kagome antiferromagnets in two and three dimensions, Phys. Rev. B 78, 094423 (2008).
  * Moessner and Chalker [1998b] R. Moessner and J. T. Chalker, Low-temperature properties of classical geometrically frustrated antiferromagnets, Phys. Rev. B 58, 12049 (1998b).
  * Sachdev [1992] S. Sachdev, Kagomé- and triangular-lattice Heisenberg antiferromagnets: Ordering from quantum fluctuations and quantum-disordered ground states with unconfined bosonic spinons, Phys. Rev. B 45, 12377 (1992).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Spin-Liquid Ground State of the $S=1/2$ Kagome Heisenberg Antiferromagnet, Science 332, 1173 (2011).
  * Han  _et al._ [2012] T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet, Nature 492, 406 (2012).
  * Taillefumier  _et al._ [2017] M. Taillefumier, O. Benton, H. Yan, L. D. C. Jaubert, and N. Shannon, Competing Spin Liquids and Hidden Spin-Nematic Order in Spin Ice with Frustrated Transverse Exchange, Phys. Rev. X 7, 041057 (2017).
  * Benton  _et al._ [2018] O. Benton, L. D. C. Jaubert, R. R. P. Singh, J. Oitmaa, and N. Shannon, Quantum Spin Ice with Frustrated Transverse Exchange: From a $\pi$-Flux Phase to a Nematic Quantum Spin Liquid, Phys. Rev. Lett. 121, 067201 (2018).
  * Schäfer  _et al._ [2022] R. Schäfer, B. Placke, O. Benton, and R. Moessner, Abundance of hard-hexagon crystals in the quantum pyrochlore antiferromagnet (2022), arXiv:2210.07235 [cond-mat.str-el] .
  * Hagymási  _et al._ [2021] I. Hagymási, R. Schäfer, R. Moessner, and D. J. Luitz, Possible Inversion Symmetry Breaking in the $S=1/2$ Pyrochlore Heisenberg Magnet, Phys. Rev. Lett. 126, 117204 (2021).
  * Astrakhantsev  _et al._ [2021] N. Astrakhantsev, T. Westerhout, A. Tiwari, K. Choo, A. Chen, M. H. Fischer, G. Carleo, and T. Neupert, Broken-Symmetry Ground States of the Heisenberg Model on the Pyrochlore Lattice, Phys. Rev. X 11, 041021 (2021).
  * Hering  _et al._ [2022] M. Hering, V. Noculak, F. Ferrari, Y. Iqbal, and J. Reuther, Dimerization tendencies of the pyrochlore Heisenberg antiferromagnet: A functional renormalization group perspective, Phys. Rev. B 105, 054426 (2022).
  * Nakamura  _et al._ [1997a] H. Nakamura, K. Yoshimoto, M. Shiga, M. Nishi, and K. Kakurai, Strong antiferromagnetic spin fluctuations and the quantum spin-liquid state in geometrically frustrated $\beta-\mathrm{Mn}$, and the transition to a spin-glass state caused by non-magnetic impurity, Journal of Physics: Condensed Matter 9, 4701 (1997a).
  * Isakov  _et al._ [2008] S. V. Isakov, J. M. Hopkinson, and H.-Y. Kee, Fate of partial order on trillium and distorted windmill lattices, Phys. Rev. B 78, 014404 (2008).
  * Canals and Lacroix [2000] B. Canals and C. Lacroix, Mean-field study of the disordered ground state in the $\beta-\mathrm{Mn}$ lattice, Phys. Rev. B 61, 11251 (2000).
  * Schiffer  _et al._ [1995] P. Schiffer, A. P. Ramirez, D. A. Huse, P. L. Gammel, U. Yaron, D. J. Bishop, and A. J. Valentino, Frustration Induced Spin Freezing in a Site-Ordered Magnet: Gadolinium Gallium Garnet, Phys. Rev. Lett. 74, 2379 (1995).
  * Nakamura  _et al._ [1997b] H. Nakamura, K. Yoshimoto, M. Shiga, M. Nishi, and K. Kakurai, Strong antiferromagnetic spin fluctuations and the quantum spin-liquid state in geometrically frustrated $\beta-\mathrm{Mn}$, and the transition to a spin-glass state caused by non-magnetic impurity, Journal of Physics: Condensed Matter 9, 4701 (1997b).
  * Okamoto  _et al._ [2007] Y. Okamoto, M. Nohara, H. Aruga-Katori, and H. Takagi, Spin-Liquid State in the $S=\frac{1}{2}$ Hyperkagome Antiferromagnet ${\mathrm{Na}}_{4}{\mathrm{Ir}}_{3}{\mathrm{O}}_{8}$, Phys. Rev. Lett. 99, 137207 (2007).
  * Paddison  _et al._ [2013] J. A. M. Paddison, J. R. Stewart, P. Manuel, P. Courtois, G. J. McIntyre, B. D. Rainford, and A. L. Goodwin, Emergent Frustration in Co-doped $\beta$-Mn, Phys. Rev. Lett. 110, 267207 (2013).
  * Chillal  _et al._ [2020] S. Chillal, Y. Iqbal, H. O. Jeschke, J. A. Rodriguez-Rivera, R. Bewley, P. Manuel, D. Khalyavin, P. Steffens, R. Thomale, A. Islam, _et al._ , Evidence for a three-dimensional quantum spin liquid in ${\mathrm{PbCuTe}}_{2}{\mathrm{O}}_{6}$, Nature communications 11, 1 (2020).
  * Koteswararao  _et al._ [2014] B. Koteswararao, R. Kumar, P. Khuntia, S. Bhowal, S. K. Panda, M. R. Rahman, A. V. Mahajan, I. Dasgupta, M. Baenitz, K. H. Kim, and F. C. Chou, Magnetic properties and heat capacity of the three-dimensional frustrated $S=\frac{1}{2}$ antiferromagnet ${\mathrm{PbCuTe}}_{2}{\mathrm{O}}_{6}$, Phys. Rev. B 90, 035141 (2014).
  * Khuntia  _et al._ [2016] P. Khuntia, F. Bert, P. Mendels, B. Koteswararao, A. V. Mahajan, M. Baenitz, F. C. Chou, C. Baines, A. Amato, and Y. Furukawa, Spin Liquid State in the 3D Frustrated Antiferromagnet ${\mathrm{PbCuTe}}_{2}{\mathrm{O}}_{6}$: NMR and Muon Spin Relaxation Studies, Phys. Rev. Lett. 116, 107203 (2016).
  * Chern and Kim [2021] L. E. Chern and Y. B. Kim, Theoretical study of quantum spin liquids in $S=\frac{1}{2}$ hyper-hyperkagome magnets: Classification, heat capacity, and dynamical spin structure factor, Phys. Rev. B 104, 094413 (2021).
  * Jin and Zhou [2020] H.-K. Jin and Y. Zhou, Classical and quantum order in hyperkagome antiferromagnets, Phys. Rev. B 101, 054408 (2020).
  * Hopkinson and Kee [2006] J. M. Hopkinson and H.-Y. Kee, Geometric frustration inherent to the trillium lattice, a sublattice of the B20 structure, Phys. Rev. B 74, 224441 (2006).
  * Luttinger and Tisza [1946] J. M. Luttinger and L. Tisza, Theory of Dipole Interaction in Crystals, Phys. Rev. 70, 954 (1946).
  * Sklan and Henley [2013a] S. R. Sklan and C. L. Henley, Nonplanar ground states of frustrated antiferromagnets on an octahedral lattice, Phys. Rev. B 88, 024407 (2013a).
  * Živković  _et al._ [2021] I. Živković, V. Favre, C. Salazar Mejia, H. O. Jeschke, A. Magrez, B. Dabholkar, V. Noculak, R. S. Freitas, M. Jeong, N. G. Hegde, L. Testa, P. Babkevich, Y. Su, P. Manuel, H. Luetkens, C. Baines, P. J. Baker, J. Wosnitza, O. Zaharko, Y. Iqbal, J. Reuther, and H. M. Rønnow, Magnetic Field Induced Quantum Spin Liquid in the Two Coupled Trillium Lattices of ${\mathrm{K}}_{2}{\mathrm{Ni}}_{2}({\mathrm{SO}}_{4}{)}_{3}$, Phys. Rev. Lett. 127, 157204 (2021).
  * Rehn  _et al._ [2017] J. Rehn, A. Sen, and R. Moessner, Fractionalized ${\mathbb{Z}}_{2}$ Classical Heisenberg Spin Liquids, Phys. Rev. Lett. 118, 047201 (2017).
  * 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 .
  * 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 .
  * Davier  _et al._ [2023] N. Davier, F. G. Albarracín, D. Rosales, and P. Pujol, A combined approach to analyze and classify families of classical spin liquids (2023), arXiv:2304.10906 .
  * Zheng and Zhang [1998] G. P. Zheng and J. X. Zhang, Thermal hysteresis scaling for first-order phase transitions, Journal of Physics: Condensed Matter 10, 275 (1998).
  * Okubo  _et al._ [2011] T. Okubo, T. H. Nguyen, and H. Kawamura, Cubic and noncubic multiple-$q$ states in the Heisenberg antiferromagnet on the pyrochlore lattice, Phys. Rev. B 84, 144432 (2011).
  * Shimokawa  _et al._ [2019] T. Shimokawa, T. Okubo, and H. Kawamura, Multiple-$q$ states of the ${J}_{1}-{J}_{2}$ classical honeycomb-lattice Heisenberg antiferromagnet under a magnetic field, Phys. Rev. B 100, 224404 (2019).
  * Brown [2006] P. Brown,  _Magnetic form factors_ (Wiley Online Library, 2006).
  * Zhang  _et al._ [2019] S. Zhang, H. J. Changlani, K. W. Plumb, O. Tchernyshyov, and R. Moessner, Dynamical Structure Factor of the Three-Dimensional Quantum Spin Liquid Candidate ${\mathrm{NaCaNi}}_{2}{\mathrm{F}}_{7}$, Phys. Rev. Lett. 122, 167203 (2019).
  * Pohle  _et al._ [2017] R. Pohle, H. Yan, and N. Shannon, How many spin liquids are there in Ca10Cr7O28? (2017), arXiv:1711.03778 [cond-mat.str-el] .
  * Balla  _et al._ [2020] P. Balla, Y. Iqbal, and K. Penc, Degenerate manifolds, helimagnets, and multi-$Q$ chiral phases in the classical Heisenberg antiferromagnet on the face-centered-cubic lattice, Phys. Rev. Research 2, 043278 (2020).
  * Mambrini  _et al._ [2006] M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, Plaquette valence-bond crystal in the frustrated Heisenberg quantum antiferromagnet on the square lattice, Phys. Rev. B 74, 144422 (2006).
  * Yamamoto  _et al._ [2014] D. Yamamoto, G. Marmorini, and I. Danshita, Quantum Phase Diagram of the Triangular-Lattice $XXZ$ Model in a Magnetic Field, Phys. Rev. Lett. 112, 127203 (2014).
  * Iqbal  _et al._ [2019] Y. Iqbal, T. Müller, P. Ghosh, M. J. P. Gingras, H. O. Jeschke, S. Rachel, J. Reuther, and R. Thomale, Quantum and Classical Phases of the Pyrochlore Heisenberg Model with Competing Interactions, Phys. Rev. X 9, 011005 (2019).
  * Hanna  _et al._ [2021] A. R. N. Hanna, A. T. M. N. Islam, R. Feyerherm, K. Siemensmeyer, K. Karmakar, S. Chillal, and B. Lake, Crystal growth, characterization, and phase transition of ${\mathrm{PbCuTe}}_{2}{\mathrm{O}}_{6}$, Phys. Rev. Materials 5, 113401 (2021).
  * Thurn  _et al._ [2021] C. Thurn, P. Eibisch, A. Ata, M. Winkler, P. Lunkenheimer, I. Kézsmárki, U. Tutsch, Y. Saito, S. Hartmann, J. Zimmermann, A. R. N. Hanna, A. T. M. N. Islam, S. Chillal, B. Lake, B. Wolf, and M. Lang, Spin liquid and ferroelectricity close to a quantum critical point in ${\mathrm{PbCuTe}}_{2}{\mathrm{O}}_{6}$, npj Quantum Materials 6, 95 (2021).
  * Pohle  _et al._ [2021] R. Pohle, H. Yan, and N. Shannon, Theory of ${\mathrm{Ca}}_{10}{\mathrm{Cr}}_{7}{\mathrm{O}}_{28}$ as a bilayer breathing-kagome magnet: Classical thermodynamics and semiclassical dynamics, Phys. Rev. B 104, 024426 (2021).
  * Wang  _et al._ [2020] T. Wang, X. Cai, K. Chen, N. V. Prokof’ev, and B. V. Svistunov, Quantum-to-classical correspondence in two-dimensional Heisenberg models, Phys. Rev. B 101, 035132 (2020).
  * Sklan and Henley [2013b] S. R. Sklan and C. L. Henley, Nonplanar ground states of frustrated antiferromagnets on an octahedral lattice, Phys. Rev. B 88, 024407 (2013b).
  * Miyatake  _et al._ [1986] Y. Miyatake, M. Yamamoto, J. J. Kim, M. Toyonaga, and O. Nagai, On the implementation of the 'heat bath' algorithms for Monte Carlo simulations of classical Heisenberg spin systems, Journal of Physics C: Solid State Physics 19, 2539 (1986).
  * Creutz [1987] M. Creutz, Overrelaxation and Monte Carlo simulation, Phys. Rev. D 36, 515 (1987).