Quantum Excitations in Quantum Spin Ice

Recent work has highlighted remarkable effects of classical thermal fluctuations in the dipolar spin ice compounds, such as"artificial magnetostatics", manifesting as Coulombic power-law spin correlations and particles behaving as diffusive"magnetic monopoles". In this paper, we address quantum spin ice, giving a unifying framework for the study of magnetism of a large class of magnetic compounds with the pyrochlore structure, and in particular discuss Yb2Ti2O7 and extract its full set of Hamiltonian parameters from high field inelastic neutron scattering experiments. We show that fluctuations in Yb2Ti2O7 are strong, and that the Hamiltonian may support a Coulombic"Quantum Spin Liquid"ground state in low field and host an unusual quantum critical point at larger fields. This appears consistent with puzzling features in prior experiments on Yb2Ti2O7. Thus Yb2Ti2O7 is the first quantum spin liquid candidate in which the Hamiltonian is quantitatively known.

duced to that of an effective spin S = 1/2 moment, with the strongest possible quantum effects expected. In this case symmetry considerations reduce the exchange constant phase space of the nearest neighbour exchange Hamiltonian to a maximum of three dimensionless parameters. 7 The compounds Yb 2 Ti 2 O 7 , Er 2 Ti 2 O 7 , Pr 2 Sn 2 O 7 1 (and possibly Tb 2 Ti 2 O 7 8 ) are of this type, and it has recently been argued that the spins in Yb 2 Ti 2 O 7 and Er 2 Ti 2 O 7 are controlled by exchange coupling rather than by the long-range dipolar interactions which dominate in spin ice. 9,10 This makes these materials beautiful examples of strongly quantum magnets on the highly frustrated pyrochlore lattice. They are also nearly ideal subjects for detailed experimental investigation, existing as they do in large high purity single crystals, and with large magnetic moments amenable to neutron scattering studies. Yb 2 Ti 2 O 7 is particularly appealing because its lowest Kramers doublet is extremely well separated from the first excited one, 11 and a very large single crystal neutron scattering data set is available, allowing us to measure the full Hamiltonian quantitatively, as we will show. Although we specialize to the latter material in the present article the theoretical considerations and parameter determination method described here will very generally apply to all pyrochlore materials where exchange interactions dominate and whose dynamics can be described by that of a single doublet.
Theoretical studies have pointed to the likelihood of unusual ground states of quantum antiferromagnets on the pyrochlore lattice. 12,13 Most exciting is the possibility of a quantum spin liquid (QSL) state, which avoids magnetic ordering and freezing even at absolute zero temperature, and whose elementary excitations carry fractional quantum numbers and are decidedly different from spin waves. 14 Intriguingly, neutron scattering measurements have reported a lack of magnetic or-dering and the absence of spin waves in Yb 2 Ti 2 O 7 at low fields. 15,16 In a recent study, sharp spin waves emerged when a magnetic field of 0.5T or larger was applied, suggesting that the system transitioned into a conventional state. 16 The possible identification of the low field state with a quantum spin liquid is tantalizing, but progress certainly requires a more detailed understanding of the spin Hamiltonian.
In this article, we present a detailed experimental and theoretical investigation of the excitation spectrum in the high field state throughout the Brillouin zone. We show that the spectrum is extremely well fit by spin wave theory, and through this fit we unambiguously extract all the microscopic exchange parameters in the spin Hamiltonian (see below). Interestingly, we find that the largest exchange interaction is of the same "frustrated ferromagnetic" Ising type as in spin ice, despite the fact that the g-tensor tends to orient magnetic moments primarily normal to the Ising axes. Moreover, spin-flip terms which induce quantum dynamics are comparable to the Ising exchange, which confirms the picture of Yb 2 Ti 2 O 7 as a strongly quantum magnet, in qualitative agreement with recent studies. 10,17,18 Strikingly, we find that the predictions of mean field theory (MFT) using these parameters disagree drastically with experiment in zero field, indicating that fluctuations strongly reduce or destroy the classically expected spin order. Taken together, these observations make Yb 2 Ti 2 O 7 a promising candidate for observation of QSL physics. The precise determination of the microscopic spin interactions sets the stage for a quantitative understanding and test of this proposal.
Time-of-flight neutron scattering measurements were performed on a 7g single crystal of Yb 2 Ti 2 O 7 , grown via the optical floating zone method. Details of the crystal growth were given elsewhere. 16,19 The neutron scattering data was collected at the Disk Chopper Spectrometer at the NIST Center for Neutron Research, using 5Å incident neutrons. This configuration allowed an energy resolution of 0.09meV. The sample environment consisted of a 11.5T vertical field magnet combined with a dilution insert that achieved a base temperature of 30mK. The scattering plane was HHL, with the field applied along the [110] vertical direction. The sample was rotated 147 degrees in 1.5 degree steps about the vertical, allowing a three dimensional data set to be acquired, i.e. two components of the wavevector Q within the scattering plane, and energy transfer. The spin excitation spectra along several high symmetry directions with the scattering plane were thereby obtained.
At 30mK, the inelastic spectrum changes qualitatively at H = 0.5T; above this field strength, resolution-limited spin wave excitations which go soft with quadratic dispersion at nuclear-allowed positions develop, indicating a transition to a field-polarized ferromagnetic state. 16 The spin wave excitations indicate that the symmetry of the underlying lattice is preserved, as is evident from gaps in the spectrum at the nuclear zone boundaries. In Figure 1 we show the spin wave dispersions along several directions in the HHL plane for H = 2T and H = 5T. These high symmetry directions are shown relative to the Brillouin zone structure within the HHL plane in Fig. 2.
We compare the experimental data to spin-wave theory. We assume nearest-neighbour exchange coupling only, as appropriate to the strongly localized f -electron states, and neglect dipolar interactions. The Hamiltonian, written in global spin coordinates, is then where J µν ij = J νµ ji is the matrix of exchange couplings between sites i and j, g µν i is the g-tensor for the spin at site i, and we takeh = 1. Symmetry allows four independent exchange constants, 7 J 1 , · · · , J 4 . To specify them, we give the exchange matrix on one pair of nearest-neighbour sites, located at positions r 0 = a 8 (1, 1, 1) and r 1 = a 8 (1, −1, −1) on a tetrahedron centered at the origin (a is the conventional cubic lattice spacing for the fcc Bravais lattice): The other exchange matrices can be obtained from this one by cubic rotations given in the Supplementary Information. The g-tensor contains two components: g z parallel to and g xy perpendicular to the local C 3 rotation axis through the Yb site.
Spin wave theory, carried out as described in the Supplementary Information, is fit to the H = 5 T, T = 30 mK measurements; the fitting procedure focuses on the dispersion relation alone, and the overall intensity of the calculated spin waves is scaled to agree with the experiment at a single wavevector and energy point. The resulting inelastic structure factor S(Q, ω) (see Supplementary   Information) is compared to both the 5T and 2T data in Figure 1. The best fit is achieved with the following exchange parameters, in meV: Here we quote rough uncertainties obtained by the visual comparison of the theoretical and experimental intensities. The fit is performed by taking the ratio of components of the g-tensor to be g xy /g z = 2.4, i.e. the ratio obtained by Ref. 11. The fit then produces g z = 1.80, in nearly perfect agreement with these studies (using the g-factor ratio of Cao et al. instead, 17 i.e. g xy /g z = 1.8, does not reproduce the data as precisely). Using these results, a high temperature expansion gives (see Supplementary Information) a theoretical Curie-Weiss temperature Θ CW = 312mK, which is comparable to but smaller than the experimentally determined values, Θ CW = 400mK 20 and 750mK. 11 The deviations may be explained by the sensitivity of the theoretical value to small changes in the g-factors and exchange parameters, and to the dependence of the experimental value on the fitting range. 17 Furthermore, and most importantly, our extracted exchange parameters correctly reproduce relative intensities as well as the shape of the spin wave dispersion for each of the five directions. Agreement is excellent for H = 2T, showing that these parameters produce a robust description of the field-induced ferromagnetic state. We note, however, that there is a significant quantitative disagreement with the exchange parameters quoted in Refs. 9, 10 (see Supplementary   Information).

