The pyrochlore S=1/2 Heisenberg antiferromagnet at finite temperature

Frustrated three dimensional quantum magnets are notoriously impervious to theoretical analysis. Here we use a combination of three computational methods to investigate the three dimensional pyrochlore $S=1/2$ quantum antiferromagnet, an archetypical frustrated magnet, at finite temperature, $T$: canonical typicality for a finite cluster of $2\times 2 \times 2$ unit cells (i.e. $32$ sites), a finite-$T$ matrix product state method on a larger cluster with $48$ sites, and the numerical linked cluster expansion (NLCE) using clusters up to $25$ lattice sites, which include non-trivial hexagonal and octagonal loops. We focus on thermodynamic properties (energy, specific heat capacity, entropy, susceptibility, magnetisation) next to the static structure factor. We find a pronounced maximum in the specific heat at $T = 0.57 J$, which is stable across finite size clusters and converged in the series expansion. This is well-separated from a residual amount of spectral weight of $0.47 k_B \ln 2$ per spin which has not been released even at $T\approx0.25 J$, the limit of convergence of our results. This is a large value compared to a number of highly frustrated models and materials, such as spin ice or the kagome $S=1/2$ Heisenberg antiferromagnet. We also find a non-monotonic dependence on $T$ of the magnetisation at low magnetic fields, reflecting the dominantly non-magnetic character of the low-energy spectral weight. A detailed comparison of our results to measurements for the $S=1$ material NaCaNi$_2$F$_7$ yields rough agreement of the functional form of the specific heat maximum, which in turn differs from the sharper maximum of the heat capacity of the spin ice material Dy$_2$Ti$_2$O$_7$, all of which are yet qualitatively distinct from conventional, unfrustrated magnets.


I. INTRODUCTION
The pyrochlore lattice, composed of corner-sharing tetrahedra, is a common motif in materials chemistry; in the context of magnetic materials, it has been prominent in a range of rare-earth 1 and spinel compounds 2,3 . Pyrochlore magnets and models have played a tremendously important role in the history of frustrated magnetism and topological condensed matter physics. One of the foundational publications, in 1956, was Anderson's identification of the classical pyrochlore Ising magnet 4 as an interesting model system. Now called spin ice 5 , this is a topological magnet exhibiting an emergent gauge field and fractionalised excitations 6 . The classical Heisenberg model, following a pioneering study by Villain 9 , turned out to be the first classical Heisenberg spin liquid 10 . This undergoes a very delicate order-by-disorder transition for large spins, as the zero-point energy induced by quantum fluctuations favours a subset of collinear states 11,12 . Beyond this (semi-)classical limit of large spin, S → ∞, little is known reliably about the properties of the pyrochlore quantum Heisenberg model. This is because the properties of the lattice conspire to frustrate not only magnetic order, but also attempts to apply standard theoretical and numerical approaches. The presence of a macroscopic number ('flat band') of gapless excitations in most bare models precludes standard perturbative schemes and mean-field theories. * dluitz@pks.mpg.de FIG. 1. Heat capacity per spin in different pyrochlore magnets, for a detailed account see Sec. IV. The red curve represents the converged part of our results from the NLCE for the S = 1/2 Heisenberg model. Results for a single tetrahedron for S = 1/2 (dashed) for S = ∞ (classical case, dots) have T scaled to match in the high-T limit. Symbols are for experiments on NaCaNi2F7 7 , a S = 1 (approximate) Heisenberg magnet, and the Ising spin ice Dy2Ti2O7 8 , both scaled in T such that their maxima coincide with that of the S = 1/2 model (with a factor ln 3/ ln 2 to account for the larger S = 1 entropy). Inset similarly shows entropy per spin.
This reflects the fact that fluctuations are typically very strong, the basic ingredient via which frustrated magnets avoid ordering. For this reason, even the relatively 'high' dimensionality, d = 3, often considered almost homologous with proximity to mean-field behaviour, is a hindrance rather than a help: the most unbiased method, exact diagonalisation, breaks down already for small linear system sizes L, as the Hilbert space dimension grows exponentially with L 3 . Similarly, DMRG-based methods are well-known to struggle beyond d = 1, while geometric frustration yields a sign problem in quantum Monte Carlo. In this sense, the pyrochlore magnet is even less tractable than the notoriously enigmatic kagome S = 1/2 Heisenberg magnet, for which decades of intensive interest have not yielded a consensus on the nature of its ground state. For the pyrochlore lattice, even a reliable ground-state energy estimate is lacking, and the proposed ground states depend strongly on the method used to study it. [13][14][15][16][17][18][19][20][21][22][23] The pyrochlore S = 1/2 Heisenberg magnet thus has many ingredients that make it one of the most likely candidates for realising new and exotic phases of matter, inter alia, quantum spin liquid states: it harbors both promise and obstacles, albeit in a somewhat imbalanced way.
In the view of all the abovementioned difficulties, the best source for information on these low temperature phases are experiments. In particular since scattering neutrons off a frustrated magnet is a priori no more involved than off an unfrustrated one, and d = 3 permits the straightforward study of bulk samples.
Up to now several compounds are known whose magnetic properties can be described by Heisenberg-like models on the pyrochlore lattice, [24][25][26] . While a good realisation of an S = 1 Heisenberg magnet is now available 7 , an isotropic nearest-neighbour-dominated S = 1/2 material still remains on the wish-list.
As experiments necessarily involve non-zero temperatures, theory is compelled to study this setting. Our study is therefore devoted to the S = 1/2 pyrochlore magnet at finite temperature. We focus on the thermodynamics -susceptibility and in particular specific heat, Fig. 1, and we also consider the spin correlations in the form of the momentum resolved static structure factor.
The remainder of this account is structured as follows. In Sec. I A we provide a summary of our results. Sec. II introduces the model and observables. The bulk of the technical advances are bundled into Sec. III, with material on canoncial typicality, the numerical linked cluster expansion, exact diagonalisation as well as DMRG. This may be skipped on first reading (as well as by the reader not interested in the underlying methodology). Our results on thermodynamics (specific heat, susceptibility in zero and nonzero fields) as well as spin correlators are presented in Sec. IV. We close with a broader discussion and an outlook in Sec. V.

Methods
Our quest to make progress has a considerable purely technical component. This involves efforts along two, principally computational, axes.
First, we devise a high-order numerical linked cluster expansion (NLCE) [27][28][29] for the pyrochlore lattice. This approach has been used before, and our contribution is to push the expansion -based on tetrahedral clusters well-suited to the corner-sharing tetrahedra of the pyrochlore lattice [30][31][32][33][34][35] -to significantly higher orders. We reach clusters of up to eight tetrahedra, involving the full solution of clusters of up to 25 spins. Going to such high order allows a significantly improved exploration of the low temperature regime, in particular permitting a more extended and controlled use of Euler transforms to extrapolate the results to low temperature.
Indeed, high expansion orders are essential to capture a range of physical processes. Concretely, up to the sixth nearest-neighbour hop, the pyrochlore (and, indeed, the kagome lattice) is equivalent to a Husimi cactus [a Cayley tree of tetrahedra (triangles)], and many series expansions are 'trivial' to high orders, e.g. with degeneracy lifting only occurring at eighth order in perturbation theory in a high-temperature kagome 36 or a strong-coupling pyrochlore expansion 37 . In a similar vein, the importance of resonance processes on more extended clusters is a recurring theme in the study of frustrated magnets. The clusters included in our expansion host not only the simplest hexagonal loop motifs but also loops of eight spins and longer decorated versions thereof (the longest loop consisting of eight tetrahedra cf. Fig. 15), which crucially encode the three dimensional structure of the lattice.
The second technical axis involves a finite-temperature DMRG analysis of the pyrochlore S = 1/2 magnet using finite clusters. Our results demonstrate that finitetemperature DMRG is a powerful approach, feasible even in this challenging three dimensional setup down to nontrivial temperatures. Here, we use a "snake" path through the lattice to map the system to one dimension with long range interactions. We take advantage of the SU(2) symmetry of the model and keep SU(2) block states up to χ = 10000 (∼ 40000 U(1) equivalent) and consider clusters up to 48 sites with periodic boundaries.
Taken together, these approaches permit us to reach converged results at temperatures down to around T = 0.25 (where the exchange constant of the Heisenberg model has been set to unity, Eq. 1) for the NLCE, and down to T = 0.6 for DMRG.

