Quantum Spin Ice Dynamics in the Dipole-Octupole Pyrochlore Magnet Ce$_2$Zr$_2$O$_7$

TL;DR

Neutron scattering on Ce$_2$Zr$_2$O$_7$ reveals a quantum spin ice state with dipole-octupole doublets, defying expectations of magnetic order and highlighting the role of disorder.

Contribution

This study provides experimental evidence of a quantum spin ice state in Ce$_2$Zr$_2$O$_7$, supporting recent theories of symmetry-enriched U(1) quantum spin liquids for dipole-octupole pyrochlores.

Findings

Discovery of a dynamic quantum spin ice state in Ce$_2$Zr$_2$O$_7$

Suppressed scattering near |Q|=0 indicating no long-range order

Disorder, especially oxidation, influences low-temperature behavior

Abstract

Neutron scattering measurements on the pyrochlore magnet Ce$_2$Zr$_2$O$_7$ reveal an unusual crystal field splitting of its lowest $J$=5/2 multiplet, such that its ground state doublet is composed of m$_J$=$\pm$3/2, giving these doublets a dipole - octupole (DO) character with local Ising anisotropy. Its magnetic susceptibility shows weak antiferromagnetic correlations with $\theta_{CW}$=-0.4(2)K, leading to a naive expectation of an All-In, All-Out ordered state at low temperatures. Instead our low energy inelastic neutron scattering measurements show a dynamic quantum spin ice state, with suppressed scattering near |$\textbf{Q}$|=0, and no long range order at low temperatures. This is consistent with recent theory predicting symmetry enriched U(1) quantum spin liquids for such DO doublets decorating the pyrochlore lattice. Finally, we show that disorder, especially oxidation of powder samples, is important in Ce$_2$Zr$_2$O$_7$ and could play an important role in the low temperature behaviour of this material.

Tables

Table 1. Table 1: Comparison between the observed (Obs 1 ,Obs 2 ) and calculated (Calc 1 ,Calc 2 ) energies and intensities of the CEF excitations in Ce 2 Zr 2 O 7 . Scenario 1 refers to the case where no vibronic bound state is present, while Scenario 2 is the case where weak inelastic scattering near 100 meV is identified as arising due to a vibronic bound state.

	Obs1	Calc1	Obs2	Calc2
E1 (meV)	55.9(2)	55.92	55.9(1)	56.06
E2 (meV)	110.5(3)	110.55	106.1(4)	105.97
I1/I2 (arb.units)	1.2(1)	0.98	0.99(15)	0.92

Table 2. Table 2: Eigenstates and eigenfunctions of the spin-orbit ground state manifold written within the | J = 5 / 2 , m J ⟩ ket 𝐽 5 2 subscript 𝑚 𝐽 |J=5/2,m_{J}\rangle basis. These eigenfunctions correspond to the scenario of no vibronic bound-state. A m J  =  ± plus-or-minus \pm 3/2 CEF ground state is stabilized for both scenarios (without and with a vibronic bound state).

E(meV)	-5/2	-3/2	-1/2	1/2	3/2	5/2
E1	0	1	0	0	0	0
E2	0	0	0	0	1	0
E3	0.725	0	0	0.688	0	0
E4	0	0	-0.688	0	0	0.725
E5	0	0	-0.725	0	0	-0.688
E6	0.688	0	0	-0.725	0	0

Equations

$${\mathcal{H}}_{CEF}=B_2^0{\hat{O}}_2^0+B_4^0{\hat{O}}_4^0+B_4^3{\hat{O}}_4^3$$

$${\chi }_{CEF}=\frac{N_A{g}_J^2{\mu }_B^{2}X}{k_{B}Z}\sum _{\alpha }(\sum _n\frac{|⟨n|J_{\alpha }|n⟩{|}^2{e}^{E_n/T}}{T}+$$

$$\sum _n\sum _{m\ne n}\frac{|⟨m|J_{\alpha }|n⟩{|}^2(e^{-E_n/T}-e^{-E_m/T})}{E_m-E_n})$$

Full text

