Hierarchy of exchange interactions in the triangular-lattice spin-liquid YbMgGaO$_{4}$

The spin-1/2 triangular lattice antiferromagnet YbMgGaO$_{4}$ has attracted recent attention as a quantum spin-liquid candidate with the possible presence of off-diagonal anisotropic exchange interactions induced by spin-orbit coupling. Whether a quantum spin-liquid is stabilized or not depends on the interplay of various exchange interactions with chemical disorder that is inherent to the layered structure of the compound. We combine time-domain terahertz spectroscopy and inelastic neutron scattering measurements in the field polarized state of YbMgGaO$_{4}$ to obtain better microscopic insights on its exchange interactions. Terahertz spectroscopy in this fashion functions as high-field electron spin resonance and probes the spin-wave excitations at the Brillouin zone center, ideally complementing neutron scattering. A global spin-wave fit to all our spectroscopic data at fields over 4T, informed by the analysis of the terahertz spectroscopy linewidths, yields stringent constraints on $g$-factors and exchange interactions. Our results paint YbMgGaO$_{4}$ as an easy-plane XXZ antiferromagnet with the combined and necessary presence of sub-leading next-nearest neighbor and weak anisotropic off-diagonal nearest-neighbor interactions. Moreover, the obtained $g$-factors are substantially different from previous reports. This works establishes the hierarchy of exchange interactions in YbMgGaO$_{4}$ from high-field data alone and thus strongly constrains possible mechanisms responsible for the observed spin-liquid phenomenology.


I. Introduction
Quantum spin-liquids (QSL) are exotic states of matter in which spins are highly correlated but remain dynamic down to zero temperature due to strong quantum fluctuations [1,2]. Many distinct QSL states have been proposed theoretically [3,4] and classified according to their non-local (topological) properties [5]. Their detection, however, remains a central challenge for condensed matter physics [6], and relies on the presence of quantum entanglement in their ground-state and fractional quasiparticles in their excitation spectra. While the former can be checked by numerics [7], the latter can be experimentally detected by thermodynamic techniques [8][9][10] or spectroscopic probes such as inelastic neutron scattering [11][12][13][14][15][16][17] and electron-spin resonance [18][19][20][21].
The newly discovered rare-earth triangular-lattice antiferromagnet YbMgGaO 4 [53,54] appears to fulfill precisely this promise. The magnetic Yb 3+ ions carry effective spin-1/2 moments in a symmetry environment allowing anisotropic exchange interactions [41,47] in the absence of antisymmetric (Dzyaloshinsky-Moriya) terms and magnetic defects, both of which are present in other two-dimensional QSL candidates such as Herbertsmithite [55][56][57]. The immediate availability of single crystals [54] uncovered a QSL phenomenology in YbMgGaO 4 characterized by the absence of spin ordering or freezing down to T = 100 mK in muon spin relaxation measurements [58], much lower than the Curie-Weiss temperature θ W ≈ −4 K, and a power-law behavior for the magnetic specific heat at low temperatures [54,59]. Perhaps the strongest evidence for a QSL in YbMgGaO 4 came from inelastic neutron scattering measurements in zero field that unraveled a broad continuum of magnetic excitations across the entire Brillouin-zone [60][61][62]. This continuum has been interpreted as composed of fractional excitations from a U(1) QSL state with spinon Fermi surface [42,44,61] or from a RVB-like state [62].
The absence of a magnetic contribution to the thermal conductivity [59], however, appears difficult to reconcile with the spinon Fermi surface interpretation. Additionally, the disordered occupancy of the inter-triangular layers by Mg 2+ and Ga 3+ ions [53] appears to affect profoundly the Yb 3+ ions with broadened crystal electricfield (CEF) levels [60,63], a distribution of g-tensors [63], and a broadened magnetic excitation spectrum at high fields [60]. The impact of disorder on the YbMgGaO 4 exchange Hamiltonian and, therefore, whether the ground state is a QSL or not, remains an outstanding issue [47]. In fact, the nature of the dominant exchange interactions in YbMgGaO 4 is also controversial. While the overall planar anisotropy is clear from susceptibility measurements [54,61], both antiferromagnetic next-nearestneighbor terms (J 2 ) [60] and nearest-neighbor anisotropic off-diagonal exchanges (so-called J ±± 1 and J z± 1 ) [43,64] have been proposed as extensions from the XXZ model. Comprehensively determining the exchange interactions in YbMgGaO 4 is of fundamental importance in deciphering the nature of its ground state.
Here, we combine time-domain THz spectroscopy (TDTS) with inelastic neutron scattering (INS) measurements in the field-polarized state of YbMgGaO 4 to extract values of the exchange interactions from a global fit to spin-wave spectra measured for different field directions and scattering wave-vectors including the Brillouinzone center (Γ-point). Previous high-field INS measurements were limited to wave-vectors around the antiferromagnetic zone boundary [60] while previous X-band electron spin resonance (ESR) measurements [54] were intrinsically limited to small fields ( 0.4 T) below the field-polarized regime. In contrast, high resolution TDTS (e.g. [65]) functioning as high field ESR on magnetic insulators allows an accurate determination of magnetic excitations in fields comparable to the saturation field in YbMgGaO 4 [54,60]. Since the wavelength of THz radiation is much greater than lattice constants, TDTS measures excitations in the first Brillouin zone with zeromomentum transfer (i.e. the q = 0 response), which is impossible with neutron scattering. Inclusion of the highfield TDTS data allows a substantial refinement of the Hamiltonian parameters compared to Ref. [60]. Moreover, an analysis of TDTS linewidths further constrains the anisotropic exchange interactions. In the context of prevailing models our results suggest that both NNN and off-diagonal anisotropic exchange interactions are present in YbMgGaO 4 -with pseudo-dipolar terms sub-leading compared with the XXZ part of the model -and that this nominal set of exchanges may lie closer to a spinliquid regime than previously thought [47,60]. When combined with the likely presence of exchange disorder, the extent to which remains to be determined, our work strongly constrains the nature of the YbMgGaO 4 ground state and points towards a more intricate underlying scenario than originally reported [61].

