The case for a U(1)_π Quantum Spin Liquid Ground State
in the Dipole-Octupole Pyrochlore Ce₂Zr₂O₇

The rare-earth pyrochlore oxides, R₂B₂O₇, where $R^{3+}$ is a trivalent rare-earth ion and $B^{4+}$ is a non-magnetic tetravalent transition-metal ion, display a wealth of both exotic and conventional magnetic ground states. Their $R^{3+}$ ions decorate a network of corner-sharing tetrahedra, one of the archetypes for geometrical frustration in three dimensions, and this crystalline architecture underlies many of their exotic properties [1]. A separation of energy scales, with crystal electric field (CEF) effects dominating over exchange interactions, often results in a well-separated CEF ground state doublet for the $R^{3+}$ ion, and interacting pseudospin-$\frac{1}{2}$ degrees of freedom at low temperatures [2, 3, 4].

It is well-appreciated that the CEF Hamiltonian determines both the size of the magnetic moment at the $R^{3+}$ site and its anisotropy, but less well-appreciated that the symmetry of the CEF ground-state can imprint itself on the exchange Hamiltonian [4, 5, 6]. The possible symmetries of the ground state doublets then lead to an important classification of the rare earth pyrochlores, which depends on how their CEF doublet transforms under time-reversal symmetry and the point group symmetry of the $R^{3+}$ site. Three classes of doublets arise, one for non-Kramers ions with an even number of electrons, and two for Kramers ions with an odd number of electrons. The non-Kramers case gives rise to a pseudospin wherein one component of the pseudospin transforms as a magnetic dipole and two transform as quadrupoles. For Kramers ions, we have the familiar case where all three components of the pseudospin in the ground state doublet transform as magnetic dipoles, as well as the more exotic one where two components transform as magnetic dipoles and one transforms as an octupole. This latter case is known to describe the CEF Kramer’s ground state of 4f¹ Ce³⁺ in Ce₂Zr₂O₇ [7, 8], a dipole-octupole (DO) ground state doublet, and also that of its sister pyrochlore, Ce₂Sn₂O₇ [9]. Fig. 1(a) pictorially displays the magnetic charge distributions associated with both magnetic dipoles and octupoles decorating the tetrahedra on part of a cubic pyrochlore lattice. As discussed above, for the dipole-octupole doublets relevant to Ce₂Zr₂O₇, a single component of the pseudospin-$\frac{1}{2}$ degree of freedom (the $y$-component) behaves as an octupole, while the $x$ and $z$ components behave as dipoles under the symmetry of the lattice and time-reversal symmetry, as schematically illustrated in Fig. 1(b).

Such DO doublets decorating pyrochlore lattices are theoretically known to allow for at least 6 distinct quantum disordered and ordered ground states, with three in each of the dipole and octupole sectors [10, 11, 12]. Recent neutron scattering measurements on single crystal Ce₂Zr₂O₇ have uncovered a signal that strongly resembles predictions for the energy-integration of emergent photon excitations in a U(1) quantum spin ice [7], while recent experiments on powder samples of Ce₂Sn₂O₇ have been interpreted in terms of a U(1) quantum spin ice ground state in the octupole sector [9].

In this paper, we present new polarized neutron diffraction and heat capacity measurements on single crystal Ce₂Zr₂O₇. The former bears both similarities and differences from that measured in the canonical dipolar spin ice compound, Ho₂Ti₂O₇, while the latter shows no sign of a thermodynamic phase transition above $T$ = 0.06 K. $C_{p}$ rises sharply at low temperatures, initially plateauing near 0.08 K, before falling off towards a high temperature zero beyond 3 K, consistent with previous measurements [8]. We have modelled the high temperature $C_{p}$, and the powder-averaged magnetic susceptibility using quantum numerical linked cluster (NLC) expansions. This allows us to estimate and constrain the parameters of the anticipated near-neighbour XYZ Hamiltonian. To the extent that interactions beyond near-neighbour do not alter ground state selection, we constrain the nature of the ground state itself, with the results indicating a U(1)_π QSL ground state is selected at low temperature.

We use the resulting near neighbour exchange parameters to calculate the equal-time spin flip (SF) and non-spin flip (NSF) structure factors in the $[HHL]$ scattering plane. This calculation resembles the new polarized neutron diffraction measurements in the SF channel from single crystal Ce₂Zr₂O₇, but cannot account for the observed zone-boundary diffuse scattering in the NSF channel. We attribute this discrepancy to interactions beyond near-neighbour in the Hamiltonian, which are expected to be small, and a full study of which is beyond the scope of our present work. The same discrepancy exists for spin-polarized neutron diffraction from Ho₂Ti₂O₇, where it was ascribed to expected long range dipolar interactions [13]. NLC calculations using the same near-neighbour exchange Hamiltonian were also carried out to 7^th order. While these agree with the 4^th order calculations above $\sim$0.5 K, they depart from the measured $C_{p}$ at lower temperatures. We interpret this as arising from the same interactions beyond near neighbour in Ce₂Zr₂O₇ that were revealed by the NSF zone boundary scattering. As these are relatively weak, they only manifest themselves at low temperatures.

