Role of defects in determining the magnetic ground state of ytterbium titanate.

Pyrochlore systems are ideally suited to the exploration of geometrical frustration in three dimensions, and their rich phenomenology encompasses topological order and fractional excitations. Classical spin ices provide the first context in which it is possible to control emergent magnetic monopoles, and anisotropic exchange leads to even richer behaviour associated with large quantum fluctuations. Whether the magnetic ground state of Yb2Ti2O7 is a quantum spin liquid or a ferromagnetic phase induced by a Higgs transition appears to be sample dependent. Here we have determined the role of structural defects on the magnetic ground state via the diffuse scattering of neutrons. We find that oxygen vacancies stabilise the spin liquid phase and the stuffing of Ti sites by Yb suppresses it. Samples in which the oxygen vacancies have been eliminated by annealing in oxygen exhibit a transition to a ferromagnetic phase, and this is the true magnetic ground state.

Introduction

Yb2Ti2O7 has attracted intense experimental and theoretical interest as a model geometrically frustrated magnet with large quantum fluctuations[1-15]. The crystalline electric field gives a Kramers doublet ground state well separated from the excited states[4], and symmetry considerations lead to an effective spin ½ Hamiltonian with four exchange constants[5]. In zero applied magnetic field there is a continuum of spin excitations as a result of quantum fluctuations[6,14,15]. In a magnetic field above B ~ 0.5 T well-defined spin waves are recovered, and it is possible to determine the exchange constants via inelastic neutron scattering[1]. It was originally proposed that there is no long-range magnetic order down to the lowest temperatures[7,8], and that the dominant exchange is a ferromagnetic Ising-like term, with a large anisotropic term, leading to the intriguing proposal of a quantum spin ice[1]. Alternatively, it has been reported that the magnetic ground state is ferromagnetic[2,9-13], and the values of the anisotropic g-tensor show that the moments are more easy-plane than easy-axis[16].

The identification of different magnetic ground states for nominally stoichiometric samples clearly points towards the importance of low levels of structural disorder. This sample dependence is strongly reinforced in measurements of the heat capacity, where the existence of a sharp peak in the millikelvin range varies dramatically for different powder and single-crystal samples[11,17,18]. To date studies of structural disorder in Yb2Ti2O7 have focussed on the effect of stuffing Yb3+ ions on Ti4+ sites[19,20]. Recent studies have investigated the effect of pressure on the ground state of nominally stoichiometric and stuffed samples[21]. The fact that pressure appears to stabilise the ferromagnetic phase has been related to the chemical pressure associated with stuffing.

We have determined the defect structures in a range of oxygen-depleted, stuffed and nominally stoichiometric single crystals using diffuse neutron scattering, which is particularly sensitive to vacancies and displacements of oxygen ions. By correlating our results with magnetic diffuse scattering we are able to determine the effect of structural disorder on the spin correlations, and to unambiguously identify the true magnetic ground state of Yb2Ti2O7.

Results

Defect structures

All of our single-crystal X-ray diffraction data for stoichiometric Yb2Ti2O7, oxygen-depleted Yb2Ti2O7− and stuffed Yb2(Ti2−Yb)O7− refine in the cubic pyrochlore structure, space group Fd-3m, see Table 1. The Yb and Ti ions are located on pyrochlore lattices, and there are two inequivalent O sites: O(1) located at the centre of the Yb tetrahedra and O(2) forming a corner-sharing lattice of octahedra surrounding Ti ions. The presence of high atomic number elements means that the sensitivity to the oxygen ions is limited, and it is not possible to distinguish between O(1) and O(2) vacancies. However, we note that oxygen depletion does not lead to an increase in lattice parameter, and this is the first difference between Yb2Ti2O7− and the related oxygen-deficient pyrochlore Y2Ti2O7− (ref. [22]). In the case of Y2Ti2O7−, the increase in the lattice parameter was associated with the expansion of the Y tetrahedra around O(1) vacancies, and this is a clear indication that the defect structure of Yb2Ti2O7− is different.

