# Low Energy Electrodynamics of Novel Spin Excitations in the Quantum Spin Ice YbTiO

In condensed matter systems, the formation of long range order (LRO) with broken symmetry is often accompanied by new types of excitations. However, in many magnetic pyrochlore oxides, geometrical frustration suppresses conventional LRO while at the same time non-trivial spin correlations are observed. For such materials, a natural question to ask then is what is the nature of the excitations in this highly correlated state without broken symmetry? Frequently the application of a symmetry breaking field can stabilize excitations whose properties still reflect certain aspects of the anomalous state without long-range order. Here we report evidence of novel magnetic excitations in the quantum spin ice material YbTiO, obtained from time-domain terahertz spectroscopy (TDTS). At large fields, both magnon and two-magnon-like excitations are observed in a 001 directed magnetic field illustrating the stabilization of a field induced LRO state. The unique ability of TDTS to measure complex response functions allows a direct study of magnetic responses in different polarization channels, revealing the existence of an unusual left-hand polarized magnon. The g-factors of these excitations are dramatically enhanced in the low-field limit, showing a cross-over of these one- and two-magnon states into features consistent with quantum string-like excitations proposed to exist in quantum spin ice in a small 001 applied field.

Geometrical frustration suppresses the formation of conventional long range order (LRO) in magnetic materials to a temperature () much lower than the Curie-Weiss temperature . This provides the possibility of realizing novel states and excitations, since non-trivial spin correlations are often observed in these materials in the range T even in the absence of LRO. Magnetic pyrochlore oxides, in which magnetic rare earth ions sit at the vertices of corner-sharing tetrahedra, provide a fascinating material system for the search of such novel quantum ground states and excitations.

In spin ice materials such as DyTiO and HoTiO, strong crystal field anisotropy forces the spins to point along the local 111 direction while long-range dipolar interaction between the spins provides an effective ferromagnetic interaction, giving a situation much resembling proton disorder in common water ice. Enforcing the ices rules in spin ice results in a ground state with each tetrahedron having the so-called “two-in, two-out” configuration. Flipping a single spin creates a pair of magnetic monopoles out of the ground state vacuum in neighboring tetrahedra, which can subsequently move away from each other through the network. In applied field, a pair of magnetic monopoles with opposite magnetic charges are connected by a string of spins aligned against the field (Fig. 1(a)). Recent neutron scattering experiments have shown evidence for the existence of thermally activated Dirac strings and magnetic monopoles in classical spin ice materials.

Recently the material YbTiO has attracted considerable attention. Although the exact form of the low temperature ground state is still under debate, the nature of the magnetic interactions place YbTiO in the quantum spin ice regime. High-field inelastic neutron scattering (INS) experiments have obtained the exchange interaction parameters at low temperature, which demonstrates that in this material the ferromagnetic Ising type exchange interaction is larger than the other terms in the Hamiltonian that lead to considerable quantum fluctuations and dynamics. From considerations of crystal-field levels, its low-energy spin sector may be reduced to that of an effective spin-1/2 moment. A number of interesting effects have been proposed for this material, including it being a quantum spin liquid, whose elementary excitations carry fractional quantum numbers. Such a state supports an emergent quantum electrodynamics with a photon mode at low energy. It was proposed recently that under weak applied 001 magnetic field (or in the presence of spontaneous magnetization), the elementary excitations of materials like YbTiO take the form of extended “quantum strings” consisting of fluctuating multiple flipped spins connecting monopole pairs. With inherent quantum dynamics, these novel excitations are extended objects rather than point particles, an interesting feature of quantum spin ice.

To explore the nature of the excitations in quantum spin ice, we performed time-domain terahertz spectroscopy (TDTS) measurements of YbTiO single crystals under magnetic field applied along the 001 direction. Spectroscopic features consistent with the picture of quantum strings are observed in the low-field regime. When the field strength increases, a crossover towards field-induced order is observed where the excitations take the form of magnons and two-magnon excitations.

## I Results

### i.1 Magneto-optical Measurement

Magneto-optical measurements are performed in a custom built time domain terahertz (THz) spectrometer. Two single crystals from the same traveling solvent floating zone boule are used in the measurement; they will be referred to as Sample A and B throughout the text. The single crystal YbTiO samples are cut and polished with their surfaces oriented along the 001 plane. The sample is hosted in a 7 T superconducting magnet. Experiments are performed in transmission geometry as illustrated in Fig. 1(b), with temperatures of the samples on the scale of the exchange interaction of YbTiO. These temperatures are about one order of magnitude higher than the value of the debated ordering temperature , while at the same time low enough that spin correlations are well developed to allow the study of quantum spin ice physics. The incident THz pulse is linearly polarized, with the wavevector k of the pulse defining the z direction. Measurements are performed in both the Faraday geometry (applied magnetic field in the z direction) as well as the Voigt geometry (field in the y direction).

In magnetic insulators, the complex transmission function we obtain is related to the complex frequency dependent susceptibility as . TDTS measures the response and as such is very useful for materials with interesting correlations. The technique operates at a frequency range most relevant to the physics of quantum spin ice materials (0.4 to 10 meV); it has the strength of high signal to noise, excellent energy resolution, and short acquisition times. Further details can be found in the Methods section.