A further consistency check is carried out via semiclassical Monte Carlo and Molecular Spin Dynamics using the best fit near-neighbour Hamiltonian. This calculation accounts for the energy dependence of the inelastic spectral weight making up the diffuse scattering at low temperatures without further adjustment of the NLC-determined near-neighbour Hamiltonian. We further show that the full $R\ln(2)$ entropy of the DO ground state doublet can be accounted for to 10 K with a smooth extrapolation of $C_{p}$ from the lowest temperature data point at $T$ = 0.06 K, to zero at $T$ = 0 K, using a theoretical form which is simultaneously consistent with both the expected behavior of a U(1) QSL at low temperature, and the high temperature limit of the NLC calculations. Interestingly, the Pauling, classical spin ice entropy $R\ln(2)$ less $\frac{R}{2}\ln(\frac{3}{2}$) is recovered from the peak in the $C_{p}$ data at $\sim$0.08 K, to 10 K.

We have carried out new polarized diffraction measurements on single crystal Ce₂Zr₂O₇ using the D7 diffractometer at the Institute Laue Langevin. This diffractometer employs a spin polarized monochromatic incident beam, which was $E_{i}$ = 3.47 meV for this experiment. This configuration effectively integrates over $\Delta E$ $\sim$0.16 meV during the course of a diffraction measurement. A single polarization direction, perpendicular to the $[HHL]$ scattering plane, was employed, and as such the spin flip (SF) and non-spin flip (NSF) diffuse scattering profiles can be independently measured. The diffuse scattering associated with these two cross sections, SF and NSF, are shown in the $[HHL]$ scattering plane for Ce₂Zr₂O₇ in Fig. 2(a) and 2(b), respectively for the temperature-difference data set $T$ = 0.045 K - $T$ = 10 K. For comparison, the corresponding SF and NSF diffuse scattering patterns as measured on single crystal Ho₂Ti₂O₇ at $T$ = 1.7 K are shown in Fig. 2(c) and 2(d), respectively [13]. These earlier spin polarized diffuse scattering measurements on Ho₂Ti₂O₇ (Ref. [13]) played a formative role in the development of classical spin ice physics, as they drew clear attention to “pinch point” scattering within the SF cross section at (0,0,2) and (1,1,1) and equivalent wavevectors, due to the presence of a classical Coulomb phase at low temperature. These measurements on Ho₂Ti₂O₇ also observed zone-boundary diffuse scattering in the NSF channel, which was later attributed to the long range dipolar interactions relevant to the large Ho³⁺ dipole moments.

The comparison between the spin polarized diffuse scattering from Ce₂Zr₂O₇ and Ho₂Ti₂O₇ in Fig. 2 is interesting both in what is similar and where the discrepancies between the two materials lie. One may note however, that the comparison is made at quite different temperatures, 0.045 K for Ce₂Zr₂O₇ but only 1.7 K for Ho₂Ti₂O₇. In fact, the large Ho³⁺ moments and effective ferromagnetic coupling cause Ho₂Ti₂O₇ to depolarize the beam at lower temperatures, whereas no such issue is present for Ce₂Zr₂O₇ due to its much smaller Ce³⁺ moments. Quasi-pinch point SF scattering is observed near (0,0,2) Bragg positions in Ce₂Zr₂O₇, but it is not as constricted as that observed at (0,0,2) in Ho₂Ti₂O₇, even though the earlier measurements in Ho₂Ti₂O₇ were taken at much higher temperature. Furthermore while diffuse SF scattering extends out in (1,1,1) and equivalent directions in a snow flake-like pattern in Ce₂Zr₂O₇, pinch points appear to be absent in these directions.

In contrast, and somewhat surprisingly, the observed NSF diffuse scattering in Ce₂Zr₂O₇ is quite similar to that seen in Ho₂Ti₂O₇. In both cases the diffuse scattering tends to follow the face-centred cubic Brillouin zone boundaries, outlined in grey in Fig. 2(b). In Ho₂Ti₂O₇, this was ascribed to interactions beyond near neighbour [13], which was not surprising, given that dipolar interactions are expected to dominate over exchange interactions even for near neighbours in Ho₂Ti₂O₇. However, the Ce³⁺ moments are $\sim$8 times smaller than those of Ho³⁺ and hence dipolar interactions are expected to be $\sim$64 times smaller in Ce₂Zr₂O₇. We revisit our new polarized neutron diffraction data in Section V where we compare the measured SF and NSF signals to NLC calculations using the near-neighbour exchange parameters yielded in this work.