Table 1

Refinement of the average structures of stuffed Yb2(Ti2−Yb)O7−, oxygen-depleted Yb2Ti2O7−, as-grown and oxygen-annealed Yb2Ti2O7 from the single-crystal X-ray structure factors measured at T ~ 300 K

	Stuffed	Depleted	As-grown	Annealed
Colour	Yellow/Brown	Black	Brown	Transparent
Space group	Fd-3m	Fd-3m	Fd-3m	Fd-3m
Lattice parameter	10.0838(14)	10.0216(16)	10.0271(2)	10.0313(17)
Yb1	1	1	1	1
Yb2	0.26(3)	0	0	0
Ti	0.74(3)	1	1	1
O(1)	1	1	1	1
O(2)	0.96	0.96(3)	1	1
x	0.3397(14)	0.3300(6)	0.3308(2)	0.3309(10)
R	4.01	1.92	4.23	2.60
R W	4.90	2.60	4.87	2.78

The fractional coordinates are given within the second origin choice of the Fd-3m space group: Yb1 (0.5,0.5,0.5), Yb2 (0,0,0), Ti (0,0,0), O(1) (0.375,0.375,0.375), O(2) (x,0.125,0.125). The goodness of fit is described by the crystallographic R-factor and the weighted R-factor, RW

The diffuse neutron scattering from oxygen-depleted Yb2Ti2O7− is very sensitive to departures from ideal stoichiometry. Figure 1 presents the defect structure, diffuse scattering in the (hk7) plane, and the calculated scattering from Yb2Ti2O7−. The diffuse scattering pattern is very different to the comparable data set from Y2Ti2O7−, which was shown to arise from the defect cluster surrounding O(1) vacancies[22]. We are able to reproduce qualitatively the main features of the diffuse scattering using a Monte Carlo simulation, where isolated O(2) vacancies and charge compensating Ti3+ ions replace nearest-neighbour Ti4+ ions[23,24]. The presence of Ti3+ ions results in displacement of nearby O(2) ions such that the Ti3+–O(2) bond length, 2.03 Å, agrees with the values for Y2Ti2O7− (ref. [22]) and Ti2O3 (ref. [25]). Thus, our single crystals grown from the melt using the optical floating zone technique in a reducing atmosphere have oxygen vacancies on the O(2) sites. Chemically, ytterbium exhibits polyvalency, existing in both Yb3+ and Yb2+ valence states, unlike yttrium which can exist only as Y3+. This is the origin of the dramatic difference in the nature of the defect structures. We note that a powder neutron diffraction study of Yb2Ti2O7 reduced at low temperature in a topotactic reaction with CaH2 also finds O(2) vacancies[26].

Fig. 1

The structure of defects in Yb2Ti2O7. a For oxygen-depleted Yb2Ti2O7− the oxygen vacancies (gold) are located on the O(2) sites (blue), and neighbouring Ti4+ (green) are replaced by charge compensating Ti3+ (purple) ions that are slightly displaced away from the vacancy. b For stuffed Yb2(Ti2−Yb)O7− neighbouring Ti4+ ions are replaced by Yb3+ (red) ions that are displaced towards the O(2) vacancy. The diffuse neutron scattering measured at T ~ 30 K in the (hk7) plane in c for oxygen depleted and d for stuffed samples is compared with the simulated scattering in e and f, calculated using the models shown in a and b, respectively. The contour diagrams are on a hot scale with an arbitrary maximum

Figure 1 also shows that the diffuse neutron scattering from stuffed Yb2(Ti2−Yb)O7− in the (hk7) plane is different. It is again possible to reproduce the scattering using a Monte Carlo simulation for O(2) vacancies, but with neighbouring Ti4+ ions replaced by Yb3+ ions, and with different displacements of surrounding ions. The Yb3+–O(2) bond length is larger than the Ti3+–O(2) bond length, as expected for the larger ionic radius of Yb3+. Recently, calculations for Yb2(Ti2−Yb)O7− using density functional theory have confirmed that O(2) vacancies are energetically stable relative to O(1) vacancies[27].