Observables
In the zero-field specific heat we resolve a pronounced maximum at T = 0.57. Crucially, at T = 0.25, where the specific heat has dropped well below its maximum value, the residual entropy is still around 0.33 k B , i.e. 47% of the value of a free spin of k B ln 2. This demonstrates the persistence of the spectral weight downshift characteristic of frustrated magnets to this case. We discuss implications of this observation in detail, in particular in comparison with the kagome magnet, as well as two simple tetrahedral models, on top of the experimental results on two pyrochlore magnets: the S = 1 Heisenberg antiferromagnet NaCaNi 2 F 7 7 , and the classic Ising spin ice Dy 2 Ti 2 O 7 8 . We find that our model at T = 0.25, Fig. 1, has a higher low-T entropy than all of these.
This entropy at T = 0.25 is in particular much greater than that proposed for singlet subspaces in resonating valence bond type effective theories. Indeed, in this regime, there is considerable admixture of triplet components in each tetrahedral wavefunction, reflecting the inability of neighbouring tetrahedra to be in singlet states simultaneously. In the magnetisation curves, this is reflected in a non-monotonic temperature dependence for fixed intermediate fields: upon cooling from the maximum magnetisation, the entropy of the magnetic excitations loses out to the singlet-dominated low-energy sector; while at high temperatures, the magnetisation assumes a conventional asymptotic 1/T behaviour. The maximum disappears at zero field, where there is no magnetisation in the absence of time-reversal symmetry breaking; and at high fields, where a conventional monotonic paramagnetic magnetisation curve is found.
For the magnetic-field dependence of the specific heat, we find a continuous drift of the location of its maximum to higher temperatures; at the same time, the amplitude of the maximum changes non-monotonically, first decreasing and then increasing again.
The spin correlators in turn, exhibit the by now familiar structure of incipient bow-ties, commonly found in various magnets on the pyrochlore lattice. These reflect the emergent gauge field and while they become arbitrarily sharp in the cases of classical magnets, their finite width indicates the presence of a nonzero net moment on the tetrahedra, on account of the abovementioned inability to have tetrahedra sharing a spin to be in a singlet configuration simultaneously.
The sum i, j in Eq. (1) runs over nearest neighbor bonds of the pyrochlore lattice. In the absence of a magnetic field h = 0, the model is SU(2) symmetric.
In this work, we focus on thermodynamic observables: the heat capacity at fixed volume C V , the magnetic susceptibility χ, the entropy S. We also consider the static spin structure factor S( Q). In the following definitions, we use the canonical ensemble averages • β = 1 Z Tr e −βH • . As usual, β = 1/T denotes the inverse temperature, and Z = Tr e −βH is the partition function. We use natural units k B = 1, = 1, J = 1.
The heat capacity is obtained either from the temperature derivative of the internal energy H β or from the fluctuations of the energy: Similarly, we obtain the magnetic susceptibility χ, defined by the change of the magnetization in z direction with respect to a change of the field h in z direction from the fluctuations of the magnetization: In the SU(2)-symmetric case, the susceptibility can be also expressed as: The thermodynamic entropy S can be calculated using the definition of the free energy: The static structure factor can be obtained from the Fourier transformation of the spin-spin correlations (the factor 4/3 stems from the normalization 1/(S(S + 1)) for spin S = 1/2): where R i denote the real-space coordinates of sites according to Eq. (2).

III. METHODS
This section is devoted to a detailed exposition of the methods used to obtain the results presented in the following section. It can safely be skipped at first reading, as well as by the reader primarily interested in the behaviour of the observables, rather than details on how they were obtained.

A. Canonical typicality
Calculating thermodynamic expectation values is possible via the density matrix ρ β = 1 Z e −βH , which can be calculated from all eigenvalues and eigenvectors of the Hamiltonian H. Due to the exponential scaling of the Hilbert space dimension with system size, this is impractical for systems with more than 25 spins (we discuss how to perform full diagonalization for such systems in Sec. III D). The concept of quantum typicality 38,39 permits a different approach, which has its foundation in Lévy's lemma. It can be summarized in the statement that for the vast majority of wavefunctions |ψ , the state |β = e −β/2H |ψ is typical [40][41][42] for the canonical ensemble.
In practise, this means that, starting from a random wavefunction |ψ , one can calculate finite temperature expectation values of observables O by The statistical error of this replacement is exponentially small in system size 43 and can be estimated (and reduced) by sampling over random (infinite temperature) initial wavefunctions |ψ . The application of e −(β/2)H corresponds to imaginary time evolution up to β/2 (the factor 1/2 stems from a symmetric splitting of the exponential) and can be carried out efficiently using Krylov space techniques 42,[44][45][46][47] , which is commonly known as the finite temperature Lanczos method 46 . The main advantage of this technique is that it can be carried out storing only the vectors of the Krylov space spanned by |ψ , while the application of the Hamiltonian to vectors during the Lanczos algorithm is either performed using a sparse matrix representation, or an on the fly generation of the corresponding matrix elements, making very large system sizes accessible which are comparable to Lanczos ground state calculations. We note that the same techniques for reducing the Hamiltonian to its symmetry sectors discussed in Sec. III D can be readily applied here.