The gold standard for determining the microscopic spin Hamiltonian of magnetic materials is inelastic neutron scattering studies of spin wave spectra. This technique can and has been successfully applied to pyrochlore magnets with pseudospin-$\frac{1}{2}$ degrees of freedom arising from well-separated ground state CEF doublets, including Yb₂Ti₂O₇ and Er₂Ti₂O₇ [14, 15, 16, 17, 18, 19]. For disordered ground states, it is necessary to perform measurements in a sufficiently strong magnetic field, so as to polarize the ground state, thus giving rise to well defined spin wave spectra. However, this is not always possible. For example, the classical spin ice ground state as appears in Ho₂Ti₂O₇, doesn’t allow transverse spin fluctuations - hence no well defined spin wave excitations are observed due to Ho³⁺’s non-Kramers CEF doublet eigenvectors [20]. No evidence for well defined spin waves has been observed to date in either zero or non-zero magnetic field in Ce₂Zr₂O₇, a likely consequence of the form of Ce³⁺’s DO CEF ground state doublet and spin Hamiltonian. Hence estimates for the microscopic spin Hamiltonian parameters for such materials can only come from sophisticated modelling of other data, such as the high temperature thermodynamic data presented here. We note that a related work has appeared coincident with this paper which performs independent modeling of heat capacity, magnetization and neutron scattering measurements in Ce₂Zr₂O₇, and reaches similar conclusions [21].

The combined analyses of the measured $C_{\mathrm{mag}}$ and $\chi$ give experimental estimates for the near neighbour exchange constants in Ce₂Zr₂O₇, yielding $\theta\sim 0$ and $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ = (0.064, 0.063, 0.011) meV or $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ = (0.063, 0.064, 0.011) meV. While neutron scattering measurements were not modelled in order to constrain the microscopic spin Hamiltonian for Ce₂Zr₂O₇, it is interesting and important to see to what extent the measured neutron scattering from Ce₂Zr₂O₇ is consistent with calculations using the near-neighbour spin Hamiltonian so derived.

To conclude, we report new spin polarized neutron diffraction and $C_{\mathrm{mag}}$ measurements on single crystal Ce₂Zr₂O₇ in zero magnetic field. Our modelling of $C_{\mathrm{mag}}$, $\chi$, and $S(\mathbf{Q},T)$ with NLC calculations provides strong constraints on the exchange terms in the microscopic near neighbour XYZ Hamiltonian. We arrive at best fit Hamiltonian parameters $\theta\sim 0$ and $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ = (0.064, 0.063, 0.011) meV or $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ = (0.063, 0.064, 0.011) meV, which indicates that a U(1)_π QSL ground state is selected near the boundary between dipolar and octupolar character.

The best-fitting exchange parameters from this work largely describe the SF neutron diffraction signal measured from single crystal Ce₂Zr₂O₇, while zone boundary scattering in the NSF channel indicates the significance of interactions beyond near-neighbour, including long-ranged dipolar interactions. The 7^th order NLC calculations for $C_{\mathrm{mag}}$ evaluated at the best fit Hamiltonian parameters do not describe the measured $C_{\mathrm{mag}}$ at the lowest temperatures, again consistent with weak interactions in Ce₂Zr₂O₇’s Hamiltonian beyond near-neighbour and beyond the scope of the present calculations.

The new $C_{\mathrm{mag}}$ data extends to temperatures as low as $T$ = 0.058 K and can be smoothly extrapolated to zero temperature using a form consistent with both the XYZ spin Hamiltonian estimated from fitting the NLC calculations to the data, and with a $T^{3}$ form for $C_{\mathrm{mag}}$ at sufficiently low temperatures, appropriate to emergent gapless photons. With such a low $T$ form for $C_{\mathrm{mag}}$ in place we show the $R\ln(2)$ entropy associated with Ce³⁺’s DO doublet ground state is recovered to 10 K. Phenomenologically, we observe that the Pauling entropy for spin ice is recovered above the onset of the $T$ $\sim$0.08 K plateau in $C_{\mathrm{mag}}$.

The powder and single crystal samples of Ce₂Zr₂O₇ used in this work were prepared and characterized as described in Ref. [7]. La₂Zr₂O₇ was synthesized in order to estimate the phonon contribution to the $C_{p}$ of Ce₂Zr₂O₇. The powder samples of La₂Zr₂O₇ measured in this work were first prepared by mixing stoichiometric amounts of La₂O₃ (Alfa Aesar 99.99%) and ZrO₂ (Alfa Aesar 99.7%). The La₂O₃ (ZrO₂) powder was precalcined (dried) at 800^∘C (200^∘C) prior to mixing. The stoichiometric mixture was pelletized and sintered in air at 1350^∘C for 36 hours, three times, with regrinding and repelletization between sinterings. Fig. 12 shows an x-ray Rietveld refinement against the $Fd\bar{3}m$ space group for a typical powder sample of La₂Zr₂O₇ synthesized for this work.

We provide further details on the results of our 4^th-order NLC calculations for $C_{\mathrm{mag}}$ with a 5% oxidation level included in the calculations. The calculated $C_{\mathrm{mag}}$ with 5% oxidation using the globally (locally) best-fitting exchange parameters that are also able to describe the measured susceptibility, A’ (B’), is shown in Fig. 13(a) (13(b)). To improve convergence of the NLC calculations, we have used the Euler transformation to the 3^rd (Euler 3) and 4^th (Euler 4) orders (see Appendix J). While the parameter sets A’ and B’ are both locally optimal, the A’ description of the $C_{\mathrm{mag}}$ data is clearly superior.