The diffuse scattering from the nominally stoichiometric sample presented in Fig. 2a resembles the scattering from the oxygen-depleted sample, but is much weaker. Annealing in oxygen removes most of the intensity of the diffuse scattering, see Fig. 2b. The 1D cuts in Fig. 2c directly compare the diffuse scattering as a function of oxygen stoichiometry. We note that the diffuse scattering remaining after annealing in oxygen in the vicinity of the {337} and {777} reflections closely resembles the scattering observed for oxygen-annealed Y2Ti2O7 (ref. [22]). We are able to demonstrate that this weak scattering arises from the inelastic scattering of phonons by examining the scattering observed in a single detector bank at a single orientation using the time-of-flight technique. Our first principles calculation of the phonon dispersion from Yb2Ti2O7 using the CASTEP code[28], and the diffuse scattering arising from one-phonon excitations in our neutron single-crystal Laue diffraction using the approach described in ref. [29], is in excellent agreement with the arc of scattering emerging from the (337) reflection, see Supplementary Figure 1. This demonstrates that annealing in oxygen completely eliminates the structural component of the diffuse scattering for Yb2Ti2O7, and proves that the dominant defects in our nominally stoichiometric samples are O(2) vacancies.

Fig. 2

Structural diffuse scattering in nominally stoichiometric Yb2Ti2O7. a The diffuse neutron scattering from the as-grown, nominally stoichiometric sample measured at T ~ 30 K in the (hk7) plane contains weaker structural scattering resembling the scattering from the oxygen-deficient sample. b The same crystal measured after annealing in oxygen, showing that most of the diffuse scattering has been eliminated. The remaining scattering is attributed to inelastic scattering from phonons, see Supplementary Figure 1. c 1D cuts along [h,4,6.8] for the oxygen-depleted, as-grown and oxygen-annealed samples. This cut was chosen to minimise inelastic scattering. The broad structural diffuse features at h ~ ±4 are clearly present for the oxygen-depleted and as-grown samples, but are greatly reduced for the oxygen-annealed sample

Magnetic ground state

The magnetic diffuse scattering from the nominally stoichiometric sample in the (hhl) plane at T ~ 50 mK shown in Fig. 3a closely resembles the features reported previously in the spin liquid phase[2,30-33]. The same crystal was measured under identical conditions after annealing in oxygen, and the diffuse magnetic scattering completely disappeared (Fig. 3b). This can be seen more clearly in the corresponding 1D cuts through these data sets in Fig. 3c, d. The only variation in intensity away from Bragg reflections in the case of the oxygen-annealed sample can be explained by attenuation in the copper mounting plate, see Supplementary Figure 2 for further details. Polarisation analysis was employed throughout, and the flipping ratio was measured for the (111) Bragg reflection. For the oxygen-annealed sample the flipping ratio drops to unity below TC ~ 425 mK, indicating a transition to a ferromagnetic phase, see Fig. 3e. In zero field, the formation of ferromagnetic domains below TC results in depolarisation of the beam and, therefore, in Fig. 3 we show only the total scattering.

Fig. 3

Spin correlations in nominally stoichiometric Yb2Ti2O7. a The unpolarised neutron scattering from the as-grown sample measured at T ~ 50 mK in the (hhl) plane, showing the characteristic rods of intensity along <111> directions. b Measurements under identical conditions from the same sample after annealing in oxygen shows the complete elimination of diffuse magnetic scattering. 1D cuts along [0.5,0.5,l] for c the as-grown and d the oxygen-annealed samples. The diffuse magnetic scattering from −1 < L < 1 observed for the as-grown sample completely disappeared for the oxygen-annealed sample. a and b are plotted with different scale bars because of different background levels, see c and d, from different volumes of copper in the beam. The dips in the scattering at L = ±1.8 in d result from the grazing-incidence/grazing-exit absorption by the copper plate, see Supplementary Figure 2c. e Temperature dependence of the integrated intensities of the (113) Bragg reflection normalised to unity at T ~ 1 K (blue), and the flipping ratio measurements (red), show that the oxygen-annealed sample enters a ferromagnetic phase below TC ~ 425 mK. Error bars are the standard deviations derived using Poisson statistics. f Magnetic scattering intensities obtained by subtracting the integrated intensities at T ~ 1.2 K from those at T ~ 50 mK. The solid lines show the comparison with the data for the collinear ferromagnet[9], the two-in two-out splayed ferromagnet[6] and the all-in all-out splayed ferromagnet[10]. The magnetic intensity at (220) rules out the collinear ferromagnet, and the absence of magnetic intensity at (002) is inconsistent with the ice-like splayed ferromagnet