The rare-earth pyrochlore oxides R2B2O7, where R3+ and B4+ consist generally of rare earth and transition-metal ions respectively, display a wealth of both exotic and conventional magnetic ground states. Their R3+ ions decorate a network of corner-sharing tetrahedra, one of the archetypes for geometrical frustration in three dimensions. Due to strong crystal electric field (CEF) effects, the nature of the magnetic interactions in such materials are strongly influenced by their single-ion physics Gardner et al. (2010); Hallas et al. (2018); Rau and Gingras (2018). A naive theoretical description of the magnetic interactions in rare-earth pyrochlores is generally performed by introducing an ad hoc effective single-ion term in addition to Heisenberg exchange interactions. For example, Heisenberg antiferromagnetism with an effective Ising anisotropy leads to non-frustrated All-In, All-Out (AIAO) magnetic order, as seen in several heavy rare earth iridate pyrochlores Sagayama et al. (2013); Disseler (2014) and illustrated in the insert to Fig.1(a). Heisenberg ferromagnetism and an effective Ising anisotropy give rise to a classical spin ice ground state Bramwell and Gingras (2001), as seen in (Ho,Dy)2Ti2O7 Harris et al. (1997); Morris et al. (2009) and illustrated as the 2I2O local structure in the inset to Fig.1(a). However, to capture all the physics that can arise at low temperatures, the magnetic interactions should be projected into pseudo-spin operators acting solely on the low energy CEF states Curnoe (2008); Onoda (2011); Ross et al. (2011); Savary et al. (2012); Rau and Gingras (2015, 2018). This procedure has been applied for example in the Yb3+ Ross et al. (2011); Robert et al. (2015); Thompson et al. (2017) and Er3+ Savary et al. (2012); Petit et al. (2014); Hallas et al. (2017); Petit et al. (2017) XY pyrochlores where CEF effects give rise to effective $S$ = 1/2 quantum degrees of freedom that interact via anisotropic exchange interactions.

More recently, it has been realized that the precise composition of the ground state crystal field doublets in rare-earth pyrochlores is crucial in determining the form of the microscopic Hamiltonian, and in itself, diversifies the possibility of quantum magnetic states Rau and Gingras (2018); Huang et al. (2014). This has been appreciated for some time in the case of non-Kramers doublets, based on magnetic ions with an even number of electrons such as the 4f2 configuration in Pr3+. Only the local $z$-component of the spin operators transforms as a dipole, with the transverse components transforming as quadrupoles Lee et al. (2012); Chen (2017); Martin et al. (2017). This restricts the form of the effective spin Hamiltonian and can stabilize quadrupolar phases that are not present in the phase diagram for dipolar doublets Onoda and Tanaka (2011); Takatsu et al. (2016). For Kramers ions with an odd number of electrons, such as 4f1 in Ce3+, 4f3 in Nd3+ and 4f5 in Sm3+, a crystal field ground state doublet with DO character can be realized where the local $z$ and $x$ components transform as a dipole, but the local $y$ component transforms as an octupole Huang et al. (2014); Li and Chen (2017); Lhotel et al. (2015); Mauws et al. (2018). After a rotation of the pseudo-spins about the $y$ axis, the DO exchange Hamiltonian on the pyrochlore lattice can be reduced to an XYZ model with three independent exchange parameters ($J_{\tilde{x}}$,$J_{\tilde{y}}$,$J_{\tilde{z}}$) Huang et al. (2014); Li and Chen (2017). This Hamiltonian allows for multiple phases to emerge such as AIAO order in the limit of isotropic antiferromagnetic interactions ($J_{\tilde{x}}\approx J_{\tilde{y}}\approx J_{\tilde{z}}$). The DO nature of the doublet allows for octupolar ordered phases Huang et al. (2014); Li and Chen (2017) and also for moment fragmentation within the ground state, as observed in Nd2Zr2O7, where static AIAO order co-exists with dynamic spin ice fluctuations Benton (2016); Petit et al. (2016a). In the limit of dominant antiferromagnetic interactions and strong easy-axis exchange anisotropy, a dipolar quantum spin ice is stabilized so long as the easy-axis is along one of the dipolar components of the DO doublet ($J_{\tilde{x}}>>J_{\tilde{z}}$,$J_{\tilde{y}}$ or $J_{\tilde{z}}>>J_{\tilde{x}}$,$J_{\tilde{y}}$). An octupolar quantum spin ice is favoured if the easy-axis is along the octupole component ($J_{\tilde{y}}>>J_{\tilde{x}}$,$J_{\tilde{z}}$) Huang et al. (2014); Li and Chen (2017).

A promising family of candidate materials for dipolar or octupolar quantum spin ice physics originating from DO doublets are the cerium pyrochlores Ce2B2O7. The Ce3+ ions in the pyrochlore Ce2Sn2O7 are believed to have a DO CEF ground state and to interact via dominant antiferromagnetic interactions, but do not magnetically order down to $T$ = 20 mK Sibille et al. (2015); Li and Chen (2017). The low energy spin dynamics of the cerium pyrochlores remains unexplored and their characterization is key in determining the nature of their possible spin liquid states. In this letter, we gain further insight into the novel magnetism of cerium pyrochlores by reporting new inelastic neutron scattering experiments on powder and single crystal samples of Ce2Zr2O7. Using high energy inelastic neutron scattering, we first confirmed the DO nature of the Ce3+ single ion ground state wave functions in Ce2Zr2O7. We also present low energy inelastic neutron scattering measurements performed on a single crystal of Ce2Zr2O7 and observe diffuse, inelastic magnetic scattering that emerges at low temperatures. The Q dependence of this diffuse scattering is consistent with a symmetry-enriched U(1) quantum spin ice state at low but finite temperatures. Furthermore, we show the quantum spin-ice correlations remain dynamic down to at least 60 mK with no sign of static magnetic order. Our results demonstrate Q signatures of a dynamic quantum spin ice ground state in Ce2Zr2O7, with associated emergent quantum electrodynamics and elementary excitations based on magnetic and electric monopoles as well as emergent photons Hermele et al. (2004); Savary and Balents (2012, 2013); Benton et al. (2012).

