Abundance of hard-hexagon crystals in the quantum pyrochlore antiferromagnet

We propose a simple family of valence-bond crystals as potential ground states of the $S=1/2$ and $S=1$ Heisenberg antiferromagnet on the pyrochlore lattice. Exponentially numerous in the linear size of the system, these can be visualized as hard-hexagon coverings, with each hexagon representing a resonating valence-bond ring. This ensemble spontaneously breaks rotation, inversion and translation symmetries. A simple, yet accurate, variational wavefunction allows a precise determination of the energy, confirmed by DMRG and numerical linked cluster expansion, and extended by an analysis of excited states. The identification of the origin of the stability indicates applicability to a broad class of frustrated lattices, which we demonstrate for the checkerboard and ruby lattices. Our work suggests a perspective on such quantum magnets, in which unfrustrated motifs are effectively uncoupled by the frustration of their interactions.

We propose a simple family of valence-bond crystals as potential ground states of the S = 1/2 and S = 1 Heisenberg antiferromagnet on the pyrochlore lattice. Exponentially numerous in the linear size of the system, these can be visualized as hard-hexagon coverings, with each hexagon representing a resonating valence-bond ring. This ensemble spontaneously breaks rotation, inversion and translation symmetries. A simple, yet accurate, variational wave function allows a precise determination of the energy, confirmed by density matrix renormalization group and numerical linked cluster expansion, and extended by an analysis of excited states. The identification of the origin of the stability indicates applicability to a broad class of frustrated lattices, which we demonstrate for the checkerboard and ruby lattices. Our work suggests a perspective on such quantum magnets, in which unfrustrated motifs are effectively uncoupled by the frustration of their interactions.
The propensity of frustrated quantum magnets to host exotic ground states makes them a rewarding target for study. For several of the most prominent models, however, it has proved difficult to come to a consensus about the exact nature of the ground state . This is particularly true for frustrated isotropic Heisenberg models with low spin S ≤ 1 in dimension d > 1, for which approximate analytical methods are generally uncontrolled and numerical studies are challenging. To improve the situation requires simple yet reliable heuristics, by which we can understand not just which state is the ground state but also whence it derives its stability.
In this Letter, we propose a family of valence-bond crystal states generated from hard (nonoverlapping) hexagon coverings of the pyrochlore lattice [43] shown in Fig. 1(a). A one-parameter variational wave function describing these states achieves an energy equal to the best-known proposals to date within numerical uncertainties. The hard-hexagon crystal states are exponentially numerous in the linear system size, demonstrating the abundance of competing low-energy states. The difficulty in arriving at a consensus over the ground state is likely in considerable part due to the presence of so many competing states with barely distinguishable energies.
The simplicity of the hard-hexagon states enables us to understand their low energy and stability, which arises from the following ingredients. The first is that the pyrochlore lattice can be decomposed into nonoverlapping hexagonal loops -the hard hexagons -which exhibit a robust finite-size gap of ∼ 0.69J. Second, the hexagons are connected by a quartet of bonds linking two pairs of antiferromagnetically correlated spins symmetrically, making the coupling doubly frustrated [ Fig. 1(d)] [44]. This latter point is perhaps conceptually the most interesting one, as it shifts the perspective from the frustration-induced degeneracy arising on a single tetrahedron to the isolation of unfrustrated (nondegenerate) geometrical motifs. Third, the kinetic energy of local defects in this background pattern is similarly suppressed. Fourth, there are no matrix elements between different hard-hexagon coverings to any finite order in perturbation theory.
From this point forward, we discuss in detail the FIG. 2. (a) Energy per site at finite temperature for the pyrochlore S = 1/2 and J = 1. The orange curves are showing a tetrahedron-based NLCE expansion up to eighth order [50], which defines an upper limit for the ground-state energy. The blue curves show the NLCE expansion based on the hexagons. Euler refers to a resummation algorithm extending the convergence down to lower temperatures [50]. (b) Different estimates for the ground-state energy per site: bare hard hexagons (ground state of H0), upper bound obtained by the converged tetrahedron expansion [50], Astrakhantsev et al. [38], Harris et al. [35], dressed hexagon state Eα (this work) [Eq. (2)], Hagymási et al. [37] and the NLCE hexagon expansion at second order for T = 0 (this work). Estimates of the ground-state energy higher than Eα the variational energy can be ruled out.
Having understood its essential ingredients, we are able to apply similar constructions to the S = 1 pyrochlore Heisenberg model as well as Heisenberg models on the ruby [ Fig. 1 (b)] and checkerboard [Fig. 1(c)] lattices. In these cases, we verify that the variational energy of these states is competitive with the ground-state energy obtained by density matrix renormalization group (DMRG) methods [45][46][47][48][49] [ Table I].
We start by covering the pyrochlore lattice with hexagons such that each site participates in exactly one hexagon. The number of ways to do this is exponentially large in the linear system size, as we discuss in more detail later. For a given covering, we then decompose the nearest-neighbor Heisenberg model into links within the hard hexagons, H 0 , and the connecting terms, V : Ground-state energy. Crucially, the decoupled hexagon Hamiltonian H 0 exhibits a large gap ∼ 0.69J, while the coupling V is effectively suppressed since it symmetrically couples pairs of neighboring-and hence antiferromagnetically correlated-spins in two hexagons via the four bonds [ Fig. 1 (5) TABLE I. Ground-state energies in units of J for different models and spin lengths. The ground-state energy is calculated using second-order NLCE, DMRG, and the variational wave function Eq. (2) with optimization of the parameter α.
The optimal values of α are given in Tab. S1 the Supplementary Material [58]. The DMRG energies for the pyrochlore lattice are from Hagymási et al. [37,51], while the DMRG results were obtained using ITensor [59]. * Note that the S = 1 pyrochlore case was obtained for 48 sites, and finite-size effects are likely to underestimate the ground-state energy.
is already not far from the pre-existing estimates of the ground-state energy of the full pyrochlore lattice. This motivates a variational approach by which we establish a strict upper bound on the ground-state energy of Eq. (1) in the thermodynamic limit which, within error bars, competes with previous extrapolations based on estimates for small clusters [35,37,38,51].
The trial wave function is constructed by dressing the ground state of H 0 , |Ψ 0 ⟩ -a simple product state. To introduce additional correlations between the hexagons and minimize the energy further, we perform imaginarytime evolution using the Hamiltonian of the tetrahedral links, V , connecting the hexagons: The variational energy per site E α , which we evaluate using an expansion in powers of α [52][53][54][55][56][57][58], exhibits a well-defined minimum at α = α 0 . The expansion is fully converged around the minimum and the truncation error of the variational energy is much smaller than the symbol size in Fig. 2 (b). Importantly, the resulting minimal energy E α0 (see Table I), due to its variational nature, rules out all larger estimates from previous studies [35,38], as indicated in Fig. 2 We further support the stability of the hard-hexagon state by a numerical linked cluster expansion (NLCE) [60,61]. This method has proven valuable in determining thermodynamic quantities at finite temperature in three-dimensional lattices, including the pyrochlore lattice [50,62,63]. The algorithm is in a spirit similar to a high-temperature expansion in the sense that it systematically includes larger clusters to obtain an estimate at a finite temperature. While most previous works on the pyrochlore exploit the tetrahedral structure (which is powerful at finite temperature), we generalized the expansion based on nonoverlapping hexagons [58]. Typically, the limit of convergence of a finite order expansion is clearly identified by the fact that successive orders rapidly diverge from each other to infinity below some temperature. This is the case for the tetrahedral expansion, as shown in Fig. 2 (orange lines). However, the situation is remarkably different for the hexagon-based expansion (blue lines in the same figure). It is clearly not converged for intermediate temperatures, 0.06J ≤ T ≤ 0.8J, but converges again for T → 0 yielding a realistic groundstate energy compatible with previous studies [35,37,38]. Strikingly, the low-temperature convergence is effectively obtained in the second order. This lends further credence to a particularly simple low-temperature state of weakly dressed hexagons. The success of the NLCE at T = 0 is reminiscent of a study of a distorted kagomé lattice [64], where a similar approach was used to support the conclusion of a dimerized ground state.
Hard-hexagon coverings and symmetry breaking. The proposed states are nonmagnetic valence-bond crystals, whose symmetry-breaking properties are determined by the underlying coverings of the lattice by hard hexagons. The pyrochlore lattice admits an exponentially large family of hard-hexagon states, all of which take the form of long-ranged ordered planes of hexagons stacked along one of the three equivalent ⟨001⟩ directions. Taking for example the covering shown in Fig. 1 (a) we can obtain a new valid covering by shifting the second plane from the bottom (in yellow-red) one unit to the right (i.e. along one of the ⟨110⟩ crystal directions). Assuming that all hard-hexagon coverings of the pyrochlore lattice can be constructed from such shifts, this yields a subextensive, yet exponentially large, number of coverings N cover = 3 × 2 4L/3 , where L is the linear system size given as the number of cubic, 16-site unit cells. We have verified numerically [58], for finite clusters up to L = 12, that these are the only coverings by using an unbiased numerical optimization algorithm.
The fact that there is an infinite family of states on the pyrochlore lattice, rather than a single hard-hexagon state, is interesting for two main reasons. First, for a single such state to be stable, it is important that a finite order of perturbation theory connect no two members of the family. This is clearly the case for the stackedlayer coverings discussed above since going from one to the other requires the translation of an entire plane of hexagons. Second, this has implications for the symmetry breaking of a state randomly selected from, or averaged over, this family. Since in a valence bond solid, two-point correlation functions decay rapidly, we instead compute a four-point correlation function which we call the dimer structure factor, or bond correlator, for an ensemble of hard-hexagon states with [001] stacking direction obtained by the numerical optimization algorithm mentioned before. The result is shown in Fig. 3. In contrast to the "usual" spin structure factor based on a two-point correlation function, the dimer structure factor is based on a four-point correlation function and carries signatures of the translational and rotational symmetry breaking [Bragg peaks in panels (a) and (b)] as well as the disorder in the stacking direction [broad features in panels (c) and (d)] and also the short-range correlations within the hexagons (low-intensity broad features rendered visible by the log scale). We note that, given the non-magnetic nature of the symmetry-breaking, probes such as Raman scattering or ultrasound measurements may be useful experimentally besides neutron scattering, which can best establish the absence of long-range magnetic order.
Robustness of the gap. Since the hard-hexagon states break only discrete symmetries, they are -if stablegenerically gapped. We analyze the excitations above our family of hexagon states to asses their (local) stability. As the candidate ground states are weakly dressed product states of single-hexagon ground states, the simplest ansatz for the excited states is given analogously by weakly dressed local triplet excitations on a single hexagon. The lowest-lying excited state on a single hexagon is given by a triplet, with a gap of ∼ 0.69J. A possibility to reduce and finally destabilize the gap would be for the localized excitations to gain kinetic energy. Numerically evaluating the excitation dispersion in three dimensions is difficult, and presently out of reach. Instead, we turn to an analytic approximation: namely the multiboson expansion, which is able to describe the hopping of local excitations above a given ground state, at first order in perturbation theory [58,65,66]. The results are shown in Fig. 4: remarkably, the lowest-lying triplet bands, |t − m ⟩ (m = −1, 0, 1) at E ∼ 0.69J, remain completely flat, while the higher energy triplets disperse (bands labeled by |k ± n ⟩ in Fig. 4) and are actually pushed below the excited singlet located at E ∼ 1.3J, which also remains flat. The flatness of the bands is due to a symmetry of the hopping matrix elements V , cf. see Eq. (1), connecting two hexagons A and B. This is illustrated in Fig. 4 (b): mirroring each hexagon individually across a plane running through the tetrahedron leaves the operator matrix element V (A, B) invariant. This implies a selection rule for the hopping matrix element ⟨s 0 , t| V |t, s 0 ⟩. Classifying all eigenstates of H 0 by the three mirror symmetries of a hexagon, it turns out that the ground state singlet |s 0 ⟩ as well as the lowestlying triplet |t − m ⟩ are simultaneous eigenstates of all three, but with eigenvalue +1 and −1 respectively. This then guarantees ⟨s 0 , t − m | V |t − m , s 0 ⟩ = 0 by symmetry, hence the flatness of the lowest band. By contrast, the two triplet pairs at E ∼ 1.5J and E ∼ 1.8J are not simultane-ous eigenstates of all three mirror operations and hence are not prevented from mixing into dispersive modes. In Fig. 4 (c), we show schematically a process with several high-energy virtual intermediate states, which would give the lowest-lying triplet kinetic energy. However, since this only appears at the fourth order in V , we do not expect this to lower the gap of the hard-hexagon state significantly. Finally, we note that the symmetry argument above does not depend on the choice of a particular tiling, but the exact shape of the dispersion of highenergy bands in Fig. 4 (a) will.
To summarize our results for the pyrochlore S = 1/2 model: we propose an exponentially large family of valence-bond crystals based on hard-hexagon coverings as candidate ground states. A variational calculation of their energy in the thermodynamic limit is within the error bars of the best numerical estimates of the groundstate energy. The stability of these states is further supported by NLCE and by a multiboson expansion which shows that the kinetic energy of excitations is suppressed. The hard-hexagon coverings break rotation, translation, and inversion symmetries of the lattice. It is important to note that the recent numerical studies that found lattice symmetry breaking [37,38] used clusters that are incompatible with the hard-hexagon states, rendering them energetically disfavorable there.
Note that the same ingredients rendering the stability of the valance bond state are found in other Heisenberg models, such as the two-dimensional ruby and checkerboard lattice for S = 1/2 and S = 1, cf. Fig. 1. We obtain analogous results, such as the flatness of the triplet band and a finite size gap, in these cases. A short discussion of spin S = 1 [67] and broken symmetries [68][69][70][71][72][73] are given in the supplement.
Discussion. Having identified a new family of energetically competitive states based on effectively decoupled close-packed motifs, we still cannot settle the question of what is the actual ground state in the pyrochlore case. On a technical level, our study underlines the need to consider possibly very large unit cells in finite-size studies. We note that such a change in perspective was already proposed in the context of the S = 3/2 spinel compound ZnCr 2 O 4 [43], where hexagonal motifs had been identified in a study of its magnetoelastic properties [74]. In this case, theoretical modeling ultimately suggested that the system is not described by independent hexagonal clusters but rather that weak further neighbor exchange can account for the salient observations [75,76]. In this work, we have demonstrated that the formation of such clusters can indeed occur, but via a different, quantum, mechanism.
The existence of potential material realizations of the pyrochlore Heisenberg antiferromagnet [77,78], and the growing capabilities of cold atom emulations [79][80][81], gives hope that some of these questions may eventually be settled by experiment. The presence of a significant gap to spin excitations above the hard hexagon states would be an important signature for both thermodynamic and spectroscopy measurements if indeed the system can equilibrate into such a state at low temperatures. The presence of sharp, gapped, S = 1 excitations, along with the presence of lattice symmetry breaking would serve to distinguish the hard hexagon states from competing quantum spin liquids.
However, due to the presence of such a large family of low energy states, as well as potential additional states beyond the ansatz considered here, it may be that the physics at the temperatures reachable in the experiment is controlled not by a single ground state but by many competing ones. Further, the eventual ground-state selection in an actual material in the presence of any near degeneracies will take place via any residual deviations from an ideal Heisenberg Hamiltonian, and may require exquisitely low temperatures. Finally, hexagonal motifs have appeared in various pyrochlore settings [74][75][76], and the mechanisms leading to their stabilization may reinforce each other in a given material.
We are very grateful to Imre Hagymási, Christopher Laumann, David J. Luitz, Frank Pollmann, Götz S. Uhrig, and Alexander Wietek for many helpful discussions on this topic. DMRG and TDVP calculations were performed using the ITensor [59] package with the global subspace expansion [82]. This work was in part sup-  The family of valence-bond crystals that we propose in the main text is applicable to Hamiltonians which can be decomposed into plaquette terms and edge terms connecting two plaquettes: In the case of the pyrochlore lattice, the "plaquettes" correspond to the hexagons and V e refer to the quartet of bonds connecting two hexagons. Another example -the checkerboard lattice -is shown on the right of Fig. S1. An obvious starting point for building a trial wave function for our family of states is the product state of singlehexagon ground states |s 0 ⟩ ) arXiv:2210.07235v3 [cond-mat.str-el] 1 Sep 2023 FIG. S1: Reduced lattice (left) for the checkerboard lattice (right). We identify the plaquettes with the nodes and the links of the tetrahedron connecting two plaquettes with the links of the reduced lattice. The same decomposition was exploited for the NLCE algorithm.
Which is an eigenstate of H 0 and the zero expectation value vanishes for the perturbation V in our case. Hence, it has an energy expectation value 47J in the pyrochlore for example). To make the energy of |Ψ 0 ⟩ competitive with values in the literature S37,S38,S70,S71 , we have to weakly dress it. This could in principle be done perturbatively. However, because of the simplicity of our proposed family of states, we are able to construct a dressed variational wave function, which has the great advantage of yielding a strict upper bound on the ground-state energy: where v is some set of variational parameters and S(v) some local operator. In this paper, we consider a very simple form S = αV with only a single variational parameter α. It can be interpreted as an imaginary-time evolution While α = 0 corresponds to just the single-hexagon product state |Ψ 0 ⟩, increasing α constitutes a trade-off between paying energy on the hexagons to gain some energy on the connecting bonds. The variational energy per site is then given by S5) where N p is the number of plaquettes and L p is the length of the single plaquette.

