Adapting Planck's route to investigate the thermodynamics of the spin-half pyrochlore Heisenberg antiferromagnet

The spin-half pyrochlore Heisenberg antiferromagnet (PHAF) is one of the most challenging problems in the field of highly frustrated quantum magnetism. Stimulated by the seminal paper of M.~Planck [M.~Planck, Verhandl. Dtsch. phys. Ges. {\bf 2}, 202-204 (1900)] we calculate thermodynamic properties of this model by interpolating between the low- and high-temperature behavior. For that we follow ideas developed in detail by B.~Bernu and G.~Misguich and use for the interpolation the entropy exploiting sum rules [the ``entropy method'' (EM)]. We complement the EM results for the specific heat, the entropy, and the susceptibility by corresponding results obtained by the finite-temperature Lanczos method (FTLM) for a finite lattice of $N=32$ sites as well as by the high-temperature expansion (HTE) data. We find that due to pronounced finite-size effects the FTLM data for $N=32$ are not representative for the infinite system below $T \approx 0.7$. A similar restriction to $T \gtrsim 0.7$ holds for the HTE designed for the infinite PHAF. By contrast, the EM provides reliable data for the whole temperature region for the infinite PHAF. We find evidence for a gapless spectrum leading to a power-law behavior of the specific heat at low $T$ and for a single maximum in $c(T)$ at $T\approx 0.25$. For the susceptibility $\chi(T)$ we find indications of a monotonous increase of $\chi$ upon decreasing of $T$ reaching $\chi_0 \approx 0.1$ at $T=0$. Moreover, the EM allows to estimate the ground-state energy to $e_0\approx -0.52$.


I. INTRODUCTION
A paradigmatic highly frustrated spin model is the pyrochlore Heisenberg antiferromagnet (PHAF). The pyrochlore lattice is built of corner-sharing tetrahedra, see Fig. 1, below. There are several compounds where the magnetic atoms reside on the sites of the pyrochlore lattice and the exchange interaction is antiferromagnetic, see, e.g., Refs. [1][2][3].
Already the classical PHAF (i.e., for spin S → ∞) exhibits interesting properties and its study is far from being trivial [4][5][6][7][8][9][10]. Thus, the ground-state manifold is highly degenerate, the model exhibits strong short-range correlations, but it does not exhibit any long-range order, and, because of the huge degeneracy of the ground state, the model is very susceptible to various perturbations.
The quantum spin S = 1/2 PHAF is even more complicated. Thus, so far no accurate values for the groundstate energy e 0 for this model are available. On the one hand, the S = 1/2 case opens the route to new quantum phases [11]. On the other hand, such powerful straightforward numerical tools like standard quantum Monte Carlo or molecular dynamics simulations are not applicable for the S = 1/2 PHAF. Moreover, several approximation methods developed for one-and two-dimensional quantum spin systems (e.g., density matrix renormalization group and tensor network methods) are very limited in three dimensions.
Theoretical studies of the quantum PHAF are mostly focused on ground-state properties, see, e.g., [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], whereas much less attention has been paid to its finitetemperature properties. One reason for that is the lack of methods to study thermodynamics of three-dimensional frustrated quantum spin systems. Among the few papers studying the thermodynamics of the S = 1/2 PHAF we mention bold diagrammatic Monte Carlo simulations (stochastic sampling of all skeleton Feynman diagrams) down to the temperature J/6 [26]. This paper reports data for the susceptibility χ(T ) but no data for the specific heat c(T ). We will refer to these data for χ(T ) in Sec. IV B. A comprehensive analysis of the spin-S J 1 − J 2 Heisenberg model by employing the pseudofermion functional renormalization group technique was presented in Ref. [11]. However, this paper does not contain data for χ(T ) and c(T ). Finally, we mention the hightemperature expansion study and the rotation-invariant Green's function study of the S = 1/2 PHAF [25,27]. In these recent papers [11,25,26] no evidence for a finitetemperature phase transition was found, i.e., the spinhalf PHAF is most likely a three-dimensional spin system without singularities in the specific heat and the susceptibility.
The goal of the present paper is to study the thermodynamics of the S = 1/2 PHAF for the whole temperature region focussing on the specific heat c(T ) and the static uniform susceptibility χ(T ), both being basic and easily accessible quantities in experimental studies of PHAF compounds. To deal with the above mentioned chal- lenges when studying the finite-temperature properties of the S = 1/2 PHAF, we follow M. Planck's ideas of his seminal paper in 1900 [28], see also Appendix A, and perform a sophisticated interpolation between the low-and high-temperature behavior of a thermodynamic potential, namely, the entropy s as a function of internal energy e. For that we exploit also sum rules valid for the specific heat as proposed by B. Bernu and G. Misguich [29,30], for details see Sec. III A. In what follows we call this approach the entropy method (EM). We complement our studies based on the EM by using the finite-temperature Lanczos method (FTLM) for a finite pyrochlore lattice of N = 32 sites and the high-temperature expansion (HTE) up to order 13.
In the present paper, we estimate the ground-state energy to e 0 ≈ −0.52 and find evidence for a gapless spectrum, i.e., for a power-law behavior of the specific heat at low temperatures, and for a single maximum in c(T ) at about 25% of the exchange coupling.
The manuscript is organized as follows. We begin with introducing the model (Sec. II) and the description of the exploited methods (Sec. III). We report our findings obtained by the FTLM (finite lattices) and by the HTE and EM (infinite lattice) in Sec. IV. We summarize our results in Sec. V.