Single crystal and powder samples of Ce2Zr2O7 have been grown using floating zone techniques and solid state synthesis. Stabilizing the Ce3+ oxidation state in Ce2Zr2O7 is not simple, and requires growth and annealing in strong reducing conditions to minimize Ce4+ Otobe et al. (2005). As discussed in the Supplemental Material (SM), this is a serious issue, especially in powder samples, where oxidization is observed to occur in powders exposed to air on a time scale on the order of minutes, complicating the exact characterization of the material’s stoichiometry. The oxidization process can be tracked through high resolution x-ray diffraction measurements of the lattice parameter, and it is much slower for single crystal samples. There we can make an estimate of the stoichiometry of the single crystal used in our experiments as Ce2Zr2O7+δ with $\delta$ $\sim$ 0.1.

We first present high energy inelastic neutron scattering measurements, which probe the degeneracy breaking of the Ce3+ spin-orbit manifold that occurs due to CEF effects. To do so, we used the SEQUOIA high resolution inelastic chopper spectrometer Granroth et al. (2010) at the Spallation Neutron Source of Oak Ridge National Laboratory and employed neutrons with incident energies (Ei) of 150 and 500 meV. The Ei = 150 meV instrument setting was chosen to resolve the crystal electrical field (CEF) states that belong to the spin-orbit ground state manifold ($J$ = 5/2). The CEF interaction lifts the Ce3+ spin-orbit ground state degeneracy into three different eigenstates that are each doubly degenerate. We also estimated a CEF Hamiltonian for Ce2Zr2O7 using a scaling procedure based on the Er3+ pyrochlore CEF scheme Gaudet et al. (2018a). This predicts two CEF excited states near 80 and 100 meV with similar inelastic neutron scattering intensity at $T$ = 5 K. As seen in Fig.1(b), this scenario is in qualitatively good agreement with our 150 meV inelastic neutron experimental spectra where two strong magnetic excitations are observed at $\sim$ 56 and $\sim$ 112 meV. The relative scattered intensity of these CEF transitions can be obtained once minor corrections for a phonon background are made. This procedure is further described in the SM and gives $I_{56meV}$/$I_{112meV}$ = 1.2(1), in good agreement with our expectations based on this scaling argument.

Additional weak inelastic scattering whose Q-dependence is inconsistent with phonons is also visible in the spectra, for example weak scattering near $\sim$ 100 meV in Fig.1 (b). It is not clear if this weak inelastic scattering is due to the influence of Ce3+ or Zr3+ in defective sites Sala et al. (2018), on residual Ce4+, or on the possible presence of hybridized phonon - crystal field excitations known as vibronic bound states, as has been recently observed in holmium and terbium pyrochlores Gaudet et al. (2018b); Ruminy et al. (2016). In any case, this unidentified contribution to the inelastic scattering yields a small fraction of the spectral weight and we conclude the features at 56 and 112 meV are the CEF excitations corresponding to the main Ce3+ site.

The details of the crystal field analysis determining the full set of eigenvalues and eigenfunctions for Ce3+ are summarized in Fig.1(c), and further discussed in the SM. The key conclusion is that the ground state Kramers doublet appropriate to Ce3+ is well separated from all excited crystal field states (with a gap of $\sim$56 meV), and is composed of pure $m_{J}$ = $\pm$ 3/2 states. A large CEF gap is consistent with the high temperature heat capacity of Ce2Zr2O7 measured in ref. Popa et al. (2008) where no Schottky anomaly is observed between 5 and 300 K. These pure $m_{J}$ = $\pm$ 3/2 states have a dipole-octupole character with a dipolar moment whose anisotropy is purely Ising and whose magnitude must be 1.286 $\mu_{B}$. This result does not originate from a fine-tuning of the CEF parameters, but is instead a property protected by the point-group symmetry of the A-site in the pyrochlore lattice.