### i.2 Magnetic Excitations in Faraday Geometry

In the Faraday geometry, the applied DC magnetic field is oriented along the THz pulse propagation direction, and the incident THz pulse is linearly polarized along the x axis of the laboratory frame. With the spins aligned and precessing around the applied field direction, one naturally expects magnetic absorption to be observed in the right circular polarized (RCP) channel, which is a direct consequence of the fundamental coupling between magnetic field and spin precession. As a result, the magnetic excitations in this geometry would best be studied in the circular polarization frame, which can be measured with a phase sensitive spectroscopic tool such as TDTS.

In Fig. 2(a) and (b) we show the transmission magnitude for the right and left circularly polarized (RCP and LCP) light as a function of frequency and field measured at 4.3 K from Sample B in the Faraday geometry (See Methods for transformation to circular base). The measurements were taken every 0.5 T from 0 to 7 T. Clear dark lines diagonally across the figures immediately show the position of the absorptions in the spectra, revealing the existence of several branches of magnetic excitations. As we show in the supplementary information, spectra taken at 1.65 K and 4.3 K are almost identical, except that the absorptions are slightly sharper and stronger at 1.65 K. Recalling that the quantity we explicitly measure (and are plotting) is related to , we observe that these excitations are generally losing spectral weight in as the field increases and they move to higher . For the RCP channel, careful examination of the data reveals a total of five branches for fields higher than 3 T. They are labeled b1, b2, b3, b5, and b6 from low to high frequency. There is a single LCP branch, b4, observed for fields higher than 4 T. This absorption in a +z magnetic field is anomalous and will be discussed in more detail below. While branches b1 to b4 have similar slope in the frequency-field plot, suggesting similar effective g-factors; branches b5 and b6 clearly have higher g-factors.

In YbTiO, there are four magnetic ions in each primitive unit cell; naively one would expect at most four spin waves modes, as has been observed in INS experiments recently. The fact that six different branches of magnetic excitations are observed reveals the unconventional nature of some of these excitations. Another interesting feature is that as the observed excitations move to low-frequencies at low-field, they become less well defined and eventually turn into a broad magnetic absorption. Similar features are also obtained with sample A, from temperatures 1.65 K upward. All the features other than b1 and b2 fade away as temperature increases above 12 K. Details of the temperature dependence as well as additional data can be found in the supplementary information.

### i.3 Magnetic Excitations in Voigt Geometry

In the Voigt geometry, the magnetic field is oriented along the y direction. As a result, the circular frame is no longer an appropriate base for the study. Symmetry considerations dictate the absence of circular effects with this geometry. Shown in Fig. 3(a) and (b) are the transmission magnitudes with electric polarization along the x direction, , measured in Sample A at 1.65 K, with Fig. 3(b) showing a more detailed measurement below 1 T, while results for are shown in Fig. 3(c). A total of four branches of magnetic excitations are observed, and labeled c1 to c4 from low to high frequency. While c1 to c3 show similar g-factors, c4 shows a higher value. Excitations c2 and c4 show a notable downward concavity at low magnetic fields. As will be discussed in detail below, we believe that these features are consistent with a picture of quantum spin ice where high-field one- and two-magnon states continuously evolve into string excitations at low field.

## Ii Discussion

### ii.1 One- and Two-Magnon Features

As shown in Fig. 2 and 3, the magnetic excitations we observe at fields greater than 3 T in both the Faraday and Voigt geometries can be categorized into two groups: b1 to b4 and c1 to c3 having a similar, lower value of the g-factor (2.96 to 3.86, values listed in Fig. 4 caption), while b5, b6 and c4 have larger g-factors (6.30 to 6.74). The frequencies of the observed magnetic excitations (with demagnetization correction applied) are summarized in Fig. 4(a).

Upon examination of the complete data set, one notices that the energies of the magnetic excitations with lower g-factors are in the same range as the magnons observed by INS in the field induced ordered state of YbTiO . Although the direction of the applied magnetic field in the INS experiment is different from the present study, this provides a useful hint in establishing the nature of those branches of excitations observed in our TDTS as magnons.

To form a better understanding of the observed magnetic excitations, we performed a spin-wave analysis with the interaction parameters obtained in Ref. 14 as the starting point. Linear spin-wave theory predicts a total of four magnon modes, given the four unique spins in the unit cell of YbTiO. Only two of the magnon modes are predicted to be visible as in the Faraday geometry, while in the Voigt geometry three modes are predicted to be visible. As shown in Fig. 4, both the energy and the polarization state of the excitations observed in the LCP in the Faraday geometry as well as the ones in the Voigt geometry match the spin-wave theory prediction very well at high fields, where the field induced LRO supports the existence of well defined magnon modes. With the high signal to noise ratio of the TDTS, the data above 3 T from these three branches allowed a further refinement of the exchange parameters in a least squares fitting procedure (see supplementary information for further details).