A. Experimental methods
The YbMgGaO 4 crystal (space group R3m) used in this work was grown by the floating-zone technique (see Supplementary Information (SI) Section 1 for details) as reported in Ref. [60] and cut with a diamond blade to present a c-axis facet, where the c-axis is orthogonal to the triangular ab plane of Yb 3+ ions. TDTS measurements were performed in a custom built polarization modulation setup with a frequency range of 0.2 to 2 THz (0.83 to 8.3 meV) (see SI Sec. 2). The complex THz transmission of a 4 × 4 × 1 mm 3 single crystal was measured down to 1.6 K (see SI Sec. 7) with external fields up to H = 6.8 T in both the Faraday (k H) and Voigt (k ⊥ H) geometries, where k is the direction of light propagation. In both cases, the THz pulse ac magnetic field, h, was applied along a * with h ⊥ H and thus the H c and H a orientations were probed in Faraday and Voigt geometries, respectively.
INS experiments were performed on the cold neutron chopper spectrometer [66] (CNCS) at the Spallation Neutron Source (SNS), and on the cold triple-axis spectrometer CTAX at the High Flux Isotope Reactor (HFIR), both at Oak Ridge National Laboratory. Experiments were performed on the same crystals as Ref. [60] with a magnetic field of H = 7.8 T along the crystal c-axis and T ≈ 0.06 K on CNCS, and at H = 10.8 T along the crystal a-axis and T ≈ 0.32 K on CTAX. The INS data taken at CNCS is the same as published in [60] but now analyzed together with the TDTS data. Given the broad spectra, even at high-fields, the energy of magnetic excitations for a given wave-vector was determined by fitting the maximum in scattering intensity (Ref. [60] and SI Sec. 3). Real (a) and imaginary (b) parts of the magnetic susceptibilityχ(ν) for right-circularly polarized light at different fields H dc c at 5 K.χ(ν) is obtained by referencing the TDTS data to the spectra at 100 K. Spectra in (a) are offset vertically by 0.02 for clarity. (c) Resonant frequency (νc) at 5 K vs magnetic field for magnetic excitation peaks in (b) in the Faraday geometry H dc c. Yellow circles represent the data. Red line is a linear best fit to extract g . The black dashed line represents the global fit to the data (see Fig 3). Error bars in (c) are smaller than the marker size.