Implications:
The excellent agreement with spin-wave theory for fields H ≥ 2T clearly indicates that the high field state is accurately modeled semi-classically, and is smoothly connected to the fully polarized limit. Theoretically, the ground state in this regime breaks no symmetries, and supports a ferromagnetic polarization along the axis of the applied field (for the 110 field used in the experiment). However, the semiclassical analysis clearly and dramatically fails at small fields, where the measurements show no signs of spontaneous long range order. 16 The classical zero field ground state for our Hamiltonian parameters has a large spontaneous polarization along the 100 axis. Extension of this analysis to a T > 0 mean-field theory wrongly predicts a continuous magnetic ordering transition at a temperature of T M F c = 3.2K (see Supplementary Information).
The experimental indications of a zero-field transition to long range order are mixed, 21, 22 but early specific heat measurements 20 found an anomaly at T c = 214mK, and Mössbauer spectroscopy 23 suggested a transition at 240mK. This temperature is approximately 14 times lower than T M F c .
If there is magnetic ordering at all, it appears to be substantially suppressed, indicating strong fluctuations -classical, quantum, or both.
The presence of strong fluctuations makes a QSL ground state plausible in low field. We now use the Hamiltonian parameters to suggest the nature of this state. To do so, we rewrite the zero field Hamiltonian in terms of spins quantized along the local C 3 axis for each site, similarly to Ref. 18: where here S µ i are local spin coordinates, , and the matrices γ ij , ζ ij consist of unimodular complex numbers (given in the Supplementary Information). From the fits in Eq. (3), we find, in meV, where the uncertainties have been estimated by treating those in Eq. (3) as Gaussian random variables. Note that the strongest interaction is J zz > 0, which precisely coincides with the nearestneighbour spin-ice model. The model with J ± and J zz only has been studied theoretically. 13,24 It does indeed support a QSL ground state, for sufficiently small J ± /J zz . For larger J ± /J zz , the ground state is instead a magnet with S ± i = 0. 24 While the actual value of J ± /J zz ≈ 0.3 would place this model in the magnetic state, 24 the J z± interaction in particular is non-negligible in Yb 2 Ti 2 O 7 , and preliminary theoretical work suggests that it tends to stabilize the QSL state.
Indeed, in perturbation theory, the leading effect of the J z± coupling is to induce in the effective A key consequence of the QSL scenario is the presence of a quantum confinement phase transition in applied field (see Fig. 3). Such a quantum phase transition is required to remove the "fractional" excitations of the QSL phase (electric and magnetic monopoles) from the spectrum in the semiclassical high field phase. Theoretically, such quantum critical points have been studied in related model Hamiltonians, and occur by a mechanism analogous to the Higgs transition in the standard model 26 or by magnetic monopole condensation. 27