The temperature dependence of the integrated intensity at the (113) Bragg reflection normalised to its intensity at T ~ 1 K is also shown in Fig. 3e. There is a clear increase in intensity indicating the onset of long-range ferromagnetic ordering below TC ~ 425 mK. We note that a single crystal of Yb2Ti2O7 cut from the same boule and annealed in oxygen under identical conditions exhibits a sharp peak in its heat capacity at TC ~ 214 mK[13]. This is consistent with the purest samples in the literature[2,19]. The different TC from neutron diffraction can readily be explained by the large hysteresis observed previously for this transition[6,9] or possibly the accidental inclusion of low-energy fluctuations. The low-temperature magnetic intensities are consistent with the results of ref. [10] where the moments have a ferromagnetic component along <100>, but are splayed towards <111> directions so that the components perpendicular to the local easy axes are of the all-in all-out type. Other nominally stoichiometric samples have adopted different collinear[9], nearly collinear[2], and two-in two-out splayed[6,12] ferromagnetic structures. The presence of strong magnetic intensity at (220) is inconsistent with the collinear ferromagnetic structure and, coupled with the absence of intensity at (002), this also rules out the two-in two-out ferromagnetic structure, see Fig. 3f.

Spin correlations

For the as-grown sample there is no depolarisation of the beam at any temperature. We present the spin-flip scattering in the (hhl) plane at T ~ 50 mK in Fig. 4a. This is a particularly clean measurement of the spin correlations, since there is no structural scattering in this polarisation channel. Comparison with the diffuse scattering intensity calculated using classical Monte Carlo and the exchange constants from ref. [30] at T ~ 450 mK (above TC) shows excellent agreement with the data (Fig. 4b). In particular, the rods of diffuse scattering along the [111] directions and the broad features near (220) and (004) are reproduced well in the calculations. For this model below TC the system orders in a ferromagnetic phase in agreement with previous MC simulations[30,34,35]. These exchange parameters place the system close to the boundary with the antiferromagnetic ψ3 phase, which has vanishing (002) scattering and strong (220) intensity, in agreement with our experiment. We note that this model is also consistent with the inelastic neutron scattering data from Yb2Ti2O7[13,30].

Fig. 4

Spin-flip scattering as a function of composition. a The spin-flip scattering from the as-grown sample measured at T ~ 50 mK in the (hhl) plane exhibits the characteristic features associated with the spin liquid phase, including the rods along [111] and broad scattering near (220) and (004). b The scattering calculated using classical Monte Carlo at T ~ 450 mK and the exchange model of ref. [30]. c The oxygen-depleted sample at T ~ 50 mK is similar to the scattering from the as-grown sample in a. d In contrast, these features are washed out for the stuffed sample at T ~ 50 mK, and instead the very broad diffuse scattering suggests mainly uncorrelated spins. The intensities are normalised to the size of crystal after background subtraction