A. Linked Cluster Theorem
The exponential form of the variational state Eq. (S4) makes it possible to develop a systematic linked cluster expansion S53,S54,S57 to evaluate the trial energy E α in Eq. (S5) in the thermodynamic limit. In the remainder of this section, we derive this expansion, explain how it is implemented and show that it reveals a pronounced minimum in the variational energy. We emphasize that the linked cluster expansion used here to calculate the variational energy is different to the numerical linked cluster expansion described in the main text and in Sec. S II, which is based on exact diagonalization and does not require a variational ansatz.
To start, note that |Ψ α ⟩ is a spin-singlet (as it should be for a SU (2) symmetric Hamiltonian), and all plaquettes, as well as the edge terms connecting two plaquettes, are equivalent, respectively. Hence, to compute the variational energy E α , it suffices to compute the energy of a single plaquette term and a single connecting term: S6) where ⟨O⟩ α = ⟨Ψ α |O|Ψ α ⟩ for some operator O and c is the number of edge terms V e attached to each plaquette. In the cases considered throughout this letter, it is given by the length of the plaquette. Following the linked cluster theorem, the expectation value of any local operator can be written in terms of connected correlation functions, also known as Ursell functions S52,S55,S56 . To obtain the expectation value of some local operator O loc we expand the exponentials in e −αV O loc e −αV and only consider terms that are connected to the local operator. Each connected term with n operators is weighted according to the nth Ursell function; the first two Ursell functions are given by The linked cluster theorem S53,S54,S57 states that non-connected terms that factorize do not contribute to the local observable. In the technical sense, connected terms are those with connected support given by the operators stemming from the exponentials. This yields In order to identify a set of edges {e 1 , . . . , e n } "connected" to O loc , it is useful to introduce the notion of the reduced lattice. The reduced lattice for the partition of plaquettes in the checkerboard lattice is shown in Fig. S1. Each vertex v in the reduced lattice corresponds to a plaquette, and two vertices are connected by an edge e = (v, v ′ ) if the two corresponding plaquettes, v and v ′ , are connected. A set of edges {e i } in this language is then "connected" to O loc if the subgraph induced by the {e i }, joint with the support of O loc , is connected. We illustrate Eq. (S9) by evaluating ⟨O loc ⟩ α to first order in α. The nominator and denominator in Eq. (S9) have the following forms: Here, N (O loc ) contains "edges" e such that V e and O loc have overlapping support. For example, a plaquette term H p in the checkerboard lattice exhibits four neighboring edges V e each containing a quartet of doubly frustrated bonds.
We have used that all operators are Hermitian and the expectation value of two operators with no common support factorizes Evaluating Eq. (S9) according to the linked cluster theorem is restricted to neighboring terms of O loc weighed with  the Ursell functions which yields Collecting the terms in from Eq. (S10), Eq. (S11) and Eq. (S13) confirms the linked clusters theorem to first order in α as stated in Eq. (S9): Analogously, it can be shown that this procedure holds for all powers of α S54 . However, the possibility to carry out the expansion straightforwardly is quickly exhausted as the complexity scales super exponentially. Therefore, in the next section we present a systematic way to determine the variational energy. Bethe (n = 6) pyrochlore ruby Bethe (n = 4) checkerboard 0. 43(3) 0.3961 (4) 0.3932 (2) 0.393 (1) 0.30529 (4) 0.30529 (4) S = 1 Bethe (n = 8) Bethe (n = 6) pyrochlore ruby Bethe (n = 4) checkerboard 0.44 (9) 0.24 (7) 0.24 (7) 0.24 (7) 0.194 (4) 0.197 (2) TABLE S1: Optimal imaginary time step α 0 to minimize the energy of the variational wavefunction |Ψ α ⟩ for different lattices and spin lengths. The table demonstrates that the underlying lattice topology has negligible influence, as the optimal parameter is mainly determined by the loop length and spin. The error is defined by half the difference between the largest and second largest order. Note that the error is larger than the optimal energy itself (cf. Tab. 1 in the main text) because the minimum is flat and derivations from the optimal α 0 have only a marginal influence.
(a) For k = 1, . . . , n we generate the sets C k containing k different edges that are, together with O loc , connected.
The expansions for H p and V e are different.
(b) Given a specific set with k edges {q 1 , . . . , q k } ∈ C k , we have to generate all possible (ordered) sequences of n edges (e 1 , . . . , e n ) such that each q i occurs at least once. The restriction that each q i has to occur at least once ensures the connectivity together with the local operator. We denote the set of all generated edge sequences by S.
(c) Each edge sequence (e 1 , . . . , e n ) ∈ S generates an operator string composed of n + 1 components where the local operator is inserted at position a + 1. This corresponds to an individual term in Eq. (S9) which is given by the (n + 1)th Ursell function and evaluates to u Pushing the expansion to high order is challenging as each individual step described above scales exponentially or even factorially with n. Luckily, in our case, the above procedure simplifies significantly since the expectation value of any single V e with respect to a product state on the two plaquettes vanishes if one of the two states is the ground state singlet This implies that the expectation value of a product of V e with respect to |Ψ 0 ⟩ vanishes unless each plaquette is touched by at least two V e . Referring to the recipe above, this means that in step (b), we only need to generate sequences such that each vertex is either part of at least two edges or lies in support of the local operator O loc .