B. Finite temperature calculations with matrix product states
Matrix-product-state (MPS) 48,49 based algorithms also provide a way to address the equilibrium thermodynamics of many-body quantum systems. [50][51][52][53] One of the most widely used methods is the purification of the finitetemperature density matrix. 50,53,54 The idea behind this approach is that one can interpret the density matrix, ρ P , as a partial trace of a Schmidt decomposition of a pure state, |Ψ , in an enlarged Hilbert space: where P and A denote the physical and auxiliary system, respectively. The state |Ψ can be easily constructed by simply creating a copy (auxiliary system) of the physical system and generating maximally entangled bonds between each physical site and its auxiliary site. This can be achieved by creating an entangler Hamiltonian whose ground state corresponds to the maximally entangled initial state. In our case this Hamiltonian simply reads: where a(i) denotes the auxiliary site belonging to site i and the sum is performed over the physical sites.
One can easily see that the density matrix calculated from this ground state, |Ψ β=0 , corresponds to infinitetemperature. Any finite-temperature density matrix can be obtained by performing an imaginary time evolution on the physical system and tracing out the auxiliary degrees of freedom: As a matter of fact, any expectation O β can be directly evaluated from |Ψ β : This provides a great advantage of this method, since thermodynamic quantities can directly be obtained by simulating the density matrix, rather than averaging over low entanglement pure states 51,52 , and therefore results are free of statistical errors. By design, this method is most efficient in one dimension. In order to use it for a three-dimensional system, one has to put a 'snake' path through the lattice sites to map the original problem to a one-dimensional equivalent one, which contains longrange couplings between the lattice sites. This is the main difficulty of this approach, since the MPSs need to encode a large amount of entanglement, i.e. if we think of the area law for a moment (valid only for ground states), the bond dimension should scale exponentially in 2 ( is the linear size of the 3D system) to accurately represent the many-body state. This obviously limits the feasible system sizes and the accessible temperatures. The presence of long-range couplings in the one-dimensional topology poses another difficulty regarding the choice of the time evolution method. 55 The timeevolving block decimation (TEBD) 56,57 is very effective if the couplings are short-ranged, otherwise one has to subsequently apply a series of swap gates to move distant sites next to each other so that the time evolving operator can be applied. This swapping procedure is extremely slow and becomes very inefficient as the bond dimension is increased. The Krylov method 58 is capable of handling long-range interactions by default, but already in the early stages of the imaginary-time evolution the Krylov vectors become strongly entangled, making the calculation unfeasible. To reach physically relevant temperatures, we demonstrate that the time-dependent variational principle (TDVP) 59,60 provides an effective way. Although it introduces another source of error by projecting the evolution vector onto the MPS manifold, this error is usually much smaller than the truncation error. At this point we also have to mention that it is not straightforward to evolve |Ψ β=0 directly with TDVP. 55 Since the initial state is essentially a product state a naive evolution with TDVP would be simply wrong due to the loss of long-range interactions already in the first projection step. To overcome this difficulty we apply the same trick that has been successfully applied in real-time evolution, 55 namely, we generate an initial state with an artificially enlarged bond dimension. This is achieved by finding the ground-state of the Hamiltonian: where the parameter a is being varied. We start with a = 1 and perform 20 sweeps then we reduce it by a factor of ten each time. During the sweeps single-site version 61 of the density-matrix renormalization group (DMRG) algorithm 49,62-64 is applied with subspace expansion and setting the truncation error to zero. Five stages are performed altogether. In the last stage we set a = 0 and perform three additional sweeps. We demonstrate that this procedure removes the above bottleneck of TDVP also for imaginary-time evolution. In addition, to encode the large amount of entanglement, the compression of the many-body states must be very efficient. To this end we exploit the SU(2) symmetry of the model and keep block states up to ∼ 10000 (∼ 40000 U(1) equivalent) to minimize the truncation error as much as possible. 65,66 It is worth mentioning that higher bond dimensions can be achieved in a ground-state search, 67 where one can take advantage of single-site DMRG as well as parallelizing the computation in real space to reduce memory usage and computation time respectively. In our case, however, the two-site TDVP update-scheme needs to be used and the serial solution of the TDVP equations is crucial.

C. Numerical linked Cluster expansion
For studying three dimensional frustrated quantum magnets most controlled algorithms are restricted to a small number of spins and thus have a hard time capturing the three dimensional structure and its correlations. Using a systematic high temperature series expansion opens up the possibility to obtain reliable results in the thermodynamic limit. The numerical linked cluster expansion (NLCE) is able to determine any extensive property P in the high temperature regime. NLCE has been applied to various geometries like the square lattice, kagome lattice or pyrochlore lattice and has provided new insights in these systems [27][28][29]31,68,69 such as the transition of different phases in quantum systems 33,35,70 or a deeper understanding of real materials 30,32,34 . Moreover, the generality allows the application of this algorithm to a variety of other systems [71][72][73][74][75] .
In general, the systematic expansion can be applied to any lattice, the crucial part is the choice of building block, which builds up the infinite lattice by translational symmetries t. All generated configurations are extended by adding the building block in each step of the expansion. There are two main aspects that need to be considered for the choice of building block. First, the number of generated clusters scales exponentially with the order of the expansion. Second, the complexity of solving these clusters scales exponentially with systems size, which limits the maximal size of practically solvable clusters. Choosing a building block with a large number of sites induces a relatively small number of clusters which are still solvable; a building block with a small number of sites induces a very large number of solvable clusters to a degree that one may not be able to reach the largest solvable cluster size. Common choices of the building block in the square lattice is a single site or a complete square. Physically motivated, most NLCE approaches in the pyrochlore lattice use a tetrahedra expansion, based on clusters of complete tetrahedra, such that no dangling spins or triangles occur. In the following, we discuss the cluster expansion for multisite unit cells and compare expansions in the pyrochlore lattice based on three different building blocks, in particular the single site expansion, the unit cell expansion and the tetrahedra expansion. The latter turns out to yield the most reliable results and represents the optimal expansion in that our results include full exact diagonalization of all clusters consisting of up to 8 tetrahedra. The largest included clusters thus consist of 25 spins 1 2 , which host crucial loops of 6 and 8 spins in the lattice.