B. Analysis of the TDTS data
In the field-polarized state in the Faraday geometry, the linearly polarized THz pulse becomes elliptically polarized as it passes through the sample due to spins precessing around the applied field direction. The complex transmission is represented by a 2 × 2 Jones matrix. Due to the three-fold rotational symmetry of the lattice, this reduces to an anti-symmetric matrix for the transmission of a linearly polarized pulse [67]. By diagonalizing this antisymmetric matrix, we can convert the transmission matrix from the linear basis into a circular basis (see SI Sec. 2) which naturally corresponds to eigenstates of the transmission in the Faraday geometry [50]. In this manner we can separate the transmission of a left (LCP) and right (RCP) circularly polarized THz pulse. Fig. 1a and Fig. 1b show such transmission magnitudes for LCP and RCP THz light, respectively, and as a function of frequency and magnetic field at a temperature of T = 5 K. The bright yellow feature in Fig. 1b shows that only one particular helicity of the light (RCP) is strongly absorbed. This indicates the presence of well defined spin-wave excitations with energy linearly dependent on the applied field. In a field polarized regime in the Faraday geometry, the direction of spin precession is determined by the orientation of the applied H dc field and it is therefore natural to expect an absorption for only one helicity and not the other. Indeed, our experiments confirm that when the polarity of the H dc field is flipped, only LCP light is absorbed (see SI Sec. 4). This confirms that the absorption line observed in Fig. 1b originates from spin-wave excitations in the field-polarized regime of YbMgGaO 4 . The observed absorption is the only feature in the data that is affected by the H dc field and temperature, suggesting the featureless background has a non-magnetic origin and is likely the low-energy tail of a crystal field level absorption (see SI Sec. 5).
The featureless nature of the background allows us to extract the frequency dependent magnetic susceptibilitỹ χ(ν) as follows. The complex transmission of a particular eigenpolarization (LCP or RCP light) is given by the μ is the complex index of refraction,˜ is the dielectric constant,μ = 1+χ and d is the sample thickness. We determineñ using the Newton-Raphson method and then isolateχ by measuring the sample at a reference temperature at which the spectrum does not show any signatures of the spin-wave absorption [51]. At this reference temperature (100 K in this case, see SI Sec. 6) χ can be taken to be zero and soñ 100K = √˜ . Thus, the low temperature magnetic susceptibility is given bỹ show the real and imaginary parts ofχ(ν), respectively, at different fields and T = 5 K in the Faraday geometry. Peaks in Imχ(ν) correspond to the spin wave excitations in the q = 0 limit. By fitting the data at each field with a Lorentzian (see SI Sec. 7) we extract the resonant energy (E) of the spin-wave absorption. The resulting E vs H dc plot is shown in Fig. 2c. The peak widths and resonant frequencies show little temperature dependence between 1.6 K and 40 K (see SI Sec. 7). Similar analysis (see SI Sec. 8) is done for TDTS in the Voigt geometry to extract the spin-wave energies for field H dc along the a-axis (reported in Fig. 3). This TDTS data is then combined with INS data to extract the exchange interactions in a global fit, as discussed below.

C. High-field spin-wave theory analysis
An effective spin-1/2 Hamiltonian relevant for YbMgGaO 4 has been given by [41,42,54,60]:  Table I). Green circles at the Γ-points are the TDTS data from where S ± = S x ±iS y , g and g ⊥ are components of the gtensor parallel and perpendicular to the c-axis, complex numbers γ ij are defined in [54] and brackets <> and represent nearest and next-nearest neighbor pairs respectively. The exchange interactions J zz and J ± originate from the standard XXZ model, with subscripts 1 and 2 indicating the nearest and next-nearest neighbor interactions respectively. We also include symmetry allowed bond-dependent interactions J ±± and J z± (also known as pseudo-dipolar interactions), which have a spin-orbit origin. The pseudo-dipolar interactions between nextnearest neighbors are neglected since they are supposed to be small [69,70] . For large applied fields and assuming a field-polarized state, the above Hamiltonian can be solved to yield the spin wave dispersions for external fields along the c-and a-axes (see SI). By setting q = 0, the spin wave excitations for the Faraday and Voigt geometry can be simplified to: According to Eq. 2, the field dependence in the Faraday geometry is particularly simple and depends only on g . From a simple linear fit to the data in Fig. 2c we obtain g = 3.81 (4). Note that the effective g-tensor here is a property of a single Yb 3+ ion and as such is independent of the exchange interactions [48]. This allows us to treat g as a fixed parameter in the subsequent global fit. In contrast, E(H x ) is not linearly related to g ⊥ (Eq. 3) and so g ⊥ must be extracted simultaneously with the exchange constants of Eq. 1. One expects the present measurements of the g-factor to be considerably more accurate than previous X-band ESR results [54] due to greater than a ten-fold increase in the magnitude of the magnetic fields. The expressions for the dispersions with q = 0, as relevant for neutron-scattering measurements, are given in Refs. [43,60] and SI Sec. 3.