Cubic Rotations and Local Bases
As described in the main text, we use the usual coordinate system for the pyrochlore lattice, with sites located on tetrahedra whose centers form a FCC lattice. We take one to be centered at the origin with its four corners at r 0 = a 8 (1, 1, 1), r 1 = a 8 (1, −1, −1), r 2 = a 8 (−1, 1, −1) and r 3 = a 8 (−1, −1, 1). The exchange matrices J ij between sites of types i and j are obtained by applying the following cubic rotations R ij to J 01 : Note J ji = J T ij .
We use the following local b i =ê i ×â i , and the 4 × 4 complex unimodular matrices for which our exchange Hamiltonian takes the form of Eq. (4).

Curie-Weiss Temperature
A high-temperature expansion of the Hamiltonian of Eq. (1) yields the O(1/T 2 ) term in the uniform susceptibility from which we extract the Curie-Weiss temperature, where k B is the Boltzmann constant. Using the formulation of Eq. (4) Θ CW takes the simpler form

Spin Wave Theory
As usual, we expand the Hamiltonian about one of the classical states using Holstein-Primakoff bosons in the spirit of large s, and keep only terms of and up to order s, which shall then be set to 1/2. We define the transverse Holstein-Primakoff bosons x a = x † a , y a = y † a , conjugate with one another on site a of the pyrochlore lattice, satisfying such that Here (v a , w a , u a ) is an orthonormal basis, chosen so that u a gives the direction of the spin in the classical ground state at site a (which we find numerically). We find that for all fields, the ground state does not enlarge the unit cell, so that there are only four distinct such bases, which we denote by a = 0, .., 3. One may choose, for example, v a = u a × (1, 1, 1)/ u a × (1, 1, 1) and Since the classical ground state does not enlarge the unit cell, we can readily proceed to Fourier space in the four site basis. Keeping only terms of order s, we arrive at the spin-wave (quadratic) Hamiltonian, where X T Y T = x 0 .. x 3 y 0 .. y 3 and D ab k =D ab cos(k·(r a −r b )), D = A, B, C, where H is the magnetic field and J ab and g a are 3 × 3 matrices with matrix elements J µν ab and g µν a , respectively. To find the modes, we resort to the path-integral formulation. The action at where ω n = 2πn β is the Matsubara frequency, where we have defined Z T = X T Y T , (1 4 is the four-by-four identity matrix). As usual, the real frequency dispersion relations ω(k) are found by solving the zero eigenvalue equations of the matrix G k +ωΓ. Here, these are equivalently the (both zero and non-zero) eigenvalues of −ΓG k .
We calculate the inelastic structure factor (to which the intensity I(k, ω) of the scattering is proportional) obtained from the moment-moment correlation function, 28 wherek is the unit vector associated with k, m µ a (k, ω) = µ B κ g µκ a S κ a (k, ω) is the moment at momentum k and real frequency ω on the a = 0, .., 3 sublattice in direction µ = x, y, z, and S κ a is, as usual, the κth coordinate of the effective spin-1/2 spin at site a. F µν (k) = δ µν − (k) µ (k) ν selects the component of the spin-spin correlations perpendicular to the scattering vector, 28 and ω) is the moment-moment correlation function that originates from the interaction between the neutrons' moment and the spins' moments, which we find to be where δ is the Dirac distribution, η µ a = κ=x,y,z g µκ a v κ a w κ a , α is the αth eigenvalue of −ΓG, and ψ R,α is its corresponding "right" eigenvector, i.e. such that −ΓG · ψ R,α = α ψ R,α . Also note that the momentum and frequency dependences are implied everywhere to be k, ω.
Now we estimate the amplitude of quantum fluctuations using our spin wave theory, by evaluating the quantum moment reduction in zero field and at zero temperature. In the Holstein-Primakoff boson language the reduction of the spin expectation value from the classical value of 1/2, averaged over the four site basis, is r = 1 4 3 a=0 0|n a |0 , where 0|n a |0 indicates a ground state quantum expectation value. Evaluating this using the path integral method, we obtain (Θ is the usual Heaviside distribution), which is a 10% reduction compared with the classical value of 1/2.