Basic recipe
NLCE generates all possible subclusters c (subject to the choice of building blocks) which are embedded in the infinite lattice structure L and contribute to the thermodynamic limit P (L)/N per site. The contribution of each cluster c is given by its weight W P (c), describing the new (i.e. not included at lower order) contribution of c to P , and its multiplicity L(c), describing the number of possible embeddings of c in L.
The generality of this idea allows the definition of various building blocks. For the pyrochlore lattice studied here, we include a) all clusters built from single lattice sites b) all clusters built from complete (tetrahedral) unit cells and c) all clusters built from complete tetrahedra.
All possible configurations of clusters, subject to the choice of building blocks, are expanded further in each step of NLCE. The initial cluster given by the building block (or multiple clusters given by inequivalent building blocks) have to respect translational symmetries t and cover the whole infinite structure L by applying these symmetries t. It is important to point out the crucial difference between translational t and non-translational symmetries s such as rotation or reflection.
Each step of the expansion (at expansion order n) generates a set of connected clusters C n . This large set can be reduced to the set of clusters which are not related by lattice symmetries S n , and subsequently to topologically distinct clusters T n of size n (number of building blocks). The set of connected clusters C n includes all possible clusters which are embedded in L that are not related by any translational symmetry t to each other. The set of clusters not related by lattice symmetries S n are given by a subset of C n ; each cluster in S n is neither related by translational t nor non-translational symmetry s to each other. Applying all non-translational symmetries s to S n generates all connected clusters C n . Moreover, the set S n can be reduced further. Even though clusters are not related by any symmetry they can exhibit the same interaction topology and hence, generate the same Hamiltonian matrix. We describe each cluster by its interaction graph G, where its nodes i ∈ N G correspond to the spins included in the cluster and its edges (i, j) ∈ E G correspond to nearest neighbor interaction terms of the Hamiltonian (1). Two clusters are topologically equivalent if there is a graph isomorphism π : N G1 → N G2 (bijective) mapping G 1 on G 2 while preserving its structure; that means if (i 1 , j 1 ) ∈ E G1 is an edge of G 1 then (π(i 1 ), π(j 1 )) ∈ E G2 needs to be an edge of G 2 . Hence, the set of topologically distinct clusters T n is a subset of S n . Using building blocks including more than one lattice site (like the tetrahedron) requires to check the topological structure on the full connectivity graph including all sites.
The multiplicity L(c) assigned to each cluster c describes the number of possible embeddings in the infinite structure L. All possible clusters (subject to the choice of building blocks) of size n are included in C n ; hence, the multiplicity of each cluster is one. Non-translational symmetries s reduce the set of connected cluster to S n and summarize multiple clusters in C n to one representative cluster c ∈ S n . The multiplicity of each cluster in S n is given by the number of non-translational symmetries s that transform the cluster to another cluster that is not related to the first one by any translation t. Again, multiple clusters in S n are summarized to one representative cluster c ∈ T n , its multiplicity given by the number of topologically equivalent clusters. Hence, the multiplicity of topologically invariant clusters is simply the sum of the multiplicities of all topologically equivalent clusters in S n . Summing all multiplicities of clusters in S n or T n equals the number of connected clusters: For clarity, we will drop the index top in what follows, i.e. L(c) = L top etc. The basic recipe to expand the clusters by one building block (n → n + 1) is equivalent for all geometries and building blocks: i) Starting from all clusters not related by lattice symmetries of size n in S n , we generate new clusters by adding a building block to every free nearest neighbor. Again, we only consider clusters that are distinguishable by translational symmetries t.
ii) Non-translational symmetries s are used to reduce the set of all newly generated expansions to create the set S n+1 . Applying all non-translational symmetries s to S n+1 generates the full set of connected clusters C n+1 with respect to translational equivalence.
iii) Clusters in S n+1 are further reduced to obtain all topologically distinct clusters in T n+1 .
Each cluster in T n contributes to the expansion according to its multiplicity L(c) and weight W P (c). The nth order of NLCE is given by: The weight assigned to a cluster c is defined with respect to all smaller subclusters s ⊂ c (subject to the choice of building blocks) which can be embedded in c; hence, it extracts contributions of c to P which are not covered by smaller clusters: In practice, the thermodynamic observable P needs therefore to be calculated for all topologically invariant clusters. The extensive property of P induces a zero weight for disconnected clusters 29 which do not have to be considered in (16). Expanding clusters only by nearest neighbors guarantees connected clusters in C n , S n and T n . The recursive definition of the weight ensures the convergence towards the infinite structure L which is the thermodynamic limit. It is not necessary to use all (or any) non-translational symmetries in the NLCE. A lower number of nontranslational symmetries increases the number of clusters |S n | such that each cluster has a lower multiplicity. In fact, it is possible to consider the identity id as the only non-translational symmetry, then C n = S n . However, checking the topological structure of these clusters generates the same set of topologically distinct clusters T n . The computational effort can increase drastically with a low number of symmetries, since the number of clusters grows exponentially.

Pyrochlore lattice and building blocks
The underlying Bravais lattice is given by a fccstructure with a tetrahedral unit cell (four sites). Hence, in order to respect the translations, all expansions we use focus on the tetrahedral unit cell and converge to the s1: id s7: 2 (0,ȳ, y) s2: thermodynamic limit per unit cell and not per site, thus accounting for an extra factor of four. Its symmetries are described by the space group F d3m (227); it contains 192 symmetry operations. However, only 12 symmetries are purely non-translational and used in the NLCE, see table I. s 2 , s 3 and s 7 , s 8 , s 9 describe three and two folded rotations, respectively. Reflections are given by s 4 , s 5 , s 6 and s 10 represents the inversion. s 11 and s 12 combine a threefold rotation with the inversion.
As discussed earlier, the choice of building block is crucial. In principle, various geometries (such as dimers, hexagons or multiple tetrahedra) can be used as building blocks as long as they respect the translations t and cover the whole lattice. Each expansion we use is embedded in the fcc structure of the pyrochlore with either equivalent (unit cell expansion) or multiple inequivalent (single site and tetrahedra expansion) building blocks.
An efficient implementation of translational symmetries is essential due to the exponentially increasing complexity. Labeling lattice sites along the translational axes automatically describes the translational symmetries t by a simple index shift. a) single site expansion: As pointed out earlier, the single site expansion generates a large number of clusters of relatively small size. The advantage of this approach is its complete generality. In contrast to the single site expansion in Bravais lattices (e.g. square or triangular lattice 28,29,68 ), the unit cell in the pyrochlore lattice consists of four sites, which need to be treated inequivalently. Hence, the starting point of the single site expansion are four sites arranged in the unit cell/tetrahedron, which covers the whole lattice by translations. Applying translation symmetries to these four sites generates the full pyrochlore lattice. Here, all symmetries in table I can be applied to find clusters which are related by lattice symmetries to reduce the complexity.
b) unit cell expansion: The unit cell expansion is related to the single site expansion in the fcc-lattice, substituting each site in the obtained clusters by the tetrahedral unit cell. However, working within the pyrochlore lattice, the symmetries in table I have to be examined more carefully: We require the symmetries to preserve the unit cell structure such that only entire unit cells are mapped to each other. Only the symmetries s 1 to s 6 in table I, which are a subset of the fcc-lattice symmetries, are unit cell conformal. As mentioned before, working with a lower number of symmetries produces the same results. Since the building block includes more than one site, the topological structure has to be compared on the level of the full connectivity graph including all lattice sites. The advantage of this approach is the consideration of larger clusters due to a much slower growth of the number of clusters with the number of unit cells. However, the connection between the unit cells are dangling bond that do not reflect the geometrical frustration. In the presence of magnetic fields, the Hamiltonian has bond and site terms. Therefore, we include a single site (yielding the simplest contribution of the site terms) as the 0th order in the expansion (m = 0) in (15), which is embedded four times in the unit cell.
c) tetrahedra expansion: The central motif of the pyrochlore lattice is the tetrahedron. An examination of the lattice shows that there are two types of tetrahedra: up pointing tetrahedra (these are the unit cells) and down pointing tetrahedra, which correspond to the interaction of each spin in the unit cell to three neighboring spins in different unit cells. Both the single site and unit cell expansion do not respect this structure, leading to dangling bonds in the case of most clusters in the single site expansions and incomplete down pointing tetrahedra in the case of the unit cell expansion.
For this reason, we use an expansion including all clusters with complete tetrahedra [30][31][32][33][34]70 . The systematic expansion focuses on an hourglass structure composed by two inequivalent building blocks of thetrahedra (up/down pointing) which are placed in the underlying fcc-lattice and expands these as described before. Our comparison of results for the heat capacity and the magnetic susceptibility demonstrates that this intuition is correct and that this expansion is indeed superior (cf. comparison in Appendix A).
It is important to note that the physical size (number of spins) of each cluster is not uniquely related to the order of the expansion, due to an overlap of up-and down-facing tetrahedra. This means that at an expansion order (given by the total number of tetrahedra) n tetra clusters of different sizes are included. Similarly to the unit cell expansion, this expansion leads to a relatively small number of large clusters. Again, we need to consider the 0th order contribution of a single site. We did not apply any symmetry due to the small number of edges of each tetrahedron and the low order of expansion and rely on the automatic identification of topologically equivalent clusters by directly comparing their interaction graphs.
In appendix B we provide a comparison of the number of clusters generated at each order in the three expansions discussed here. We list the number of connected |C n |, not related by lattice symmetries |S n | (if present) and topologically distinct clusters |T n | in tables II, III and IV. The visualization of these results is shown in fig-ure 12, clearly showing that the number of clusters with at most n sites is smallest in the tetrahedra expansion, leading to a tractable number at the edge of feasibility of full diagonalization (the maximal symmetry block dimension of the Hamiltonian is 228592) for clusters with 25 sites (full Hiblert space 3.355 · 10 7 )). Additionally, all topologically distinct clusters of the tetrahedra expansion can be found in the appendix B.