For the oxygen-annealed sample above TC the neutron beam remains polarised. The magnetic scattering in the spin-flip channel resembles the as-grown sample at a comparable temperature, see Supplementary Figure 3. In the case of the oxygen-depleted sample, there is no depolarisation of the beam down to base temperature, and Fig. 4c shows that the diffuse scattering in the (hhl) plane at T ~ 50 mK closely resembles the as-grown sample. Hence the presence of isolated O(2) vacancies appears to suppress the transition to the ferromagnetic phase and preserve the spin liquid behaviour down to the lowest temperature studied. The transition to the ferromagnetic phase is also suppressed for the stuffed sample. However, in this case the spin correlations are qualitatively different, see Fig. 4d. The rods of scattering along [111] directions are heavily suppressed, and the uniform diffuse scattering with a magnetic form factor instead points towards mainly uncorrelated spins. Hence the presence of additional magnetic Yb3+ ions on Ti4+ sites leads to additional exchange pathways and this suppresses the characteristic spin correlations associated with Yb2Ti2O7 in the spin liquid phase.

Discussion

The observation of rods of scattering has been interpreted as a consequence of the proximity of competing antiferromagnetic ψ2 and ψ3 states[34,35]. It is interesting to observe that, at the classical Monte Carlo level, one of the most discriminating features between the different sets of exchange parameters proposed in the literature is the (220) scattering intensity, see Supplementary Figure 4. This is an antiferromagnetic feature that could be interpreted as either all-in all-out correlations in the transverse components of the incipient splayed ferromagnetic order, or as a competition between a collinear ferromagnet and the ψ2 or ψ3 antiferromagnetic phases. A comparison between quantum and classical numerical methods applied to Yb2Ti2O7 was presented in ref. [34], and the numerical linked cluster (NLC) expansion results differ from the classical Monte Carlo ones most strikingly by the change in (220) intensity. It would be interesting to determine the origin of the (220) intensity in NLC—whether it is quantum effects or the perturbative approximation favouring antiferromagnetic correlations—as it could become a smoking gun for quantum mechanical effects in Yb2Ti2O7.

The sensitivity of the magnetic ground state to modest hydrostatic pressures is also consistent with the proposed proximity to a phase boundary. It was recently reported that a nominally stoichiometric sample that did not order long range at ambient pressure, adopted a two-in two-out splayed ferromagnetic structure at P ~ 11 kbar[21]. It was suggested that the stabilisation of the ferromagnetic phase through the application of hydrostatic pressure implies that the observation of a ferromagnetic ground state arises when stuffing leads to chemical pressure. It is then surprising that the deliberately stuffed sample does not exhibit a ferromagnetic ground state at any pressure. We demonstrate that stuffing has the opposite effect of washing out the characteristic spin correlations in the spin liquid phase.

In summary, we have shown that low levels of intrinsic defects have a crucial role in determining the magnetic ground state of geometrically frustrated systems with enhanced quantum fluctuations. Single-crystal diffuse neutron scattering is an extremely sensitive probe of defect structures, and should be employed more widely in studies of model magnetic systems, particularly candidate quantum spin liquids. In the case of Yb2Ti2O7 we identify the introduction of isolated O(2) vacancies as a mechanism to stabilise the spin liquid, and annealing in oxygen as a means to obtain ultra-pure samples. Understanding the gapless continuum of quantum excitations at low temperature in the presence of competing classical ordered phases remains an outstanding challenge[15]. More generally, understanding defect structures in pyrochlores is becoming increasingly important since they can lead to spin glasses[36], topological spin glasses[37], and they may induce entirely novel quantum spin liquid phases[38-40].

Methods

Sample preparation

A single crystal of Yb2Ti2O7 was grown at the Clarendon Laboratory by the floating zone technique[41], and this as-grown nominally stoichiometric sample was dark brown. After studying the structure and magnetism, this crystal was annealed in O2 at a flow rate of 50 ml min−1 and a temperature of 1200 °C for 2 days, producing a transparent sample. Thermogravimetric analysis reveals a change in the oxygen stoichiometry of 0.38 ± 0.14%. The oxygen-depleted sample was grown and annealed in flowing mixed gas of hydrogen and argon, and this lead to a black crystal. The stuffed sample was grown by increasing the ratio of Yb2O3 to TiO2 in the starting materials, and in this case the crystal was yellow/brown.

Structure determination