The inset to Fig. 13(b) shows the two locally optimal regions of parameter space for the 4^th-order NLC calculations of $C_{\mathrm{mag}}$, for both a 0% oxidation level and a 5% oxidation level, defined by $\log\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{C\mathrm{mag}}<2.7$ for the purposes of the visualization. From the similarity of these regions and their local minima (A and A’, B and B’), we conclude that the results of the NLC calculations for $C_{\mathrm{mag}}$ are robust to the sample oxidation up to oxidation levels of at least 5%. In Table 1 we summarize the results of our NLC fittings to $C_{\mathrm{mag}}$ and list the best-fitting exchange parameters corresponding to each fitting.

In this section we discuss the results of the NLC-fitting to the magnetic susceptibility measured from a powder sample of Ce₂Zr₂O₇. The magnetic susceptibility is dependent on the values of $J_{\tilde{x}}$, $J_{\tilde{y}}$, $J_{\tilde{z}}$, and $\theta$. The exchange parameters $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$, are given by some permutation of $(J_{a},J_{b},J_{c})$. We allow $\theta$ to vary in the range from 0 to $\pi$/4. This is enough to cover all distinguishable scenarios, since changing the sign of $\theta$ does not affect any quantity considered here, and shifting $\theta$ to $\theta+\pi$/2 is the same as reversing the sign of $\theta$ and swapping the values of $J_{\tilde{x}}$ and $J_{\tilde{z}}$, which is already covered by considering all six permutations of exchange parameters.

NLC calculations up to 2^nd order were performed to compute the powder-averaged magnetic susceptibility and to compare the calculations to the corresponding measurement on Ce₂Zr₂O₇. A constant term was added to the NLC calculations to account for the effect of mixing in higher crystal field levels due to an applied magnetic field. This term is calculated from the low-temperature limit of single ion susceptibility using the crystal-field scheme of Ce³⁺ in Ce₂Zr₂O₇ reported in Ref. [7]. The level of sample oxidation for the measured powder sample had an upper limit of $\sim$14%. This upper limit was estimated from fits to the single ion susceptibility at high temperature using the crystal-field scheme of Ce³⁺ in Ce₂Zr₂O₇ reported in Ref. [7]. Accordingly, a 14% oxidation level is included in our NLC calculations of the magnetic susceptibility.

NLC calculations of the magnetic susceptibility were performed for parameter sets throughout the A and B regions identified by the $C_{\mathrm{mag}}$ fittings. The calculations were compared with experimental data between 1 K and 10 K. Fig. 5 shows the regions of the phase diagram for which it is possible to obtain simultaneous agreement with $C_{\mathrm{mag}}$ and $\chi$, for $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ equal to the different permutations of $(J_{a},J_{b},J_{c})$. We define these regions by the simultaneous satisfaction of $\log\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{C\mathrm{mag}}<2.7$ and $\log\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{\chi}<12.1$. The goodness-of-fit measure $\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{C\mathrm{mag}}$ is defined in the main text, and $\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{\chi}\propto\sum\frac{(\chi^{\mathrm{NLC}}(T_{\mathrm{exp}})-\chi^{\mathrm{exp}}(T_{\mathrm{exp}}))^{2}}{\epsilon(T_{\mathrm{exp}})^{2}}$, where the sum is over measured temperatures $T_{\mathrm{exp}}$ between 1 K and 10 K and $\epsilon(T_{\mathrm{exp}})$ is the experimental uncertainty on the magnetic susceptibility at temperature $T_{exp}$. We allow $\theta$ to vary in the range from 0 to $\pi$/4 in finding the best agreement with the susceptibility data for each permutation. The relatively small experimental uncertainties on the magnetic susceptibility contribute to the larger upper limit for $\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{\chi}$ in comparison to the upper limit used for $\langle\frac{\delta^{2}}{\epsilon^{2}}\rangle_{C\mathrm{mag}}$.

In this section we discuss the analysis of our elastic neutron scattering data, measured on an annealed powder sample of Ce₂Zr₂O₇ and used to place an upper limit of $\mu_{\mathrm{ordered}}$ $\leq 0.04$ $\mu_{\mathrm{B}}$ on the ordered moment corresponding to any All-in, All-out (AIAO) dipole order in Ce₂Zr₂O₇’s magnetic ground state. The strongest magnetic Bragg peaks associated with AIAO order are expected to reside at the $\textbf{Q}=(2,2,0)$ and $\textbf{Q}=(1,1,3)$ positions of reciprocal space. Accordingly, we can examine the temperature dependence of the scattered intensity at these locations in order to look for increases of intensity with decreasing temperature, which would signal the onset of a magnetic Bragg peak and associated magnetic order. As shown in Fig. 14(a), no such increase in intensity is detected upon lowering temperature.

In Fig. 14(b), we show the measured intensity around the $\textbf{Q}=(2,2,0)$ (left) and $\textbf{Q}=(1,1,3)$ (right) positions at $T$ = 0.06 K (blue) and as averaged over the temperatures $T$ = 0.25 K, $T$ = 0.5 K, $T$ = 0.75 K, $T$ = 1 K, and $T$ = 1.5 K (red). Gaussian fits to the peak at each of the locations are also shown for each temperature (or temperature-average) and the area underneath of these Gaussian peaks was used in order to determine the corresponding integrated intensity for each peak. From these values for the integrated intensity, we place an upper limit of $\mu_{\mathrm{ordered}}$ $\leq 0.04$ $\mu_{\mathrm{B}}$ for the Ce³⁺ ordered moment corresponding to any AIAO magnetic dipole ordering in Ce₂Zr₂O₇.