Resummation algorithms
Correlations are increasingly long ranged as temperature is lowered. Hence, contributions of larger cluster have to be taken into account to converge to the thermodynamic limit. Accessing these orders is limited due to the exponentially increasing complexity regarding the number of clusters or Hilbert space dimension.
Resummation algorithms rely on a systematic usage of lower orders of the series and are most effective if many terms are included. They are guaranteed to converge to the limiting value of the underlying series, or not at all. In this work, we use the Euler transform of our NLCE data for the expansion up to n = 8 tetrahedra and compare the highest order Euler transform of the Euler transform up to expansion orders n = 7 (and also n = 6), to ensure convergence of our results.
Euler's transformation is particularly useful for alternating series, which are transformed according to 76 where ∆ s is the s fold application of the forward difference operator ∆, defined by ∆u n := u n+1 − u n . Euler's method can be derived by repeated application of summation by parts 77 .
In the present work, we find that the Euler transformation of our NLCE series up to n = 8 tetrahedra (i.e. containing 8 terms) yields a significant improvement of convergence at low temperatures. We always compare the Euler transformation of the first n = 7 and n = 8 terms of the series to ensure that the results are indeed converged, yielding reliable results down to T ≈ 0.25 for thermodynamic properties of the S = 1/2 pyrochlore antiferromagnet as shown in detail in what follows.

Cluster symmetries from graph automorphisms
The numerical linked cluster expansion (NLCE) expresses thermodynamic observables as a series expansion in terms of the exact solution of a large number of finite size clusters. Since the number of clusters grows factorially with the number of constituents, it is useful to use an automatic strategy to identify cluster symmetries, which are used for block diagonalizing the Hamiltonian. Here, we provide a practical description of the method with only minimal reference to graph and group theory to make it accessible. A similar pedagogical description for the exploitation of translational symmetry can be found in Ref. [78]. In all of the following discussion, we use the computational S z basis, in which each basis state |σ 1 , σ 2 , . . . is an eigenstate of all local S z i operators and labelled by their eigenvalues σ i .
As described in the NLCE part, we identify a finite cluster with its interaction graph G, where its nodes i ∈ N G correspond to the spins and its edges (i, j) ∈ E G correspond to nearest neighbor interaction terms of the Hamiltonian (1). Hence, the Hamiltonian is defined by: The sum (i, j) ∈ E G runs over all edges (i, j) of the graph G, the sum i ∈ N G runs over all sites. We notice that any automorphism of the graph G leaves the Hamiltonian invariant, since an automorphism is a permutation π of nodes which maps the graph onto itself, such that for any edge (i, j) ∈ E G , the mapped edge has to be in G as well: (π(i), π(j)) ∈ E G . This means that graph automorphisms are symmetries of the Hamiltonian and commute with it, which implies that we can simultaneously diagonalize the automorphism and the Hamiltonian.
The matrix representation A of a graph automorphism A (which necessarily is a permutation π A of graph nodes) can be obtained by noticing that any basis state is transformed as A is a permutation matrix on the set of basis states and A T A = 1, i.e. A is orthogonal with eigenvalues on the complex unit circle. Since A is a permutation matrix, it is idempotent with a certain order N A < |N G |: A N A = 1. Therefore, its eigenvalues are given by the N A -th roots of unity: e i2πn A /N A , with n A ∈ {0, . . . , N A − 1}.
This means, we can block diagonalize H into n A blocks, labelled by the eigenvalues e i2πn A /N A of A. We denote the order o(π) of an automorphism π by the minimal number n ∈ N such that π n = 1.

Identification of largest commuting automorphism subgroup
Typical interaction graphs G have a large number of (independent) graph automorphisms. However, typically not all of them commute! Since only commuting graph automorphisms can be used together to further reduce the block size of the Hamiltonian, we want to find the largest abelian (commuting) subgroup of the complete automorphism group of G.
In order to do so, we start by generating the automorphism group of the graph. In the next step, we check for each pair of automorphisms if they commute. This can be interpreted by a new graph C, in which each automorphism of G corresponds to a node and two nodes are connected if and only if the corresponding automorphisms commute. What we are looking for is a subgroup U ⊂ N C in which each node in U is connected to all other nodes of the subgroup. Hence, each element in U commutes with each other. This is called a clique in graph theory. Finding the largest abelian subgroup of the automorphism group is therefore identical to finding the largest clique in C. In general, the largest clique is not uniquely determined.
After determining the largest clique U of automorphisms, we need to identify a minimal set H of independent generators of the abelian subgroup. Each element of the clique U has to be generated uniquely by elements of H. Assume, H includes m := |H| elements h ∈ H of order o(h), then there is exactly one tuple (n h1 , . . . , n hm ) of integers for each element u of the subgroup, which generates u from the corresponding integer powers of the generators: ∀u ∈ U : ∃! (n h1 , . . . , n hm ) ∈ N m with 0 ≤ n hi < o(h i ), The product in equation (20) refers to the composition of permutations. The ordering is arbitrary since all generators commute. Each element in U is represented uniquely by a multiindex (n h1 , ..., n hm ) ∈ N m where each number is smaller than the corresponding order of the generator. This multiindex determines a phase in the symmetrized basis (cf. Sec. III D 3). The bijective mapping from U to N m , respecting the order o(h i ) of the generators h i , induces the following relation between the cardinality of the abelian subgroup |U| and the orders of its generators: In practice, we create the minimal set of generators by starting with the elements exhibiting the highest order in the commuting subgroup U. First, one element with the highest order will be added in H. A new element g is added to H if it does not violate the bijective mapping described through (20). That is, all possible compositions g n •h, for n = 0, ..., o(g)−1 and h ∈ H, generate uniquely a new element in the subgroup U. The generating set H is not uniquely determined. In Fig. 2, we show an illustration for an example interaction graph (hourglass composed of two corner sharing tetrahedra, corresponding to the n = 2 cluster in NLCE) along with its 71 nontrivial automorphisms and their commutation relations. The largest clique for this example has 8 nontrivial automorphisms with two generators A and B [e.g. given by the node permutations π A = (2, 0, 1, 3, 5, 6, 4), π B = (2, 0, 1, 3, 6, 4, 5)], generating independent C 3 rotations of the upper and lower tetrahedron. In the present work, we rely on established algorithms for the identification of graph iso-and automorphisms bundled in the package nauty 79 and used the clique maximization algorithm described in Ref. 80 .