A. Fits to spin-waves energies
We now turn to our central result, which is to refine the parameters of YbMgGaO 4 by combining all the data at hand. In Fig. 3, we show the spin wave dispersion ob-tained through INS [ Fig. 3(a-d)] as well as the TDTS data in Voigt geometry [ Fig. 3(e)]. The entire dataset is fit simultaneously to the spin-wave dispersions in the field-saturated state obtained from the Hamiltonian in Eq. 1. There are six target parameters for this global fit, which are J zz 1 , J ± 1 , J ±± 1 , |J z± 1 |, the ratio J 2 /J 1 and g ⊥ with a fixed g = 3.81. Note that only the magnitude of J z± 1 can be determined from the spin-wave dispersions (see SI Sec. 9). The signs of the other J terms obtained from the fit is verified by starting from both negative and positive initial values. As spin-space anisotropy is primarily a property of the effective spin-1/2 doublet of Yb 3+ , we adopt the same overall XXZ exchange anisotropy for both nearest neighbor and nextnearest neighbor interactions, i.e. J zz 2 /J ± 2 = J zz 1 /J ± 1 . This reasonable assumption helps in reducing the size of the parameter space for the global fit.
We obtain an excellent fit to the data using the above model [see lines in Fig. 3 and spread in Fig. 4] and report our global fitting parameters in Tab. I (Model C). The χ 2 goodness of fit and correlation plots of these parameters are shown in the SI Sec. 10. We note that the above fit is performed with the constraint |J zz 1 | > |J z± 1 |. Without this, we obtain |J z± 1 | = 0.45 which we regard as unphysical as it is nearly three times larger than J zz 1 and also considerably outside the bounds set by the linewidth analysis performed below.
In previous works, two different sets of parameters have been proposed for the exchange interactions of YbMgGaO 4 . The first set [60], model A in Tab. I, includes both nearest and next-nearest neighbor interactions but sets J z± 1 = 0 based on previous X-band ESR results [54]. The second set [43], model B in Tab. I, includes only nearest-neighbor interactions, with the J 2 = 0 constraint argued as a consequence of the localized nature of Yb 3+ 4f electrons [43,64]. Both models used the g-factors determined from low frequency X-band ESR. Model C is the new best fit from our spin-wave analysis. To compare these models, we plot in Fig. 4 the experimental vs. calculated spin-wave energies for each of model A, B and C. Clearly, the points obtained for our model C lie significantly closer to the E exp = E calc line than models A and B, highlighting the advantage of our enriched datasets over previous works.
Our global analysis yields finite NNN interactions with J 2 /J 1 = 0.18(7) and a slightly sub-dominant pseudodipolar J ±± 1 ≈ 0.07 meV interaction. The large fitting error bars on some parameters of the model reveal the high degree of correlation between J 2 and J ±± 1 and the very weak sensitivity of our fit to |J z± 1 |. This is apparent in the χ 2 goodness of fit 2D plots (see SI Sec. 10) in the parameter spaces of J ±± 1 versus |J z± 1 |, and J ±± 1 versus J 2 /J 1 , respectively. The poor sensitivity to |J z± 1 | is also clear from the analytical expression for the fieldpolarized spin-wave dispersion with field along the a-axis (see SI Sec. 9).
Given the correlation between J ±± 1 and J 2 , it is natural to analyze how our results change by enforcing the  Tab. I). This model gives only a slightly worse fit to the spinwave data than model C, as can been seen from a plot of the experimental against calculated spin-wave energies in Fig. 4. While a unique set of exchange parameters cannot be obtained from the above analysis, a definitive hierarchy of interactions nevertheless emerges. It yields easy-plane XXZ terms for YbMgGaO 4 with J zz 1 = 0.15(1) meV, ∆ = J zz 1 /2J ± 1 = 0.9(1), a subleading pseudo-dipolar term J ±± 1 /J zz 1 ≤ 1 and a small NNN exchange J zz 2 = 0.03(1) meV.