For each selected Bragg peak position Q, an upper limit is calculated in accordance with the equation,

where ${I}^{\mathrm{exp}}_{\mathrm{mag}}$ and ${I}^{\mathrm{exp}}_{\mathrm{nuc}}$ are the measured magnetic and nuclear contributions to the integrated Bragg intensity, respectively. $F(\textbf{Q})$ is the nuclear structure factor and $\mathbf{F^{\mathrm{mag}}_{\perp}}(\textbf{Q})/\mu$ is the component of the magnetic structure factor that is perpendicular to Q, after dividing out the magnitude of the ordered moment ($\mu$) from the calculation [27].

We have used 3^rd order NLC calculations to compute the energy-integrated scattering signals corresponding to a polarized neutron scattering experiment with sample alignment in the $[HHL]$ scattering plane. The exchange parameters are set to $\theta=0$ and $(J_{\tilde{x}},J_{\tilde{y}},J_{\tilde{z}})$ = (0.064, 0.063, 0.011) meV for the calculation and we perform the NLC-3 calculation with $T$ = 0.5 K as that is the lowest temperature for which the NLC expansion converges. Specifically, we compute,

and,

where $S^{\mathrm{SF}}(\mathbf{Q}))$ ($S^{\mathrm{NSF}}(\mathbf{Q})$) denotes the energy-integrated structure factor for SF (NSF) scattering. $N$ is the number of spins in the lattice, $f(|{\bf Q}|)$ is the magnetic form factor for Ce³⁺ (calculated using the analytical approximation in Ref. [28]), $\hat{{\bf n}}$ is the neutron polarisation direction, $\hat{{\bf z}}_{i}$ is the local anisotropy axis for the site $i$ and,

We compute $S^{\mathrm{SF}}(\mathbf{Q}))$ and $S^{\mathrm{NSF}}(\mathbf{Q})$ at $T$ = 0.5 K and in each case we subtract the corresponding calculation at $T$ = 10 K for better comparison with the temperature-subtracted experimental data. In Fig. 8(a, b) (Fig. 8(c, d)) of the main text, we compare the NLC-calculated $S^{\mathrm{SF}}(\mathbf{Q})$ ($S^{\mathrm{NSF}}(\mathbf{Q})$) to polarized neutron scattering measurements performed on an annealed $\sim$1.5 g single crystal sample of Ce₂Zr₂O₇ using the D7 diffractometer at the Institut Laue-Langevin with an incident energy of $E_{i}$ = 3.47 meV and a dilution refrigerator sample environment. The sample was aligned in a copper sample holder in the $[HHL]$ scattering plane with the uniaxial polarization direction perpendicular to the $[HHL]$ plane, and the sample was rotated in 0.5^∘ steps over a total of 250^∘. The data is subsequently folded into a single quadrant of the $[HHL]$ plane and further symmetrized. We have further discussed this symmetrization process in the supplemental material of Ref. [7]. For each data set, we reduce the data in a manner that avoids adding artefacts arising from the subtraction of strong nuclear Bragg peaks. Allowed nuclear Bragg peaks are located at $\mathbf{Q}=(1,1,1),(2,2,2),(1,1,3),(2,2,0),(0,0,4)$, and symmetrically equivalent locations. The intensity at each Bragg peak location is masked in performing the temperature subtraction, and we then show the intensity at these masked Bragg peak locations as the average intensity of the surrounding points in reciprocal space.

The microscopic spin Hamiltonian parameters A and B can be employed in 3^rd order NLC calculations to calculate equal-time (energy-integrated) structure factors, and these can be compared to inelastic neutron scattering measurements on Ce₂Zr₂O₇. The energy-integrated structure factor is:

where $i,j,$ are sublattice indices and $f(|{\bf Q}|)$ is the magnetic form factor for Ce³⁺. Averaging over directions at fixed magnitude $|{\bf Q}|=Q$, gives the powder structure factor $S(Q)$. We integrate over $Q=[0.46,0.93]~{}\text{\AA}^{-1}$ and the result represents the energy-integrated neutron scattering response of a powder sample integrated over that momentum range, as a function of $T$. The structure factor at $T$ = 9.6 K is subtracted to replicate the background subtraction used in the experiment. Lines in Fig. 10 of the main text represent the powder integrated equal-time structure factor calculated in 3^rd order NLC using parameter sets A with $\theta=0$, and B with $\theta=0.561$ radians.

We compare the NLC calculations to the experimentally measured neutron scattering response of a powder sample, integrated over the energy range $[-0.2,0.4]$ meV. This integration range is chosen to enclose the low-lying excitations of the system while avoiding unnecessary contamination to the temperature-subtracted signal, which often becomes more prevalent at higher energies. To further reduce the effect of noise on the experimental data, we integrate in momentum-transfer over the range $\lvert\textbf{Q}\rvert=[0.46,0.93]~{}\text{\AA}^{-1}$. This integration range is chosen to avoid nuclear Bragg peaks while still enclosing the dominant portion of the measured magnetic signal. We find that parameter sets from region A correctly predict the observed increase in scattering over the range $\lvert\textbf{Q}\rvert=[0.46,0.93]~{}\text{\AA}^{-1}$ with decreasing temperature, while parameter sets from region B do not, as shown in Fig. 10 of the main text.