Symmetrized basis
Once we have obtained a set of independent generators H and the multiindices referring to the largest abelian subgroup U of the graph automorphism group, we can proceed with the block diagonalization of the Hamiltonian. A subgroup of size |U| induces the same number of blocks, each is uniquely identified by m (number of minimal generators) quantum numbers given by multiindices described in (20). Each index refers to the phase of the eigenvalue of the corresponding generator; the commutation relations of the generators allows the simultaneous diagonalization of all generators and the Hamiltonian. Each computational basis state | σ = |σ 1 , . . . σ N has to be replaced by a symmetrized state induced by the quantum numbers ε = (ε 1 , ..., ε m ) ∈ N m which is given by where N ε σ is the normalization constant of the state, n u = (n u 1 , ..., n u m ) ∈ N m is the multiindex referring to u ∈ U defined by (20), and o(h i ) is the order of the generator h i ∈ H. The eigenvalue of each generator h i ∈ H is obtained (using the properties of the generators and the complex roots of unity) from It is important to note that multiple unsymmetric basis states | σ typically generate the same (apart from a phase) symmetric state. This is in fact true for any basis state which is in the set which we call the "family of the state | σ " generated by the commuting subgroup U of the graph automorphism group. Therefore, each symmetric basis state has to be added only once to the symmetric basis. This is typically ensured by using a parent state of the family, for example the state | σ with the lowest binary representation; this parent state is denoted by p(F σ ). Crucially, some unsymmetric basis states | σ do not generate any symmetric state in a given symmetry sector. This happens if the basis state is incompatible with the symmetry sector and the sum of phase factors cancels, leading an unnormalizable state. An example is the state | ↑↑↑ . . . ↑ . For any graph and any automorphism u, it is mapped to itself. Therefore, its family is F ↑↑...↑ = {| ↑↑↑ . . . ↑ }. In the symmetry sector 0 m := (0, . . . , 0), we obtain | ↑↑ . . . ↑; 0 m = | ↑↑ . . . ↑ . In all other sectors, however, we get | ↑↑ . . . ↑; n = 0 m = 0, i.e. this state does not appear in other sectors. As in the general example, this mechanism leads to an imbalance of the size of sectors, the sector 0 m always being the largest.
It is crucial to ensure the correct normalization of the symmetric basis states | σ and | σ introduced in Eq. (22). We require: i.e. states are orthogonal if they correspond to different parents (and therefore families), or if they correspond to different symmetry sectors n. We note that in addition to the graph automorphisms, the Heisenberg model we study has additional spin symmetries. Here, we exploit the conservation of the total z component of the spin, because [ i S z i , H] = 0. Since all computational basis states we use are already eigenstates of the total z component i S z i , and since i S z i also commutes with the graph automorphisms, this symmetry is trivial to exploit: A simple reordering of basis states by their z magnetization makes the Hamiltonian block diagonal. Additionally, we exploit the spin inversion symmetry, [Q, H] = 0 with Q := i S x i in the sector m z = 0. The spin inversion is also used to deduce the results for the −S z sector from an already solved S z sector.

Hamiltonian submatrix in symmetry sectors
Before we can fully diagonalize the Hamiltonian, we need to represent each block of H (labelled by the quantum numbers ε) in the symmetric basis. For simplicity, we only focus on the graph automorphisms and ignore symmetries defined by S z and Q (which are a trivial extension). Hence, we need to construct the matrix elements σ; ε |H| σ ; ε ; note that by construction the interblock matrix elements ε = ε are zero.
Let us apply the Hamiltonian to a symmetrized basis state, exploiting the fact that H commutes with u ∈ U: The Hamiltonian is expressed as a sum of non-branching terms: H = b h b ; note that the permutations u ∈ U commute with each single term h b . The operators h b can be divided into diagonal operators (which do not change the state | σ ) and off-diagonal operators, which yield a different state | σ b together with the corresponding matrix element by is in general not a parent state. It is however related to its parent state by a symmetry operation | p b = u 0 | σ b = p(F σ b ). Note that u 0 is not determined uniquely; multiple permutations can fulfill this mapping. Also, the new state is not necessarily a valid state in the symmetry sector ε, in which case it will be canceled by later terms. The parent state is assigned to an index of the basis; the referring matrix element can be calculated as follows: Terms in Eq. (27) only have a non-zero contribution if and only if | σ b is mapped by u u to its parent state | p b . The matrix elements of the Hamiltonian in the symmetrized basis are therefore given by the unsymmetrized matrix elements h σ σ b b multiplied by symmetry sector dependent phase factors and normalization constants. We note in passing that the sums over phase factors can yield the normalization constant N ε p b 78 .

IV. RESULTS
Using the combination of the methods described in Sec. III, we address thermodynamic properties of the pyrochlore quantum Heisenberg antiferromagnet. We start by considering the SU(2) symmetric case without an applied magnetic field and present our results for the heat capacity (IV A), the magnetic susceptibility (IV B), and the thermodynamic entropy. We compare the results for different orders in the numerical linked cluster expansion (NLCE) and their Euler transform to show that the high temperature regime down to T ≈ 0.25 is converged to the thermodynamic limit. We furthermore compare the NLCE results to the solution of finite size clusters obtained using canonical typicality and finite temperature DMRG.
We also present results obtained from DMRG for the static spin structure factor at finite temperature (IV C), as well as NLCE results for the heat capacity and the magnetization at finite applied magnetic field h (IV D). To accelerate the convergence of the NLCE series at lower temperatures, we apply the Euler transformation up to order n in the series (cf. III C 3). We have furthermore calculated the specific heat capacity for a larger cluster with N = 48 sites and periodic boundary conditions using our SU(2) symmetric finite temperature DMRG method with finite bond dimensions ranging from χ = 2000 to χ = 10000. These results were obtained from a numerical derivative of the (spline interpolated) energy as a function of inverse temperature. For T > 2, the results are converged with bond dimension and agree with the NLCE and finite N = 32 cluster results (bottom panel). At lower temperatures, the dependence of the results on the bond dimension becomes significant and we extrapolate to infinite bond dimensions using a quadratic polynomial in 1/χ yielding a very good match with the N = 32 and NLCE results within the accessible temperature range (T > 0.62). The errorbars indicate the distance of the extrapolated value from the largest bond dimension (χ = 10000).

A. Heat capacity in zero field
The specific heat capacity quantifies the change of the internal energy as a function of temperature and is directly accessible in experiment. Our results and a comparison to experimental results are summarised in Fig. 1.
Starting at high temperatures, T 1, the heat capacity decays as 1/T 2 , as required in the leading order high temperature expansion. In the regime down to T = 2, all orders of the NLCE with n > 5 agree with each other, and also with the finite size N = 32 result from typicality. The Euler transform of our NLCE results agrees between orders n = 7 and n = 8 down to T ≈ 0.25, which we take to indicate that the series is converged over this range.
Crucially, this allows us to resolve unambiguously the maximum of the heat capacity located at T ≈ 0.57. In the proximity of the maximum, the results for the finite size cluster N = 32 deviate significantly from the NLCE, which indicates that in the regime T < 2 the correlations beyond the size of the N = 32 cluster start playing a discernible role. From the strong dependence of the DMRG results on the bond dimension around the location of the maximum of the heat capacity, as well as from the visible discrepancy of the specific heat for finite size clusters compared to the converged NLCE results close to the maximum, we conclude that the system enters a nontrivial quantum regime at temperatures T ≈ 1, where it exhibits entanglement beyond what is representable faithfully by χ = 10000 matrix product states.   The heat capacity of the N = 32 cluster exhibits a second maximum at low temperatures, similarly to what was observed previously in a different pyrochlore model 81 . However due to the divergence from converged results of the NLCE in this regime, which is not subject to finitesize effects of this kind, we conclude that this feature is likely not representative for the thermodynamic limit.
Indeed, the converged part of our NLCE data shows a rapid decrease of the heat capacity as T is lowered from the maximum. Therefore, if there is an additional feature, it must be well separated from the maximum we have found.
In order to gain further insight into the low-T regime, we have calculated the thermodynamic entropy as a func-tion of temperature, with S(∞) = ln 2. A direct calculation of the entropy per site in NLCE using Eq. (6) shown in Fig. 4 indeed agrees with this temperature integral down to the lowest temperatures for which our (Euler transformed) NLCE is converged. Interestingly, we find that just over half the total entropy is released down to T ≈ 0.25, where the entropy is S/N ≈ 0.33 ≈ 0.47 ln 2. This in turn means that the spectral weight below the maximum is huge, and there is plenty of scope for further interesting behaviour, the nature of which we are unfortunately unable to determine from our approach. [In several classical models 6 , as well as some fine-tuned quantum models [82][83][84] , not all the entropy is released even at T = 0, but in real systems, the third law of thermodynamics stipulates that this is not what actually happens. For instance, for the finite size N = 32 cluster, a steep decrease at low T 0.1 is associated with its low-T peak in the specific heat.] We will return to more detailed comparisons of this behaviour with other models and experimental systems in the final discussion, Sec. V.