B. Further constraints: TDTS lineshapes
To further constrain these parameters, we require additional information. For that, we analyzed the linewidths of the spin-wave peaks in the TDTS data of Fig. 2. As noted above, TDTS functions here as high-field ESR [56]. ESR linewidths are strongly sensitive to the magnitude of the off-diagonal anisotropic interactions J ±± Here we use this formalism to show that the observed spin-wave linewidths cannot be reconciled with the large anisotropic exchange terms suggested by model B* and therefore J 2 must be non-zero. In Fig. 5, we plot the deviation between the calculated and the experimental spin-wave resonance linewidths in both the Faraday ∆ f and Voigt ∆ v geometry as a function of |J z± 1 | and |J ±± 1 |. We analyze the case of J 2 = 0 that is relevant for model B*. We define a function , that is minimal along the yellow contour in Fig. 5, for which a rough analytical description is |J z± 1 | 2 + |J ±± 1 | 2 ≈ 0.04(1) meV. This calculation assumes that all the broadening of the THz data comes from exchange anisotropy with no distribution of g-factors beyond the values determined above. In fact, a relatively large spread in the latter values, ∆g /ḡ ≈ 0.3 and ∆g ⊥ /ḡ ⊥ ≈ 0.1, has been demonstrated by a careful analysis of CEF excitations linewidths [63] and must also play a role in the lineshapes observed here. Our analysis of the high-field TDTS lineshapes therefore provides an upper bound on any off-diagonal anisotropic exchanges |J z± 1 | and |J ±± 1 |. The parameters determined in model B* are incompatible with this bound and so model B* is ruled out as a realistic description of the exchange parameters of YbMgGaO 4 .
Taken together, our results strongly favor an easyplane XXZ scenario for YbMgGaO 4 with the combination of finite J 2 and relatively small pseudo-dipolar exchanges. Thus, the refinement of model C with a range of best possible parameter values we found for YbMgGaO 4 can be summarized as a model C*: (7), |J ±± 1 | 2 + |J z± 1 | 2 0.05 meV. The above analysis assumes that finite J 2 will also give an additional contribution to the linewidth and hence the bounds on |J ±± 1 | and |J z± 1 | for J 2 = 0 represent the maximum values that these parameters can take. We note that the bound determined for |J z± 1 | justifies the earlier constraint of |J zz 1 | > |J z± 1 | in determining the parameters for model C in Table I.

C. Discussion
Experiments using other techniques, such as unpolarized [60] and polarized [64] neutron diffraction in zero field and low-temperatures may help to constrain the values further. Indeed, classical Monte-Carlo simulations for the instantaneous spin structure-factor S(Q) reveal that either large J ±± 1 or relatively small J 2 ∼ 0.2J 1 yield correlations peaked at the M-point of the triangular-lattice Brillouin zone [60]. Maximal correlations at the M-point have been reported in all neutron scattering investigations to date [60][61][62]64]. In this regard, we note that Monte-Carlo simulations with a large J ±± 1 yield a strong modulation in S(Q) from the first to the second Brillouin zone [60]. To the best of our knowledge, such a modulation has not yet been observed experimentally, which appears consistent with the relatively small J ±± 1 term indicated by our spin-wave resonance linewidth analysis. In general, the sensitivity of S(Q) to spin anisotropy arises because neutrons scatter only from spin components perpendicular to Q. However, this projection has not always been fully included in theoretical calculations of the scattering pattern, which may help to account for the somewhat larger values of J ±± 1 proposed in Refs. [43,64]. Future magnetization studies at very low temperatures will be another important test for the presence of further-neighbor or anisotropic exchange interactions in YbMgGaO 4 . The presence of multiple phase transitions and magnetization plateaus in a moderate magnetic field [73] is expected in the case of a XXZ triangularlattice antiferromagnet [74,75]. This may also help to constrain the possible values of J 2 [76] and anisotropic exchanges, although exchange disorder may be very efficient at suppressing these features [77]. Due to their non-spin-conserving nature, dominant pseudo-dipolar interactions should manifest as an asymptotic approach to saturated magnetization, as recently observed in the candidate Kitaev spin-liquid α-RuCl 3 [78]. Finally, we note that the role and extent of CEF and exchange disorder remains an outstanding issue in YbMgGaO 4 with possible impact ranging from the mimicry of a spin-liquid [47] to disorder-induced entanglement [79,80].