Figure 1(a) shows the inverse magnetic susceptibility of a 107 mg powder sample of Ce2Zr2O7 measured with a Quantum Design MPMS magnetometer equipped with a 3He insert. The dc magnetic susceptibility measurements were made with applied magnetic fields of $\sim$0.001 T, $\sim$0.01 T, and $\sim$0.1 T, and no significant field dependence was observed. The main panel of Fig.1(a) shows the high temperature susceptibility of Ce2Zr2O7 and reveals strong non-linearity. Assuming a dilution of the Ce3+ moments by non-magnetic Ce4+ ions at the $\sim$8 % level, the Van-Vleck susceptibility calculated with the CEF Hamiltonian of Ce2Zr2O7 reproduces the high-temperature susceptibility data well and yields an antiferromagnetic Curie constant of - 0.4(2) K. We expect conventional and unfrustrated AIAO order in Ce2Zr2O7 based on the effective antiferromagnetic interactions and the Ising anisotropy associated with its magnetism. However, our magnetic susceptibility measurements (inset of Fig.1(a)) as well as both powder and single crystal neutron diffraction experiments show no indication of long range magnetic order down to $T$ = 0.06 K. In particular and as shown in the SM, no new Bragg scattering or enhancement of the Bragg scattering associated with any k = 0 magnetic structure is observed, including at those wave vectors characteristic of the AIAO, $\Gamma_{3}$ structure. Ce2Zr2O7 therefore remains disordered to $T$ = 0.06 K, our lowest temperature measured. The lack of magnetic order in Ce2Zr2O7 likely indicates a strong easy-axis anisotropy in its exchange Hamiltonian.

We examined the low-temperature spin dynamics in Ce2Zr2O7 using the low energy DCS neutron chopper spectrometer at NCNR with Ei = 3.27 meV incident neutrons giving an energy resolution of $\sim$0.09 meV at the elastic line. One experiment was performed on a $\sim$ 6 gram powder sample and a second one was performed on a $\sim$ 5 gram single crystal, which was mounted with its [HHL] plane coincident with the horizontal plane of the spectrometer. Figure 2(a) shows the DCS measurements on our powder, where the integration in $|$Q$|$ is 0.35 $\text{\AA}^{-1}$ to 0.85 $\text{\AA}^{-1}$. This integration in momentum transfer $|$Q$|$ corresponds to integrating over the $|$Q$|$ = $|$(001)$|$ position ($\sim$ 0.59 $\text{\AA}^{-1}$), where quantum spin ice correlations are expected to be strongest Benton et al. (2012). A build up of inelastic spectral weight below $\sim$ 0.4 meV is observed on decreasing the temperature.

Low energy inelastic neutron scattering from our single crystal is shown in Fig.2(b,c,d) and Fig.3(a). All this data was acquired using the same Ei = 3.27 meV instrument configuration of DCS, and Fig.2(b,c,d) shows powder-averaged single crystal data. Figure 2(b) shows the full powder-averaged spectrum at $T$ = 0.06 K with a $T~{}=~{}2~{}K$ data set subtracted from it. This result shows enhanced inelastic scattering at low temperature, which peaks up at $|$Q$|$ $\sim$ 0.59 $\text{\AA}^{-1}$, that is the magnitude of the Q = (001) position in reciprocal space. This is explicitly shown via the $|$Q$|$-cut of the data presented in Fig.2(c). Importantly, Fig.2(c) shows no enhancement of the low energy inelastic scattering around $|$Q$|$ = 0, consistent with expectations for a U(1) quantum spin ice. Finally, Fig.2(d) shows energy cuts through the full difference spectrum shown in (b), taken by integrating in $|$Q$|$ from 0.35 to 0.85 $\text{\AA}^{-1}$, so around $|$Q$|$ = $|$(001)$|$, and also well away from $|$Q$|$ = $|$(001)$|$, integrating from 1.3 to 1.8 $\text{\AA}^{-1}$. This clearly shows the quantum spin ice correlations to be dynamic in nature, characterized by an energy less than $\sim$ 0.15 meV.

With the energy range of the dynamic quantum spin ice correlations identified, we can look explicitly at this scattering from the single crystal, but now comparing Q maps of these correlations to the expectations of both classical, near-neighbor spin ice (without dipolar interactions) and a U(1) quantum spin ice. Fig.3(a) shows $T$ = 0.06 K - $T$ = 2 K data integrated between 0 and 0.15 meV, folded into a single quadrant of the [HHL] map and further symmetrized. The details of this data symmetrization are in the SM. For reference, a theoretical simulation of the structure factor expected for classical near-neighbor spin ice Benton et al. (2012) is shown in Fig.3(b), and that for a U(1) quantum spin ice at low but finite temperature Benton et al. (2012) is shown in Fig.3(c). While these theoretical predictions have similarities, the structure factor for U(1) quantum spin ice has minima in intensity near Q = 0, while the intensity of the structure factor is maximal there for classical near-neighbor spin ice.