Fig. 15(a) (15(b), 15(c), 15(d), 15(e), 15(f), 15(g)) shows the temperature-difference neutron scattering spectra measured from an annealed powder sample of Ce₂Zr₂O₇ for a $T$ = 0.06 K (0.25 K, 0.5 K, 0.75 K, 1 K, 1.5 K, 3 K) data-set with a $T$ = 9.6 K data-set used a background. The data in Fig. 15(a) (15(c), 15(g)) is also shown in Fig. 9(a) (9(b), 9(c)) with different $\lvert\textbf{Q}\rvert$-range and $\lvert\textbf{Q}\rvert$, $E$ pixel size. These data sets were used to compute the measured integrated-intensity over the energy range $[-0.2,0.4]$ meV and the momentum range $[0.46,0.93]~{}\text{\AA}^{-1}$, forming the data points shown in Fig. 10 of the main text.

The new measured $C_{\mathrm{mag}}$ data also allows for a comparison with the temperature-dependent inelastic neutron scattering signal measured on an annealed powder sample of Ce₂Zr₂O₇ [7]. Specifically, we compare the $C_{\mathrm{mag}}$ data with the imaginary part of the dynamic spin susceptibility, $\chi^{\prime\prime}(\textbf{Q},E)$, calculated from our previously reported neutron data. As shown in Fig. 15(a-g), a signal with the approximate energy range $E=[0,0.15]$ meV is seen to onset in the inelastic neutron scattering spectra with decreasing temperature. The dominant intensity within this signal was used to calculate $\langle\chi^{\prime\prime}(\textbf{Q},E)\rangle$ for each temperature, giving rise to the data points shown in Fig. 16 and allowing us to further examine the temperature dependence of the measured neutron scattering signal. $\chi^{\prime\prime}(\textbf{Q},E)$ was calculated via $\chi^{\prime\prime}(\textbf{Q},E)=S_{0}(\textbf{Q},E,T)(1-e^{-E/k_{\mathrm{B}}T})$, where $S_{0}(\textbf{Q},E,T)$ = $S(\textbf{Q},E,T)-S(\textbf{Q},E,$ $T$ = 9.6 K). This subtraction is used to isolate the magnetic contribution to the measured neutron scattering spectra, and assumes that $\chi^{\prime\prime}(\textbf{Q},E)=0$ at $T$ = 9.6 K. This was used to calculate the average of $\chi^{\prime\prime}(\textbf{Q},E)$ over $\lvert\textbf{Q}\rvert=[0.46,0.93]~{}\text{\AA}^{-1}$, $E=[0,0.15]$ meV, denoted as $\langle\chi^{\prime\prime}(\textbf{Q},E)\rangle$. As shown in Fig. 16, the temperature onset of $\langle\chi^{\prime\prime}(\textbf{Q},E)\rangle$ coincides well with that of the broad-hump in $C_{\mathrm{mag}}$, and $\langle\chi^{\prime\prime}(\textbf{Q},E)\rangle$ continues to grow, separating from $C_{\mathrm{mag}}$, below $T$ $\sim$0.3 K.

Recent theory work on the XYZ model Hamiltonian with $J_{\tilde{x}}$ = $J_{\tilde{y}}$ (which is a relevant approximation for the best-fitting exchange parameters found in this work) has predicted that a U(1) quantum spin ice ground state can be realized upon decreasing temperature through a classical spin ice regime [23, 29, 12]. Furthermore, these works predict that a broad hump in $C_{\mathrm{mag}}$ onsets slowly on entrance into the classical spin ice regime upon decreasing temperature. This prediction is consistent with the coincidence of the temperature onsets of $\langle\chi^{\prime\prime}(\textbf{Q},E)\rangle$ and $C_{\mathrm{mag}}$ shown in Fig. 16.

Here we discuss the semiclassical molecular dynamics calculations of $S(\lvert\mathbf{Q}\rvert,E,T)$ that lead to the calculated spectra shown in Fig. 9(d-i) of the main text. First, classical Monte Carlo simulations were performed using the best fit A exchange parameters, to obtain an ensemble of spin configurations sampled at temperature $T$. We then use these configurations as initial configurations and solve the semiclassical Landau-Lifshitz equation $\frac{d}{dt}\mathbf{S}_{i}$ = - $\mathbf{S}_{i}\times\mathbf{h}_{i}$, where $\mathbf{h}_{i}$ is the effective magnetic field on the spin $\mathbf{S}_{i}$. The dynamical structure factor is obtained as the time-space Fourier transform of the time-evolved magnetic moments, averaged over the ensemble of initial states.