One surprising aspect of the data was the clear observation of a LCP branch b4 when the magnetic field was applied in the +z direction (Fig. 2b). As mentioned above, at first sight this seems to be in contradiction to the natural coupling between magnetic field and spin moments. Previously such LCP spin precession has only been proposed for metallic ferromagnets with strong spin-orbit interaction. This interesting observation finds a natural explanation in the context of the spin-wave theory. With the easy plane anisotropy of the Yb ion, the equilibrium direction of the spin is tilted away from the 001 axis, moreover, the spin precession is slightly elliptical around the applied field in the 001 direction. Thus, the precession motion has both LCP and RCP components when viewed along the field direction. This is schematically illustrated in Fig. 4(c) and (d). For the highest-energy magnon mode, which corresponds to the b4/c3 branches, the phase relationship between the spins in the unit cell is such that the RCP components cancel exactly, leaving only the LCP component in the Faraday geometry.

Another interesting feature of the data is the observation of excitations with clearly higher g-factors than the magnons. A closer examination of the data reveals that excitations b5, b6 and c4 have a g-factor about twice the value of the magnon branches. Interestingly they also appear in the middle of the two-magnon continuum band obtained from the spin-wave calculation, as shown in Fig. 4(a). These suggest that, while the branches b1 to b4 and c1 to c3 can be understood as magnons, b5, b6 and c4 are two-magnon excitations. The fact that they appear in the middle of the calculated two-magnon continuum suggests they are unlikely to be simple two-magnon bound states. Thus it remains to be understood how these excitations can be stable and have well defined energies. Moreover our spin-wave calculations show that these excitations do not match any peaks in the joint two-magnon density of states. Note that the two-magnon continuum band intersects with the single magnon modes at low fields. This gives the possibility of decay for the higher energy single magnons at low fields.

Also at odds with the spin-wave calculation is the observation of four single magnon-like branches in the Faraday geometry, as only the highest and the lowest modes are expected to be visible. In comparing the geometries, b1, b2, and c1 all appear in the energy range of the calculated lowest magnon mode. Our data suggest that b1 and b2 are distinct features (more clearly shown in Fig. SI1 and SI5), although they are difficult to distinguish due to the strong absorption in this frequency range. They may derive from the same lowest-energy magnon state, but with some unexplained splitting. Moreover, it is difficult to understand b3 as the second or third highest magnons that have picked up a finite intensity by some means (from misalignment of the applied field for instance). Its energy is below the c2 excitation in the Voigt geometry, which is straightforwardly assigned to the second highest magnon mode. These discrepancies which appear only in the Faraday geometry are as of yet unexplained. Further details of the spin-wave analysis of the four magnon modes are given in the supplementary information.

### ii.2 The Crossover to Quantum Strings

As field decreases, the match between the spin-wave calculation and the experiments becomes progressively worse (Fig. 4(a)). At the temperature range of the study, spin correlations are well developed and at low fields the material approaches the nominal quantum spin ice regime. Thus it is not surprising to see the spin-wave theory fails at low fields. The spectroscopic features in this region reveal the nature of the excitations unique to the quantum spin ice state. Among those features, the most prominent one is the change of slope of the magnetic excitations at low fields. We believe this large concave downward non-linearity is a general signature of quantum mixing of states with different numbers of spin flips. When the field is lowered from the high-field regime (where the excitations should be understood as spin flips of a quantized number e.g. magnons) the excitations of the system are now expected to be admixtures of chains of flipped spins with different lengths, as proposed in the quantum string picture and also observed in spin chains . In this scenario one expects that as the field diminishes, strings get longer. As a result their g-factor increases and the slope of the curve goes up. The strong enhancement of these g-factors below 2 T clearly observed for branches c2 and c4 as field decreases provides strong support for this picture, as shown in Fig. 4(b). This strong increase of g-factors is evidence of significant admixing of states of many flipped spins in this part of the phase diagram.

The quantum string picture was derived in the low-field limit and requires the breaking of time-reversal symmetry either by a 001 directed field or a spontaneous magnetization in this direction. It was proposed that a hierarchy of states would be observed corresponding to strings of different lengths. However, in real materials with multiple spins in the unit cell and finite non-Ising exchange terms, such a hierarchy of states might be difficult to realize explicitly due to the decay channels opening up with overlapping bands of multi-string and multi-magnon continuums (as implied in Fig. 4(a)). In the present study high-field two-magnon excitations are observed to evolve continuously into string states at low fields, although their stability remains to be understood. In any case, one may still expect that the lowest magnon features observed at high fields evolve into string-like states at low fields. Our observations are consistent with such a picture that the lowest elementary excitations of a quantum spin ice in small magnetic fields are quantum strings connecting weakly bound monopoles. As discussed above, the application of a symmetry breaking field can stabilize excitations whose properties reflect certain aspects of an anomalous state without long-range order. In the present case, the observation of a crossover from magnons to strings as the field is decreased reflects the presence of deconfined monopoles in the zero-field quantum spin ice state.

## Iii Methods

### iii.1 Sample Fabrication and Characterization