IV. Conclusion
In conclusion, combining time-domain terahertz spectroscopy and neutron scattering in the high-field regime of YbMgGaO 4 yields the strongest constraints to date on possible mechanisms for the observed spin-liquid phenomenology. We note that a QSL regime is predicted for the spin-1/2 J 1 -J 2 Heisenberg triangular-lattice antiferromagnet for 0.06 J 2 /J 1 0.19 [36][37][38]. While a previous analysis obtained J 2 /J 1 = 0.22 [60], the present work yields a slightly smaller J 2 /J 1 = 0.18 (7). Additionally, the ratio J zz 1 /J ± 1 = 1.75 (5) suggests that YbMgGaO 4 is more spin-isotropic than previously thought (J zz 1 /J ± 1 = 1.16 [60]) with J zz 1 /J ± 1 = 2 and 0 the Heisenberg and XY limits, respectively. We note, however, that subleading pseudo-dipolar interactions are also necessary to best explain our data. When combined with the likely presence of exchange disorder due to Mg 2+ and Ga 3+ disorder, this makes the spin-liquid mimicry [47] or the J 1 -J 2 quantum spin-liquid [36][37][38] mechanisms serious contenders to the various scenarios proposed to explain the physics of YbMgGaO 4 thus far [44,61,62].

Sample Preparation
Polycrystalline samples of YbMgGaO 4 were synthesized by a solid state method. Stoichiometric ratios of Yb 2 O 3 , MgO, and Ga 2 O 3 fine powder were carefully ground and reacted at a temperature of 1450 • C for 3 days with several intermediate grindings. Single-crystal samples of YbMgGaO 4 were grown using the optical floatingzone method under a 5 atm oxygen atmosphere. The best single crystals were obtained with a pulling speed of 1.5 mm/h, and showed [001] surfaces after several hours of growth.

Time-domain THz (TDTS) setup and analysis
Time-domain THz spectroscopy was performed in a home built setup at Johns Hopkins University. THz pulses, with a bandwidth between 0.2 to 2 THz, were generated by a photoconductive antenna (Auston switch) -emitter -upon illumination by an infrared laser and then detected by another Auston switch (receiver). The sample was mounted on a 4 mm aperture. The electric field profiles of the THz pulses transmitted through the sample and an identical bare aperture were recorded as a function of time and then converted to the frequency domain by Fast Fourier Transforms (FFTs). By dividing the FFTs of the sample and aperture scans, we can obtain the complex transmission of the sample.
In linear basis, the electric field vectors are represented as E i and E t for the incident and transmitted THz pulses respectively. Similarly, we define E i,cir and E t,cir in the circular basis. Λ represents the transformation from circular to linear basis and we have E i = ΛE i,cir and E t = ΛE t,cir where Λ = 1 T cir = Λ −1 T Λ is the transmission matrix in circular basis. The transmission matrix of a sample with three-fold rotational symmetry in Faraday geometry is fully antisymmetric [1], i.e. T xx = T yy and T xy = −T yx . Therefore, the transmission matrix in the circular basis is diag- where T r and T l denote transmission matrix of right and left polarized light, respectively. With the polarization modulation technique as discussed in [2], T xx and T xy can be obtained simultaneously with a fast rotator. This allows us to determine two the complex transmission coefficients T xx and T xy simultaneously and convert these to T r and T l .

Fitting INS data
The inelastic neutron scattering lineshapes obtained on CTAX are presented in Fig. S1 with the scattering wave-vector spanning the MΓ 2 direction with a magnetic field of 10.8 T applied along the a-axis of the crystal. The position of the most intense peak is indicated by vertical mark.