C. Results
The final results for the variational energy E α [Eq. (S5)] on the pyrochlore and checkerboard lattice are shown in Fig. S2. We compute the variational energy using the linked cluster expansion derived above up to order six (seven) for the pyrochlore (checkerboard) lattice. The order refers to the number of edge terms V e included in the expansion. As a function of α, it shows a well-pronounced local minimum at α 0 . The optimal imaginary time steps for the different models considered here are listed in Table S1. Note that in the relevant range of α, orders larger than four (five) are barely distinguishable by the eye in Fig. S2. At the minimum, the error in the pyrochlore lattice estimated by the difference between order five and six is of order 10 −4 J reflecting the convergence of the algorithm. We define the error of the expansion by the difference between the largest and second-largest order at the optimal value α 0 . This is shown in Tab. 1 in the main text. For the checkerboard lattice, we also perform the time evolution on a large finite cluster of size N = 100 with periodic boundary conditions using the time-dependent variational principle (TDVP) to check consistency with our linked cluster expansion, finding good agreement (see also Sec. S III). Smaller clusters additionally allow the application of an exact finite-temperature Lanczos calculation S42 starting from an initial plaquette product state where 2α serves as the inverse temperature and V as the Hamiltonian. Implementing this on a finite cluster of the checkerboard lattice with N = 36, we find the same qualitative dependence of E α on α. However, the result does not agree quantitatively with the other two methods, which we attribute to the presence of strong finite-size effects.
The coefficients of the expansion of the variational energy in powers of α are listed in Table S2. We list not only the pyrochlore and checkerboard but also results for the ruby lattice as well for different Bethe lattices. By "Bethe (n = 4, 6, 8)" we here denote that the reduced lattice is a Bethe lattice. The number n given in parenthesis is the size of the unfrustrated loops constituting H 0 . These loops are then connected via a quartet of bonds similar to the pyrochlore, checkerboard, and ruby case, but in such a way that the resulting reduced lattice is a regular tree of degree n.
We note that for some of the models (marked by an asterix in table Table S2), the coefficients of ⟨V e ⟩ α at the largest order in the expansion are computed only approximately. Each order in α contains contributions of operator strings with support on clusters of at least one (two) and at most n + 1 (n + 2) plaquettes for O loc = H p (O loc = V e ). Due to a large number of arising terms we omitted, for some coefficients of O loc = V e , the contributions of the largest clusters. We expect this to not make any difference in practice for the following reasons. (i) Most importantly, the expansion is already converged at much lower orders, and (ii) also the contribution ⟨V e ⟩ α is small compared to that of ⟨H p ⟩ α . (iii) The omitted terms are only a small fraction of those contributing at the given order, and they also typically exhibit small weight. (iv) The omitted cluster occurs at high order, which is suppressed for α 0 ≈ 0.4. All this is also evidenced by the excellent quantitative agreement between the numerical imaginary-time evolution and the linked cluster expansion.
Remarkably, the lattice geometry and dimensionality have little impact on the numerical values of the coefficients in the linked cluster expansion as long as the plaquette length stays fixed. This is true even at orders where the lattices are topologically distinct (e.g. orders larger than three for the ruby and pyrochlore lattices). This again validates the (local) stability of the hard-hexagon crystals, as it implies that even closed loops of plaquettes, contributing at higher orders, are negligible.
For some of the models, the largest order of the expansion of E weak in α is incomplete and we mark the corresponding coefficient by an asterix * . However, these coefficients only miss a small number of contributions from some of the largest clusters, which we expect to be negligible.