Clearly, the measured dynamic S(Q) shows a qualitatively stronger resemblance to the expectations of the symmetry enriched U(1) quantum spin ice Hermele et al. (2004); Savary and Balents (2012, 2013); Benton et al. (2012). The quantum spin ice ground-state is one of various spin liquids that are supported by a model of well isolated DO CEF doublets on the pyrochlore lattice Huang et al. (2014); Li and Chen (2017). A similar dynamic S(Q) is expected in the case of classical dipolar spin ice (here dipolar refers to long range dipolar interactions between magnetic dipoles, as opposed to quadrupoles, octupoles etc), which also shows the suppression of diffuse scattering near $|$Q$|$ = 0 den Hertog and Gingras (2000); Bramwell et al. (2001). Although a definitive conclusion can only be reached once a full spin Hamiltonian is parametrized, the Ce3+ ions in Ce2Zr2O7 have a moment of 1.286 $\mu_{B}$. These moments are roughly a factor of 8 smaller than those associated with Ho3+ or Dy3+ in the classical dipolar spin ices Ho2Ti2O7 and Dy2Ti2O7. The resulting long range dipole terms are expected to be $\sim$ 64 times weaker in Ce2Zr2O7, making such a scenario unlikely. This suggests the spin-ice correlations in Ce2Zr2O7 originate from quantum effects. An octupolar ordered state is also consistent with the lack of magnetic dipole order in Ce2Zr2O7. However, the neutron scattering spectra associated with such an octupolar ordered phase has yet to be calculated, thus we cannot compare our data in Fig.3(a) to it.

The effect of disorder in Ce2Zr2O7 is still an open question as we are aware that our single crystals have some low levels of oxidation. Furthermore, stuffing Ross et al. (2012); Koohpayeh et al. (2014); Baroudi et al. (2015); Arpino et al. (2017) (site-mixing) is expecting to be important in Ce2Zr2O7, because both undesired Ce4+ and Zr3+ ions are chemically stable. It is known that small amounts of disorder can have a drastic impact on the physical properties of frustrated pyrochlore magnets Ross et al. (2012); Koohpayeh et al. (2014); Arpino et al. (2017). It will then be important to further optimize the growth procedure and annealing techniques of Ce2Zr2O7. However, we believe that our inelastic neutron scattering results rule out the scenario of a sensitive AIAO order. Indeed, the conventional impact of quenched disorder on a pyrochlore antiferromagnet would be spin glass physics with diffuse scattering peaked for Q’s corresponding to the Bragg positions of the AIAO state. Here, we observe strong diffuse scattering at Q = (001), which is not only strictly zero for an AIAO state, but also forbidden for all k = 0 long-range ordered magnetic structures allowed by symmetry of the pyrochlore lattice. We thus conclude that our work demonstrates Ce2Zr2O7 to be one of a very few candidates for quantum spin ice physics. Other candidates for quantum spin ice physics are based on Pr3+ and Tb3+ pyrochlores Zhou et al. (2008); Fritsch et al. (2013); Kimura et al. (2013); Sibille et al. (2016, 2018). However, in contrast to Pr3+ and Tb3+, Ce3+ is a Kramers ion and its magnetism is thus further protected against disorder, which in and of itself, can drive a spin liquid state for non-Kramers doublets Savary and Balents (2017); Wen et al. (2017); Martin et al. (2017); Benton (2018). Furthermore, Tb3+ and Pr3+ pyrochlores display low lying CEF field states, which complicate their theoretical understanding due to multipolar interactions Molavian et al. (2007); Takatsu et al. (2016); Petit et al. (2016b). For all these reasons, the cerium pyrochlores are an excellent theoretical and experimental template to investigate quantum spin ice physics.

I Supplemental Material:

I.1 Powder synthesis and Single crystal growth

The polycrystalline and powder samples of Ce2Zr2O7 used in this work were first prepared by arc melting stoichiometric amounts of CeO2, Zr, and ZrO2 in an argon atmosphere, followed by a regrinding and firing at 1000∘C for two days in flowing hydrogen. Prior to any measurements, the powder samples were further annealed in flowing hydrogen at 1000∘C for several hours. An x-ray refinement against the $Fd\overline{3}m$ space group is shown in Fig.S1(a) for one of the Ce2Zr2O7 powder samples, synthesized and annealed using this protocol. Previous studies have characterized the level of oxidation of Ce2Zr2O7+δ powder samples using the value of the lattice parameter, $a$, refined from x-ray diffraction measurements. For example, the values $a$ $\sim$ 10.735 Å has been refined for $\delta$ = 0 and $a$ $\sim$ 10.66 Å has been refined for $\delta$ = 0.5. A linear relationship has been observed between these limits Otobe et al. (2005); Thomson et al. (1999, 1996). Our materials synthesis and powder x-ray diffraction is consistent with these previous works, and, as shown in Fig.S1(a), our stoichiometric powder sample has a lattice parameter of $a$ = 10.735(5) Å. As we will discuss, we observe the oxidation of the powder sample left exposed to air to begin to occur on the order of minutes. There is an obvious colour change which occurs between the stoichiometric (annealed) and oxidized samples, wherein a light green/yellow sample is obtained for the $\delta$ = 0 annealed sample, and a black powder is obtained for the $\delta$ = 0.5 oxidized sample.