The molecular dynamics solution computes the classical dynamics. That is, it treats the spins as classical magnetic moments precessing in their local field. To compare this to the (quantum) experiment or a theoretical method such as linear spin wave theory, one has to re-scale the classical calculation. This is because the classical dynamical structure factor is symmetric with respect to neutron energy-transfer $E$, and it vanishes as $T$ approaches zero for all $E>0$. Neither of these is the case for the dynamical structure factor of the quantum system. Another, more quantitative, way to think about this is via the fluctuation-dissipation theorem by comparing the version for classical and quantum systems [30]. In particular, for a classical system we get $(\beta E)S_{\mathrm{classical}}(\textbf{Q},E,T)=\chi^{\prime\prime}(\textbf{Q},E,T)$ while for the quantum system it reads $(1-e^{-\beta E})S_{\mathrm{quantum}}(\textbf{Q},E,T)=\chi^{\prime\prime}(\textbf{Q},E,T)$, where $\beta=1/(k_{\mathrm{B}}T)$. It is then reasonable to equate the imaginary part of the susceptibility, $\chi^{\prime\prime}(\textbf{Q},E,T)$, as this quantity is real and symmetric for both the classical and the quantum system. Furthermore, as shown in Ref. [31], $\chi_{\mathrm{quantum}}^{\prime\prime}$ = $\chi_{\mathrm{classical}}^{\prime\prime}$ within linear spin wave theory. Using the quantum and classical fluctuation dissipation theorem for the respective sides then yields,

which is what we use to estimate the dynamical structure factor of the (quantum) experiment using our classical simulation. The dynamical structure factor is then powder-averaged to obtain $S_{\mathrm{quantum}}(\lvert\textbf{Q}\rvert,E,T)$, and convolved with the experimental resolution. In Fig. 9(d-i) of the main text, we show the calculated powder-averaged dynamical structure factor at 0.06 K, 0.5 K, and 3 K, with the powder-averaged dynamical structure factor at $T$ = 9.6 K subtracted from the result.

Note that Eq. (H1) accounts for detailed balance, $S_{\mathrm{quantum}}(\mathbf{Q},-E,T)$ = $e^{-\beta E}S_{\mathrm{quantum}}(\mathbf{Q},E,T)$, since $S_{\mathrm{classical}}(\mathbf{Q},E,T)=S_{\mathrm{classical}}(\mathbf{Q},-E,T)$. Zhang et al. (Ref. [31]) derive the conversion factor $\beta E$ by comparing the classical spin wave theory at finite temperature with the quantum spin wave theory at zero temperature. It is thus valid in the case $\beta E\gg 1$, which is well fulfilled in their case, but not applicable to a large part of our energy and temperature range. However, note that our factor $\frac{\beta E}{1-e^{-\beta E}}$ reduces to $\beta E$ for $\beta E\gg 1$, so our calculation is entirely consistent with this argument.

Heat capacity measurements were performed on our single crystal Ce₂Zr₂O₇ sample, along with a polycrystalline sample of La₂Zr₂O₇, which is used as a 4f⁰ analogue of Ce₂Zr₂O₇. Heat capacity measurements on a polished single crystal of Ce₂Zr₂O₇ (smooth-surfaced pressed powder pellet of La₂Zr₂O₇) were carried out on a Quantum Design PPMS down to $T$ = 0.058 K ($T$ = 2.5 K) using the conventional quasi-adiabatic thermal relaxation technique. The heat capacity of La₂Zr₂O₇ is very small at $\sim$2.5 K, and there was no need to pursue measurements at lower temperatures.

We provide further details on the analysis of $C_{\mathrm{mag}}$’s approach to zero at $T$ = 0 K. Fig. 17(a) shows the results of fitting simple cubic and exponential extrapolations to the measured $C_{\mathrm{mag}}$ data, as well as the low-temperature extrapolation to $C_{\mathrm{mag}}$ which is based upon an interpolation between the results of NLC calculations at $T$ $>$ $\sim$0.5 K and a $T^{3}$ low temperature form appropriate to emergent photons in a U(1) QSL. These extrapolations are also shown in Fig. 3 and Fig. 11(a) of the main text, respectively. We label the latter extrapolation as “interpolation” in Fig. 17 and the following discussion. This interpolation method is introduced in Ref. [26] and discussed for the current context below.

The interpolation method first involves performing a high temperature expansion of the magnetic heat capacity, $C_{\mathrm{mag}}(T)$, corresponding to the XYZ Hamiltonian and the A set exchange parameters, and then turning this into an expansion for the entropy density as a function of energy density, $s(e)$, around $e=0$. If $C_{\mathrm{mag}}(T)$ $\propto$ $T^{3}$ at low temperature, then for $e$ close to the ground state energy density $e_{0}$, $s(e)\propto(e-e_{0})^{3/4}$. A Padé approximant is used to interpolate between those two limits, to obtain $s(e)$ over the region $e$ = $[e_{0},0]$, which can then be converted to $C_{\mathrm{mag}}(T)$ over the range $T$ = $[0,\infty]$. This approach requires an estimate of the ground state energy per site, $e_{0}$. We treat this estimate as an adjustable parameter and set $e_{0}=-0.385J_{a}$ for best agreement with experiment, which is in a physically plausible range.