Faraday transmission for negative fields
To confirm the nature of the observed YbMgGaO 4 excitations in the Faraday geometry, we performed TDTS for both positive and negative magnetic field polarities. Here positive field refers to the field direction described in Fig. 1 of the main text. There the absorption appeared for only right circularly polarized (RCP) light and not for left circularly polarized (LCP) light. By reversing the field polarity, LCP light is now absorbed and not RCP light, as shown in Fig. S2. Compare this plot to Fig. 1 in the main text.

FTIR spectra
We also performed Fourier transform infrared spectroscopy (FTIR) on the YbMgGaO 4 single crystal at T = 5 K. The reflectivity of the sample was measured using a commercial FTIR spectrometer in the far and mid infrared spectral ranges, i.e. 50 cm −1 to 7000 cm −1 (1.5 THz to 250 THz). The complex optical conductivity is calculated from the reflectivity spectrum using Kramers-Kronig constrained variation dielectric function fitting [3]. The real part of the optical conductivity σ 1 is shown in Fig. S3. The peaks in σ 1 match with the crystal field excitations observed by neutron scattering at 9.2 THz (38 meV) and 17.7 THz (61 meV) and 23.5 THz (97 meV). The low energy tail of the lowest energy crystal field excitation gives rise to the background observed in the TDTS measurements.

TDTS -Temperature dependence
In the main text we isolate the magnetic susceptibilitỹ χ(ν) by measuring the sample at a reference temperature at which the spectrum does not show any signatures of the spin-wave absorption. We show the measured trans- mission spectrum in TDTS for the YbMgGaO 4 sample at different temperatures between T = 5 K and T = 100 K at H = 6.8 T in the Faraday geometry (Fig. S4). The observed spin-wave absorption is the only feature in the data that is affected by the static field and temperature, suggesting that the featureless background has a nonmagnetic origin. We do not observe any spin-wave absorption at T = 100 K which allows us to use this as the reference temperature. We also measured YbMgGaO 4 at T = 1.6 K but with a thinner sample because the original sample is too thick for Terahertz light to pass. The peak positions don't shift as we lower the temperature to T = 1.6 K, hence the g-factor doesn't change either (Fig. S5).

Fitting TDTS data
The spin-wave peak at each field and temperature in Imχ(ν) was fitted to a Lorentzian: to extract the resonant frequency ν c and the FWHM width ∆. A representative fit is show in Fig. S6a. The  FIG. S5. (a) The imaginary part of the susceptibility Imχ(ν) at +6.8 T for RCP light in the Faraday geometry taken at 1.6 K and 5 K. (b) Spin wave excitation peaks versus magnetic fields at 1.6 K and 5 K, showing the same g-factor. value of ∆ as a function of temperature and field is shown in Fig. S6b and c, respectively. We find that ∆ is mostly independent of temperature and field below ∼ 40 K.

TDTS data in the Voigt geometry
The TDTS spectra taken in the Voigt geometry is shown in Fig. S7. TDTS transmission specta were measured at T = 5 K for various fields where the THz pulse H dc ⊥ c geometry. The spectra at each field are referenced to the spectra taken at 100 K for the same geometry. The resulting imaginary part of the magnetic susceptibility Imχ(ν) as calculated in the linear basis is FIG. S7. The imaginary part of the susceptibility Imχ(ν) for linearly polarized light at different fields H dc a at 5 K as obtained by referencing the TDTS data to the spectra at 100 K.
Recall that model B* gives only a slightly worse fit to the spin-wave data than model C, as can been seen in Fig. 4 of the main text. To distinguish between the two, we determine the discrepancy between the calculated and the experimental linewidths in both the Faraday ∆ f and Voigt ∆ v geometry as follows: we define a function R p = Fig. 5 in the main text shows R p , calculated using the expressions for ∆, M 2 and M 4 , for the case of model B* (J 2 = 0) as a function of |J z± 1 | and |J ±± 1 |. R p is minimized along the yellow contour, and hence within a model where J 2 = 0 one expects the anisotropic exchange parameters to fall in this yellow region. However, note that the model with J 2 = 0 only fits the neutron and THz data in the global fit if J ±± 1 is about 0.13(2) meV (model B* in Table I of the main text). But if J ±± 1 was as large as 0.13(2) meV with J 2 = 0, then a much broader linewidth in the THz data would have been observed. This gives further evidence that model B*, with no next-nearest neighbor interaction is inappropriate.