S II. Numerical Linked Cluster Expansion
The numerical linked cluster expansion (NLCE) S60,S61 was originally developed to compute extensive thermodynamic observables at finite temperature in the thermodynamic limit. Starting from a pre-defined unit, e.g. a hexagon or tetrahedron, the lattice is built up systematically by expanding the clusters unit by unit. Each generated cluster c has a weighted contribution W [⟨O c ⟩ β ] to the observable ⟨O⟩ β at some inverse temperature β. The weight is obtained by subtracting the weights of smaller connected clusters c ′ that are contained in the cluster c ′ ⊂ c. Each cluster has L c topologically equivalent instances arising during the expansion.
The numerical determination of the observable for a single cluster does not allow for approximate methods such as quantum typicality S42 or finite-temperature DMRG S50 as it is extremely sensitive to errors. Hence, as far as finitetemperature properties are concerned, observables have to be evaluated using full ED. In this letter, however, we also use the NLCE to compute zero temperature properties. Since at zero temperature, only the ground state contributes, we can use the Lanczos algorithm (which is numerically exact for T = 0) to compute the expansion to much larger orders (six instead of three) than would be possible for an expansion at finite temperature.
The NLCE algorithm is in a spirit similar to a high-temperature expansion, where larger clusters contribute at higher orders in β. Physically, this is motivated by a growing correlation length at lower temperatures. Once the correlation length is larger than the largest cluster included in the expansion, consecutive orders will usually diverge away from each other immediately, signaling the failure of convergence. This is seen in Fig. 2 (a) of the main text for the expansion based on tetrahedra (blue), which has proven extremely powerful as it converges down to temperatures T ≈ 0.2J, inducing a relevant upper bound on the ground-state energy. The behavior is quite different if the expansion is instead based on hexagons, shown as orange curves in the same figure. While convergence is lost around T ≈ 1J, the individual curves do not diverge to infinity and the algorithm seems to converge again for T → 0. Remarkably, the convergence at T = 0 is reached at the second order in the NLCE expansion. Hence, the weight of larger clusters and their contribution is very small. This can be observed in Fig. S3 (a), which shows the absolute value of the weight of the energy at T = 0, |W [⟨E c ⟩ ∞ ] |, for each cluster from order two to five. The fact that the individual orders do not diverge to infinity at intermediate temperature suggests that the low energy spectrum for the clusters is similar. Other literature results S35,S37,S38 and DMRG calculations, cf. Sec. S III, for the checkerboard lattice as well as the introduced wavefunction propose a ground-state energy that is slightly higher than the result from the NLCE calculation.
As stated before, the loss of convergence is associated with the presence of correlations beyond the largest clusters considered. In light of this, the fact that the series seems to converge at second order in the limit T → 0 suggests a simple structure of the ground states of considered clusters. In particular, in the second order, one considers exactly two coupled hexagons. This, when put on the pyrochlore lattice given a hard-hexagon covering, can capture only a uniform energy density on all hexagons as well as a (different) uniform energy density on the bonds connecting any two hexagons. This is fully consistent with the weakly dressed valence-bond crystals proposed in this letter. Note that even though the convergence is basically achieved at order two, we do compute the expansion up to order six. Higher orders include clusters which in principle could break the symmetries of the hard-hexagon states, but computing their ground states shows that they do not. Instead, the ground states of each individual cluster realize a hard-hexagon crystal, described by the variational wave function introduced in Sec. S I to remarkable accuracy, as shown in Fig. S3 for the clusters up to order five. For each cluster, we optimize the hard-hexagon wavefunction and find the optimal value of the parameter α using imaginary-time evolution. In each case, the resulting energy agrees closely with the exact ground-state energy of the cluster. The optimal values of α hardly vary from cluster to cluster and agree closely with the value of α 0 obtained in the thermodynamic limit. We do not include the 283 clusters appearing at order six in the figure for clarity. The final NLCE results for the energy at zero temperature at orders two and six are NLCE 2 = −0.4917J and NLCE 6 = −0.4919J indicating the little influence of larger clusters. The convergence in the second order is further reflected by the small error in Tab. I of the main text for the different models considered in this work. It is obtained by the difference between the second and the third order.