Single crystals of Ce2Zr2O7 were obtained using optical floating zone growth, with Xe lamps, from polycrystalline feed stock. The powder rods were first prepared through a solid state reaction using high purity materials. Stoichiometric mixtures of Ce2O3 (99.995$\%$) and ZrN (99.5$\%$) were mixed in a ball mill and pressed into rods. The rods were heated in air to 900∘C in a covered alumina crucible for 5 h. The solidified rods were then re-ground in a ball mill, repressed into rods 8 cm and 6 cm in length to be used as the feed and seed rods for the optical floating zone growth. These rods were then heated to 1550∘C for 3 h in an atmosphere containing a ratio of 95/5 argon to hydrogen. During the actual optical floating zone growth, we employed a growth rate of 2.5 mm/hour while counter-rotating feed and seed rods at 15 rpm in an argon atmosphere containing 5$\%$ hydrogen. Further annealing of the single crystals for 36 hours at 1200∘C was also performed prior to all experiments. A photograph of a typical single crystal of Ce2Zr2O7 obtained through this protocol is shown in Fig.S1(b). A broken surface is shown and one can see that the inside of the single crystal piece is bright yellow while the surface exposed to air is black, indicating some remaining oxidation at the surface after the annealing process. Refinement of a yellow single crystal piece from the unexposed inner layer of the crystal is shown in Fig.1(c). Such x-ray diffraction measurements on crushed single crystals give a lattice parameter refinement of $a$ = 10.718(5) Å for the bulk inner layer of the crystal, as shown in Fig.1(c), and this can be use to estimate the oxidation state of our annealed single crystal sample, Ce2Zr2O7+δ, to be $\delta\sim$ 0.1.

Finally, we characterized the time dependence of Ce2Zr2O7+δ by collecting several powder diffraction patterns from a broken, ceramic rod of material, for differing exposure times in air following a 1200∘C, 6 hour hydrogen annealing protocol. The lattice parameter extracted from x-ray diffraction measurements is shown as a function of time in Fig.S1(d) for this polycrystalline sample, whose oxidation rate should be between the fast rate of the powder samples of Ce2Zr2O7+δ and slow rate of the relatively-well behaved single crystal samples of Ce2Zr2O7+δ. The inset shows x-ray diffraction scans of the Q = (222) Bragg peak as a function of scattering angle and exposure time. Not only do the Bragg peaks shift higher in 2$\theta$ with exposure (sample oxidation), they also significantly broaden and eventually split into multiple distinguishable peaks. This broadening and splitting corresponds to a distribution of oxidation in the sample, and to the eventual formation of separate majority and minority phases with different levels of oxidation, such as the bulk inner and thin outer layers of differing color in the case of single crystal Ce2Zr2O7+δ. The x-ray intensity at the Q = (222) Bragg position for each data set was fit to Lorentzian peak shapes according to the number of peaks present around that Bragg position. The average lattice constant for each data set with multiple distinguishable peaks ($>$ 29.5 hours) was determined by using a weighted average of the lattice constants determined from these peaks, with the integrated intensity of the peaks used as the weights for this average. Figure S1(d) should then give the time dependence of the volume-averaged lattice constant in this polycrystalline sample.

I.2 Determination of the CEF Hamiltonian

As discussed in the main manuscript, the CEF eigenstates for Ce3+ in Ce2Zr2O7 were expressed in the $|J=5/2,m_{J}\rangle$ manifold and fitted using the low energy CEF transitions observed in the Ei = 150 meV neutron scattering spectra. We refined a CEF Hamiltonian for Ce2Zr2O7 by constraining the relative scattered intensities and energies of the CEF transitions observed in our Ei = 150 meV spectra (Fig.1(a) of the main letter). The following CEF Hamiltonian was used in this work:

[TABLE]

The protocol and the expression of the CEF Hamiltonian in terms of Stevens operators used in this work is identical to ref. Gaudet et al. (2015, 2018a) and we refer the reader to these works for further details. It is worthwhile mentioning that the CEF Hamiltonian appropriate for the R3+(D3d) site of the pyrochlore lattice includes six Stevens operators (B${}_{2}^{0}$, B${}_{4}^{0}$, B${}_{4}^{3}$, B${}_{6}^{0}$, B${}_{6}^{3}$, B${}_{6}^{6}$), however the three B${}_{6}^{n}$ terms are zero for Ce3+.