Further Details of Experimental Data and Fits
The inelastic neutron scattering data (rows 1 and 3 in Figure 1) contains several features worth commenting on further. First, the darkest blue areas do not contain any data, either for kinematic reasons near Q = (0, 0, 0), or because of the finite angular extend of the scan. Second, at E = 0 one observes intensity which is due to coherent and incoherent elastic scattering from the sample, and hence is more intense than the inelastic features by up to several orders of magnitude, thus appearing red (off scale). Third, near the (0, 0, 0) position there is higher background leading to unphysical intensities due to contamination from the un-scattered incident beam. This is observable in the HHH and HH0 data sets near the zero position.
The data sets were collected by counting for 8 minutes per angular rotation of the sample. The total time per magnetic field setting was about 12 hours.
Supplementary Figure 1 shows the dispersions obtained from the fitting procedure over-plotted on the data. These fits were accomplished by digitizing the shape of the dispersion from the experimental data, and performing a least squares minimization routine to match it. Intensities were calculated based on the spin wave theory using the extracted exchange parameters, and were not fit to the data. The 5D parameter space (four exchange plus one g-tensor parameter) was explored using a uniform search technique with the resulting excellent description obtained ( Figure   1).
Supplementary Figure 1: Dispersions (white curves) obtained from the fitting procedure overplotted on the data at H = 5 T.

Features of the Spin Wave Spectrum
The spectra of Figure 1  around Yb 2 Ti 2 O 7 and for H = 5 T, we find its energy to be, numerically: where the J α 's must be input in meV, and which for our fit gives E Yb 2 Ti 2 O 7 2d flat ≈ 1.45 meV. Note that the energy of this feature is most sensitive to J z± .
Moreover, we observe the following trends in the region around Yb 2 Ti 2 O 7 : • as |J z± | increases, all bands go up in energy (especially the two-dimensional flat one), • increases in J ± and J ±± seem to have more or less the same effect: the bands get closer, and this happens in particular because the energy of the top bands decreases.

Mean Field Theory Calculation
The Curie-Weiss mean field Hamiltonian obtained from Eq. (1) takes the form where H µ is the magnetic field in the µ direction, The ground state does not enlarge the unit cell at zero temperature, and we assume this is still the case at non-zero temperature. Thus, we define, for every sublattice a and every axis µ the average magnetization which is the same for every tetrahedron t, which can be either "up" or "down". We arrive at the twelve consistency equations where h eff a = 2 b m b · J ba − µ B H · g a ; the free energy per site is