The average structures were determined by single-crystal X-ray diffraction using a Molybdenum source Oxford Diffraction diffractometer at Royal Holloway. A large CCD detector captures full reciprocal space coverage to a real space resolution of 0.6 Å. 3D profile analysis of each reflection is performed using the CrysAlisPro software[42] and refinements of the average structure, including anisotropic thermal parameters, are performed using the Jana2006 software[43]. All refinements had R-factors less than 5%, exhibiting excellent fits to the data. Full details are available in the CIF files. The defect structures were studied by diffuse neutron scattering using the single-crystal diffractometer (SXD) at the ISIS pulsed neutron source at the Rutherford Appleton Laboratory. SXD combines the white-beam Laue technique with area detectors covering a solid angle of 2π steradians, allowing comprehensive diffraction and diffuse scattering data sets to be collected. Samples were mounted on aluminium pins and cooled to 30 K using a closed-cycle helium refrigerator in order to minimise the phonon contribution to the diffuse scattering. A typical data set required four orientations to be collected for 4 h per orientation. Data were corrected for incident flux using a null scattering V/Nb sphere. These data were then combined to a volume of reciprocal space and sliced to obtain single planar and linear cuts.

For the Monte Carlo code used to simulate the structural diffuse scattering, a crystal comprising 64 × 64 × 64 unit cells was generated from the average structure obtained from the refinement of the diffraction data. From a statistical perspective, the use of a large supercell helps us to average over the disorder in the system (self-average), and suppresses the background noise. O(2) ions are removed at random until we obtain the depletion concentration. Large displacements of neighbouring Ti3+ or Yb3+ ions are introduced by hand for the oxygen-depleted and stuffed samples, respectively. The distortion of the surrounding lattice is simulated using the balls and springs model in which hard spheres are connected to neighbouring ions by springs, and the simulation randomly displaces ions in order to minimise the elastic energy[23,24]. We were not able to reproduce the observed diffuse scattering using simulations with other point defects, such as O(1) vacancies[22]. We were further able to rule out antiphase domain boundaries[44], since all of the pyrochlore Bragg reflections were sharp, and static disorder from the flexibility of the corner-sharing network[45], since structural diffuse scattering is absent in our stoichiometric sample.

Lattice dynamics

The phonon dispersion of Y2Ti2O7 was measured using the MERLIN spectrometer at the ISIS pulsed neutron source at the Rutherford Appleton Laboratory[46]. An 8 g single crystal was mounted on an aluminium plate with the (hhl) horizontal scattering plane and the sample was cooled to T ~ 30 K using a closed-cycle cryostat. MERLIN was set up using the G chopper running at 450 Hz, giving incident energies Ei ~ 120, 52, and 29 meV. The sample was rotated over 180° in 0.5° steps with data collected at each angle. Each data file was reduced using Mantid[47] and files were combined using Horace[48].

Calculations of the titanate phonon dispersion have previously been shown to be extremely challenging[49]. However, the zirconate pyrochlores avoid many of these issues[50] so to produce the computational phonon dispersion of Yb2Ti2O7 and Y2Ti2O7 we have initially computed the dispersion for La2Zr2O7 (LZO) and mass substituted the La and Zr sites. The LZO dispersion was computed using density functional theory within the plane-wave pseudopotential approach as implemented in the CASTEP code[28]. The local density approximation was used with the default CASTEP ultrasoft pseudopotentials (version C9)[51]. The primitive electronic Brillouin zone was sampled with a Monkhorst-Pack grid[52] of 4 by 4 by 4 points with a plane-wave cutoff of 600 eV. The lattice and atomic positions were relaxed such that the residual stresses and forces were less than 0.01 GPa and 0.005 eV/Å, respectively, with the quasi-Newton method[53]. The phonon dispersion was then computed using the finite-displacement, supercell method[54] with a single cubic unit cell as the supercell. Mass substitution was performed as a post processing step using the CASTEP utilities. The time-of-flight one-phonon neutron scattering intensity was calculated for the particular sample orientation and scattering geometry on SXD[29]. Note that the inelastic scattering is calculated for energies above 1 meV in order to avoid singularities and, therefore, the calculations do not work precisely at Bragg reflections. For the heavily defective samples, the inelastic scattering is negligible in comparison to the structural diffuse scattering. However, for the pristine, oxygen-annealed sample, the inelastic scattering can be clearly identified.