S III. Matrix-Product State Approaches
Approaches based on matrix-product states (MPS) S45-S48 are the method of choice for strongly correlated onedimensional systems. Being a one-dimensional technique, generalizations to higher dimensions are obtained by linearizing the system using a one-dimensional path -a "snake" -that traverses the system. The class of models proposed here is well suited for the latter flavor by choosing the snake path to traverse the plaquettes covering the lattice. We implement such a "snake" MPS using the ITensor S59 package. First, we use it to perform imaginary-time evolution on large finite clusters of the checkerboard and ruby lattices and use this to corroborate the correctness The difference between the exact ground-state energy of each cluster and the variational energy obtained by optimizing |Ψ α0 ⟩ = e −α0V |Ψ 0 ⟩ on each cluster. The two agree very closely, suggesting that the dressed hexagon state is stable on every individual cluster contributing to the NLCE expansion. (c) The optimal variational parameter, α 0 for each cluster, which varies very little between clusters, always agrees closely with the value we obtain in the thermodynamic limit [Sec. S I]. The x-axis labels different clusters organized by (order,label). of our determination of the variational energy of the dressed hard-hexagon state (see Sec. S I). Second, we also use DMRG on finite clusters of the checkerboard and ruby lattices to study the structure of their local excitations. The result for the ruby lattice, which is equivalent to the pyrochlore up to the second order in perturbation theory and at larger orders has more closed loops, also yields a heuristic lower bound on the gap of the pyrochlore lattice.