The approach based on $s(e)$ is generally better behaved than performing the interpolation on $C_{\mathrm{mag}}(T)$ directly, and it obeys the physical constraints on the total energy and entropy $\int_{0}^{\infty}C_{\mathrm{mag}}(T)dT=-e_{0}$, $\int_{0}^{\infty}\frac{C_{\mathrm{mag}}(T)}{T}dT=$ $R\ln(2)$, respectively, by construction. The choice of Padé approximant $P(m,n)$ is constrained to $m+n\leq k$, where $k$ is the maximum order obtained for the high temperature expansion of $C_{\mathrm{mag}}(T)$. In our case $k=13$, and we take the approximant $P(7,6)$, again guided by best agreement with experiment. The estimate of $e_{0}$ and the choice of $m,n$ are the only adjustable parameters in the comparison, with the exchange parameters equal to the set A parameters (see main text or Tab. 1). The comparison between theory and experiment is good, particularly for the entropy curve $S_{\mathrm{mag}}(T)$, when one considers that the experimental entropy is missing $\sim$5% of the expected $R\ln(2)$, due to Ce⁴⁺ substitution which is not incorporated in the interpolation calculation. This demonstrates that the observed $C_{\mathrm{mag}}(T)$ can be consistent with a smooth crossover to a $T^{3}$ form, even though we do not reach the $T^{3}$ regime with the present experimental data.

As shown in Fig. 17(a), the cubic extrapolation cannot be made to connect smoothly to the data at the lowest temperature data points, while the best-fitting exponential extrapolation and the interpolation both meet the data in a smooth manner. The inset to Fig. 17(a) shows the best-fitting simple exponential extrapolation when locking the gap energy to the values $\Delta$ = 20, 30, 40, and 50 mK, and a gap of $\Delta$ = 35(5) mK results from such a naive analysis.

Using each of these extrapolations for $C_{\mathrm{mag}}$ in order to describe the data below the lowest temperature data point, we calculate the entropy recovered via $S_{\mathrm{mag}}=\int_{0}^{T}\frac{C_{\mathrm{mag}}}{T}dT$ and show the results in Fig. 17(b). As shown in Fig. 17(b), the best-fitting cubic extrapolation grossly underestimates the $R\ln(2)$ entropy associated with the CEF ground state doublet, while the exponential extrapolation and the interpolation both saturate to $R\ln(2)$ within the $\sim$5% tolerance associated with the sample oxidation.

We use the NLC method to calculate thermodynamic quantities throughout this work. The method is described in Refs. [32, 33, 34] (for example). Extensive quantities per site $\frac{\langle\mathcal{O}\rangle}{N}$ are represented as sums over contributions from clusters $c$:

where $M_{c}$ is the cluster multiplicity, defined as the number of times $c$ can be embedded in the lattice, per site $N$. $W_{c}$ is the cluster weight:

where $\langle\mathcal{O}\rangle_{c}$ is the expectation value of the quantity $\mathcal{O}$ taken from exact diagonalization on cluster $c$ with open boundary conditions. The second term in Eq. (J2) is a sum over the weights of all subclusters of $c$. The sum in Eq. (J1) is arranged in order of increasing cluster size. At high temperatures, terms from larger clusters vanish faster with increasing temperature and the series converges in the same manner as high temperature expansion. At sufficiently high temperature, one can then justify truncating the sum at finite cluster size.

We employ a series of clusters starting with a single site and then all further clusters are constructed from full tetrahedra. The $n^{\mathrm{th}}$ order of the expansion incorporates clusters of size up to $n$ tetrahedra. We denote the $n^{\mathrm{th}}$ order calculation as “NLC$n$”. For the heat capacity we have performed calculations up to 4^th order (NLC4). For the A parameter set we additionally performed NLC calculations of $C_{\mathrm{mag}}$ up to 7^th order (see inset of Fig. 6). The methodology for these 7^th order calculations is described in Ref. [35]. For $S(Q)$ and $S(\mathbf{Q})$ we have performed calculations up to 3^rd order (NLC3). For the susceptibility we have performed calculations up to 2^nd order (NLC2).

For example, to estimate $S(\mathbf{Q})$ using Eq. (F1) and the NLC method, we define for each cluster $c$ entering the expansion, the extensive quantities:

The NLC estimate of $S({\bf Q})$ is then:

where in this case (3^rd order NLC) we truncate the sum at a maximum cluster size of three tetrahedra. $M_{c}$ are the cluster multiplicities and $W_{c}({\bf Q})$ are the cluster weights

where the sum on the RHS is over subclusters of $c$.

To improve convergence of the $C_{\mathrm{mag}}$ calculations, we have used Euler transformation to the 3^rd and 4^th orders [34]. The Euler transformed results at 3^rd and 4^th order are:

and

where $\langle\mathcal{O}\rangle_{\sf NLCn}$ is the estimate of $\langle\mathcal{O}\rangle$ up to $n^{\mathrm{th}}$ order in NLC.

For the susceptibility calculations we included a population of 14% vacancies in the calculation, with disorder averaging. The disorder average can be taken as order-by-order in NLC. Since vacancy disorder is binary, the disorder average can be done exactly [33]. We have also performed heat capacity calculations with 5% vacancy disorder, as a point of comparison to the calculations with the clean model. The fits of these calculations to the experimental data produce very similar results to those found for the clean model, as shown in Fig. 13.