Magnetic scattering experiment

The magnetic diffuse scattering was measured using the D7 diffuse scattering spectrometer at the Institut Laue-Langevin in Grenoble. A pyrolytic graphite monochromator was employed to select cold neutrons of wavelength λ ~ 4.8 Å, and a beryllium filter was used to remove higher harmonics. Single crystals were mounted on a copper base plate inside a cryomagnet with the [00l] axis perpendicular to the plane of the plate, rotated so that the (hhl) crystallographic plane coincided with the horizontal scattering plane. A dilution refrigerator was used for temperature control. The instrument was operated in the diffraction configuration which integrates the dynamical response and, therefore, this technique measures the instantaneous correlations. Uniaxial polarisation analysis was employed in order to separate the magnetic and structural correlations. The polarisation direction was fixed normal to the scattering plane and this coincides with the [1-10] crystallographic direction. A polarising super-mirror (bender) and a Mezei flipper were inserted on the incoming neutron beam to select neutrons with a given spin. After the sample, the final polarisation was analysed using an array of polarising benders in front of the helium detectors. For each scattering vector, Q, the intensity was measured in both non-spin-flip (structural + magnetic correlations) and spin-flip (magnetic correlations only) channels. In the case of a ferromagnet in zero field, the formation of domains results in the depolarisation of the neutron beam, but in this case the non-spin-flip and spin-flip intensities can be combined to give the total (unpolarised) cross section. The integrated intensities of the Bragg reflections were obtained by summing the data above the background level in the vicinity of the reciprocal lattice points, and these data were normalised to data integrated in the same way well above TC.

Magnetic scattering theory

In our work, we chose to focus on classical simulations of the spin system to compare to experiments. This choice was made for two reasons: (i) tuning the parameters in these simulations, one can obtain an excellent quantitative agreement with experiments; (ii) the ability to include quantum fluctuations is currently limited to the mean field random phase approximation and the perturbative NLC expansion techniques, where it may be unclear whether the differences with respect to the classical results are due to quantum effects or to the approximations involved. We performed classical Monte Carlo simulations of nearest-neighbour exchange Yb2Ti2O7 to compare with the experimental neutron scattering results. We neglected dipolar interactions as they are about 100 times weaker than the exchange coupling, and we are interested primarily in temperatures above the transition. We restricted our modelling to the crystal-field ground-state doublet, since neither temperatures nor neutron energies in the experiments are high enough to involve states above the doublet gap (~ 900 K[4,16]). The resulting model comprises classical Heisenberg spins S = ½ on the sites of a pyrochlore lattice, interacting via the generic exchange couplings J1, J2, J3, J4 allowed by the lattice symmetries[1,33]. Using single spin-flip updates, which are sufficient to ensure thermalisation above the critical temperature, and accounting for the single-ion anisotropy via an appropriate g-tensor, we computed the neutron structure factor of the system in the (hhk) plane. The SF and NSF components were computed with respect to the relevant neutron polarisation direction [1,−1,0]. We note that the numerical results can equivalently be interpreted as the thermodynamic average of instantaneous (t = 0) neutron scattering structure factor, or as the time-integrated (ω = 0) structure factor within the crystal-field ground-state doublet. We used the canonical 16-spin cubic unit cell, and system sizes L = 20 (main paper) and L = 16 (Supplementary Figure 4). At every temperature step, the autocorrelation function was allowed to drop twice to 0.01 to ensure equilibration before lowering the temperature further. Upon reaching the target temperature, an instance of the structure factor was computed, and the entire process was then repeated from high temperature (with a different random number seed) to generate statistically independent instances. All the figures in the paper were averaged over 1000 instances.