B. Magnetic susceptibility at h = 0
We consider the magnetic susceptibility χ/N 85 per lattice site in Fig. 5 as a function of temperature. As in the case of the heat capacity, we perform NLCE calculations up to clusters of 8 tetrahedra and apply the Euler transform to these results. Fig. 5 shows the raw NLCE results at order n = 8 to be converged down to temperatures of about T ≈ 0.8. The Euler transform improves the convergence of the series significantly, again down to a temperature of T ≈ 0.25.
At high temperatures T 1, the susceptibility obtained by the various methods agrees: typicality for N = 32 site cluster, the N = 48 cluster result obtained in finite temperature DMRG, extrapolated to infinite bond dimension using a quadratic polynomial in inverse bond dimension. At T ≈ 0.6, the finite size magnetic susceptibility exhibits a pronounced maximum and decays rapidly to zero at lower temperatures. The Euler transform for the largest NLCE order clearly reveals a decrease of the magnetic susceptibility after a maximum at T = 0.54 in the thermodynamic limit.
We note that the magnetic susceptibility for very large system sizes was previously calculated in diagrammatic Quantum Monte Carlo simulations in Ref. 86, corresponding to the black crosses with errorbars in Figs. 5. These results agree with our NLCE and finite size results at temperatures above the maximum of the susceptibility. At low temperatures, however, they suggest a steady increase or plateau of the susceptibility instead of the maximum that we find. It would be desirable to push both our and the diagrammatic Monte Carlo method to higher orders in order to see which of the two apparently irreconcilable behaviours is the correct one.
C. Static spin structure factor at zero field The static spin structure factor quantifies the spin correlation patterns present at a given T : Here, we use finite temperature DMRG on the two clusters with 32 sites and full cubic symmetry and 48 sites with reduced symmetry, to investigate the static spin structure factor at finite temperature. Fig. 6 shows the result in the 32 site cluster (top rows) and for the 48 site cluster (bottom rows), in the the (h, h, l) plane (i.e. Q x = Q y , rows 1, 3) and the (h, l, 0) plane (Q z = 0, rows 2, 4) of momentum space, extended over several Brillouin zones for different temperatures (columns). Both figures show the clear emergence of a correlation structure already at high T . This becomes more pronounced as T is lowered, without acquiring much additional structure: certainly, as expected for a highly frustrated magnet, no sharp Bragg peaks appear, which would have been indicative of magnetic ordering.  Indeed, with decreasing temperature, the weight in the center of the Brillouin zone ( Q = 0) decreases and moves to the boundary of the extended Brillouin zone, the most visible location of increasing intensity being at (0, 0, ±4π).
Due to the low resolution in k-space our finite system sizes up to N = 48 do not permit to investigate more closely how sharp the pinch points become, since the associated lengths at low temperatures may become longer than our linear cluster sizes. However, a very rough idea can be gleaned by considering the T -dependence of the energy. The reason for this is that the energy simply encodes the total spin of each tetrahedron, The pinch-points become infinitely sharp in the limit of vanishing tetrahedral magnetic moment, i∈tet S i = 0. 87 This condition is known to be met in classical Ising (where it is known as the ice rule 4,6,89 ), XY and Heisenberg models, but it cannot hold in the quantum realm, since a spin which is part of two tetrahedra cannot enter singlet bonds in both of them simultaneously. The deviation of the energy from the minmal value for is thence a proxy for the possible sharpness of the pinch-points. The total spin of a tetrahedron is plotted in Fig. 7 from our various methods. Eyeballing an extrapolation of this curve to T = 0, one finds i∈tet S i ≈ 1, which confirms that frustration precludes that all tetrahedra are in singlet states and hence a finite width of the pinch points; for an estimate of the pinch-point width based on PF-FRG, see Ref. 23.

D. Heat capacity and entropy in a magnetic field
We complement the above discussion by a further analysis of the behavior of the pyrochlore Heisenberg magnet in the presence of a finite field. For readability of our figures, we only show the converged part of eighth order NLCE results using the tetrahedra expansion. As before, we use the agreement of eighth and seventh order Euler transform as convergence criterion. Fig. 8 shows the heat capacity and entropy per site as a function of temperature for a range of different fields h applied in the [001] direction. We observe a shift of the maximum of the specific heat to higher temperatures at  strong fields, as well as a non-monotonic change of the height of the maximum.
An overall upward shift of the weight is not particularly surprising given the presence of an additional term in the Hamiltonian. Indeed, at high fields, the curve resembles an unspectacular paramagnet.
However, the structure of the low-energy spectral weight and its rearrangement at intermediate fields is complex. The weight at lowest energies is in large part due to non-magnetic states which are not favoured by the magnetic field. At the same time, the more numerous states with nonzero magnetisation spread out as the field is applied. The failure of the Euler transform to converge to similarly low temperatures at intermediate fields is presumably due to a more complex behaviour of the specific heat in this intermediate field regime.

E. Magnetization at h > 0
We finally consider the effect of a magnetic field on the magnetization per site m z /N . At zero field, there is no net magnetization in the absence of spontaneous symmetry breaking. Fig. 9 shows converged NLCE results (with the highest order Euler transform) for the finite temperature magnetization (solid lines) in comparison with the result for a single tetrahedron (dashed lines). At high temperatures, these results agree and yield a Curie law. At intermediate temperatures, a difference due to the finite size of the tetrahedron is noticable, with the selection of different total magnetization groundstates in the low-T limit, with m z /N being 0, 0.25, 0.5 respectively depending on the field 90 .
Interestingly, a nonmonotonic dependence of the magnetisation on T can be observed at low fields, Fig. 9, both in the NLCE results (emphasized in the inset) as well as -more visibly -in the single tetrahedron case: the magnetisation vanishes at both low T , because the weak-field low-energy states are dominantly non-magnetic; and at high T , for the usual entropic reasons. At intermediate T , by contrast, the large weight of magnetic states oriented by the fields dominates, whence the maximum.

