Dipolar Spin Ice Regime Proximate to an All-In-All-Out Néel Ground State in the Dipolar-Octupolar Pyrochlore Ce 2 Sn 2 O 7

The dipolar-octupolar (DO) pyrochlores, R 2 M 2 O 7 (R = Ce, Sm, Nd), are key players in the search for realizable novel quantum spin liquid (QSL) states as a large parameter space within the DO pyrochlore phase diagram is theorized to host QSL states of both dipolar and octupolar nature. New single crystals and powders of Ce 2 Sn 2 O 7 , synthesized by hydrothermal techniques, present an opportunity for a new characterization of the exchange parameters in Ce 2 Sn 2 O 7 using the near-neighbor XY Z model Hamiltonian associated with DO pyrochlores. Utilizing quantum numerical linked cluster expansion ﬁts to heat capacity and magnetic susceptibility measurements, and classical Monte Carlo calculations to the diffuse neutron diffraction of the new hydrothermally grown Ce 2 Sn 2 O 7 samples, we place Ce 2 Sn 2 O 7 ’s ground state within the ordered dipolar all-in-all-out (AIAO) Néel phase, with quantum Monte Carlo calculations showing a transition to long-range order at temperatures below those accessed experimentally. Indeed, our new neutron diffraction measurements on the hydrothermally grown Ce 2 Sn 2 O 7 powders show a broad signal at low scattering wave vectors, reminiscent of a dipolar spin ice, in striking contrast from previous powder neutron diffraction on samples grown from solid-state synthesis, which found diffuse scattering at high scattering wave vectors associated with magnetic octupoles and suggested an octupolar quantum spin ice state. We conclude that new hydrothermally grown Ce 2 Sn 2 O 7 samples host a ﬁnite-temperature proximate dipolar spin ice phase, above the expected transition to AIAO Néel order.