Single crystal samples of YbTiO were grown at McMaster University with the optical floating zone technique. Details of the sample growth are discussed elsewhere. Samples A and B used for the THz measurements are two single crystals with their largest surface oriented parallel to the 001 plane. Sample A has a circular cross section with a diameter of 4.2 mm, with thickness of 0.725 mm, while Sample B has a rectangular cross section that measures 5.46 by 2.37 mm, with a thickness of 0.648 mm. Both samples are transparent. Laue x-ray diffraction was performed to assure the samples used do not contain domains of misaligned microstructure or crystallites.

The low temperature heat capacity of Sample B has been measured from 150 mK to 800 mK, and shows the same behavior as the non-annealed crystal reported in Ref. 17. As shown in Ref. 33 and 35, at temperatures of our current study, the exact nature of the low temperature state is not relevant, as the entire phase diagram is expected to be in the thermal spin liquid state at temperatures on the scale of . Also, as reported in Ref. 17, 20, and 21, sample variations are minimum for temperatures above 1 K, thus the possible sample variation does not affect our conclusion in this paper.

### iii.2 Time Domain Terahertz Spectrometer Set Up

A home-built time-domain terahertz spectroscopy was set up at Johns Hopkins University with two dipole switches as THz emitter and receiver. The sample was hosted in a cryogen free split coil superconducting magnet. The magnet provides magnetic fields up to 7 T, and a base temperature down to 1.6 K for the measurements. The magnet has four windows allowing for optical access, and can be switched between Faraday and Voigt geometries. The single-crystal samples were mounted on a metal aperture. The electric field profiles of the terahertz pulses transmitted through the sample and an identical bare aperture were recorded as a function of delay stage position, and then converted to time traces. FFTs of the sample and aperture scans gave the complex transmitted and incident pulse spectra, from which the complex transmission coefficients were obtained.

### iii.3 Polarization Modulation Technique in Terahertz Polarimeter

As shown in the previous sections, the simultaneous measurement of two orthogonal components of the transmitted THz pulse is critical for this experiment. This is realized by a terahertz polarimeter with polarization modulation technique based on a rotating polarizer developed recently. The main set up contains three polarizers, one rotates at angular velocity that modulates the polarization state of the THz pulse to be measured (shown in Fig. 1(b)); another one downstream oriented with the transmission axis parallel to the x axis. Another polarizer upstream creates a linearly polarized incident pulse before the sample.

The effect of the two polarizers in the polarimeter on an arbitrary electric field vector can be expressed as a product of two matrices:

(1) |

where:

(2) |

straightforward algebra gives:

(3) |

Here =t+ is the angle of the rotating polarizer with respect to the x axis, is the initial angle. The signal detected after the polarimeter is then fed into a lock-in amplifier with a reference signal of frequency 2. In this way, the in-phase and out-of-phase readings of the lock-in amplifier measure and simultaneously, with an appropriate choice of . Note that it is important that the sense of rotation (direction of the angular velocity ) of the rotating polarizer be the same as the THz pulse propagation direction; otherwise instead of , will be measured, which would give a reversed definition of right/left circular polarization. See Ref. 25 for further details on this technique.

### iii.4 Transforming into Circular Base in the Faraday Geometry

Following the previous section, when and of the transmitted pulse are measured while retaining their phases, it is possible to transform the frame of reference into circular polarized frame, thus providing an opportunity to study the magnetic excitations in their natural circularly polarized basis in the Faraday geometry. With , standing for the electric field component in RCP and LCP channels, we have:

(4) |

More details of this transformation of bases can be found in the supplementary information.

## References