A. Imaginary-Time Evolution
We use the global subspace expansion for the time-dependent variational principle (TDVP) S82 to perform imaginarytime evolution starting from the single-plaquette ground state |Ψ 0 ⟩. The main result is shown in Fig. S2b, which compares the variational energy E α , obtained from imaginary-time evolution on the checkerboard lattice using (i) exact diagonalization on a cluster of N = 36 sites with periodic boundary conditions, (ii) the linked cluster theorem expansion described in Sec. S I, and (iii) TDVP on a finite cluster of N = 100 sites with periodic boundary conditions. For the ED result, the dependence of the variational energy on α agrees qualitatively but not quantitatively with the other two results. We attribute this to pronounced finite-size effects. In contrast, the result from linked cluster expansion and from TDVP are in excellent quantitative agreement. This confirms the validity of the expansion, which we also use for the pyrochlore lattice where full TDVP on large clusters in unattainable.  TABLE S3: DMRG data for checkerboard and ruby lattice (J = 1). The number of plaquettes is L 2 such that the clusters contain N = 4L 2 and N = 6L 2 spin-1/2 or spin-1 sites with periodic boundary conditions, respectively. The one-dimensional "snake" path is placed within the plaquettes and breaks the inversion symmetry by choosing one over the other possible tiling in the checkerboard lattice. All results are obtained for a bond dimension of χ = 4000, and the maximal truncated weight is of the order 10 −5 . ∆ mz refers to the gap obtained in the magnetization sectors m z = 0, ±1. The derived variational energy obtained from the cluster expansion is E α0 = −0.51344 for the checkerboard and E α0 = −0.48947 for the ruby lattice. Note that the variational energy for the ruby lattice is slightly lower than the DMRG calculation.