I. INTRODUCTION
In f -electron systems, even apparently simple pseudospin-1/2 magnets can exhibit intriguing single-ion symmetry properties: beyond the "usual" case of three dipolar spin components, some components may either be time-reversal symmetric quadrupoles, or transform as octupoles.These multipolar moments, especially when combined with magnetic frustration, can lead to a wealth of exotic phenomena.Of particular interest are "hidden" orders 1-4 elusive to common experimental probes, and long-range entangled quantum spin liquids (QSL).
The theoretical description of these systems naturally goes along with sets of new model Hamiltonians, which have indeed been found to host interesting phases, although the balance between the competing interactions which stabilize these phases can be delicate.In fact, frustrated materials commonly find themselves near phase boundaries between different orders, or proximate to exotic phases [5][6][7][8][9][10][11] .The possibility of tuning across these phase boundaries through chemical pressure, applied external pressure, or disorder makes the f -electron materials attractive systems for revealing unconventional physics.
To advance our understanding of these multipolar systems, it is therefore imperative to develop a good microscopic description of relevant real materials, with respect to both the chemical properties of the samples under investigation and the appropriate Hamiltonian for their theoretical analysis.
One ideal set of materials to develop such understanding is the dipolar-octupolar (DO) pyrochlores, R 2 M 2 O 7 , where R 3+ ={Ce, Sm, Nd} and M 4+ can be a range of transition metal ions, as they host a wealth of exotic phenomena due to their canonical frustrated geometry, range of interactions, single-ion anisotropies, and the ability to map to a simple pseudospin-1/2 system.In these 4f rare-earth magnets, the hierarchy of energy scales is headed by Coulomb interactions and spin-orbit coupling, followed by crystal electric field (CEF) effects from the surrounding charged environment, followed by intersite exchange interactions.It is the CEF that breaks the (2J + 1) degenerate free-ion ground state and leads to a ground state doublet that is classified by how it transforms under the point group symmetry of the rare-earth site (D 3d ) and timereversal symmetry 12 .When the ground state CEF doublet is described by a pseudospin-1/2 in which different components transform as magnetic dipoles and octupoles, respectively, it is referred to as a dipole-octupole (DO) doublet 13 .
The symmetry of the DO pyrochlores constrains the nearestneighbor Hamiltonian to take the form of a simple XY Z Hamiltonian, where S α are the pseudospin components in the local α = x, ỹ, z coordinate frame (related to the local x, y, z coordinate frame, by a rotation of θ about the y axis), g z is the g factor, h is the magnetic field, and ẑi is the local anisotropy axis.The ground state phase diagram for the DO pyrochlores has been investigated in detail [14][15][16][17] , and the best current understanding is presented in Figure 1.It is clear and encouraging that this class of materials has huge potential to realize QSL states of both dipolar and octupolar character, as they occupy a large region of parameter space, alongside dipolar and octupolar -1 -0.5 0 0.5 1 1.5 (J x -J y )/(J x +J y ) J z /(J  16 , which find four U(1) QSL phases, labeled QSL-0 or π (dip or oct), alongside dipolar and octupolar AIAO Néel phases.The dip or oct distinction refers to whether the dominant pseudospin component transforms as a dipole or an octupole, while the 0 or π refers to the flux through the hexagonal plaquettes within the pyrochlore structure.The strong potential for realizable QSL materials is a huge motivation to further explore the DO pyrochlores in detail.
all-in-all-out (AIAO) Néel phases.Within this exciting class of materials, the 4f 1 Ce 3+ family in particular has attracted intrigue as experimentally realizable QSL candidates [18][19][20][21][22] .Experimental work on Ce 2 Zr 2 O 7 , including extensive neutron scattering, heat capacity, and magnetic susceptibility, combined with theoretical calculations have constrained Ce 2 Zr 2 O 7 to lie within the U(1) π QSL state at low temperatures, on a boundary between octupolar and dipolar character, although the exact nature of the ground state is under debate [22][23][24] .Similarly, much recent interest was generated toward Ce 2 Sn 2 O 7 when neutron diffraction measurements on powder samples grown by solid-state synthesis found broad diffuse scattering at high scattering wave vectors (|Q| > 5 Å −1 ), attributed to scattering from an octupolar spin ice state or from higher order multipoles 19,25 .
Here we present an unexpected and remarkable result obtained by comparison of detailed quantum numerical linked cluster expansion (NLCE) calculations and Monte Carlo simulations to new Ce 2 Sn 2 O 7 single crystal and polycrystalline samples grown by hydrothermal synthesis techniques.Our work provides a constrained parametrization of a model Hamiltonian that allows the results for Ce 2 Sn 2 O 7 presented here to be placed on a phase diagram and contrasted with previous measurements of Ce 2 Sn 2 O 7 , as well as the related Ce 2 Zr 2 O 7 .Contrary to the previously reported octupolar spin ice state, we report that a dipolar all-in-all-out ground state is appropriate for Ce 2 Sn 2 O 7 , with a spin ice regime above the AIAO ordering.New diffuse neutron scattering measurements show a broad signal at low scattering wave vectors only (peaked near |Q| ∼ 0.6 Å −1 ), and no evidence of higher-Q scattering indicative of octupolar physics, consistent with the newly constrained ground state.

II. SAMPLE SYNTHESIS AND QUALITY
Ground state selection in several pyrochlore magnets is known to be sensitive to disorder.Perhaps the best-known example is associated with Yb 2 Ti 2 O 7 [26][27][28][29][30][31] .Ce-based pyrochlores are likely candidates for such sensitivity because, as with other light rare-earth ions, Ce is stable in both the 3+ and 4+ oxidation states.Gaudet et al. showed that powder samples of Ce 2 Zr 2 O 7 grown by solid-state synthesis oxidize in air in a matter of hours 20 .Unlike the rest of the magnetic rare-earth series, however, this type of defect is relatively benign, as Ce 4+ is 4f 0 and consequently nonmagnetic.Nevertheless, a detailed characterization and understanding of Ce-based pyrochlores grown by different techniques is clearly a crucial open issue for this branch of topological materials physics.Therefore, we briefly outline the hydrothermal synthesis method and characterization of the sample quality.

A. Hydrothermal synthesis of Ce2Sn2O7
Single crystal and powder samples of Ce 2 Sn 2 O 7 were grown via a hydrothermal synthesis method 32 [Figure 2 (a) and (b)].The advantage of hydrothermal synthesis, from a sample quality standpoint, is the significantly lower temperature required (700 • C), as compared with solid-state synthesis (typically > 1000 • C), avoiding the temperatures at which tin oxide becomes volatile, thereby ideally minimizing oxidation.Approximately 13 g of the pyrochlore stannate powder was prepared for neutron measurements using a mixture of 30 − 50% SnO and 70 − 50% SnO 2 in order to reduce any residual Ce 4+ to Ce 3+ , thus minimizing external impurities.
We expect negligible stuffing of the A-site Ce 3+ on the Bsite Sn 4+ due to the large size difference (1.01 versus 0.69 Å), unlike Yb and Ti (0.87 versus 0.61 Å) in Yb 2 Ti 2 O 7 27 .If present, such stuffing would create a corresponding defect in the oxygen sublattice, leading to a noticeable change in color.As shown in Figure 2 (a) and (b), the single crystal and polycrystalline samples are bright yellow with a green tint, and we observe no evidence of oxidation effects over time while in air, in stark contrast to the sister compound Ce 2 Zr 2 O 7 which turns black on the order of hours.Furthermore, we have collected temperature-dependent, full-sphere single crystal x-ray datasets at 100, 200, and 300 K.The oxygen atom positions were refined on lower symmetry sites and we see no flattening or movement of the thermal anisotropic displacement parameters, suggesting an absence of defects in the oxygen sublattice at any measurable concentration.The refined lattice parameter of 10.6464(4) Å, within 0.01 Å of Tolla et al. 33 and Sibille et al. 19 , also suggests negligible oxidation.The remaining question is whether there is evidence of Ce 4+ stuffing on Sn 4+ sites, which is possible given the comparable ionic radii (0.87 versus 0.69 Å).Our x-ray refinement places an upper bound of 3% for such B-site stuffing of Ce 4+ .

B. Room temperature powder neutron diffraction (PND) and analysis
Room temperature powder neutron diffraction (PND) measurements were performed on beam lines HB-2A at the High Flux Isotope Reactor (HFIR) and NOMAD at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (ORNL) 34 .Room temperature PND on HB-2A [Figure 2 (c utilized a Ge(115) monochromator with an incident wavelength of 1.54 Å.The resulting powder neutron diffraction pattern was refined using the FULLPROF software suite which implements the Rietveld refinement method 35 .We find a lattice parameter of 10.64542(5) Å, in agreement with our single crystal x-ray diffraction.A small nonmagnetic impurity of CeO 2 is present, and makes up ∼ 3% of the powder material by weight.The total scattering measurement on NOMAD, with a Q max = 40 Å −1 , allows an atomic pair distribution function (PDF) analysis to be carried out to further constrain the local structure of the hydrothermally grown sample [Figure 2  (d)].The PDF analysis gives an atom-by-atom histogram of all of the pair-pair correlations in a material as a function of real-space distance r.The data were processed through the NOMAD autoreduction software and we utilized PDFGUI 36 to perform the PDF refinement.This analysis resulted in a similar lattice parameter as was found using HB-2A, and there was no evidence of local structure distortions.

III. ESTIMATING THE EXCHANGE PARAMETERS IN THE SPIN HAMILTONIAN FOR CE2SN2O7
A. Heat capacity and susceptibility Heat capacity measurements were performed on a 0.26 ± 0.03 mg single crystal of Ce 2 Sn 2 O 7 (Figure 3) from T = 1 K down to T = 0.05 K using a Quantum Design PPMS with dilution insert.Similar low temperature measurements were carried out on hydrothermally grown powder samples, with additional measurements taken above 1.9 K.The large uncertainty (∼ 10%) at low temperatures is attributed to the extremely small sample mass, a contribution that is typically ignored, and is determined by multiple successive weight measurements.Note that this uncertainty is associated with the overall scale of the heat capacity; the relative uncertainty, data point to data point at different temperatures, is much smaller.
Our results above 0.5 K are in good agreement with the literature 18,19 .We find no evidence for a transition to long-range FIG.3: Single crystal (T< 1 K) and powder (T> 1 K) heat capacity on hydrothermally synthesized Ce2Sn2O7 compared to previously published polycrystalline data by Sibille et al. 19 as well as Ce2Zr2O7 from Smith et al. 22 .Note that the data presented in this work are the total heat capacity, while the comparative works present the isolated magnetic contribution.However, the phonon contribution below 10 K is negligible and thus does not affect the conclusions of this work.We find good agreement with Sibille et al., and the shift of the high temperature tail to lower temperatures when compared to Ce2Zr2O7 suggests a smaller scale of the intersite interaction strengths.
order down to the lowest measured temperatures, and instead observe a broad peak in the heat capacity, the shape of which is similar to that of the sister compound Ce 2 Zr 2 O 7 21,22 .Importantly, we note the high temperature tail of the heat capacity peak is shifted to lower temperature compared to Ce 2 Zr 2 O 7 , indicating that Ce 2 Sn 2 O 7 hosts comparatively smaller intersite interaction strengths, as we will show in Section III C. Susceptibility measurements were also carried out using a Quantum Design MPMS3 with iHelium3 insert.Powder and single crystal samples were mounted to straws using grease and measured with an applied field of H = 0.01 T (field aligned along [111] for single crystals).Results are shown in Figure 7 (a).

B. Low temperature neutron diffraction
We collected low temperature PND on the WAND 2 beam line at HFIR (ORNL) using a Ge(113) monochromator for an incident wavelength of 1.488 Å (E i ∼ 37 meV), and a closedcycle refrigerator.These measurements on the hydrothermally grown Ce 2 Sn 2 O 7 samples were performed as a function of temperature from a base temperature of T = 0.3 K, well below the temperature at which octupolar correlations are expected ∼ 1 K, up to 10 K.The high flux and 3 He position-sensitive detector makes WAND 2 an ideal instrument for probing weak, diffuse neutron diffraction signals, as may be expected from octupolar scattering at high Q.The diffuse scattering signal is isolated by performing a high temperature subtraction of the paramagnetic background 10 K dataset.This protocol allowed us to put the net diffuse neutron scattering on an absolute scale, to compare with earlier diffuse scattering measurements per- formed on Ce 2 Sn 2 O 7 powder samples grown by solid-state synthesis.However, neither work is normalized using a vanadium standard and are therefore susceptible to uncertainties associated with absorption effects.This has no bearing on the difference in the Q dependence of the diffuse neutron scattering signals.
Intriguingly, our new neutron diffraction results show the diffuse scattering to peak at small |Q| ∼ 0.6 Å −1 , and no diffuse scattering above background was observed beyond ∼ 2 Å −1 [Figure 4 (a)].The temperature dependence of the diffuse scattering from the new hydrothermally grown Ce 2 Sn 2 O 7 powder sample is shown in Figure 4 (b).These data show that the low-Q diffuse scattering onsets near 1.5 K where we also see the initial increase in the heat capacity.This diffuse scattering was further corroborated by similar low temperature measurements on beam line HB-2A at T = 0.3 K and T = 4 K.
This new diffuse neutron scattering from Ce 2 Sn 2 O 7 strongly resembles that from the canonical dipolar spin ice pyrochlore Ho 2 Ti 2 O 7 .We compare the diffuse scattering results from these two pyrochlore powders in Figure 4 (a) (see Appendix B for experimental details).In both cases, this diffuse scattering is inelastic scattering that has been integrated in the diffraction measurement, as the incident energy is much larger than the bandwidth of inelastic scattering known to be relevant to Ce

C. Numerical linked cluster expansions and Monte Carlo calculations and their comparison to experiment
To understand the underlying behavior of Ce 2 Sn 2 O 7 , we invoke theoretical calculations that allow a parametrization of the spin Hamiltonian given in Eq. (1).We utilize a quantum numerical linked cluster expansion of the specific heat and susceptibility as well as classical, and wherever possible quantum, Monte Carlo simulations of the powder-averaged neutron structure factor to determine the coefficients of the nearestneighbor exchange interactions and the corresponding ground state, as well as finite-temperature properties of Ce 2 Sn 2 O 7 .At finite temperatures, NLCE provides an opportunity to calculate observables in the thermodynamic limit without parameter restrictions plagued by other numerical methods (see App. D), enabling parametrization throughout the whole DO pyrochlore phase space.
The procedure is similar to the one used to constrain Ce 2 Zr 2 O 7 to lie within the U(1) π QSL ground state 22 and is described in detail in Appendix D as well as briefly out- All data have been converted to absolute units, and paramagnetic background has been removed.We compare the Ce2Sn2O7 diffuse scattering to high-Q diffuse scattering found by Sibille et al. 19 from samples grown by solid-state synthesis.We do not see evidence of octupolar scattering at high Q, and instead find subtle diffuse scattering akin to a dipolar spin ice [see Figure 4 (a)].
lined here.First, we rewrite Eq. (1) in terms of J a , J b , and J c , where {a, b, c} are simply a permutation of {x, ỹ, z} such that We compute the specific heat for the full range of possible values of J b /J a and J c /J a , optimizing the energy scale, J a , in each case using the high temperature tail of the experimental data C(T ).An additional fitting parameter, is included in the fits to the heat capacity data to account for the systematic mass uncertainty from the extremely small samples.Using the estimated energy scales, we then compute the susceptibility χ(T ) and structure factor, S(|q|, T = 0.3 K) − S(|q|, T = 10 K), for each value of J a,b,c , each value of θ, and each permutation of exchange parameters ({x, ỹ, z} = {a, b, c}, {b, a, c}, {b, c, a}) (see Appendix D), optimizing the value of the g factor g z for each case.Finally, we minimize the total goodness of fit χ 2 to find the set of exchange parameters and permutation that best agrees with the experimental data.
Our best fit is achieved with the following exchange parameters in meV: This parametrization gives excellent agreement with the experimental heat capacity data down to 0.1 K [Figure 6  Importantly, these parameters place the ground state of Ce 2 Sn 2 O 7 within the ordered dipolar AIAO region of the phase diagram, proximal to the U(1) 0 QSL phase, as seen in Figure 8.This suggests a phase transition to AIAO order will occur at low temperatures, although we emphasize that we observe no direct evidence for this in experiment.
As our best fit estimates for the terms in the spin Hamiltonian place Ce 2 Sn 2 O 7 in a nonfrustrated region of the DO phase diagram, quantum Monte Carlo (QMC) calculations can be carried out without a sign problem using these parame-Experiment (rebinned) Theory 1.39 The measured Ce2Sn2O7 heat capacity, along with appropriate best fits to the NLCE theory up to order 7.The experimental heat capacity is corrected by the fitted parameter εc to account for the systematic mass uncertainty.Inset: the full QMC calculated heat capacity beyond experimentally accessible temperatures using parameters near the best fit parameters (see Appendix G).This finds a sharp anomaly signifying an expected phase transition to the AIAO state near 0.04 K. (b) The measured S(|q|) for hydrothermally grown Ce2Sn2O7 along with the theoretical calculation using best fit parameters.This shows both Cp and S(|q|) can be well described by calculations based on Eq. (1), and allow us to determine our estimate for the terms in the exchange Hamiltonian (except for θ).
ters.Our QMC simulation using the parameters in Eq. (2) predicts a phase transition temperature just within the temperature range of our experiments (see Appendix G).However, the transition temperature is very sensitive to small variations of the Hamiltonian parameters in this region of the phase diagram.For example, a small modification of the parameters to (J x, J ỹ , J z ) = (0.0457, −0.0014, −0.0110) meV produces a transition temperature below the experimental range and therefore consistent with experiment, without significantly reducing the agreement with other quantities.A comparison of the experimental data with the QMC results for the slightly modified parameter set is shown in the inset of Figure 6 (a).
It is worth noting that for the dipolar, classical spin ices, Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 , classical Monte Carlo simulations with realistic interaction parameters show a phase transition at very low temperatures from a classical, disordered dipolar spin ice state to an ordered spin ice state 38 .However, the temperature at which this predicted phase transition occurs is also very sensitive to the precise Hamiltonian parameters employed.This predicted phase transition has yet to be observed experimentally in either Ho 2 Ti 2 O 7 or Dy 2 Ti 2 O 7 , despite numerous attempts.

IV. DISCUSSION
While the protocol for determining Eq. ( 2) is already outlined, it is worthwhile to describe how this combination of heat capacity, magnetic susceptibility and diffuse neutron scattering measurements strongly constrains the best fit Hamiltonian parameters.First, the high temperature tail of the heat capacity constrains the overall scale of exchange interactions, while the height and shape of the smooth peak of the heat capacity constrain the ratios.In particular, a strongly negative value of J ± leads to a rather lower heat capacity peak, as observed in Ce 2 Zr 2 O 7 .On the other hand, a sufficiently large positive J ± leads to a sharp phase transition without any preceding smooth peak at higher temperature.The presence of smooth peak with height 2.0 J K −1 mol −1 , thus forces a relatively small value of |J ± |/J a .The shape of the curve on the approach to the peak, is further constraining, and narrows the region of good fits to the region around our optimal parameter set.
The permutation of parameters {J a , J b , J c } → {J x, J ỹ , J z } is not fixed by the heat capacity but can be fixed using neutron scattering and susceptibility data.In particular, the welldeveloped, spin-ice-like scattering at low Q could only be reproduced in our simulations by allowing the largest exchange parameter to be one of the dipolar couplings J x or J z .When all three datasets (heat capacity, neutron scattering, and suscep-  19,22,23,39 . Here the phase diagram is shown without the reduction of parameter space used in Figure 8 so that dipolar and octupolar phases are distinguished.Phase boundaries are based on the calculations from Benton 16 .The dashed lines represent the classical phase boundary of the AIAO order.The parameter estimate for Ce2Sn2O7 from the present work falls within the AIAO phase, in a region where the ordering is classically unstable but stabilized by quantum fluctuations for S = 1/2.tibility) are taken together, the fits are overconstrained rather than underconstrained, with a relatively small region of parameter space enabling acceptable agreement with all three.
With the NLCE calculations having established a dipolar AIAO ground state for Ce 2 Sn 2 O 7 , we interpret the measured diffuse neutron diffraction at T = 0.3 K as a spin ice regime above the AIAO transition.Further, it is interesting to note that the θ > 0 term in the best fit Hamiltonian implies spin waves should be visible in nonzero magnetic field.Inelastic neutron scattering measurements could therefore corroborate the fitted Hamiltonian parameters, but are not presently straightforward due to the small sizes of available single crystals.Nonetheless, this would be an important future step after sample size optimization is developed.
Figure 9 shows the ground state phase diagram of the pyrochlore XY Z Hamiltonian without the reduction of parameter space by permutations, in contrast to Figure 8 22,23 .This surprising degree of sensitivity of the exchange parameters to both the transition metal M and, at least potentially, to disorder, implies a high degree of tunability in Ce pyrochlores which could be used to further explore the full phase diagram.
The phase diagram in Figure 9 includes some regions where quantum fluctuations substantially renormalize the phase boundaries from what one would obtain in a classical treatment of Eq. ( 1).This is illustrated by the white dashed lines in Figure 9, which show the classical phase boundary of the AIAO antiferromagnetic phase.According to our parametrization, Ce 2 Sn 2 O 7 falls in a regime where, classically, one would expect a disordered spin ice ground state, but where quantum fluctuations stabilize antiferromagnetic AIAO order.This then rationalizes our observation of a dipolar spin ice regime at temperatures above an expected phase transition to AIAO order.
Ce 2 Sn 2 O 7 thus appears as an example of a material where strong quantum effects favor ordering over disorder.Such order by disorder is a known phenomenon in frustrated magnets, but Ce 2 Sn 2 O 7 presents a striking example both because it is far from the regime where quantum effects can be treated perturbatively and because the resulting order is not selected from among the classical manifold of ground states.

V. SUMMARY AND CONCLUSIONS
In summary, we report new thermodynamic, magnetic, and diffuse neutron scattering data from powder and small single crystals of Ce 2 Sn 2 O 7 grown by new hydrothermal techniques.Modeling of this experimental data using NLCE, classical Monte Carlo techniques and (where possible) quantum Monte Carlo techniques allow us to place Ce 2 Sn 2 O 7 in the region of the general XY Z DO phase diagram where a dipolar AIAO ground state is selected.
Diffuse neutron scattering measurements performed on the new samples of Ce 2 Sn 2 O 7 grown by hydrothermal techniques possess a similar Q dependence to that established in the canonical dipolar spin ice Ho 2 Ti 2 O 7 .These new measurements are distinctly different from the diffuse scattering observed at high Q and previously ascribed to octupolar correlations in Ce 2 Sn 2 O 7 samples grown from solid-state synthesis.The origin of the difference between these new diffuse neutron results and those reported earlier is not understood.On the one hand, we may expect different levels, and perhaps types, of weak disorder present in the same material synthesized using different techniques.Nevertheless, while the samples themselves could be different, the stark difference in low temperature diffuse neutron scattering would imply an extreme sensitivity to small levels of weak disorder.
Our results, however, are internally consistent.Modeling of our comprehensive experimental datasets places our new sample of Ce 2 Sn 2 O 7 in a region of the phase diagram where a dipolar AIAO ground state is selected, and the octupolar QSL states are not nearby.Our new diffuse neutron scattering measurements are consistent with this, and can themselves be described by this same theory.That is, for this dipolar AIAO region of the DO phase diagram, strong diffuse scattering at large Q is not expected.All of this supports our conclusion that for these new hydrothermally grown Ce 2 Sn 2 O 7 samples, we observe a proximate dipolar spin ice regime at finite temperature above a predicted dipolar AIAO ground state.0.5 0.5 0.5 0.0005(1) 0.0005(1) 0.0005(1) -0.00031(10) -0.00031(10) -0.00031(10) Sn (16c) 0 0 0 0.00019(10) 0.00019(10) 0.00019(10) -0.00002(9) -0.00002(9) -0.00002( 9) O (48f ) 0.3307(1) 0.125 0.125 0.00076 (10)  3% of the sample by weight.Importantly, the addition of an oxygen atom at the 8a crystallographic site (the oxygen vacancy site in a perfect pyrochlore) did not result in an improved fit and the occupancy refinement oscillated around zero.Given the resolution of the measurement, we can obtain an upper bound of ∼ 2% oxygen occupancy.Refined parameters are therefore found to be in agreement with PND from Sibille et al. 19 taken on the HRPT beam line at SINQ, Paul-Scherrer Institute.
Room temperature time-of-flight total scattering was performed on the NOMAD beam line, taking advantage of the mail-in program, in order to perform an atomic pair distribution function analysis.This analysis technique is sensitive to local (short-range) disorder effects and is complementary to traditional crystallographic analyses that probe the long-range average structure.
The returned data was carefully background subtracted and given in S(|q|) and G(r) for a variety of Q max cutoff values, processed by the NOMAD auto-reduction software.A Q max = 40 Å −1 was chosen to obtain sufficient real-space resolution while minimizing the high-frequency noise from the high-Q data.The PDFGUI software was used to perform least-squared refinements of the modeled data with fitted lattice parameters, thermal ADPs, and the 48f oxygen position.The refined pattern agrees well with the experimental data and fitted lattice parameter of 10.6476 Å, and the lack of spurious peaks at expected vacancy positions suggests no local distortions are present [Figure 2 (d)].The residual found is comparable to that found from fitting a silicone standard under the same instrumental conditions.

Appendix B: Ho2Ti2O7 measurements and Data treatment
Powder neutron diffraction was performed on canonical dipolar spin ice Ho 2 Ti 2 O 7 at the HB-2A beam line (HFIR, ORNL) using an incident wavelength of 1.54Å using the Ge(115) monochromator.Powder samples were prepared using the conventional solid-state synthesis method 40,41 and subsequently loaded into a standard aluminum can as a pressed pellet.Diffraction patterns were measured as a function of temperature from 0.3 to 50 K, both well below and above the temperature at which spin ice correlations emerge, around 2 K.
In both low temperature PND measurements of Ce 2 Sn 2 O 7 and Ho 2 Ti 2 O 7 , high temperature background subtraction was utilized to isolate the pure magnetic signal at low temperatures.In the case of Ce 2 Sn 2 O 7 , Sibille et al. find a temperatureindependent, featureless paramagnetic scattering above 2 K, justifying our choice of using 10 K as paramagnetic background.Ho 2 Ti 2 O 7 however, has been shown to posses shortrange correlations up to nearly 20 K 42 , and thus 50 K was chosen for the paramagnetic background.The results of this subtraction are shown in Fig. 4; and we note that the nuclear Bragg peaks have been masked for clarity in both materials.

Appendix C: Powder heat capacity measurements
For powder heat capacity measurements above 1.9 K, a pellet was pressed from the sample powder, while measurements below 1.9 K required a pellet made using an equal-mass combination of the sample powder and silver powder to enhance the thermal conductivity and sample coupling.The contribution of the silver to the specific heat was subtracted using in-house data for the specific heat capacity of silver.The relative contribution of silver decreases from 8% at 1.9 K, to 1% at 1 K, and drops below 0.1% at 0.4 K.After this careful background subtraction, the powder heat capacity data were found to be in good agreement with the single crystal measurements.

Appendix D: Numerical linked cluster expansion
Obtaining thermodynamic observables in three-dimensional frustrated materials, such as dipolar-octupolar pyrochlores, is extremely challenging as different numerical methods face different limitations.
First, the exponentially growing complexity limits exact diagonalization techniques such as finite-temperature Lanczos to small systems 43 .Second, finite-temperature methods based on tensor networks or matrix product states beyond one dimension fail to handle the entanglement and tend to develop biases toward ordered states 44,45 .Lastly, Monte Carlo approaches like stochastic series expansion (SSE), cf.Appendix G, are restricted to unfrustrated parameter regimes due to the infamous sign problem 46,47 .
At finite temperature, the numerical linked cluster expansion 48 opens up the possibility to calculate observables in the thermodynamic limit without any parameter restriction.Starting from a single tetrahedron, the algorithm generates clusters of finite size by building up the whole lattice in a systematic way.These clusters are then solved using exact diagonalization.The generic XY Z Hamiltonian in Eq. (1) does not preserve any continuous spin symmetry which makes numerical methods extremely challenging.Thus, we rely only on spatial symmetries which are determined and optimized in an automatic way 49 .The tetrahedral expansion of the pyrochlore lattice exhibits a large number of independent spatial symmetries.The region of convergence (T 0.1 K) is comparable to temperatures used in the experiment and therefore allows a systematic analysis in the full parameter regime.A similar procedure was followed in our previous analysis of Ce 2 Zr 2 O 7 22 .

Specific heat
Many thermodynamic observables like the specific heat or entropy depend only on the spectrum of the Hamiltonian.This can be exploited to reduce the possible parameter sets (J x , J y , J z ).First, the overall energy scale α can be neglected in numerical simulations as it simply adjusts the temperature unit of the specific heat curve C(T ).Adjusting the energy scale by α (E i → E i = αE i ) induces a shift of α −1 in the specific heat: Thus, the three-dimensional parameter space is reduced to two parameters.Second, (J x , J y , J z ) can be permuted by unitary rotations in spin space and the permutations do not affect the spectrum of the Hamiltonian for h = 0. Therefore, we can reduce the possible parameter sets to the characteristic triangular shape of possible exchange parameters shown in Figure 8.Given a specific parameter set, we calculate the specific heat curve in sixth order in a tetrahedral expansion including nontrivial hexagonal loops and compare it to the experiment in Figure 6 (a).The convergence depends strongly on the parameters.The NLCE data of the specific heat and the error are obtained by the largest two orders: The error associated with a given parameter set is determined as follows.First, the high temperature tail (1.9 K < T < 4 K) of the experimental data which exhibits a small error is used to fix the overall energy scale α such that the exchange parameters are scaled to Kelvin.Second, the remaining data points (T 0 < T < 1.9 K) are used to estimate an error for the specific parameter set for different T 0 = 0.1, 0.125, 0.15 K.The experiment involves a large systematic error δ sys (T ) at low temperature due to the uncertainty of the sample mass.To incorporate the systematic error, we introduce a fitting parameter ε 2 is minimized.The goodness-of-fit for the specific heat data is given by where δ unsys is the unsystematic error.The derived goodness of fit for each parameter set of the exchange constants is shown in Figure 10.The best approximation (red cross) defines all three exchange parameters (J a , J b , J c ).

Magnetic susceptibility
Determining the susceptibility in an external field is not as simple as for the specific heat.We generalized the tetrahedral NLCE from Schäfer et al. 49 to include the site-dependent Zeeman energies due to the direction of the applied magnetic field.The anisotropy requires calculation of all three independent permutations: In addition to the permutation, we performed the calculation for different values of θ ∈ [0, π/2).Since we vary θ which mixes S x and S z , we do not need to separately account for parameter permutations which swap J x and J z .We use the same overall energy scale α extracted from the high temperature tail of the specific heat (1.9 K < T < 4 K) such that the exchange parameters are in units of kelvin.Then, we apply a field of |h| = 0.01 T in the (1,1,1) direction such that the easy axis on one site of the tetrahedral unit cell is aligned with the field yielding h • z i = |h|, while the other sites yield h • z i = − 1 3 |h|.We scaled the external field by k −1 B to match the unit of the exchange parameters.We can evaluate the total magnetization using NLCE at finite temperature via where e is aligned with the field with |e| = 1.The susceptibility is defined by where c = 9.65 emu T 2 mol Ce −1 meV −1 converts the susceptibility into the correct units.The susceptibility measurements were performed at higher temperature than the measurements of the specific heat such that fourth order NLCE is sufficient.Since the NLCE data are completely converged for the region of interest, we define the susceptibility and its error the following way: To begin with we set g z = 2.57 in accordance with expectations for the CEF ground state wave function.We include a scaling factor, ε χ ∈ [0.5, 1.5], to allow a variable g factor for the effective system which minimizes the error: Note that ε χ scales the (simulated) NLCE data of the susceptibility while ε C scales the (experimental) systematic error of the specific heat.Since χ NLCE ∝ g 2 z , ε χ scales the final result for the g factor as g z → √ ε χ g z .we also rebin the experimental data into the same shells and define , (F9) where λ is an additional scale parameter that is minimized for each dataset as a consistency check and q min = 0.66 Å −1 and q max = 5 Å −1 .Since we are comparing absolute differential cross sections, λ for good fits should be close to one.Indeed we find that for the best parameter estimate [Eq.(2)], λ = 1.39.Note that the exact value of this scale factor λ is quite sensitive to the mixing angle θ.To reproduce the exact value of λ = 1.39 as a best fit, one has to choose the exact best fit value of θ = 0.185π ≈ 0.19π.We show the goodness of fit, optimized over the three inequivalent permutations of exchange parameters and all values of the mixing angle θ in Figure 11.Remarkably, for all good fits, the optimal permutation of exchange parameters is (J x, J ỹ , J z ) = (J a , J b , J c ), and the optimal value of the mixing angle is θ ≈ 0.2π.Also the optimal scale λ is of order one for all good fits.Hence, while the fit to the structure factor does constrain the exchange permutation and the value of the mixing angle quite reliably, it does not, by itself, strongly constrain the value of exchange parameters.FIG.
12: Comparison between specific heat of Ce2Sn2O7 and as computed for different parameter sets using quantum Monte Carlo simulations.The inset shows the goodness of fit obtained from SSE for parameters around the best high temperature fit.The white dashed line indicates the phase boundary between the ordered and the liquid phase 47 .Importantly, even when considering all experimentally accessible temperatures the best fits lie in the ordered phase.
Appendix G: Quantum Monte Carlo simulations Since the estimated parameters of Ce 2 Sn 2 O 7 [Eq.( 2)] fall into the region in parameter space that allows unbiased quantum Monte Carlo simulations, we implement the stochastic series expansion 52 to compute the specific heat for the best parameter fit down to lower temperatures than possible using NLCE.This also allows us to estimate the location of the phase transition that would be expected for these parameters at low temperature 16,47 .A quantitative comparison between experimental data and the SSE simulation is shown in Figure 12.Importantly, the best fit to the high temperature data [Eq.(2)] does not compare well to the experimental data at low temperatures because the sample would be expected to undergo a phase transition at the lowest measured temperatures already.However, due to the proximity to the phase boundary (indicated by a white dashed line in the inset of Figure 12), the transition temperature is expected to be very sensititve to the exchange parameters.This is exemplified already by considering the parameters (J x, J ỹ , J z ) = (0.0457, −0.0014, −0.0110) meV, (G1) which differ from the best high temperature fit by only ∼ 10 −3 meV, but this difference suffices to push the critical temperature the lowest experimentally measured temperature.Indeed, performing a SSE parameter sweep in the neighborhood of the best high temperature fit, we find the above parameters to be the best fit when considering the full temperature range of the experimental data.We also stress that even when considering this full temperature range, the best fits all lie within the ordered phase, although very close to the phase boundary to the liquid phase.We show the full goodness of fit χ 2 c of the specific heat when compared to the SSE simulation result for a limited parameter region in the inset of Figure 12, where and, as for the NLCE fit, c is fitted for each parameter set to minimize χ 2 c .

FIG. 1 :
FIG.1: Ground state phase diagram of the dipolar-octupolar exchange Hamiltonian.Phase boundaries are based on the calculations from Benton16 , which find four U(1) QSL phases, labeled QSL-0 or π (dip or oct), alongside dipolar and octupolar AIAO Néel phases.The dip or oct distinction refers to whether the dominant pseudospin component transforms as a dipole or an octupole, while the 0 or π refers to the flux through the hexagonal plaquettes within the pyrochlore structure.The strong potential for realizable QSL materials is a huge motivation to further explore the DO pyrochlores in detail.

FIG. 2 :
FIG. 2: Representative (a) single crystal and (b) powder Ce2Sn2O7 samples, both grown using hydrothermal techniques, used for heat capacity and magnetization measurements.Samples remained a bright yellow-green in air, with no visible signs of oxidation effects over time.(c) Rietveld refinement of room temperature powder neutron diffraction on beam line HB-2A.Data (blue points), calculated pattern (red line), difference pattern (black line) and Bragg positions of Ce2Sn2O7 (red tick marks) are shown.Refined parameters can be found in Table I.(d) Pair distribution function (PDF) analysis of Ce2Sn2O7 from total neutron scattering data taken on the NOMAD beam line at room temperature.The fitted curve (red) shows good agreement with the data (blue), finding no evidence of local distortions in the sample.

FIG. 4 :
FIG. 4: (a) High temperature subtracted powder diffuse neutron scattering of the known dipolar spin ice compound Ho2Ti2O7 compared with newly synthesized Ce2Sn2O7 in absolute units.Both materials show a broad diffuse peak at |Q| ∼ 0.6 Å −1 .(b) Temperature dependence of the diffuse scattering of Ce2Sn2O7 with paramagnetic background removed.The intensity increase near |Q| ∼ 0.6 Å −1 is seen to onset around 1.5 K.

2
Sn 2 O 7 and Ho 2 Ti 2 O 7 19,37 .The dipole moments associated with Ho 3+ in Ho 2 Ti 2 O 7 are about 8 times larger than those associated with Ce 3+ in Ce 2 Sn 2 O 7 ; hence the overall intensity of the diffuse dipolar spin ice diffraction pattern is much stronger in Ho 2 Ti 2 O 7 compared with the new Ce 2 Sn 2 O 7 sample.Nonetheless, the qualitative similarity between the two sets of diffuse scattering is clear.The diffuse scattering of both Ce 2 Sn 2 O 7 samples grown by hydrothermal and solid-state synthesis techniques is directly compared in Figure 5.The two works strongly contrast, with the earlier Ce 2 Sn 2 O 7 measurements showing a much increased diffuse scattering beyond |Q| ∼ 3 Å −1 , peaking at |Q| ∼ 8 Å −1 with nearly a factor of 4 times the intensity than observed at low Q.That is, it is hard to imagine that the diffuse scattering at low temperature could be more different between the new hydrothermally grown Ce 2 Sn 2 O 7 powder and the previously studied Ce 2 Sn 2 O 7 powder sample grown by solid-state synthesis.

FIG. 5 :
FIG.5: Diffuse neutron scattering from hydrothermally grown Ce2Sn2O7 at 0.3 K taken on WAND 2 is shown as the blue circles.All data have been converted to absolute units, and paramagnetic background has been removed.We compare the Ce2Sn2O7 diffuse scattering to high-Q diffuse scattering found by Sibille et al.19 from samples grown by solid-state synthesis.We do not see evidence of octupolar scattering at high Q, and instead find subtle diffuse scattering akin to a dipolar spin ice [see Figure4 (a)].
(a)], including the fit to the systematic mass uncertainty ε c .The fit to the magnetic susceptibility strongly constrains the value of mixing angle θ.While two permutations ({a, b, c} and {b, c, a}) result in minima in the goodness of fit [Figure 7 (b)], only the {a, b, c} permutation in Eq. (2) produced a reasonable S(|q|) fit [Figure 6 (b)].
FIG. 6: (a)The measured Ce2Sn2O7 heat capacity, along with appropriate best fits to the NLCE theory up to order 7.The experimental heat capacity is corrected by the fitted parameter εc to account for the systematic mass uncertainty.Inset: the full QMC calculated heat capacity beyond experimentally accessible temperatures using parameters near the best fit parameters (see Appendix G).This finds a sharp anomaly signifying an expected phase transition to the AIAO state near 0.04 K. (b) The measured S(|q|) for hydrothermally grown Ce2Sn2O7 along with the theoretical calculation using best fit parameters.This shows both Cp and S(|q|) can be well described by calculations based on Eq. (1), and allow us to determine our estimate for the terms in the exchange Hamiltonian (except for θ).

FIG. 7 :FIG. 8 :
FIG. 7: (a) Measured magnetic susceptibility χ(T ) with an applied field strength of 0.01 T, and NLCE fits that establish θ ∼ 0.2π in the {a, b, c} permutation.(b) The goodness of fit to the susceptibility data is shown as a function of θ, and we find a strong minimum (i.e.strong constraint) near θ ∼ 0.2π.While a better goodness of fit is found for the {b, c, a} permutation (which would constrain θ ∼ 0.4π), this permutation does not result in a reasonable S(|q|) fit.

4 FIG. 10 :
FIG.10:The of fit from comparing the calculated specific heat to the experimental data.The best parameter estimate is marked as a red cross.

FIG. 11 :
FIG.11: Right: Goodness of fit from comparing the calculated structure factor S(|q|) from powder of Ce2Sn2O7 to the experimental data.The best fit from fitting both structure factor and specific heat is also indicated.For each pair (J±, J±±) we show the goodness of fit optimized over both the three permutations of exchange parameters and mixing angles θ.Center (right): value of θ (exchange permutation) with best goodness of fit for each pair (J±, J±±).

TABLE I :
Refined room temperature parameters for Ce2Sn2O7 from neutron diffraction on HB-2A using the FULLPROF Rietveld refinement.