- (1) Lacroix, C., Mendels, P., Mila, F. (Eds.), Introduction to Frustrated Magnetism, Springer Series in Solid-State Sciences, Vol. 164, Springer, 2011.
- (2) Balents, L. Spin liquids in frustrated magnets. Nature, 464, 199-208 (2010).
- (3) Gardner, J.S., Gingras, M.J.P. & Greedan, J.E. Magnetic pyrochlore oxides. Rev. Mod. Phys., 82, 53-107 (2010).
- (4) Bramwell, S.T. & Gingras, M.J.P. Spin ice state in frustrated magnetic pyrochlore materials. Science, 294, 1495-1501 (2001).
- (5) Bramwell, S.T., Gingras, M.J.P. & Holdsworth, P.C.W. chapter on Spin Ice in Frustrated Spin Systems, (ed. Diep, H., World Scientific, 2004).
- (6) Ryzhkin, I.A. Magnetic relaxation in rare-earth pyrochlores. J. Exp. Theor. Phys, 101, 481 -â 486 (2005).
- (7) Castelnovo, C., Moessner,R. & Sondhi, S.L. Magnetic monopoles in spin ice. Nature 451, 42 – 45 (2008).
- (8) Fennell, T. et al. Magnetic Coulomb phase in the spin ice HoTiO. Science 326, 415 – 417 (2009).
- (9) Morris, D.J.P. et al. Dirac strings and magnetic monopoles in the spin ice DyTiO. Science 326, 411 – 414 (2009).
- (10) Hodges, J.A. et al. The crystal field and exchange interactions in YbTiO. J. Phys.: Condens. Matter 13, 9301 – 9310 (2001).
- (11) Malkin, B.Z. et al. Optical spectroscopy of YbTiO and YTiO: Yb and crystal-field parameters in rare-earth titanate pyrochlores. Phys. Rev. B 70, 075112 (2004).
- (12) Ross, K.A. et al. Two-dimensional kagome correlations and field induced order in the ferromagnetic xy pyrochlore YbTiO. Phys. Rev. Lett. 103, 227202 (2009).
- (13) Thompson, J.D. et al. Rods of neutron scattering intensity in YbTiO: copelling evidence for significant anisotropic exchange in a magnetic pyrochlore oxide. Phys. Rev. Lett. 106, 187202 (2011).
- (14) Ross, K.A., Savary, L., Gaulin, B.D. &. Balents, L. Quantum excitations in quantum spin ice. Phys. Rev. X 1, 021002 (2011).
- (15) Ross, K.A. et al. Dimensional evolution of spin correlations in the magnetic pyrochlore YbTiO. Phys. Rev. B 84, 174442 (2011).
- (16) Chang, L.J. et al. Higgs transition from a magnetic Coulomb liquid to a ferromagnet in YbTiO. Nat. Comms. 3, 992 (2012).
- (17) Ross, K.A. et al. Lightly stuffed pyrochlore structure of single-crystalline YbTiO grown by the optical floating zone technique. Phys. Rev. B 86, 174424 (2012).
- (18) Hayre, N.R. et al. Thermodynamic properties of YbTiO pyrochlore as a function of temperature and magnetic field: Validation of a quantum spin ice exchange Hamiltonian. Phys. Rev. B 87, 184423 (2013).
- (19) D’Ortenzio, R.M. et al. Unconventional magnetic ground state in YbTiO. Phys. Rev. B 88, 134428 (2013).
- (20) Lhotel, E. et al. Mapping the first-order magnetic transition in YbTiO. arXiv:1405.2704 (2014).
- (21) Chang, L.-J. et al. Static magnetic moments revealed by muon spin relaxation and thermodynamic measurements in the quantum spin ice YbTiO. Phys. Rev. B 89, 184416 (2014).
- (22) Savary, L. & Balents, L. Coulombic quantum liquids in spin-1/2 pyrochlores. Phys. Rev. Lett. 108, 037202 (2012).
- (23) Shannon, N. et al. Quantum ice: a quantum Monte Carlo study. Phys. Rev. Lett. 108, 067204 (2012).
- (24) Benton, O., Sikora, O. & Shannon, N. Seeing the light: Experimental signature of emergent electromagnetism in a quantum spin ice. Phys. Rev. B 86, 075154 (2012).
- (25) Wan, Y. & Tchernyshyov, O. Quantum strings in quantum spin ice. Phys. Rev. Lett. 108, 247210 (2012).
- (26) Kozuki, K. Nagashima, T. & Hangyo, M. Measurement of electron paramagnetic resonance using terahertz time-domain spectroscopy. Opt. Express 19, 24950-24956 (2011).
- (27) Morris, C.M., Valdes Aguilar, R., Ghosh, A., Koohpayeh, S. M., Krizan, J., Cava, R. J., Tchernyshyov, O, McQueen, T. M., & Armitage, N. P. A hierarchy of bound states in the 1D ferromagnetic Ising chain CoNbO investigated by high resolution time-domain terahertz spectroscopy. Phys. Rev. Lett. 112, 137403 (2014).
- (28) Morris, C.M., Valdes Aguilar, R., Stier, A.V. & Armitage, N.P. Polarization modulation time-domain terahertz polarimetry. Opt. Express 20, 12303-12317 (2012).
- (29) Menzel, C., Rockstuhl, C & Lederer, F. Advanced Jones calculus for the classification of periodic metamaterials. Phys. Rev. A 82, 053811 (2010).
- (30) van Mechelen, J.L.M., van der Marel, D., Crassee, I & Kolodiazhnyi, T. Spin resonance in EuTiO. Phys. Rev. Lett. 106, 217601 (2011).
- (31) N.P. Armitage, arXiv:1311.5556, submitted Phys. Rev. B (2014).
- (32) Onoda, M., Mishchenko, A. S., & Nagaosa, N. Left-Handed Spin Wave Excitation in Ferromagnet. J. Phys. Soc. Jpn. 77 013702 (2008).
- (33) Applegate, R., Hayre, N. R.,Singh, R. R. P., Lin, T., Day, A. G. R., Gingras, M. J. P., Vindication of YbTiO as a Model Exchange Quantum Spin Ice. Phys. Rev. Lett. 109, 097205 (2012).
- (34) Torrance, JR., J.B., Tinkham, M., Magnon bound states in anisotropic linear chains. Phys. Rev. 187, 587 (1969).
- (35) Savary, L., Balents, L., Spin liquid regimes at nonzero temperature in quantum spin ice. Phys. Rev. B 87, 205130 (2013).

