Constraining the Parameter Space of a Quantum Spin Liquid Candidate in Applied Field with Iterative Optimization

The quantum spin liquid (QSL) state is an exotic state of matter featuring a high degree of entanglement and lack of long-range magnetic order in the zero-temperature limit. The triangular antiferromagnet YbMgGaO4 is a candidate QSL host, and precise determination of the Hamiltonian parameters is critical to understanding the nature of the possible ground states. However, the presence of chemical disorder has made directly measuring these parameters challenging. Here we report neutron scattering and magnetic susceptibility measurements covering a broad range of applied magnetic field at low temperature. Our data shows a field-induced crossover in YbMgGaO4, which we reproduce with complementary classical Monte Carlo and Density Matrix Renormalization Group simulations. Neutron scattering data above and below the crossover reveal a shift in scattering intensity from M to K points and, collectively, our measurements provide essential characteristics of the phase crossover that we employ to strictly constrain proposed magnetic Hamiltonian parameters despite the chemical disorder. Constrained exchange parameters further suggest the material's proximity to the QSL state in the clean limit. More broadly, our approach demonstrates a means of pursuing QSL candidates where Hamiltonian parameters might otherwise be obscured by disorder.


I. INTRODUCTION
The search for an experimental realization of the quantum spin liquid (QSL) state has been central to condensed matter physics for decades [1,2]. Anderson first proposed triangular antiferromagnets to exhibit the QSL state via the resonating-valence-bond (RVB) [3], and the hunt for promising triangular lattice candidates systems is still ongoing [4][5][6][7][8][9][10][11][12], though QSLs have been proposed for a variety of other models and systems, such as the kagome (corner-sharing triangular), honeycomb, and pyrochlores [1,2,13]. Indeed, recent years have offered an explosion in possible theoretical models for QSLs, and a wide range of new candidate material hosts, but few new experimental methods to detect them [14]. With the lack of a "smoking gun," a variety of complimentary experimental methods, such as magnetic susceptibility and neutron scattering, are needed to establish the possible presence of the QSL state. Equally important, numerical techniques such as DMRG are required to validate interpretation of the data as evidence of a QSL. Thus the parameters describing the magnetic Hamiltonian must be determined before the question of whether or not a material hosts the QSL ground state can be answered. For very recent examples, consider work on the exactly-solvable Kitaev system α-RuCl 3 [15], or on the hyperkagome PbCuTe 2 O 6 [16], or in the study of the triangular antiferromagnet YbMgGaO 4 [17], the focus of this work.
Efforts to understand the impact of chemical disorder due to site-mixing between non-magnetic Mg 2+ and Ga 3+ have enriched the discussions [18,19,21,22,[40][41][42]. Inelastic neutron scattering measurements suggested that crystalline electric field (CEF) levels are broadened by a distribution of ytterbium-oxygen bond distances due to the Mg 2+ /Ga 3+ site mixing, leading to a distribution of effective spin-half g-factors and consequently broadened low-energy magnetic excitations in the polarized state [42]. Theoretical work suggests that, barring disorder in the charge environment from Mg 2+ /Ga 3+ site mixing, YbMgGaO 4 should have a collinear/stripe ground state [19,24], while other calculations have suggested that the system could exhibit either a striped or 120 • ordered state [32]. While ac susceptibility measurements down to 60 mK imply spin freezing [36], inelastic neutron scattering suggests roughly 16% of spins frozen [43], and dc susceptibility measurements provide evidence of dynamic spins down to 40 mK with only approximately 8% frozen [37]. The low fraction of frozen spins measured in Ref. 43 has previously been attributed to a high degree of frustration [36]. Heat capacity measurements [18] and µSR studies [45] at 60 mK similarly oppose a spin freezing scenario. Very recently, µSR studies were carried out down to temperatures as low as 22 mK without observing ordered or disordered static magnetism [47], further opposing interpretations of the ground state as a spin glass.
Precisely measured exchange parameters are needed to evaluate the possibility of a QSL state. Conventionally, in systems with complex Hamiltonians, exchange parameters are extracted from fits of the single-magnon dispersion obtained from inelastic neutron scattering. Quenched disorder presents a significant obstacle to this endeavor, introducing broadening observed in the dispersion in YbMgGaO 4 [43]. Though a timed-resolved THz spectroscopy study significantly refined constraints [17], difficulties are compounded by the distribution of effective g-factors [25]. Furthermore, the bond-dependent J z± exchange parameter, important in the ground states at fields below saturation, is absent from the analytical expression of the dispersion in the polarized state obtained from the linear spin wave approximation [25]. In order to properly evaluate proposed Hamiltonian parameters, systematic measurements at low and intermediate field are required. More generally, magnetic field is an effective and extensively used tuning parameter in the study of QSL candidates [51], such as to suppress a competing ordered state [52], or to suppress spin fluctuations and measure spin waves in high field [43], and a number of recent triangular antiferromagnetic materials show interesting phase transitions or features in applied field [12,53]. A detailed study of YbMgGaO 4 's field dependence, including at low and intermediate fields, is a natural next step to understand the ground states of the system, and to provide further data for comparison to proposed models.
To that end, this work presents new data revealing the field evolution of YbMgGaO 4 , including a magnetic phase crossover, and by establishing relevant criteria which we employ to systematically constrain the possible Hamiltonian parameters. We conducted measurements of high-quality, single crystalline YbMgGaO 4 using the ultra-sensitive tunnel diode oscillator technique [54] (TDO), magnetization via SQUID, cantilever torque magnetometry, and diffuse neutron scattering. As a function of temperature and applied field, varying sample orientations using TDO, SQUID, and cantilever torque magnetometry measurements are used to show anisotropy of the bulk response, while diffuse neutron scattering provides the finite Q dependence of spin correlations. All provide clear evidence of a field-induced phase crossover, which must be reproduced by any set of Hamiltonian parameters describing the magnetic ground states and phase diagram of YbMgGaO 4 . We compare our experimental data to complementary Monte Carlo (MC) and Density Matrix Renormalization Group (DMRG) calculations, where the results demonstrate that despite the diverse range of exchange parameters proposed for this system [17,25,38,39,43], only strictly constrained values can reproduce the field-induced phase crossover observed. We further employ a methodical and systematic approach to reduce the possible volume in the 7-dimensional parameter space. These results further elucidate the effects of disordered exchange interactions and g-factors and suggest a magnetic fieldinduced phase transition in the disorder-free limit. Furthermore, the strictly constrained parameter set suggests that YbMgGaO 4 's ground state is proximal to the QSL state in the disorder-free limit predicted by [23], as will be discussed below.

A. Sample Synthesis and Characterization
Powder of YbMgGaO 4 (structure shown in Supplementary Figure 1) was produced by finely grinding mixed powders of Yb 2 O 3 , MgO, and Ga 2 O 3 , and reacting at 1150 • C in a box furnace. The product was ground again to a fine powder, compressed hydrostatically into a rod, and then sintered at 1500 • C in a vertical Bridgman furnace. Finally, large single crystals of YbMgGaO 4 were grown using the optical floating zone method (Supplementary Figure 2). A typical crystal was grown in an O 2 atmosphere at 1 MPa, with an initial growth speed of 20 mm/hr, and upon stabilization of the liquid zone, 4 mm/hr until finished.
We confirmed the powder and crystal phase at each step of the synthesis using ground powders and the Panalytical X'Pert PRO MRD HR XRD System (using Cu K-α 1.5418 nm X-rays) and Rietfeld refinement via FULLPROF [55]. Single-crystal quality and alignment was determined using Laue X-ray diffraction (Supplementary Figure 3). Preliminary susceptibility and measurements were carried out on a powder sample using a Quantum Design MPMS XL7 SQUID magnetometer down to 1.8 K, and confirmed previous measurements of the Curie-Weiss temperature Θ ≈ −4 K.

B. High Sensitivity Magnetization Measurements
High-sensitivity measurements of magnetization were achieved with the complimentary tunnel diode oscillator (TDO) technique ( Fig. 1) and torque magnetometry. In a TDO measurement, a tunnel diode is biased to operate in the "negative resistance" region of the IV-curve. This provides power that maintains the resonance of a LC-circuit at a frequency range between 10 and 50 MHz. An approximately cylindrical single-crystal sample with dimensions of 2 mm in length and 1 mm in diameter was placed inside a detection coil, with the c axis of the sample aligned with the coil axis (Supplementary Figure  5). Together, they form the inductive component of the LC circuit. Changes in sample magnetization induce a change in the inductance, hence a shift in the resonance frequency. Highly sensitive measurements in changes of magnetic moments 10 −12 e.m.u., therefore, are enabled by the ability of measuring the resonance frequency to a high precision [56]. The magnetization and susceptibility results collected at temperatures down to T = 1.8 K (Supplementary Figure 4) were consistent with earlier reported results [39,43,57].
Magnetization was also directly measured via an inhouse Cryogenic S700X SQUID magnetometer in temperatures down to 280 to 300 mK using a Helium 3 probe. A 1.3 mg sample was mounted on a silver straw with vacuum grease in H c and H ⊥ c orientations. See Fig.  1d and 1e.

C. Neutron Scattering
Neutron scattering data was collected at the CORELLI spectrometer at Spallation Neutron Source, Oak Ridge National Laboratory [58]. CORELLI is a quasi-Laue TOF instrument equipped with a large 2D detector, with a -20 • to +150 • in-plane coverage and ±28.5 • out-of-plane coverage. The incident neutron energy was between 10 meV and 200 meV. A superconducting magnet was used to provide a vertical magnetic field up to 5 T, which reduced the effective out-of-plane coverage to ± 8 • . An 896.10 mg single crystal was mounted on a Cu plate in a dilution refrigerator. The sample was aligned with the (h, k, 0) plane horizontal and the magnetic field along the [0,0,l] direction. Neutron-absorbing Cd was used to shield the sample holder to reduce the background scattering. Experiments were conducted with applied fields at the base temperature 130 mK by rotating the crystal through 180 • in 3 • steps, and then at 20 K in the same fields for background subtraction. The data were reduced using Mantid for the Lorentz and spectrum corrections [59].
To ensure our observations utilizing total scattering mode were consistent with measurement at the best possible energy resolution of the instrument, we also measured 130 mK and 20 K at 0, 3, and 5 T with the correlation chopper on (see supplementary figure 8). We further filtered the TOF data to provide an estimated energy resolution of 0.2 meV and identified features completely consistent with our measurement in total scattering mode.
To account for the temperature factor in our background subtraction in total-scattering mode, we compared the ratio of integrated intensities of a small region in reciprocal space at 5 T for both temperatures. The region was bounded by −0.9 < h < −0.75, 0.75 < k < 1.1, and −0.5 < l < 0.5 (the l range is consistent with our BZ edge and planar integrations). We used the second BZ and 5 T data to avoid skewing the integration at low temperature due to diffuse magnetic scattering. This integration approximates the ratio of Bose population factors -our scaling factor was 1.133 ± 0.005. To improve statistics, we used the 6-fold symmetry (in accordance to the symmetry of the l = 0 plane) to add our data 6 times at subsequent 60 • angles and average. All analysis and visualization were performed using Mantid and Python.

III. RESULTS
We performed ultra-highly sensitive magnetic susceptibility measurement on YbMgGaO4 using TDO technique (Supplementary Figure 5) up to fields exceeding saturation and at temperatures down to 24 mK ( Fig. 1). Comparison to magnetization measured at 300 mK in SQUID (Fig. 1e, f) and cantilever torque magnetometry (Supplementary Figures 5 and 6) results confirms the appearance of an anomaly near 2 T. The anomaly weakens with increasing temperature, beginning to soften after 190 mK and vanishing near 4 K (Fig. 1b,c). The greater prominence observed in both techniques with applied field perpendicular (versus parallel) to the c-axis (Fig. 1b) is consistent with easy-plane anisotropy in agreement with earlier experiments [39]. The anomaly is likely a remnant of the quantum mechanical 1/3 magnetization plateau [60], as will be discussed below.
To gain insight into the changing magnetic structure corresponding to the anomaly detected in the magnetometry measurements, and to exploit the advantages offered by neutron scattering in revealing spin correlations at finite Q, we employed diffuse neutron scattering in total scattering mode at CORELLI, an instrument optimized for measuring short-range correlations. Data was collected with the sample c axis parallel to the applied field direction (Fig. 2a). 20 K background signal was subtracted from 130 mK data for integer fields from 0 to 5 T, isolating the magnetic contribution. To improve statistics, we averaged intensity about the l-axis by applying allowed symmetry operations (see Supplementary Materials for details). Integration over −0.5 < l < 0.5 in reciprocal space for zero field shows diffuse magnetic scattering centered at high-symmetry M points, indicative of short-range correlations and consistent with previous neutron measurements [38,43,57]. However, we discovered that with increasing applied field the spectral weight at the M points subsides. At µ 0 H = 1 T the spectral weight is flat along the zone edge, while at µ 0 H = 2 T most of the diffuse spectral weight is distributed at the K points, coinciding with the feature in our magnetic susceptibility data (Fig. 1). Importantly, we note that the intensity at the K point below 2 T is mostly independent of the applied field. When field is further  increased above 2 T, the scattering intensity at K points diminishes as the system approaches the polarized state. The field dependence for 2 to 5 T is consistent with the gradual saturation of the curve observed in our magnetic susceptibility data. Integration of the spectral weight in a narrow rectangular volume along the first Brillouin Zone (BZ) edge shows the shift in spectral weight from the M to K points (Fig. 2c). This shift in intensity indicates a change in short-range correlations and underlying spin structure, and likely relates to a crossover from a stripetype state to 120 • -type state in the clean limit. As field is further increased, scattering intensity along the zone edge is reduced as seen in scattering mode using a correlation chopper and time-offlight (TOF) filtering to achieve an energy resolution of E ≈ 0.2 meV (see Supplementary Figure 8).

I n t e g r a t e d I n t e n s i t y (
To shed light on the nature of the observed phenomena, we start by considering a classical limit of the spin Hamiltonian proposed in previous works [61]: with phase factors γ ij = 1, e i2π/3 , e −i2π/3 for each of the three principal directions of the triangular lattice. The bond notation ij ( ij ) indicates that the corresponding sum runs over pairs of nearest (second-nearest) neighbors. The corresponding exchange nearest-neighbor and second nearest-neighbor exchange tensors are,
In a first attempt, we manually optimized the Hamiltonian parameters starting with a potential set given in the Ref. 25 and exploring the parameter space to capture details of the field dependent behavior and single magnon dispersion. We considered the set of Hamiltonian parameters, J ± = 0.066J (1) zz and J (1) zz = 0.164 meV, that reproduce the phenomenology better (featured as the magenta star in Figure 4). For the average g-factor, we used g = 3.7227 [38].
To quantify the uncertainty of the proposed model Hamiltonian for a clean version of YbMgGaO4, we applied the iterative optimization procedure explained in Ref. 62. The experimental input includes the magnon dispersion at 7.8 T measured with Inelastic Neutron Scattering (INS) [43] and the field-induced crossover revealed by our magnetic diffuse scattering measurements. The challenge resides in the combination of a highdimensional (d=7) Hamiltonian space (H includes seven independent parameters J (1) ± , J (2) zz and g ), and the significant amount of disorder present in YbMgGaO 4 . Given that both sources of complexity are present in many other frustrated materials with large spin-orbit coupling, it is important to develop new protocols to extract models from data and simultaneously estimate their uncertainty.
We then implemented an optimization protocol in three steps. First, we introduced the the cost-function: to fit the measured magnon dispersion at µ 0 H = 7.8 T with the analytical expression obtained from linear spin wave theory (see Ref. [43]) along two reciprocal space pathways (see Supplementary Figure 9). As explained in Ref. 62, for each iteration we used random samples over the whole Hamiltonian space to build a low-cost estimator of χ 2 IN S . We then usedχ 2 IN S to evaluate the next set of Hamiltonian parameters uniformly distributed over the Hamiltonian space and subjected to the constraintχ 2 IN S < c. The cutoff c is lowered after each iteration. The last iteration is reached when c reaches its final value, c f inal , for which the calculated dispersion agrees with the INS data within the experimental uncertainty. We note that experimental uncertainty is higher than the instrument resolution because the chemical disorder broadens the single-magnon lineshape. The manifold of model Hamiltonians that fit the measured spin-wave dispersion is indicated by the blue contours in Fig. 3.
Second, we performed classical Monte-Carlo simulations (see Methods: Simulation) on the disorder-free model within the aforementioned manifold. For a given Hamiltonian space sample, we computed the dynamical spin structure factor: (3) where Q is the wavevector in the scattering process, α, β = x, y, z are cartesian coordinates indicating initial and final spin polarization of the neutron, F (Q) is the magnetic form factor and S αβ (Q, ω) is the Fourier transform of the two-point spin correlation function. Note that for a classical spin model, S(Q, ω = 0) S(Q) ≡ ∞ 0 S(Q, ω)dω at low enough temperature. The MC simulations are performed at T = 130 mK for three values of the magnetic field (0 T, 1 T and 2 T). We then applied the constraint that the M-peak must be dominant below µ 0 H = 1 T, while the K-peak should become dominant at µ 0 H = 2 T (H c axis). The resulting submanifold of possible model Hamiltonians is delimited by the green color contour in Fig. 3. Finally, to further reduce the model uncertainty, we introduced another cost-function, where I K is the integrated intensity under the K-peak, and applied the same optimization scheme to the MC calculations. We chose the K-peak intensity because our MC simulations of the disordered model indicate that it is less-susceptible to chemical disorder in comparison to the M -peak intensity (see Supplementary Figure 14).
The region constrained by the condition χ 2 combined = χ 2 IN S + χ 2 k < c k is shown by the magenta color contours in Fig. 3. Due to our limited knowledge of the chemical disorder in this material, it is not possible to set a definite value for c k f inal . Consequently, we are presenting three solid, dashed and dotted contour lines, corresponding to higher to lower values of c k f inal , that reveal the landscape of χ 2 combined . We note that the set of parameters that we found by manually optimizing the Hamiltonian parameters, indicated with a star in Fig. 3, lies within the range of possible solutions.
As we discuss below, the disorder present in YbMgGaO 4 introduces a significant variation in the value of these Hamiltonian parameters. To demonstrate this statement, we introduce chemical disorder in the exchange interactions (J-disorder) and in the g-factors (gdisorder) and assume that both have a uniform distribution of predefined width around the aforementioned mean values, g i = g av + ∆ i (σ g ), and J ij = J av (1 + ∆ ij (σ J )), where σ g = 1.2 and σ J = 0.5. The static magnetic structure factor S(Q) is calculated by a standard Metropolis sampling algorithm while the dynamical structure factor S(Q, ω) is computed by Landau-Lifshitz dynamics [63]. Both types of structure factors shown throughout this paper are calculated by averaging over 60 independent sets of 2000 configuration samples from a standard Metropolis sampling algorithm on a 48 × 48 supercell (2304 spins), followed by a slow annealing process starting from a random disordered spin configuration.
In Fig. 2 we show the comparison between experimental data and the theoretical simulations throughout the field-induced phase crossover, where main aspects of the experimentally observed phenomenology can be reproduced. This agreement includes the evolution of the diffuse scattering and the non-linearity observed in the magnetization curve at finite fields. Interestingly, only a strictly selective set of parameters can reproduce the experimental results for all fields and capture the crossover (see Fig. 5 and discussion below). This is so despite the broadening of the observed magnon dispersion in the polarized state that is induced by the chemical disorder (see Supplementary Figure 17).
To account for the effects of quantum fluctuations, we have computed the M(H) curve with the DMRG method on the S = 1/2 version of H (without disorder) for field directions parallel to the c-axis and to the ab-plane. The resulting M(H) curve exhibits a characteristic plateau at 1/3 of the saturation value [60,64], which is bigger for H ⊥ c because of the easy-plane anisotropy (Fig. 3). This plateau phase is a true quantum mechanical signature where quantum fluctuations favor collinear configurations. The effect of quantum fluctuations can then be reproduced by adding an effective bi-quadratic interac- where H is the classical spin Hamiltonian given in Eq.
(1) and λ 1 = 0.3J (1) ± and λ 2 = 0.03J (1) ± are effective biquadratic couplings for next and second nearest-neighbor spins, respectively. The values of λ 1 and λ 2 have been obtained by comparing the M (H) curve obtained with DMRG for the S = 1/2 version of H and the M (H) curve obtained from a classical MC simulation of H 1 . To compare the calculations and TDO magnetic susceptibility results, we compare the change in resonant frequency (directly proportional to dM/dH measured with SQUIDsee Fig. 1d,e) then integrate with respect to applied field to obtain a curve proportional to magnetization (Fig.  3c), revealing a distinct non-linearity. This non-linearity can also be seen in the low-temperature calculation result (Fig. 3a), where disorder has been included. These results suggest that the suppression of the feature in M (H) curves reported in earlier work [39,43,57] could be due to thermal fluctuations (since these were measured from 1.7 K to 2 K). Furthermore, we note that the lower 0.5 K curve measured from a powder sample as reported in Ref. 18 captures the field induced anomaly and agrees with our results at comparable temperatures. Lower temperature and more sensitive techniques are required to detect the anomaly with better resolution.
A comparison of the calculated S(Q) in the absence of disorder at T = 0 (Supplementary Figure 18) to T = 130 mK (Supplementary Figures 19) shows that thermal fluctuations are enough to disrupt long-range order, owing to the high degree of frustration. Furthermore, the calculated S(Q) in the presence of disorder as T approaches zero also shows disruption of long-range order (Supplementary Figure 21). These findings indicate that the lack of magnetic Bragg peaks and consequently the broad continuum of magnetic excitations can likely be attributed to thermal fluctuations or chemical disorder. In other words, even if it were possible to produce cleaner samples of YbMgGaO 4 , our classical analysis indicates that lower temperatures would be necessary to measure the presence of long range order. While this observation does not preclude the possibility of a spin liquid state induced by quantum fluctuations, it indicates that the lack of Bragg peaks down to T = 130 mK cannot be regarded as a strong indicator of spin liquid behavior.
Recent theoretical studies of the effect of disorder on the J 1 − J 2 Heisenberg triangular antiferromagnets show that disorder may contribute more intensity along the BZ edge (in zero field) [22]. Likewise, as noted above, chemical disorder is essential to reproduce the broadening of the magnon peaks (Supplementary Figure 17) at µ 0 H = 7.8 T as also observed in experiment [43]. We therefore suggest that both thermal fluctuations and chemical disorder contribute to the observed lack of a long-range order in YbMgGaO 4 . In particular, we note that as QSL candidates in general are highly frustrated, they may be similarly susceptible to the effects of thermal fluctuations (even at very low temperatures) and thus these effects should be taken into consideration in future studies of this compound and other highly frustrated systems with relatively small Curie-Weiss temperatures.
The above observations leave us with the following question: how can we use the available experimental information to diagnose quantum spin liquid behavior in the clean limit? One possible approach is to extract a model Hamiltonian for the clean limit of YbMgGaO 4 and compare the result against the quantum phase diagram that is obtained from numerical techniques, such as DMRG [23], that can account for non-perturbative effects of quantum fluctuations. As we will briefly elaborate upon below, when we extend our work to optimize a case with isotropic J 2 interactions, we find parameters that lie in the spin liquid region indicated in Fig. 4 of Ref. 23).
FIG. 5. Comparison of field dependence of S(Q) calculated for parameters from four unique models, including those found in this study and as suggested by constraints determined in Ref. 17. S(Q) calculated for integer fields 0-5 T at 130 mK. Note that the field-induced transition visible in the diffuse scattering data is absent for the parameters other than those used for this study, suggesting stringent constraints on the possible models for the system.

IV. DISCUSSION
Definitively proving the existence of the quantum spin liquid ground states in candidate systems remains a critical experimental challenge [2]. The lack of long-range order in QSL states, coupled with the quenched disorder of many real experimental systems like YbMgGaO 4 , suggests that future studies will benefit greatly from an improved understanding of the role disorder can play in QSL phenomena. Indeed, disorder has been identified as a contributing factor for several proposed QSL candidates [66,67]. However, it is important to recognize that a good characterization of the kind of disorder that is present in a specific compound can be extremely challenging. If the minimal spin Hamiltonian of a given material has multiple parameters (seven for the case of YbMgGaO 4 ), extracting an accurate model Hamiltonian is already very challenging in the clean limit.
An alternative approach is to study whether a system with quenched disorder has a QSL state in the disorderfree limit to determine if it is worth investigating ways of reducing the level of chemical disorder in a specific material. The present work proceeds accordingly: the presence and unique characteristics of a phase crossover have been used to narrowly constrain the possible range of magnetic Hamiltonian parameters, so that the Hamiltonian in the clean-limit can be established and the ef-fects of disorder interpreted with that in mind. We emphasize that despite the high quality of past studies, more information about the low and intermediate field behavior are critical to map out the low temperature magnetic phase diagram. Our observation of the crossover in YbMgGaO 4 yields a lower-field reference point for comparison to measurements of the polarized state. Furthermore, our results show a remarkably high sensitivity of the field-induced crossover to Hamiltonian parameters. As an example, we consider recent timedomain THz spectroscopy (TDTS) measurements employed in conjunction with inelastic neutron scattering to measure the disorder-averaged exchange interactions in YbMgGaO 4 [17]. This work found a slightly larger gfactor (g = 3.81) than previously reported. We applied the exchange parameters and g-factor proposed in model C of Ref. 17 as well as a few variants of this model. Remarkably, we find that these parameters do not reproduce either the diffuse scattering, or the field-induced crossover we observe (Fig. 5) -this is despite the close agreement between the S(Q, ω) spin waves at higher field calculated with both sets of parameters (Supplementary Figure 22). The sensitivity of the phase crossover to the parameters is exemplified by Fig. 3 b, where only a small discrepancy in the value of J ±± contributes to the absence of the crossover in applied field (see Fig. 5). However, we note that optimizing the g-factor improves the quality of ours fits and yields a value of g = 3.84, in good agreement with Ref. 17 (see supplementary figure  13).
Conservatively, our results indicate an intriguing competition of factors. The spin liquid signatures observed in YbMgGaO 4 can arise from high frustration working in concert with thermal fluctuations even at very low temperatures, and these effects are likely enhanced by the inherent chemical disorder of the system. The diffuse scattering is greatly enhanced by finite temperature, owing to the very high frustration (see supplementary figure 21). Intriguingly, the scattering intensity at the M points is particularly sensitive to the quenched disorder when compared to the K points (see supplementary fig.  14). This can likely be understood by the greater degeneracy of possible orientations of stripe order-type phases on the triangular lattice, which generally can have up to 3 orientations. At the same time, the presence of the field-induced crossover also reaffirms quantum mechanical fluctuations as opposed to spin freezing at temperatures as low as 24 mK. Along with the revised constraints on exchange parameters, these results invite further inquiries into the ground state of YbMgGaO 4 .
The central aim of studies of YbMgGaO 4 and QSL candidates more generally is to determine the nature of the ground state and to confirm or reject a QSL interpretation. In this work, we have considered a more general Hamiltonian than the one considered in previous works by allowing the second nearest neighbor interaction, J (2) , to be anisotropic. The optimal solution for an isotropic (Heisenberg) second nearest neighbor interaction J (2) is J  figure 15). It is interesting to note that this optimal solution belongs to a spin-liquid phase of the quantum phase diagram reported in Fig. 4 of Ref. 23. Given the high frustration and consequential sensitivity to thermal fluctuations, as well as the presence of quenched disorder, more studies will be needed to disentangle these various contributions to spin liquid features from a possible QSL state.
The effects of disorder on possible QSL states merits further exploration, and a few recent studies provide intriguing comparisons to YbMgGaO 4 . In stark contrast to YbMgGaO 4 , where the magnetic ground state defies long range ordering even as the field is raised above the phase crossover (below the polarized state), the closely related triangular antiferromagnet QSL candidate NaYbO 2 antiferromagnetically orders at a critical applied field [5,68]. On the other hand, spin-liquid-like states, as have been suggested for YbMgGaO 4 [22], have been proposed as a consequence of cation mixing in the perovskite Sr 2 Cu(Te 0.5 W 0.5 )O 6 [69].
Recently, a few studies have likewise used our data to place constraints on the possible parameters [34,35]. We note the conclusions compare favorably despite alterna-tive theoretical techniques.

V. CONCLUSION
In summary, our magnetometry and neutron diffraction experiments show a field-induced crossover, which we compare to DMRG and Monte Carlo simulations to place the strictest systematically developed constraints to date on the possible Hamiltonians describing the system. The presence of this crossover acts as a litmus test for the possible Hamiltonians for the system, and this approach can be extended to other systems where disorder makes the typical means of measuring exchange parameters more challenging. The range of data used illustrates the importance of establishing sufficient experimental evidence to determine appropriate model parameters in systems with complex Hamiltonians and especially in the presence of quenched disorder. Furthermore, our results invite further investigation into the competition of high frustration, quenched disorder, and quantum and thermal fluctuations whose significance persists to the lowest investigated temperatures. Finally, our constraints suggest a proximity to the possible QSL state in the disorderfree limit for the case of isotropic next-nearest neighbor interactions.

VI. ACKNOWLEDGMENTS
We are grateful to N. Peter Armitage and Alexander L. Chernyshev for their thoughtful comments and discussion. We are grateful to Andrei T. Savici for help in analyzing the neutron scattering data. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by the National Science