V. DISCUSSION AND EXPERIMENT
Having laid out the results, we now place them in the broader context of other highly frustrated model systems on one hand, and experiments on magnetic materials on the other. Here, we focus on the specific heat, not only because it is a quantity which is readily available across the board; but also because it allows for a relatively straightforward comparison between different platforms thanks to the integral 1 We have collated the data for a number of models and materials in Fig. 1. There, we compare (i) our converged results with what is found for single tetrahedra with (ii) S = 1/2 or (iii) in the classical limit, S = ∞; as well as experiments on (iv) the Ising spin ice pyrochlore magnet Dy 2 Ti 2 O 7 and (v) the S = 1 Heisenberg antiferromagnet NaCaNi 2 F 7 , rescaled by ln 2/ ln 3 to take into account the greater spin length. The T -axes of the experimental data, (iv,v), have been scaled so that the temperature of their respective maxima coincide with the one of our data; while the single tetrahedron results (ii,iii) were scaled to agree in the asymptotic limit of high-T .
All of these have in common a considerable spectral weight downshift -at T = 0.25, all of them exhibit a significant residual entropy, see inset of Fig. 1. There are, however, considerable differences of detail (leaving aside case (iii) on account of the unbounded classical entropy). The single S = 1/2 tetrahedron, (ii), releases its entropy more swiftly at low-T than our NLCE results on account of its singlet gap coupled with a small residual entropy S(0) = 1 4 ln 2. The spin ice experiment, (iv), with a slightly higher residual entropy of around S p = 1 2 ln 3 2 , in fact releases its entropy even more swiftly, with the peak in the specific heat peak being the most narrow on the low-T side.
The NaCaNi 2 F 7 experiment, (v), shows an initial high-T release of the entropy remarkably close to that of our S = 1/2 results. However, already above the peak in C V , the release in NaCaNi 2 F 7 is comparatively considerably greater, meaning that the spectral weight downshift in our results is stronger.
Indeed, the breadth of the peak in C V we find for the S = 1/2 Heisenberg model is broader not only than all the cases (ii-v), but also than the other paradigmatic highly frustrated S = 1/2 Heisenberg model, that on the kagome lattice. By comparison with the high order series expansion results in Ref. 91, we find that in the isotropic S = 1/2 Heisenberg antiferromagnet on the kagome lattice, the residual entropy of 0.475k B ln 2 is reached already at a higher temperature of T ≈ 0.30, whereas in the pyrochlore lattice our results suggest that this entropy is retained at a lower temperature of T ≈ 0.254, corresponding to the larger spectral downshift in pyrochlore.
It should be emphasized that this is not at all what would obviously have been expected. Generally, low dimensionality is considered to favour spectral weight downshift, as encoded e.g. by the Mermin-Wagner theorem. Also, in the Ising setting, triangular motifs are considerably more frustrated -S kagome = 1 24π 2 2π 0 dxdy ln [21 − 4 (cos x + cos y + cos(x + y))] ≈ 0.50183 92 for the kagome Ising magnet, much larger than in the pyrochlore case S Pauling ≈ 1 2 ln 3 2 ≈ 0.2027 4,93 . This implies that there is huge scope for unusual behaviour of this model at low-T . Alas, our results provide little indication of the detailed nature of the low-energy space of states. Indeed, many proposals for the behaviour of this magnet have been made, and it is hard to choose between them based on presently available information, as there is not even compelling evidence in favour of a particular physical picture. The concurrent lack of a pristine experimental realisation goes a long way towards explaining the divergence of theoretical predictions [15][16][17]23,84,94 so that different methods arguably come up with the conclusion most suited to them.
We are therefore left with the twin higher-level insights, namely that the pyrochlore S = 1/2 Heisenberg magnet is at least as frustrated, and arguably interesting, as the one on the kagome lattice; and that it is at least as intractable. We hope that future work will be able to build on the advances reported in this work. And, of course, that the low-T regime will become accessible experimentally in a suitable magnetic material.
Note added: As we were concluding this work, a preprint 95 appeared which also studied the thermodynamic properties of the pyrochlore S = 1/2 Heisenberg model using a combination of methods including canonical typicality, high temperature series expansion and the entropy method. It placed particular emphasis on extrapolation schemes in order to access the low-T regime. These constraints should respect the underlying physical properties of the system as much as possible to get a rapidly converging series expansion. We show a comparison of three different NLCE expansions for the pyrochlore lattice in Figs. 10 and 11: the single site expansion, the unit cell expansion (i.e. clusters on the fcc lattice, decorated by the tetrahedral unit cell) and the tetrahedral expansion used in the main text.
All expansions converge to the same curve at high temperatures, however the highest order single site expansion does not reach temperatures low enough to resolve the maximum of the specific heat, even after the Euler transform is applied. The unit cell expansion includes in principle much larger clusters (up to 6 unit cells correspond to 24 site clusters), however it is barely possible to obtained converged results for lower temperatures than in the single site expansion, while the tetrahedral expansion yields convergence to much lower temperatures and a clear resolution of the specific heat capacity maximum. There are several reasons for this superior behavior: First, the central motif of the pyrochlore lattice is the tetrahedron. It is crucial to avoid dangling bonds and incomplete tetrahedra, which are present in both the unit cell and single site expansion. Secondly, due to the construction of unit cell clusters, the clusters used here do not include similarly large loops as in the tetrahedral expansion, which are crucial at low temperatures. Thirdly, and this is purely technical, due to unfavorable symmetry properties of the n = 6 unit cell clusters, we were unable to solve all 100 topologically distinct clusters with full diagonalization, since the largest remaining symmetry blocks were too large. Finally, the Euler transform for the unit cell expansion can only rely on 5 complete expansion orders, which is far less than in the case of the single site and tetrahedral expansion.
Therefore, we conclude that the tetrahedral expansion is far superior to any other approach in the pyrochlore lattice and allows to access the lowest temperatures.  of sites, unit cells, or tetrahedra respectively. The numbers listed correspond to the total number of clusters |C n |, the number of clusters which are not identical under application of non-translational symmetries |S n | and the number of topologically distinct clusters |T n |, which is computationally relevant since this is the number of clusters for which the Hamiltonian has to be diagonalized. The growth of the number of clusters with order n is depicted in Fig. 12.
It should be noted that it is computationally challenging to check the topological equivalence of two clusters, since their interaction graphs need to be checked for an isomorphism, which is an NP hard problem. Therefore, e.g. the reduction of 184140685 symmetrically distinct clusters at n = 14 in the unit cell clusters to 8002 topologically distinct clusters is already difficult and limits severely the access to higher orders.
For the tetrahedral expansion, there is only one topologically distinct cluster for the orders 1 to 3 and two clusters with four tetrahedra. At the highest order we could reach, there are 24 clusters composed of eight tetrahedra, which have either 24 or 25 spins, just at the limit of what can be solved with full exact diagonalization using all symmetries of the clusters. We show all topologically distinct clusters included in the NLCE up to 8 tetrahedra in Figs. 14 and 15.

Appendix C: Finite size clusters
In the present study, we consider two finite size clusters, the first being standard 32 site cluster consisting of two unit cells in direction a 1 , a 2 , a 3 studied e.g. in Refs. 81,95 , which we depicted in the inset of Fig. 3. The second cluster we study is the 48 site cluster shown in Fig. 13 FIG. 13. 48 site cluster used in our finite temperature DMRG simulations. We use periodic boundary conditions (periodic bonds not shown for clarity).