B. DMRG results
In addition to the imaginary-time evolution, we also utilize DMRG. We compute estimates of the energy and the real-space correlations of the ground state and the first excited state in the m z = −1, 0, +1 magnetization sectors. We consider finite clusters of both the checkerboard and the ruby lattice, with the linear system size L defined such that the total number of sites in a cluster is 4L 2 and 6L 2 , respectively. In all cases, we construct the snake path such that it passes through the plaquettes and adjacent plaquettes in real space is traversed by it with a maximal plaquette distance of 2L in the one-dimensional topology. A naive linearization of a cluster with periodic boundaries would yield a maximal distance of L 2 .
Results are shown in Table S3. The ground-state energy per site is denoted by E 0 /N . We carefully examine the ground-state energy following a similar procedure as carried out in Ref. S37 , which serves as an estimate for the pyrochlore S = 1/2 case. We first extrapolate to infinite bond dimension and then to the thermodynamic limit as illustrated in Fig. S4. To obtain an energy estimate for a cluster of finite size N , we extrapolate the energy in 1 χ where is χ is the bond dimension. The error is defined by the difference between the last data point (largest bond dimension) and the extrapolated limit E N 0 . The energy estimates obtained for finite-size clusters are shown in the insets of Fig. S4 on a 1 N scale. As done in Ref. S37 , we use a quadratic fit in 1 N 2 to determine the energy in the thermodynamic limit that takes the uncertainty of the finite clusters into account. The confidence interval is obtained by the difference between the largest system size and the extrapolated result E 0 . We find the ground-state estimate is in excellent agreement for all lattices and spin lengths with the variational estimates E α (cf. in Tab. 1 of the main text). This applies not only for the ground state but also for the reals-space energy distribution on the lattice, that is, the energy densities on strong bonds (E strong = ⟨H p ⟩ α0 /L p ) and weak bonds (E weak = ⟨V e ⟩ α0 /4). The individual contribution of weak and strong bonds to the energy are summarized in Table S4.
Motivated by these excellent results for the ground state, we now turn to a careful study of the lowest-lying excitations. We initially prepare the lowest lying excitation on a single plaquette which is given by the triplet t − mz for the respective magnetization sectors m z = 0, ±1 on top of the of a product state of single plaquette ground states |Ψ 0 ⟩. The resulting triplet gap is 1J for an isolated square and ∼ 0.69J for an isolated hexagon. DMRG preserves the total magnetization such that the lowest lying excitation in the m z = 1 sector is obtained by minimizing the energy without constraints. In contrast, the m z = 0 excitation has to be orthogonalized to the plaquette ground state. For both models, we found a lower gap in the m z = 0 sector compared to m z = ±1. This suggests that the lowest lying excitation in both cases is a singlet. This is somewhat surprising, given that the lowest lying excitation on both, a FIG . S4: Extrapolation of the ground-state energy for the checkerboard and ruby lattice for S = 1/2 and S = 1 using DMRG data at finite bond dimension. We first extrapolate to infinite bond dimension 1/χ → 0 and then to the thermodynamic limit 1/N → 0. We use a linear fit in 1 χ to extrapolate to infinite the bond dimension and a quadratic fit for the thermodynamic limit 1 N 2 to be consistent with the pyrochlore results from Ref. S37 .
single square as well as a single hexagon, is a triplet. Possible scenarios include breaking plaquettes into different (larger) unfrustrated motifs or strongly dressed localized singlets excitation. Fig. S5 shows the real space correlation of the individual states for the checkerboard lattice with L = 5 and the ruby lattice with L = 4. The left columns show the respective ground states followed by the lowest-lying excitations in the m z = 0 and m z = ±1 sectors. In the case of the checkerboard lattice, the lowest lying excitation in the m z = 0 sector indeed forms a larger motif [bottom center of Fig. S5a]. In contrast, the m z = 1 excitation appears localized on a single plaquette for both considered models, consistent with a weakly dressed, localized triplet excitation. The initial triplet gaps (1J and ∼ 0.69J) are reduced significantly by an increasing contribution of the neighboring weak bonds.
Our results on the ruby lattice suggest a triplet gap of size ∆ mz=±1 ∼ 0.36J. We note that this is similar in value to estimates on the pyrochlore lattice using DMRG calculations on a N = 64 cluster, yielding ∼ 0.42J S41 and variational Monte Carlo yielding ∼ 0.40J S38 . and t − ±1 , in the m z = 0 and m z = ±1 sector initially. The lower rows show the ground state and the respective excitations in the different magnetization sectors of the full system. The data were obtained for a bond dimension χ = 4000 using U (1) DMRG, and the truncation error is of order 10 −5 . As DMRG does not conserve the total spin and the m z = 0 gap differs from the m z = ±1 gap, the observed m z = 0 excitation is a singlet state.

Pyrochlore
Ruby Checkerboard  S4: Contribution of strong and weak bonds to the energy (S = 1/2, J = 1) of the hard-hexagon state, obtained by different methods. The first row refers to the imaginary-time evolution described in Sec. S I and the second row presents DMRG results which are obtained within this work, cf. Sec. S III, for the ruby and checkerboard lattice (bond dimension χ = 4000 and truncated error is of order 10 −6 ) and by Ref. S37 for the pyrochlore lattice. Row three and four refer to the first and second-order NLCE expansion. The last row presents ED results for a finite system with open boundary conditions (N = 42 and N = 20) where a center plaquette of length n is coupled to n other plaquettes (N = (n + 1)n) via the doubly frustrated quartets. The outer plaquettes exhibit open boundaries and the strong and weak bonds are extracted from the center plaquette. This reflects the robustness of the plaquette crystal.

A. Finding Hexagon Coverings by Simulated Annealing
A hard-hexagon covering of the pyrochlore lattice is a choice of non-overlapping hexagons such that each vertex is part of exactly one hexagon. In all such coverings, one example is shown in Fig. 1 (a) (main text), the hexagons are arranged within one of the {001} planes which are then stacked along the third direction. Each plane can be translated independently by either one of the ⟨110⟩ or ⟨110⟩ directions while preserving the hard-covering constraint. This yields an exponential, but subextensive number of coverings N cover = 3 × 2 4L/3 where L is the linear system size in units of cubic, 16-site, unit cells. As the main text explains, these coverings break rotation, translation, and inversion symmetry. The translation symmetry is particularly striking as its unusual periodicity (three) means that most finite clusters used in numerical studies of the pyrochlore Heisenberg model to date are incommensurate with any hard-hexagon covering.
We have verified that the planar coverings as discussed above and in the main text are the only possible coverings by classical Monte-Carlo simulation of a hard-hexagon model where the sum is over sites j and n hex (j) is the number of occupied hexagons that the site is part of. The zero-energy states of this model are exactly the hard-hexagon coverings of the pyrochlore lattice. We performed 10 4 simulated annealing runs on a system with 6×6×6 unit cells (that is N v = 3456 sites) and found all 68 non-symmetry-equivalent hexagon coverings of the expected form and nothing else (note that 10 3 runs already suffice to find all non-symmetry equivalent coverings). We take this as strong evidence that indeed all possible coverings of the pyrochlore lattice with hard hexagons take the form shown in Fig. 1 (a). For L = 12 (that is N v = 27648 sites), we are not able to find all possible symmetry-equivalent planar states due to their exponentially large number. However, we did perform 10 4 simulated annealing runs and verified that all coverings found are of the expected, planar form.