The B${}_{2}^{0}$, B${}_{4}^{0}$ and B${}_{4}^{3}$ terms of the CEF Hamiltonian were refined using the observed energies and scattered neutron intensities of the magnetic excitations near 58 and 110 meV. It was possible to extract this information using a constant $|$Q$|$-cut from the 150 meV data. We obtain this by integrating the 150 meV spectra (Fig.1(a) of the main letter) between $|$Q$|$ = 4.5 and 5.5 Å*-1*. The resulting integration is shown in Fig.S2(a). The two CEF transitions, so-identified by their $|$Q$|$-dependence, are indicated by the black arrows in Fig.S2(a), while other features observed in Fig.S2(a) are likely optical phonon excitations. The CEF transitions as well as the phonon excitations were fitted using a Lorentzian function in energy. The resulting fit of the 150 meV data set is plotted in Fig.S2(a). From this fit, we determined that the energies of the CEF levels within the spin-orbit $J$-manifold are E1 = 55.9(1) and E2 = 110.5(1) meV. The relative intensity of the CEF transitions between the CEF ground state to the first (I1) and the second (I2) excited state is found to be I1/I2 = 1.2(1). Using these three constraints (E1, E2, I1/I2), we refined B${}_{2}^{0}$ = 0.835 meV, B${}_{4}^{0}$ = 0.299 meV and B${}_{4}^{3}$ = 2.875 meV. Comparison between the calculation and data is shown in Table.S1. The calculated eigenstates and eigenfunctions using this CEF Hamiltonian are reported in Table.S2.

In the main letter, we also reported a weaker magnetic excitation near 100 meV, so just below E2. It is possible that magneto-elastic coupling between the CEF ground state and the 2nd excited CEF state (E2) is sufficiently strong to drive a vibronic bound state, which would effectively split the single-ion CEF excitation at E2 into two different excitations Gaudet et al. (2018b). Here, these two excitations would correspond to the one at 110 meV and the one at 100 meV. Thus, within the scenario of a vibronic bound state as the origin of the inelastic scattering at 100 meV, we underestimated I2. We can compensate for this by adding the integrated intensity of the 100 meV transition. By doing so, a new optimization of the CEF Hamiltonian could be refined, leading to B${}_{2}^{0}$ = 0.455 meV, B${}_{4}^{0}$ = 0.295 meV and B${}_{4}^{3}$ = 2.582 meV. Table.S1 also shows the comparison between the data and the calculation within such a scenario. This scenario in which a vibronic bound state is responsible for the weak inelastic intensity near 100 meV also leads to a CEF ground state that is a pure mJ = $\pm$3/2 state. Thus the conclusion that the CEF ground state doublet for Ce2Zr2O7 is a pure mJ = $\pm$3/2 doublet with dipolar octupolar character is robust.

We also present the Ei = 500 meV inelastic neutron scattering spectra of Ce2Zr2O7 in Fig.S2(b). This inelastic scattering data shows a clear magnetic feature near 275 meV that corresponds to a transition originating from the CEF ground state to the highest $J$-manifold ($J$ = 7/2). Two additional, albeit weaker, excitations can be identified by taking appropriate $|$Q$|$ cuts of this data and comparing to the measured empty can, background scattering. This is shown in Fig.S2(c), where a $|$Q$|$ integration from $|$Q$|$ = 7 to 10 Å*-1* is performed. This reveals three CEF transitions $\sim$ 260, 310 and possibly 340 meV. The location in energy of these highest $J$ = 7/2 states is consistent with estimates of the spin-orbit coupling strength ($\lambda$) for Ce3+ Freeman and Watson (1962); Blume et al. (1964).

Finally, the determination of our CEF Hamiltonian appropriate for Ce2Zr2O7 was further validated via calculation of its Van-Vleck susceptibility that can be compared with the temperature dependence of its measured magnetic susceptibility. The following Van-Vleck susceptibility term (ref.VanVleck (1932)) was computed for a powder sample:

[TABLE]

where $\alpha=x,y,z$, NA is the Avogadro constant, $g_{J}$ is the Landé $g$-factor, $k_{B}$ is the Boltzman constant, $\mu_{B}$ is the Bohr magneton and $Z=\sum_{n}e^{-E_{n}/T}$ is the partition function. The factor $X$ was used to parametrize the dilution of the Ce3+ ions into non-magnetic Ce4+ and was refined to 0.92(2), which is in good agreement with our estimate of the oxidation level in our samples using the refined lattice parameters (see section 1 of the SM).

I.3 Details on the lack of magnetic order in Ce2Zr2O7

An orange cryostat with a dilution refrigerator insert was used for both of our DCS experiments. The powder sample was wrapped in of a copper foil inside a copper can sealed under 10 atm of He. For the single crystal experiment, the sample was aligned and mounted on a copper mount. For the data collection of our single crystal DCS experiment, the single crystal was rotated for a total of 270∘ with steps of 0.5 ∘ and we counted for 5 minutes per angle at $T$ = 0.06 K and 3 minutes per angle at $T$ = 2 K. The data was analyzed and plotted using DAVE Azuah et al. (2009).