A. Entropy method (EM)
In accordance with M. Planck's strategy to derive the energy distribution of the black-body radiation [28,40], the EM is an interpolation scheme that combines presumed knowledge on high-and low-temperature properties and, in addition, exploits sum rules for the specific heat c(T ) in a clever way. The EM as used in the present paper was introduced in 2001 by B. Bernu and G. Misguich [29]. The method has been afterwards used, modified, and extended in Refs. [30,[41][42][43][44]. Below we explain briefly this procedure for self consistency.
Within the framework of the EM, we use the microcanonical ensemble working with the entropy per site s as a function of the energy per site e, s(e), in the whole (finite) range of energies. The temperature T and the specific heat per site c are given by the formulas where the prime denotes the derivative with respect to e. These equations form a parametric representation of the dependence c(T ). Knowing the high-temperature series for c(T ) up to nth order, c(T ) = n i=2 d i β i + O(β n+1 ) (d 1 = 0), β = 1/T , we immediately get the series for s(e) around the maximal energy e ∞ = 0 up to the same order n, s(e)| e→e∞=0 → ln 2 + n i=2 a i e i , (3.2) where the coefficients a i are known functions of the coefficients d i , see Appendix A of Ref. [29]. The behavior of s(e) as e approaches the (minimal) ground-state energy e 0 (i.e., as the temperature approaches 0) is also supposed to be known. It is, if c(T ) vanishes as T α when T → 0 (gapless excitations) and if c(T ) vanishes as T −α exp(−∆/T ), α = 2, when T → 0 (gapped excitations). Therefore we proceed differently in the gapless case and in the gapped case. Here it is assumed that e 0 and α are known (gapless case) or e 0 is known and α = 2 (gapped case).
In the gapless case, we introduce the auxiliary function [30] G(e) = (s(e)) The prefactor A in the power-law decay of the specific heat c(T ) for T → 0, c(T ) → AT α , is given by In the gapped case, we introduce the auxiliary function From the technical point of view, before performing the integration in the right-hand side of Eq. (3.11) one may perform the partial fraction expansion of the integrand which is obviously a rational function. The excitation gap ∆ in the decay of the specific heat c(T ) for T → 0, c(T ) ∝ T −2 exp(−∆/T ), is given by .
Until now we considered the EM for zero magnetic field h = 0. Of course, for non-zero h the thermodynamic functions depend on h, i.e., the entropy is now s(e, h). The magnetization per site m and the uniform susceptibility per site χ are given by the formulas [44] where the last equation implies that h is infinitesimally small. Clearly, the HTE coefficients for the specific heat are also changed. Simple algebra yields . . . , n; (3.14) we use here the high-temperature series for the static uniform susceptibility χ(T ) = n i=1 c i β i + O(β n+1 ), β = 1/T . The expression (3.2) for the series of s is valid, however, the coefficients a i are now known functions of the coefficients d i , c i , and h. For the gapless case all reasonings in Eqs. (3.3), (3.5) to (3.8) hold with the only difference that the ground-state energy now is e 0 − χ 0 h 2 /2, where χ 0 ≡ χ(T = 0) is the ground-state susceptibility which is assumed to be known. The approximate entropy in Eq. (3.7) now also depends on h, i.e., s app (e, h). For the case of gapped magnetic excitations the ground-state energy remains unchanged, because χ 0 = 0, and therefore all the equations (3.4), (3.9) to (3.12) are valid. Again, the approximate entropy in Eq. (3.11) now also depends on h, i.e., s app (e, h).
In summary, knowing the high-temperature series of c(T ) and χ(T ) together with (i) the ground-state energy e 0 , the exponent α, and the ground-state susceptibility χ 0 for the gapless case or (ii) only the ground-state energy e 0 for the gapped case, we obtain c(T ) and χ(T ) at all temperatures. For that, we use s app (e, h) which yields the specific heat c(T ) by Eq. (3.1) and the susceptibility χ(T ) by Eqs. (3.13) and (3.1).
Based on previous experience with the EM [29,30,[41][42][43][44], we use the following strategy: We discard those Padé approximants in G app (e), Eqs. (3.6) and (3.10), which give unphysical solutions; the remaining ones are called "physical". Moreover, we focus on those input parameter sets for which interpolations based on different Padé approximants lead to data sets for c(T ) and χ(T ) being quite close to each other. For further details about the EM in the context of the S = 1/2 PHAF see Sec. IV B.

B. Finite-temperature Lanczos method (FTLM)
The FTLM is an efficient and very accurate approximation to calculate thermodynamic quantities of quantum spin systems on finite lattices of N sites at arbitrary temperatures. It is an unbiased numerical approach, where thermodynamic quantities such as the specific heat and the susceptibility are determined using trace estimators [45][46][47][48][49][50][51][52][53]. The key element is the approximation of the partition function Z using a Monte-Carlo like representation of Z, i.e., the sum over a complete set of 2 N basis vectors present in Z is replaced by a much smaller sum over R random vectors | ν for each subspace H(γ) of the Hilbert space, where except the conservation of total S z we also use the lattice symmetries of the Hamiltonian to decompose the full Hilbert space into mutually orthogonal subspaces labeled by γ. The exponential of the Hamiltonian is then approximated by its spectral representation in a Krylov space spanned by the N L Lanczos vectors starting from the respective random vector | ν . The FTLM representation of the partition function finally reads where | n(ν) is the nth eigenvector of H in the Krylov space with the corresponding energy ǫ (ν) n . To perform the symmetry-decomposed numerical Lanczos calculations we use J. Schulenburg's spinpack code [54,55].

C. High-temperature expansion (HTE)
The HTE is a universal approach to discuss the thermodynamics of spin systems [56]. In the present study we use the Magdeburg HTE code developed mainly by A. Lohmann [27,57,58] (which is freely available at http://www.uni-magdeburg.de/jschulen/HTE/) in an extended version up to 13th order, see Appendix B. With this tool, we compute the series of the specific heat To extend the region of validity of the "raw" HTE series we may use Padé approximants [m, n] = P m (β)/Q n (β), where P m (β) and Q n (β) are polynomials in β of order m and n, respectively. The coefficients of the polynomials P m (β) and Q n (β) are determined by the condition that the expansion of [m, n] has to agree with the initial power series up to order O(β m+n ).

A. Finite lattices
In Ref. [24] finite lattices of N = 28 and N = 36 sites are used to discuss ground-state properties. These lattices are built by stacked alternating triangular and kagome layers imposing periodic boundary conditions within the layers, but open boundary conditions perpendicular to them. We have calculated the HTE series for these finite lattices. We compare these finite-lattice series with the corresponding HTE series of the infinite pyrochlore lattice to judge the finite lattices. We found that for c(T ) all HTE coefficients are different. For the susceptibility χ(T ) only the lowest-oder term coincides, i.e., the agreement is only marginally better. This drastic difference between the finite lattices and the infinite pyrochlore lattice can be attributed to the edge spins stemming from the imposed open boundary conditions. Thus, we conclude that the finite lattices of N = 28 and N = 36 used in Ref. [24] are not appropriate to discuss the thermodynamics of the PHAF. However, we note that they can be useful to discuss the ground-state properties, e.g., spin-spin correlations when considering spins away from the edge spins.
A more suitable finite lattice is the one with N = 32 sites imposing periodic boundary conditions in all directions. This lattice contains eight face-centered-cubic cells, i.e., the edge vectors go along the face-centeredcubic basis vectors and have twice the length of these. For this lattice the HTE series for c (χ) coincides up to 3rd (4th) order with that of the infinite lattice.
The FTLM is the adequate approach to study the finite S = 1/2 PHAF of N = 32 sites. In Fig. 2 we show data for c(T ) (top) and χ(T ) (bottom) over a wide temperature range using a logarithmic T scale. The specific heat exhibits the typical main maximum at T = 0.53 and, in addition, two low-T maxima at T = 0.012 and at T = 0.117. While the maximum at T = 0.012 is certainly a finite-size effect, one can speculate that the other low-T maximum at T = 0.117 signals an extra lowenergy scale set by low-lying singlets (see the density of states shown in the inset of the middle panel of Fig. 3) that might be also relevant for the infinite system. Such a feature has been observed in low-dimensional highly frustrated quantum magnets, e.g., the spin-half kagome Heisenberg antiferromagnet (HAF), where the existence of such an extra low-T peak is a subject of a long-standing and ongoing debate [30,52,[59][60][61][62][63][64]. However, in the three-dimensional PHAF the finite-size effects are undoubtedly stronger than in the two-dimensional kagome HAF. Thus, to conclude a double-peak structure in c(T ) from our FTLM data is inappropriate. For the static uniform susceptibility χ(T ) the low-lying singlets are not relevant and χ(T ) does not show extra-peaks except the well-pronounced maximum that is typical for finite spin systems with χ 0 ≡ χ(T = 0) = 0. Again, this behavior might be not representative for the infinite system, particularly, in case that χ 0 > 0 for N → ∞. To quantify the temperature region where the finite N = 32 lattice may be representative for the infinite lattice we compare in Figs. 4 and 5 several Padé approximants of the HTE series of the finite and the infinite lattices. Obviously, the data for N = 32 and N = ∞ coincide only down to T ∼ 0.7, thus, indicating that finite  pyrochlore lattices accessible by FTLM are not suitable to discuss the thermodynamics of the spin-half PHAF below this temperature.
Nevertheless, the finite-size data for the PHAF are useful to demonstrate frustration effects. For that we may compare the PHAF with the S = 1/2 HAF on the simple-cubic lattice, because both have six nearestneighbors, i.e., a simple mean-field decoupling of the Heisenberg Hamiltonian would yield identical thermodynamics. However, the simple-cubic HAF exhibits a finite-temperature phase transition to Néel order, but the PHAF does not order. Thus, we compare both models on finite lattices of N = 32 sites, where the simplecubic finite-temperature phase transition is irrelevant, see Fig. 3, where we compare FTLM data of the specific heat, the entropy, and the susceptibility using a linear T scale for N = 32. The tremendous influence of frustration is visible at all temperature scales. In particular, the spectrum at low energies in the frustrated system is much denser than that of the unfrustrated one [see the density of states (histogram, ∆E = 0.02) in the inset in the middle panel of Fig. 3], thus leading to the drastic differences at low T , see the upper and middle panels of Fig. 3. Remarkably, the noticeable differences in all quantities are present at pretty high temperatures. Only, beyond T 3 the corresponding curves approach each other. A striking effect of frustration is also the shift of the maximum in χ(T ) to lower temperatures, see the lower panel of Fig. 3.

B. Infinite lattice
Let us now move to a detailed investigation of the infinite PHAF by using the EM interpolation scheme, see Sec. III A. Using our Magdeburg HTE code [27,57,58] we have created the HTE series for c(T ) and χ(T ) up to order 13 (see Appendix B) that provides the hightemperature input for the EM. Since the EM finally uses Padé approximants of a power series of s(e) derived from the initial HTE series, the HTE input determines the highest order of the Padé approximants of the EM interpolation scheme.
As a low-temperature input for the EM we need the ground-state energy e 0 . There is a large variety of reported values for e 0 of the spin-half PHAF ranging from e 0 = −0.57 to e 0 = −0.45, namely, e 0 = −0.572 [65], −0.56 [15], −0.49 [12,16], −0.482081 [24], −0.466971 [24], −0.459 [22], −0.457804 [13], −0.45093 [25], −0.4473 [23], i.e., accurate values for e 0 are missing. We also need the low-temperature law for c(T ), where we have to distinguish between gapped and gapless behavior (cf. Sec. III A). Finally, in case of gapless excitations the exponent α of the power law is an input parameter, and for the susceptibility χ 0 ≡ χ(T = 0) is required as input. Fortunately, for gapped excitations the value of the gap is not needed as an input, it is rather an output of the EM. Moreover, we have χ 0 ≡ χ(T = 0) = 0 in this case.
We begin with the specific heat c(T ) for which the EM is better justified (two sum rules are exploited). Actually, the low-temperature law for c(T ) is not known for the S = 1/2 PHAF. Advantageously, the previous ample of experience with the EM [29,30,[41][42][43][44] provides valuable hints to overcome this difficulty. Thus, in case that these input data are too far from the true (but possibly unknown) data one gets inconsistent or unphysical results. Hence, to get physical (i.e., pole free) Padé approximants requires reasonable values for e 0 and reasonable assumptions on the low-T behavior of c(T ). Moreover, getting a large number n P of similar Padé approximants for a certain input data set indicates the physical relevance of this set. On the other hand, the appearance of significant differences between the various Padé approximants is a criterion to discard the corresponding input set. The successfulness of this strategy has been demonstrated for several examples where excellent reference data are available, in particular, for the S = 1/2 XY , HAF, and Ising chains [29,42], for the S = 1 HAF (Haldane) chain [29], for the S = 1/2 square-lattice and  triangular-lattice Heisenberg ferro-and antiferromagnets [29] or for the S = 1/2 kagome-lattice HAFs [30,42,44].
By combining different assumptions on the groundstate energy and low-T behavior of c(T ) of the PHAF we have generated a large set of temperature profiles for the specific heat. In the next step, we use the guidelines described above to evaluate the used input data, this way obtaining definite conclusions on their relevance. In sum, a crucial criterion is that a particular input data set leads to a close bundle of temperature profiles obtained by many physical Padé approximants. Since, the high-temperature part is per se identical this criterion concerns the temperature region T 0.5. In what follows, we will present in the main text only data for the most relevant input sets, whereas presentation of some other illustrative results, only briefly mentioned in the main text, are transferred to Appendix C.
As a first result, we found that the assumption of gapped excitations is not favorable for the following rea-  Fig. 16 in Appendix C), and there is a decent number n P of Padé approximants yielding similar c(T ) profiles, see Fig. 7. For e 0 less than −0.521 the number of physical Padé approximants noticeably decreases. Although the results for c(T ) do not entirely discard a gapped ground state, further EM analysis of χ(T ) under this assumption leads to disagreement with diagrammatic Monte Carlo simulations of Ref. [26] at T < 0.7, see Fig. 13 below. We may consider these findings for c(T ) and χ(T ) as an indication to favor a gapless ground state. We mention that the Green's function results indicate gapless magnetic excitations [25], and, as mentioned already above, the data for χ(T ) given in [26] down to T = J/6 seem also to be in favor of gapless excitations.
We focus now on the gapless case with a power-law decay of the specific heat c(T ) = AT α as T → 0. Since we do not know the exponent α, we study different values α = 1, 3/2, 2, 5/2, 3. As for the gapped case, we again varied e 0 in the region −0.57 . . . − 0.45. We observed that only in a much smaller region around e 0 = −0.52 reasonable results can be obtained (see also the discussion of Fig. 8, below). Thus, in what follows we consider preferably the region −0.53 . . . − 0.50 in more detail.
First we consider the prefactor A app that is given within the EM by Eq. (3.8). According to above outlined criteria for physically relevant EM outcomes the values of A app obtained by different Padé approximants must be very close to each other. From Fig. 8 (see also   [12,13,15,16,[22][23][24][25]65]. We consider the coincidence of A values for different Padé approximants as a necessary criterion to figure out regions of relevant values for e 0 and α. Since A (together with α) determines the c(T ) profile at sufficiently low T , we can get additional reliability by examining the region around the main ("high"-temperature) maximum. For that we compare the position T max and the height c(T max ) of this maximum obtained from different Padé approximants in dependence on e 0 within the regions guided by the previous inspection of A for α = 3/2, 2, 5/2 in Fig. 9 (Fig. 18 in Appendix C reports such data for a wider region of e 0 including also α = 1 and 3). For each α we find pretty small regions of e 0 which yield almost identical T max and c(T max ), and, consistently, this region fits well to the region obtained by inspection of A. For example, for α = 2 all Padé approximants give almost identical T max and c(T max ), cf. Fig. 9 (red symbols), if e 0 is taken within the region −0.522 . . . − 0.519, which is in excellent agreement with that obtained from the prefactor A, see above. Note, however, that for α = 1 and 3 the diversity of T max and c(T max ) is noticeably larger than for α = 3/2, 2, 5/2, cf. Fig. 18 in Appendix C, indicating that the exponents α = 1 and 3 are less favorable.
Finally, after the specification of the ground-state energy values as outlined above, we present in Fig. 10 the full c(T ) curves obtained by the EM for α = 3/2, 2, 5/2 and a few related optimal values of e 0 . [Corresponding curves for α = 1 and α = 3 are shown in Fig. 19 in Appendix C. Moreover, Fig. 20 in Appendix C provides additional results of c(T ) comparing data for all α = 1, 3/2, 2, 5/2, 3 for various values of e 0 .] As can be seen in Fig. 10, there is a quite large number of Padé approximants (see the parameter n P in brackets behind energy values) yielding very similar temperature profiles c(T ). Thus, for α = 2 we show in the middle panel of Fig. 10  Finally, we compare our EM results for c(T ) with data calculated by HTE (without subsequent EM interpolation) and by the Green's function approach [25], cf. Fig. 11. The Green's function results deviate from the EM results already below T ∼ 1, whereas the HTE data deviate below T ∼ 0.7.
The EM straightforwardly also provides the temperature profile of the entropy s(T ), see Fig. 12. Since the finite-temperature phase transition present in the simplecubic HAF does not influence s(T ) as much as c(T ), we compare the data for the PHAF with corresponding ones for the simple-cubic HAF taken from Ref. [66]. We also show HTE data. Similar as already found for the finite system, cf. the middle panel of Fig. 3, the frustration leads to a much faster increase of s at low temperatures for the PHAF. Thus, at T = 0.5 the entropy already amounts to more than 50% of its maximal value ln 2 ≈ 0.69. (Note that in Fig. 22 in Appendix C we present a direct comparison of the gapped and gapless temperature profiles of s.) We consider now the static uniform susceptibility χ calculated by the EM as described in Sec. III A. According to the results of the rotation-invariant Green's function method, there is χ 0 ≡ χ(T = 0) > 0. A non-zero χ 0 may be also expected from the diagrammatic Monte Carlo simulations [26]. Nevertheless, we will not exclude from the beginning χ 0 = 0, i.e., a non-zero spin gap. Although the above discussed EM data for c(T ) are in fa- Figure 11. Specific heat c(T ) of the S = 1/2 PHAF: Comparison of our EM data (e0 = −0.522, α = 5/2; e0 = −0.520, α = 2; e0 = −0.517, α = 3/2) with results obtained by the rotation-invariant Green's function method [25] (dashed blue) and by the HTE (thin solid curves; we show the same HTE data for N = ∞ as shown in Fig. 4). Figure 12. EM data for the entropy for the gapless case with e0 = −0.522, α = 5/2; e0 = −0.520, α = 2; e0 = −0.517, α = 3/2. We also show quantum Monte-Carlo data for the simple-cubic HAF [66] as well as HTE data (thin solid lines), where we show the same Padé approximants as used in Fig. 4 for c(T ). (Note the the HTE data are very close to the EM data.) vor of a gapless spectrum, the specific heat profiles were not fully conclusive to entirely reject the gapped spectrum. Moreover, one could have gapless singlet (i.e., nonmagnetic) excitations but gapped triplet (i.e., magnetic) excitations. Thus we studied the case with exponential decay of χ(T ) and c(T ) as T → 0 (i.e., gapped singlet and triplet excitations), see Fig. 13, as well as the case with exponential decay of χ(T ) and power-law decay of c(T ) as T → 0, (i.e., gapless singlet and gapped triplet excitations), see Fig. 14, where we consider those values for e 0 and α which are previously used to get c(T ) (see also Figs. 23 and 24 in Appendix C). From these figures, it is obvious that the susceptibility profiles based on gapped magnetic excitations are not compatible with the diagrammatic Monte Carlo data of Ref. [26] in the temperature region below 0.7, thus providing further evidence against a gapped spectrum. On the other hand, the EM results for gapless excitations with α = 3/2, 2, 5/2 and e 0 = −0.522 . . . − 0.516 and with non-zero χ 0 as shown in Fig. 14 fit much better to the data of Ref. [26].
As for the specific heat, we finally compare our EM results for χ(T ) with data calculated by HTE (without subsequent EM interpolation), by the Green's function approach [25], cf. Fig. 15, where we also include the data of the diagrammatic Monte Carlo simulations [26]. The Green's function results deviate already below T ∼ 1.5, and the HTE coincides down to T ∼ 0.7, whereas χ(T ) profiles with χ 0 = 0.1 are in excellent agreement with the data of Ref. [26].

V. CONCLUSIONS AND SUMMARY
We have studied the specific heat c(T ), the entropy s(T ), and the static uniform susceptibility χ(T ) of the spin-half pyrochlore Heisenberg antiferromagnet (PHAF) on a finite lattice of N = 32 sites using the finitetemperature Lanczos method (FTLM) and on the infinite lattice using the high-temperature expansion (HTE) up to order 13 and a sophisticated interpolation between the low-and high-temperature behavior of the thermodynamic potential entropy s as a function of internal energy e [the entropy method (EM)].
We found that finite lattices of such size are not appropriate to get reasonable results below T ∼ 0.7, but they might be useful to get a general impression on the strong frustration effects present in the PHAF by comparison of the HAF on finite pyrochlore and simple-cubic lattices. A similar limitation to temperatures above T ∼ 0.7 is valid for the HTE even if the range of validity of the hightemperature series is extended by Padé approximants. Only the EM interpolation is suitable to overcome these limitation and to provide sound data for the whole temperature range.
Our main findings for the specific heat c(T ) are as follows. (i) Contrary to the two-dimensional kagome HAF, we do not find hints neither for an extra low-T peak nor an extra shoulder below the main maximum. However, the absence of an extra low-T feature goes hand in hand with a significant shift of the single maximum towards Figure 15. Susceptibility χ(T ) of the S = 1/2 PHAF: Comparison of our EM data (e0 = −0.522, α = 5/2; e0 = −0.520, α = 2; e0 = −0.517, α = 3/2) with results obtained by the rotation-invariant Green's function method [25] (dashed blue), by diagrammatic Monte Carlo [26], and by the HTE (thin solid curves; we show the same HTE data for N = ∞ as shown in Fig. 5).
T ∼ 0.25, which is much lower than for the kagome HAF, where the main maximum is at T max /J = 0.67 [44,52]. This conclusion is robust, i.e., it is obtained not only for gapless excitations for all reasonable exponents α of the low-temperature power law of c(T ), but holds also for gapped excitations. (ii) A gapless spectrum is more favorable than a gapped one, i.e., most likely there is power-law low-T behavior of c(T ). Although best results are for an exponent α = 2, other exponents (α = 1, 3/2, 5/2, 3) cannot be excluded. (iii) We predict a ground-state energy e 0 ≈ −0.52 [67].
Our EM data for the susceptibility χ(T ) in comparison with data obtained by diagrammatic Monte Carlo [26] provide further evidence for a gapless spectrum with a ground-state energy e 0 ≈ −0.52. The temperature profile of χ most likely does not show a maximum, rather there is a monotonous increase of χ upon decreasing of T reaching a zero-temperature value of χ 0 ≈ 0.1.
Recently, we have learned that R. Schäfer et al. [68] also examined thermodynamics of the quantum PHAF, however, using for that numerical linked-cluster expansions. In his revolutionary paper in 1900 [28] (see also the Nobel Prize address [Notes 7, 12, and 13 at the end of the Nobel Prize address "The origin and development of the quantum theory" by Max Planck delivered before the Royal Swedish Academy of Sciences at Stockholm, 2 June, 1920] [40]), M. Planck investigated the energy distribution of electromagnetic radiation emitted by a black body in thermal equilibrium. For that he considered the entropy of the equilibrium radiation S [69] in relation with its energy U , or more accurately, the second derivative d 2 S/dU 2 . According to Wien's law it is which interpolates between both limiting cases. By integrating we get which yields Planck's formula These arguments can be formulated within the setup of the entropy method to find the specific heat c(T ) of a (Bose) oscillator, which represents the electromagnetic radiation with the frequency ν = 2πω. Now we know that b = ω and c = 1. Taking into account the zeropoint energy we have to replace U → e = U + ω/2. Then the Planck's interpolation formula (A.3) for the auxiliary function G(e) = (s(e)) ′′ reads: , e ≪ ω,  In this appendix we collect more figures presenting EM data for the specific heat, the entropy, and the susceptibility which are briefly discussed but not shown as figures in Sec. IV B.
In Fig. 16 we show some results related to Fig. 6 using a smaller region of e 0 .
In Fig. 17 we show some results related to Fig. 8 using a smaller region of e 0 . The region of e 0 in which various Padé approximants yield almost the same prefactor A app (3.8) is different for α = 3/2, 2, 5/2. Clearly, the values of α and A are linked. For example, after assuming α = 2 for e 0 = −0.520 . . . − 0.519 the EM prediction for the specific heat as T → 0 reads: c(T ) = A app T 2 with A app = 29 . . . 31.
In Fig. 18 we show some results related to Fig. 9 using a wider region of e 0 and including exponents α = 1 and α = 3. The position and the height of the maximum of the specific heat, as they follow from raw HTE series ex-  Fig. 20 provides temperature profiles of the specific heat which complement those shown in Fig. 10 and Fig. 19: We show c(T ) for e 0 = −0.521, e 0 = −0.519, e 0 = −0.517, and e 0 = −0.515 (from bottom to top) and compare data for α = 1 (cyan), α = 3/2 (green), α = 2 (red), α = 5/2 (magenta), and α = 3 (blue). The shown temperature profiles allow to estimate how close to each other are the various EM data based on different Padé approximants.
In Fig. 21 we collect the best (i.e., for such a value of e 0 which gives the largest number of almost coinciding resulting curves) EM predictions for c(T ) for the gapped and the gapless spectrum. The best EM results for the temperature dependence of the entropy are shown in Fig. 22. Note, all curves for each color are indistinguishable in this figure.  In Fig. 24 we collect the best (i.e., for such a value of e 0 which gives the largest number of almost coinciding resulting curves) EM predictions for χ(T ) under the gapped assumption and the gapless assumption with several values of α. Note, χ(T ) as it follows from the assumption about gapless singlet excitations but χ 0 = 0 (all colors except light brown) deviates stronger from the diagrammatic Monte Carlo result than χ(T ) as it follows from the assumption about gapped excitations (light brown). However, the agreement of different seven χ(T ) curves for the latter case is rather poor.