## Appendix A addendum

This work at JHU was supported by the Gordon and Betty Moore Foundation through Grant GBMF2628 to NPA and the DOE through DE-FG02-08ER46544. The crystal growth work at McMaster was supported by NSERC. We would like to thank C. Broholm and J. Deisenhofer for helpful conversations.

Competing Interests: The authors declare that they have no competing financial interests.

Correspondence: Correspondence and requests for materials should be addressed to NPA (email:npa@pha.jhu.edu).

## Appendix B Author Contribution

LP performed the THz experiments and analysis with assistance from CMM. SKK, AG, and OT performed the theoretical calculations. KAR, EK, BDG, provided the high quality single crystals. SMK cut and aligned the crystals. NPA conceived and directed the project. All authors contributed to discussions on data analysis and writing of the manuscript.

.

## Appendix C Supplementary Information

### c.1 Determination of Transmission Matrix in Linear Base

Here we provide the detailed Jones matrix analysis of the time domain terahertz spectroscopy measurement, and the analysis of the transmission matrix in the linear and circular bases. With the direction of the THz pulse propagation (wave vector k direction) set as the z axis, the electric field vector at any point along the optical path can be written as a column vector with complex entries.

The transmission of the THz pulse through the sample can be viewed as a transformation of the incident vector by a 2 by 2 complex transmission matrix:

In the Faraday geometry, since the sample is oriented with surface parallel to the 001 plane, the transmission matrix should preserve the 4-fold rotational symmetry of the cubic crystal along the surface normal. Requiring the transmission matrix to remain unchanged after a rotation of 90 along z direction shows that the transmission matrix must be antisymmetric, that is: , and . Now the transmission matrix in the linear laboratory frame (frame xyz) has only two independent component, and can be determined by:

(5) |

(6) |

With the x axis set in the direction of the polarization of the incident THz pulse, would be identically zero.

In a real time-domain terahertz spectrometer, a weak non-zero can be present due to small imperfections in alignment that give ellipticities in the incoming beam. For instance they may come from the modification of the incident pulse polarization due to the optical components along the beam path. With the full Jones matrix analysis presented above, a clean separation and calculation of and can be obtained.

### c.2 Determination of Transmission Matrix in Circular Base

As discussed in the main text, when magnetic excitations are present in the Faraday geometry, circular base is a more natural choice to analyze the data. Following the discussion in Refs. 26, 28 we present in detail the determination of transmission matrix in the circular base.

Suppose in the linear base, the electric field vectors are written as E and E for the incident and transmitted pulses; while in the circular base we have E and E. If we write the transformation matrix from the circular to linear base as , then we have:

in the linear base with the transmission matrix T, we have

then, follows from those equations, we have

(7) |

clearly now T is the new transmission matrix for the circular base.

With , denoting the electric field component in the right/left circular polarization channel, we have:

(8) |

combined with:

gives the transformation matrix:

(9) |

and the transmission matrix in the circular frame:

(10) |

This matrix is diagonal in the bases of right/left circular polarized light. The diagonal elements gives and , which can be calculated once and are known. Another way of obtaining and is by decomposing the incident and transmitted THz pulse into right/left circular bases with the transformation matrix , then taking the complex ratio of the corresponding spectra. Those two methods are equivalent.

In the Faraday geometry, the magnetic field is oriented in the same direction as the k propagation direction of the THz pulse (z direction in Fig. 1(b)). With the polarization modulation technique, the two orthogonal components (along x and y axes) of the transmitted THz pulse can be obtained in a single measurement. This allows the simultaneous determination of two complex transmission coefficients and . As mentioned above, for the Faraday geometry, a 4-fold rotational symmetry is expected along the surface normal direction with the FCC symmetry of the material, which requires the transmission matrix to have only two independent components:

(11) |

Thus the measurement of and completely determines the transmission matrix in the linear basis.

Following the discussion above, the transmission matrix can be diagonalized in the circular basis:

(12) |

where and are the complex transmission coefficient of the right/left circular polarized light.

An example of the transformation of basis is shown in Fig. SI1(a) and (b). In Fig. SI1(a) we plot the magnitude and from 0.1 to 0.7 THz measured from Sample A at 1.65 K with an applied field of 7 T. Fig. SI1(b) shows the same set of data in the circular base. Dashed lines in the figure mark the position of spectroscopic features. It is immediately evident that the circular base presents the data in a more straightforward fashion. It is clear that the pairs of dips and peaks in and corresponds to absorptions in the right/left circular channel, seen as the dips in and , respectively. The fact that those absorptions show up in the spectra of one type of circular polarized light only clearly demonstrates their nature as magnetic excitations. Another important point is that for a spectra as rich as the ones shown here, analyzing the data in an incorrect frame might lead to the misidentification of excitations. We should emphasize that the use of the time domain technique is vital for the types of analysis reported in this paper, as it allows the determination of both amplitude as well as phase of the terahertz pulses. The transformation between linear and circular bases is possible only when the phase information is retained.