B. Dimer Structure Factor
One way do diagnose such a hard-hexagon valence-bond-crystal state would be the unusual translational symmetry breaking with period three. Since two-spin correlators vanish beyond single hexagons, one needs higher-weight correlation functions. A simple, transitional invariant quantity to diagnose such a symmetry breaking would be the dimer-dimer structure factor where ⟨ij⟩ and ⟨kl⟩ are nearest-neighbor pairs on the lattice. We show this quantity in Fig. 3 of the main text. To compute it in practice, we approximate the four-spin correlation function in Eq. (S18) as That is we include the dominant short-range correlations exactly while approximating the long-range correlations as the Fourier transform of the energy density. This is readily computed for a specific hexagon covering and assuming self-averaging, we then obtain the plots in Fig. 3 of the main text by computing Eq. (S19) for all coverings of a L = 6 (N v = 3456) pyrochlore cluster with stacking direction [001] and taking the average.
To include longer-ranged correlations is possible, but including for example also neighboring hexagons and bonds connecting them exactly only change S dimer (q) by about 5%. In fact, the Bragg peaks and even the high-intensity broad features in Fig. 3 of the main text are already resolved without even treating correlations within one hexagon exactly. This is needed only to resolve the small-intensity broad features between resonances. . The difference to Ref. S43 comes from the fact that we consider S = 1/2 rather than classical vector spins and also we also do not multiply by any magnetic form factor.

C. Structure factor
For completeness, we show in Fig. S6 the spin-spin structure factor of the undressed hard-hexagon state |Ψ 0 ⟩. In the undressed case, all inter-hexagon correlation functions vanish and we only have to average the structure factor of a single hexagon over the four possible orientations that it can have in the pyrochlore lattice.

S V. Kinetic Energy from Multiboson Theory
In the following, we describe the multiboson theory approach S65,S66 used to compute the kinetic energy of triplet excitations up to first order in V , as shown in Fig. 4 (a) of the main text.
We start from a general Hamiltonian as decomposed into a term H 0 constituting the individual hexagons and a term V , which couples these hexagons via tetrahedra. To fix notation, we write explicitly where we have chosen to split up the site index r into the index of the hexagon h the site belongs to and a site index within that hexagon j ∈ {1, . . . 6}, that is r := (h, j). The notation t = ⟨h, h ′ ⟩ implies that the tetrahedron t couples hexagons h and h ′ .
We now proceed similarly to the bond-operator formalism of Ref. S65 to rewrite the spin operator in the basis of single-hexagon eigenstates. We denote the state n ∈ {0, . . . , 63} on hexagon h as |h, n⟩ and define a set of bosons a h,n such that |h, n⟩ = a † h,n |vacuum⟩ , S α h,j = S α h,j nm |h, n⟩⟨h, m| , = S α j nm a † h,n a h,m , where S α j nm = ⟨n| S α j |m⟩ = ⟨h, n| S α h,j |h, m⟩ ∀h (S25) and α = x, y, z. We have used that this matrix element does not depend on the hexagon h but only on the site index i, j. Then, while H 0 is quadratic in the a n , the coupling V consists of four-boson terms. The bosons are constrained such that there must be exactly one boson on each hexagon mn h,m = 1 ∀ h The idea of the linear multiboson theory is that, similar in spirit to linear spin wave theory, one can derive a quadratic approximation to H, assuming the excitation density is small. That is to say, we assume a † h,m a h,m =n h,m ≪ 1 for m > 0.
Eq. (S26) then implies thatn h,0 ≈ 1. The desired quadratic approximation is then obtained by neglecting all terms that are not at least quadratic in the ground state Boson operators and setting all quadratic terms to be 1. Formally, this can be done by substituting the ground state Boson operators with a number V = t=⟨h,h ′ ⟩ n,m ⟨ij⟩∈t α (S α i ) 0n S α j 0m a h,n a h ′ ,m + (S α i ) n0 S α j 0m a † h,n a h ′ ,m + (S α i ) 0n S α j m0 a h,n a † h ′ ,m + (S α i ) n0 S α j 0m a † h,n a † h ′ ,m where the ω n are just the eigenenergies of a single hexagon. Hence, we have obtained from H a standard (anomalous) bosonic hopping problem that is readily solved by Fourier transform.

S VI. Application to other models
We briefly discuss the application of the variational wavefunction to other models experiencing the same heuristic of unfrustrated loops with highly frustrated couplings. In particular, the same physics is observed in the ruby and checkerboard lattice [ Fig. 1 in the main text], which have been discussed throughout the supplementary material. The ground state energies for each case obtained by our variational construction are compared with DMRG and second-order NLCE in Tab. I of the main text.
To begin with, we consider the S = 1 version of the pyrochlore Heisenberg model. In this case, the entire argument carries over from the S = 1/2 case, with the additional twist that for integer spins the gap on an even-length loop is not an effect of its finite length S67 . We obtain a variational energy slightly above a recent estimate from a DMRG study on a cluster of 48 sites S51 , but given the likely considerable finite-size effects on the DMRG result, this constitutes a competitive estimate.
For the ruby lattice, there is a unique hard-hexagon covering, preserving all lattice symmetries. The hard-hexagon state on this model is thus a featureless quantum paramagnet, comparable in character to the ground state of the Shastry-Sutherland model S68 .
Finally, for the checkerboard lattice, the shortest unfrustrated loops that tile the lattice are squares. There are two possible hard-square coverings, and thus the hard-square state is a plaquette crystal breaking a Z 2 symmetry. This is consistent with previous results for the ground state order on the checkerboard lattice S69-S73 and so, in this case, our approach confirms the established consensus.
In all cases, the multiboson theory predicts gapped flat bands for the low-lying triplet excitations. In the case of the checkerboard and ruby lattices, we can use DMRG S45-S48 , cf. Sec. S III, to conduct an additional test of the robustness of the gap S49 . The knowledge of the (possible) ground state, together with the irrelevance of kinetic energy, enables us to use DMRG to reliably estimate the robustness of the gap to the local dressing of the lowest-lying excitations. We find that the gaps remain robust, with ∆ ≈ 0.47J for the checkerboard lattice and ∆ ≈ 0.27J for the ruby lattice. The latter is particularly interesting since it is up to second order in perturbation theory equivalent to the pyrochlore lattice. Even at higher order, the ruby lattice has more closed loops than the pyrochlore, and hence its gap serves as a heuristic lower bound for that of the pyrochlore.
Note that this construction is not straightforwardly generalized to arbitrary frustrated models. For example, while the kagomé lattice famously exhibits the hexagonal loop motif, it neither allows a hard-hexagon covering nor is the coupling between the hexagons doubly frustrated.