Perturbation Theory
We show that the U (1) QSL described by Hermele et al. 13,24 is stable to the addition of the terms in the symmetry-obtained Hamiltonian Eq. (4), provided all coupling constants are small with respect to the Ising exchange parameter J zz . In other words, we show that the U (1) QSL exists in a finite region of parameter space, specifically where J ± , J z± and J ±± are small with respect to J zz .
In that limit, one may apply perturbation theory. When J zz > 0, the ground state manifold of the unperturbed Hamiltonian is the extensively degenerate "two-in-two-out" manifold. When J z± = J ±± = 0, Eq. (4) can be mapped exactly onto the Hamiltonian of Ref. 13, where the first non-vanishing and non-constant term in perturbation theory (above the "two-in-two-out" manifold) was shown to be third order in J ± /J zz : where K = Refs. 13, 24, 29, this ring Hamiltonian has as its ground state a U (1) QSL, whose low energy physics is described as the Coulomb phase of a U (1) gauge theory. This phase is locally stable to all perturbations 13 in three dimensions, which is enough already to guarantee the persistence of the QSL state when the other exchange couplings are sufficiently small, i.e. when the induced terms in the effective Hamiltonian are much smaller than the ring coupling K.
We can, however, go further and consider these effects explicitly in the perturbative limit, which extends the discussion to the case when J z± , J ±± J zz but with no particular assumptions placed upon the magnitude of the induced terms in the effective Hamiltonian relative to K. When J z± or J ±± are non-zero, other "non-ring" effective Hamiltonians are allowed. In particular, the J z± term gives rise to an effective third neighbour ferromagnetic Ising Hamiltonian: where J (3) = S z i S z j . A blue sphere represents an "in" spin while a black sphere represents an "out" one. Because two of the four chain types contain non-alternating spins, there are no flippable hexagons in any of the six (equivalent) ground states.
contain flippable hexagons represent an energy cost. We can therefore consider, H eff 3rd Ising as an analog of the Rokhsar-Kivelson term introduced by Hermele et al., 13, 25 H RK = V N f , where N f is the operator that counts the number of flippable hexagons. This term actually stabilizes the QSL state. 13,29 In particular, the QSL phase grows in stability as V is increased from zero up to the point V = K, beyond which (for V > K), the system undergoes a first order transition to a degenerate set of classical unflippable states, which includes the six ordered 100 ferromagnetic ground states described above.
In the case relevant here, when J (3) /K is sufficiently large, the system must undergo a transition to the unflippable 100 ground states. Since this model and the RK one of Ref. 13 agree on the phases both when J (3) K and when J (3) K, it is natural to expect that the intervening phase diagram coincides in the two models as well. Therefore we expect that the QSL state is maximally stable when J (3) /K takes some value of O(1). (For the values of the exchange constants given by our fits we find J (3) /K = J 2 z± Jzz 4J 3 ± = 6.2, but we caution that this is probably outside the perturbative regime). We note in passing that the J (3) exchange can also be expressed in terms of purely plaquette interactions, which might allow further analytical connection to the RK theory.
We will not, however, pursue this further here.
Inclusion of the other coupling J ±± , higher order effects, and cross terms amongst the exchange couplings does not lead to any new effects. Indeed, all the associated terms in the effective Hamiltonian assume a ferromagnetic or ring form, and can be subsumed in the above couplings.
8 Comparison between our Exchange Constants and those of Thompson et al. 10 The correspondence between our effective spin-1/2 operators S i in Eq. (1) and the full 7/2-angularmomentum operators J i used by Thompson et al. in Ref. 10 is given by projecting the full angular momentum into the ground state Kramer's doublet where P 1/2 is the projection operator to the ground state Kramer's doublet, and g J = 8/7 is the Landé factor.
We use g xy /g z = 2.4 and g z = 1.79 for concreteness, but the results do not depend too much upon the details of this choice within the range of parameters found in the literature. With this different from those found in our fits, Eq. (3) and Eq. (5), and indeed yield spin wave spectra in strong disagreement with experiment. In principle, one resolution of the difference could be that the exchange couplings are actually strongly temperature dependent, and distinctly different in the temperature range studied in Ref. 10 and at the low temperatures studied here. Such a change in exchange parameters could conceivable occur if the transition observed at 214-240mK had a substantial structural component. However, we do not have any a priori reason to suspect this. In any case, it would be very interesting to resolve the differences between the two exchange models.