From the above discussion, it is evident that by analyzing the transmission in the circular bases, we can separate contributions from different channels in the magnetic transitions. This provides a huge advantage when trying to understand the nature of those excitations, as we demonstrated in the main text of this paper. Again, measurement with a phase insensitive technique does not allow such a transformation of bases.

### c.3 Temperature dependence of the spectra in the Faraday geometry

We present in Fig. SI2(a) to (h) the frequency and field dependence of the transmission magnitude of the right and left circular polarized light at temperatures from 1.65 K to 13.3 K, measured from Sample A. As Sample A is thicker than Sample B, data from A show stronger absorption for b1 and b2. However, the lower transmission from A also lowers the contrast between the magnetic excitation and the background, thus diminishes the visibility of b5 and b6. Several features are evident from the data: 1) as temperature increases, all absorptions widen in frequency and lose intensities; 2) the slopes of those excitations have no obvious temperature dependence in the frequency - field plot; 3) b3, b5 and b6 disappear at a temperature around 8 K; 4) b4 disappears at temperature between 8 and 13 K; 5) b1 and b2 widen and merge together as temperature increases; 6) the spectra at 1.65 K and 4.3 K appear almost identical, except that at 1.65 K the absorptions are slightly sharper and stronger; 7) the higher frequency branches b5, b6 broaden significantly at low field; and 8) for the frequency range below 0.25 THz, the application of magnetic field increases the transmission in both the right and left channel, an effect which reduces as temperature is increased. To demonstrate the temperature dependence of the magnetic excitations in detail, measurements on Sample A are shown as a function of temperature while holding the applied field constant at 7 T. Fig. SI2(i) shows the ellipticity of the transmitted pulse at 7 T as calculated from , from 1.65 K to 40 K. The results show similar trends as those observed from the plots shown in Fig. 3(a) to (h). The temperature dependence is consistent with the accepted scale of the exchange constants; Fig. SI2(i) shows that the single ion regime is reached above 12 K ( ).

Excitations b1 and b2 are the strongest, and stay separated until around 25 K, above which temperature the magnetic absorption becomes considerably weaker. On the other hand, b3 and b4 disappear at temperatures as low as 8 and 10 K, respectively. The resonant frequency of all the excitations show insignificant temperature dependence; except for b3and b4 which showed a slight softening before disappearing.

As one can see from the main text, systematic errors in TDTS show up as horizontal lines in the intensity plots of the field dependent transmission data (e.g. Fig. SI2 (a) to (h)). Those type of systematic modulations of the data can be significantly reduced by plotting the ellipticity instead of the transmission amplitude. From the definition , it is clear that the common error in the measurement will be canceled out in computing . Thus it is a very useful quantity to look at when studying the subtle trends of the data, as shown in Fig. SI3.

### c.4 Field-Induced Transparency

The field induced transparency feature presents itself in the form of increased transmission in the frequency ranges away from magnetic excitations as an external magnetic field is applied (most noticeably below 0.25 THz). The observation of such effects indicates the existence of anomalous magnetic absorptions at low fields, which in our case, hints at the existence of string-like excitations.

In Fig. SI4, we present the data taken from sample A at 1.65 K in the Voigt geometry, plotted in the form of ratio of absorption under a fixed magnetic field () with that obtained from zero field scans (). Plotting the data in this fashion would reduce the spectroscopic features due to systematic noise in the measurement, thus makes it easier to observe the subtle trends. Since () is proportional to the absorption, a value larger than unity of this quantity would mean at that particular frequency, the transmission in field is lower than zero field. On the other hand, when the absorption ratio is lower than unity, there is a field induced transparency, as clearly shown in the data in Fig. SI4. In this case, the sample becomes more transparent as the field is increased and magnetic absorptions move to higher energy. It is also interesting to notice that, in the two branches visible in Fig. SI4 (c2 and c4), the width of c2 is relatively constant, while c4 broadens significantly as field is decreased. The original transmission spectra as a function of frequency at 1.65 K at different fields for the t in the Voigt geometry and RCP in the Faraday geometry are shown in Fig. SI5.

To understand the low frequency absorption better, we look at the temperature dependence in Fig. SI6, with data taken from the Faraday geometry. Shown in the figure is the absorption ratio for RCP and LCP between scans taken at 7 T and zero field at temperatures from 1.65 K to 20.4 K. It is clear from the data that the field induced transparency reduces as temperature increases, and is almost completely absent for temperature higher than 20 K, about the same temperature thermal fluctuations destroy spin correlations.

### c.5 Spin Wave Calculation

The spin wave calculations take into account the Zeeman term with an anisotropic g tensor and four possible bilinear spin interactions as discussed in Ref. 14. The classical spin wave modes are obtained in the linear spin wave theory by expanding the Hamiltonian in the large limit about one of the classical ground states, which is numerically determined. We use the interaction parameters as obtained from the recent fits to INS experiment on YbTiO as the starting values. From the detailed field dependent scans and high signal to noise in TDTS, we can refine the values for the exchange parameters. When fitting to the magnetic absorptions in the Voigt geometry above 3 T with the method of least squares, we find a best fit in units of meV of , , , and within the parameter search intervals which are twice as large as the uncertainty intervals in Ref. 14, which correspond to , , and .