We demonstrate the evidence for no magnetic order in Ce2Zr2O7 to temperatures as low as $T$ = 0.06 K using our single crystal DCS experiment. We do this by noting that no new Bragg peaks are observed at $T$ = 0.06 K compared with $T$ = 2 K, and by explicitly isolating the elastic scattering for each k = 0 Bragg peak accessible in our experiment. This elastic scattering is shown in Fig.S3 at both $T$ = 0.06 K and $T$ = 2 K for the ${\bf Q}$ = (220),(113),(111),(222) and (002) positions of reciprocal space, with $E_{i}$ = 3.27 meV incident neutrons and energy integration in the range from -0.1 to 0.1 meV. None of these selected k = 0 positions in reciprocal space show any significant changes in intensity with temperature as can be seen by the subtraction of the $T$ = 2 K data set from the $T$ = 0.06 K data set in each plot of Fig.S3. The elastic cuts through the ${\bf Q}$ = (220),(113) and (111) positions were taken along the (HH0) direction with integration in (00L) from L = -0.1 to 0.1 r.l.u., 2.9 to 3.1 r.l.u., and 0.9 to 1.1 r.l.u. respectively and the elastic cuts through the ${\bf Q}$ = (222) and (002) positions were taken along the (00L) direction with integration in (HH0) from H = 1.9 to 2.1 r.l.u. and -0.1 to 0.1 r.l.u. respectively. Figure S3(f) shows a table qualitatively outlining which of these positions is expected to show magnetic intensity for each k = 0 magnetic structure permitted by the pyrochlore lattice Hallas et al. (2016). A green check mark indicates the presence of magnetic intensity at that location due to the corresponding ordered structure, while the red symbol indicates that the corresponding structure results in no magnetic intensity at that location. For example, the $\Gamma_{3}$ structure corresponding to AIAO order generates magnetic intensity at the (220) and (113) positions only. As can be seen from this Table, the fact that we have measured no magnetic intensity at each of the (220), (111), (113), (222) and (002) positions, signifies a lack of ${\bf k}$ = 0 magnetic order in Ce2Zr2O7 down to $T$ = 0.06 K. Furthermore, the whole elastic ${\bf Q}$ map within the (HHL) plane can be plotted for 0.06 K and this is shown in Fig.S3(g). Only the expected structural Bragg peaks are visible in Fig.S3(g), which confirms the lack of ${\bf k}$=0 magnetic order.

We note that our powder and single crystal DCS experiments that both show no magnetic order in Ce2Zr2O7 are further corroborated with cold neutron triple-axis experiment using the SPINS instrument, also at NIST. For this experiment, the sample was also mounted on a copper mount and aligned within the (HHL) plane. An orange cryostat with an He3 insert was used for this experiment. Figure S4 shows $\theta$-2$\theta$ scans for several Bragg peaks at both 5.2 K and 0.3 K collected with neutron incident energies of 5 meV. It is clear from the temperature difference data that our subtraction of both temperature data sets that our triple-axis experiment also does not support the existence of ${\bf k}$ = 0 magnetic order in Ce2Zr2O7 down to 0.3 K.

Finally, low temperature magnetic susceptibility measurements of Ce2Zr2O7 were performed in several fields with both field cooled (FC) and zero field cooled (ZFC) protocol. As seen in Fig.S5, our low temperature susceptibility measurements do not reveal any sign of magnetic ordering or spin freezing down to $\sim$ 0.5 K. Furthermore, the general behaviour of the susceptibility is not impacted by the field value used for these measurements.

I.4 Details of the symmetrization process in our low energy inelastic neutron scattering (DCS) experiment on a single crystal

In this section, we outline the symmetrization process used to average the Q-dependence of the inelastic temperature difference data from an annealed single crystal sample of Ce2Zr2O7 using the DCS instrument, shown in Fig.3 of the main manuscript. Figure S6(a) displays the subtraction of a $T$ = 2 K data set from a $T$ = 0.06 K data set, integrated in energy over the range from 0 to 0.15 meV , with $E_{i}$ = 3.27 meV incident neutrons. This is the unsymmetrized data that was used to produce the symmetrized data shown in Fig.3(a) of the main letter. The intermediate steps in this symmetric averaging are shown in Fig.S6(b) and Fig.S6(c) where the data has been folded upon itself and averaged with respect to the (00L) and (HH0) lines of symmetry. The fully symmetrized data is displayed in Fig.S6(d) to further illustrate the process. Only one quadrant of this data set is truly independent; the other three quadrants are required to impose the symmetry of the reciprocal lattice. No other constraints are imposed on the data. This data set corresponds to symmetric averaging with respect to both the (00L) and (HH0) lines of symmetry. This symmetrization helps to bring out relatively weak signals, as such signals are measured at more than one equivalent Q position, and averaged over. Non-symmetric features in the scattering, resulting from, for example self-absorption, are also averaged over and hence diminished in the process.