There are a total of four non-degenerate spin wave modes, but two of these are predicted to have vanishing spectral weight as in the Faraday geometry. With 4 spins in the unit cell, there are 8 transverse spin components. Their harmonic interactions have the symmetry of the dihedral group , including two and three proper rotations (group ) times the inversion. These 8 transverse spin components contain all four 1-dimensional irreducible representations (, , , ) and two copies of the two-dimensional irreducible representations . Since the magnetic field components and transform as , only the modes are expected to be visible in optics in the Faraday geometry. The visible doublets have the highest and lowest frequencies. The modes invisible in the Faraday geometry are + and + with the former slightly higher than the latter for coupling constants of Ref. 14. The splitting of the visible modes is field-independent and equal to = 0.62 meV = 0.15 THz with the same coupling constants. The visibility of these modes can be understood heuristically in that one can see that in all four modes individual spins precess elliptically around the applied 001 field, with a large right-hand component and a small left-hand component. In the lowest-frequency mode the right-hand components add up in phase. In the highest-frequency mode they cancel out, leaving a small left-hand component. In the intermediate modes right- and left-hand components cancel out, leaving only a longitudinal oscillation.

The spin wave calculation results bear some resemblance and some notable difference with our TDTS results. We see a lower multiplet of 4 states that presumably corresponds to magnon-like (e.g. strings of length one) excitations. Note that all of them exhibit about the same -factor of approximately 3.2 as predicted. The spin-wave result also explains the otherwise anomalous appearance of the LCP polarized absorption in a 001 directed magnetic field. But there are also notable differences. In the lowest multiplet of states we see all 4 instead of 2 excitations predicted by theory. We also see two excitations ( and ) in another group of states at higher energy with approximately the twice the -factor of the lower multiplet. If the lower multiplet can be understood as spin-wave like excitations, then we can understand these higher energy excitations in the high field regime as two-magnon states or strings of approximate length two.

### c.6 Two-Magnon Calculation

As an attempt to understand the origin of the two-magnon peaks, we calculated the spectral function of the two-magnon continuum. To do this, we rewrote the Hamiltonian in a set of locally rotated axes such that at a given site the spin in the classical ground state points along the local z-direction. From this one can readily write an effective Hamiltonian in the subspace of one and/or two magnons at zero total momentum. has nearest-neighbor magnon hopping (deriving from terms like and in the full Hamiltonian), nearest-neighbor magnon interaction (deriving from the terms) and single-magnon to two-magnon tunelling terms (deriving from terms like and ). The spectral function of this Hamiltonian (with the optimal parameters obtained from linear spin wave theory mentioned in Section E) was obtained using Lanczos diagonalization. Convergence of the spectral function was checked by varying both the system size (achieved at ) and the number of Lanczos iterations (achieved at ).

At a given external magnetic field, four spectral functions were obtained corresponding to the RCP and the LCP channels in Faraday geometry and the X and the Y channels in Voigt geometry. Each spectral function has sharp peaks corresponding to the single-magnon excitations followed by a two-magnon continuum. The single-magnon energies thus obtained were found to differ somewhat from the classical spin-wave results for the lowest branch. This difference was larger at low fields and became smaller as the field was increased ( at 3 T to at 7 T). The peaks in the two-magnon continuum region were not found to agree well with the experimental peaks (except perhaps for the second peak in the RCP channel).

A comparison between the calculations and the experimental data is shown for one representative field in Fig. SI7. Here the parameter being plotted is the imaginary part of the dynamic susceptibility . The experimental spectra are obtained from /, with d being the thickness of the sample. Spectra obtained this way represent the dissipative part of material parameter, which in the particular case studied here, is proportional to . Note at frequencies above 0.7 THz, absorption due to optical phonons become substantial, and causes additional contribution to the experimental spectra. As shown in the figure, the calculation provides qualitative description of the data.

We also examined the variation in the single-magnon energies at as one moves away from the optimal point in the four dimensional parameter space, which is achieved from fitting described in Section F. To do this, we obtained the matrix that connects a small change in the coupling constants to the corresponding change in the single magnon energies (ordered such that is the highest mode) at a given external field and then looked at its eigenvectors and eigenvalues. This matrix has an eigenvector with a very small eigenvalue , which implies that the single magnon energies are not very sensitive to small changes of the coupling constants in that particular direction. This direction translates to . This result was also found to be quite independent of the magnitude of the applied external field. Thus, in the present calculation, the value of the (or equivalently, the ) parameter is not very well determined from fitting the single-magnon energies to spin wave calculations.

We wanted to check whether the insensitivity of the magnon energies along the direction listed above persists at finite momenta also. But the full Lanczos diagonalization was carried out only at . So, here we will just quote what linear spin wave theory says at finite momenta. The linear model was found to predict a greater shift in the single-magnon energies with a change in . More specifically, the vector at was found to be . The components of the corresponding vectors at and differed by at most from these values, indicating that using the full magnon dispersions will not significantly improve the uncertainty in .