Phase Diagram of YbZnGaO4 in Applied Magnetic Field

Recently, Yb-based triangular lattice antiferromagnets have garnered significant interest as possible quantum spin liquid candidates. One example is YbMgGaO4, which showed many promising spin liquid features, but also possesses a high degree of disorder owing to site-mixing between the non-magnetic cations. To further elucidate the role of chemical disorder and to explore the phase diagram of these materials in applied field, we present neutron scattering and sensitive magnetometry measurements of the closely related compound, YbZnGaO4. Our results suggest a difference in magnetic anisotropy between the two compounds, and we use key observations of the magnetic phase crossover to motivate an exploration of the field- and exchange parameter-dependent phase diagram, providing an expanded view of the available magnetic states in applied field. This enriched map of the phase space serves as a basis to restrict the values of parameters describing the magnetic Hamiltonian with broad application to recently discovered related materials.


I. INTRODUCTION
Recent years have seen renewed interest in triangular-lattice antiferromagnets featuring anisotropic interactions and other traits conducive to exotic quantum ground states, particularly in the hunt for experimental realizations of quantum spin liquids. Mapping the phase diagrams of these materials is thus of paramount importance, as the variation in magnetic anisotropy, relative exchange coupling strengths, and corresponding magnetic Hamiltonians offered by different materials are fundamentally important to the pursuit of such states.
One system that sparked an explosion of experimental and theoretical interest is the triangular-lattice antiferromagnet YbMgGaO 4 , which has led to a broad research effort that currently involves a number of related compounds [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. There are two important physical aspects of this group of materials that have become clear due to recent insights. Because of the localized nature of the f -orbitals, the ranges of the effective spin-spin interactions in these insulating materials are strongly limited, implying that all of the Kramers-ion-based materials are expected to be closely described by the same nearest-neighbor anisotropic-exchange model with the parameters that are permitted by triangular-lattice symmetry [17][18][19][20]. This reasonably compact model should provide a consistent interpretation of the current and future experiments and give important new insights into fundamental properties of these materials [21].
The second generic feature is the strong effect of disorder observed in some of the well-studied representatives of these compounds. It has been argued theoretically that even a benign form of bond disorder necessarily generates perturbations that are relevant in the Imry-Ma sense [22][23][24], making a consideration of the defects an inevitable and essential part of a realistic description of most anisotropic-exchange magnets. Empirically, many of the newly synthesized materials seem to show no magnetic ordering [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Then, given the strong disorder effects, it is a question whether the non-magnetic ground states of these materials are due to a genuine quantumdisordered spin-liquid (SL) phase, or due to a scenario similar to the disorder-induced "spin-liquid mimicry," suggested for YbMgGaO 4 [23].
There is also an intriguing broader question on whether the disorder-induced spin-liquid-like behavior retains any of the unique and desired properties of the intrinsic spin liquids, thus potentially turning disorder into a feature rather than an obstacle [25][26][27]. A counterintuitive example of the role of disorder in a related material is the case of NaYbO 2 , where introducing Na + site vacancies leads to an antiferromagnetic transition at a few Kelvin and thus disorder supports an ordered phase [28]. Further understanding and insight into the role of disorder are needed to make progress in this direction.
One of the persistent issues in the studies of the anisotropicexchange magnets in general and of the rare-earth family in particular is the identification of their model parameters [29]. In the case of the disorder-induced pseudo-SL state, this problem is further aggravated, as it is not clear what state the disorderfree system would have assumed. In this work, we propose that the experimental and theoretical investigations of the fieldinduced phases offers a powerful instrument to significantly narrow the allowed parameter space to a region that is consistent with the material's phenomenology.
In YbMgGaO 4 and YbZnGaO 4 , the source of disorder is due to the R3m symmetry, leading to fifty-fifty site mixing of the non-magnetic cations Mg 2+ or Zn 2+ with the Ga 3+ . Efforts to determine the exchange parameters for placing YbMgGaO 4 in proposed phase diagrams and otherwise comparing to the numerous theoretical investigations to affirm or deny a QSL state were obstructed by the various broadening effects and consequentially enhanced uncertainty in the measurements. Several studies concentrated their efforts on further refining measurements of the exchange parameters [30,31].
In this work, we provide a detailed study of the field-induced effects and characteristics of YbZnGaO 4 and offer comparison with the results in YbMgGaO 4 [31]. Informed by measurements from high resolution magnetometry and a variety of neutron scattering techniques, we put forward a theoretical analysis of the structure of their field-induced phase diagram and propose a parameter region for both materials that is compatible with our empirical findings.

A. Preliminary characterization
Susceptibility measurements on a 1.81 mg single crystal sample of YbZnGaO 4 conducted using an in-house Cryogenic S700X SQUID magnetometer (with 3He probe) yield low-temperature fits suggesting Curie-Weiss temperatures of Θ CW = −2.67 K and Θ CW = −2.62 K for parallel and perpendicular to the sample c axis, respectively (see Supplementary Figure 4). These values are more isotropic than those reported for YbMgGaO 4 or for YbZnGaO 4 in earlier works though YbZnGaO 4 was also more isotropic than YbMgGaO 4 in those measurements) [32].
We emphasize that while the difference in the Curie-Weiss temperature for the in-plane vs out-of-plane field has initially suggested a rather strong easy-plane character of YbMgGaO 4 [33], the subsequent spectroscopic studies have hinted at a rather moderate XXZ anisotropy, yielding a nearly Heisenberg value of ∆ = 0.8-0.9 [30]. In the case of YbZnGaO 4 , Curie-Weiss temperatures for the in-plane and out-of-plane fields from our measurements and those of earlier works [32] are much closer to the Heisenberg limit. Given the trend, this indicates that the anisotropy in YbZnGaO 4 may, in fact, be of the easy-axis type. A more direct demonstration of that is offered by the results that are discussed in Sec. II.B.
A unique opportunity afforded by a close comparison between YbMgGaO 4 and YbZnGaO 4 is the qualitative contrast of the effect of the cation substitution on the disorder between the two materials. One obvious consideration is the differing ionic radii of the Mg 2+ and Zn 2+ , which are approximately 72 and 74 pm respectively (a 2.7% difference). This small difference yields slightly smaller lattice parameters for YbMgGaO 4 when compared to YbZnGaO 4 , as shown by comparison of parameters given by references [34] and [32], respectively, which, in turn, may yield marginally stronger exchange for YbMgGaO 4 , also explaining the smaller Curie-Weiss temperatures and lower field onset for the anomaly discussed in Sec. II.B for YbZnGaO 4 . This is consistent with the observation that the related compound NaYbO 2 , with an in-plane lattice parameter smaller by only about 1.8% from YbMgGaO 4 , shows a significantly larger Curie-Weiss temperature (Θ CW = −10.38 K [12] as opposed to −4 K for YbMgGaO 4 ) [31]. Furthermore, both Zn 2+ and Ga 3+ have d 10 electronic configurations, whereas Mg 2+ is p 6 . While the displacement of Yb 3+ can still be expected based on the charge difference between the cations, leading to the observed broadening in inelastic neutron scattering studies of the single magnon dispersion and crystal electric field levels [35], the local environment may be more homogeneous for YbZnGaO 4 , and may be related to the difference in anisotropic response under field. However, further studies will be required to compare the effective role of disorder between the two systems.

B. Tunnel diode oscillator technique and SQUID magnetization
The first indications of the magnetic transitions in YbZnGaO 4 , which are similar to the ones we have previously identified in YbMgGaO 4 [31], were found using the tunnel diode oscillator (TDO) technique (see methods). As the applied magnetic field is changed, the sample magnetization M(H) with respect to field is altered, thus changing the inductance of the coil and the measured resonant frequency of the circuit, yielding a signal proportional to χ(H). From Figure 1 (a), where the change in resonant frequency is plotted versus the applied field, a clear nonlinearity is apparent beginning just below 1 T. This nonlinearity persists to at least 2 T for the field parallel to sample c axis. Upon raising the temperature, the feature is completely suppressed at about 4 K, affirming its magnetic origin. This behavior is also consistent with a similar feature in YbMgGaO 4 [31]. As the sample orientation is rotated with respect to the field, the anisotropic response is made apparent, with the feature encompassing a broader range of the field for H ⊥ c. This feature is confirmed by more conventional magnetization and susceptibility measurements carried out using the in-house SQUID magnetometer with a 3He probe in temperatures down to 300 mK. The same distinct plateau-like feature is apparent in the χ(H) in both TDO and SQUID measurements, as is clear from Figures 1 (d) and (e). Integrating the change in frequency with respect to applied field yields a curve consistent with magnetization as measured in SQUID (see Figure 1 (c)). We note that the anomaly in YbZnGaO 4 measured via TDO and SQUID magnetization occurs at a slightly lower field compared to the analogous feature measured in YbMgGaO 4 [31] (see Supplementary Figure 5).
The features in the magnetization derivative and tunnel diode oscillator technique (TDO) measurements in both in-plane and out-of-plane fields are reminiscent of the plateau-like behavior that is expected in the canonical Heisenberg or XXZ nearest-neighbor triangular-lattice magnets [36]. Importantly, the TDO and magnetization measurements in YbZnGaO 4 and YbMgGaO 4 , together with the Curie-Weiss extrapolations from the susceptibility mentioned above, suggest an important distinction between the two materials. The effects of the in-plane and the out-of-plane field directions show different behavior  [32] as a function of the field orientation -compare Figure 1(b) of this work and Figure 1(c) in Ref. [31]. Specifically, the relative severity of the anomaly in TDO data in YbZnGaO 4 in H c are resemblant of that in H ⊥ c data of YbMgGaO 4 , and vice versa. This suggests that YbZnGaO 4 and YbMgGaO 4 correspond to different types of the XXZ anisotropy, the easy-axis and the easy-plane, respectively.

C. ac susceptibility and disorder
From our measurements of ac susceptibility of YbZnGaO 4 (see Figure 1 (f)), we find critical temperatures corresponding to the characteristic cusps of the ac susceptibility occur at substantially lower temperatures than previously measured [32]. Indeed, our characteristic temperatures are around 20% − 30% lower for all comparable frequencies -see Figure 1 (g). Consequently, for our measurement of ∆P = ∆T f T f ∆ log(ω) , a quantitative measure of the freezing temperatures per decade of characteristic temperatures, we find a value of ∆P = 0.139(4), substantially larger than the ∆P = 0.053(4) of previous work. This value of ∆P is typically associated with superparamagnetic behavior [37] as opposed to spin glasses. That being said, insulating spin glasses typically show greater frequency dependence [37]. If this ∆P is interpreted as indicative of superparamagentic behavior, it may be understood as a consequence of many microscopic domains with insignificant cooperative freezing. This phenomenon is likely due to disorder, but further suppression of the freezing temperature could be attributed to the high degree of frustration as well as disorder. The percentage of frozen spins is estimated to be approximately 16% (see Supplementary Figure 9), comparable to the previous study [32] and to the case of YbMgGaO 4 [33].
General questions about the effects of disorder in spin systems persist, especially in light of its aforementioned potential relationship to QSL states for some materials [25][26][27]. The origin and effects of disorder related to the observed spin-liquid features of YbMgGaO 4 have been addressed earlier [23,35].

D. Neutron scattering
For this work, diffuse neutron scattering data were collected at CORELLI at Oak Ridge National Laboratory in total scattering mode (see Figure 2). The sample was aligned with the [h, k, 0] scattering plane, and applied field along the sample c axis. Data were collected at 130 mK for 0, 1, 1.5, 2, 3, 4, and 5 T. Color maps of the neutron scattering intensity after subtracting the 20 K as background reveal the evolution of the magnetic structure, which is qualitatively comparable to the YbMgGaO 4 data. With no applied field, the intensity largely resides at the high-symmetry M points on the edges of the Brillouin zone (see also Supplementary Figure 6, left panel). At 1 T the scattering intensity is almost completely uniform along the zone edge, while at 1.5 T the intensity is found predominantly at the high-symmetry K points. Intensity at the M points is well established to correspond to stripe-ordered states in long-range ordered triangular systems, while intensity at the K point is generally suggestive of 120 • -type ordering or other three-sublattice states [20]. As in the case of YbMgGaO 4 , this migration of the scattering intensity with applied field corresponds to the anomaly observed in the magnetometry data. The changes in intensity were further confirmed by measurements with the triple-axis spectrometer at SPINS (see Supplementary Figure 7). We further measured inelastic neutron scattering (INS) from a YbZnGaO 4 single crystal sample at CNCS [38] in applied field(see Figure 3) at Oak Ridge National Laboratory. Here we again see a clear evolution of the intensity as a function of energy and Q with increasing applied field. We note that at low field the intensity is notably concentrated on the M points. The intensity has no clear dispersion up to about 3 T, but instead consists of a broad continuum, similar to the 0 T with weakening intensity as the spins are presumably canted further from the scattering plane with increased field. At 4 T, the scattering remains broad in energy, but a dispersion is faintly visible at low energy along the zone edge. This dispersion has minima at the zone edges and rises as it approaches the Γ points. The shape of this dispersion closely resembles what is measured in the polarized state measured at 8 T, indicating that the system is approaching polarization. Measurements at DCS [39] (see Supplementary Figure 6) and SPINS (see Supplementary Figure 7) at the National Institute of Standards and Technology confirm the features shown in Figure 3 in a diversity of backgrounds and instrument setups. At 8 T the system is nearly completely polarized, and a clear dispersion is evident. As in the case of YbMgGaO 4 [33], the dispersion is broadened, likely due to the disorder and the resulting distribution of exchange parameters and g-factors.
Additional measurements to characterize YbZnGaO 4 's response to applied field were conducted using polarized neu-trons at BT7 [40] at the National Institute of Standards and Technology, with a vertical guide field set up (see Supplementary Figure 8). After subtracting background measurements (40 K) from base (0.3 K) and correcting for the polarization rate, comparison of the spin flip (SF, measuring the in-plane component) shows a stronger in-plane component along the zone edge at 0 T, compared to 2 T, with particularly high intensity in the vicinity of the high-symmetry M point, likely affirming increased canting of the spins with increasing field, but also possibly indicating a reduced spin component parallel to the zone edge (pointing to nearest neighbors in real space). This greater intensity near the zone edge for the SF scattering at 0 T can also be seen in the orthogonal cut from M to Γ.

E. Model and zero-field phases
The interplay of the crystal field and spin-orbit coupling on the magnetic moment of the Kramers ion results in the splitting of its levels into a well-separated doublet structure built from a mix of various spin and orbital states. The exchange interactions of the lowest doublets (pseudo-spins- 1 2 ) are constrained only by the discrete lattice symmetry. For the triangular lattice of YbMgGaO 4 and YbZnGaO 4 , this general anisotropic-exchange nearest-neighbor model is given by where bond anglesφ α are that of the primitive vectors of the lattice with the x axis,φ α = {0, 2π/3, −2π/3}. The J ±± and J z± bond-dependent terms are due to the strong spin orbit coupling [17].
There are three main ordered states in the antiferromagnetic limit of the XXZ interaction: a coplanar three-sublattice 120 • state, which corresponds to the ordering vector at K points, a collinear stripe-x two-sublattice state with spins pointing along the nearest-neighbor bonds and the ordering vector at M points, and a second collinear stripe-yz state with the same ordering vector but spins tilted out of the lattice plane and perpendicular to the nearest-neighbor bonds.
We should note that the phase diagram in Figure 4 is simplified. The simplification comes from only taking the single-Q spiral ansatz that does not include more complicated multi-Q states [42,44]. Moreover, the quantum version of the phase diagram also has a spin-liquid phase [19,20], which is located along the tricritical boundary between stripe and 120 • states for a limited range of the XXZ anisotropy near the Heisenberg limit.

F. Exploring the XXZ parameter space
One of the puzzling features observed in some of the first experiments in YbMgGaO 4 [34] was an indication of the fieldinduced phase crossovers, seen in magnetic susceptibility. Recent susceptibility and TDO measurements of YbMgGaO 4 have further supported these observations [31]. Here we present a variety of measurements on high quality single crystal samples to show that very similar features, indicating fieldinduced crossovers for both H c and H ⊥ c also occur in YbZnGaO 4 .
The neutron scattering measurements in the out-of-plane magnetic field H c [31], also indicated a field-induced crossover and brought another piece of evidence to light. Neutron diffraction has showed that the field-induced crossover is accompanied by a shift of magnetic intensity from the Mpoints at the lower fields to K-points at the higher fields. In an ordered state, such an intensity shift would correspond to a transition from a four-sublattice to a three-sublattice state. Our key finding is that this feature alone allows one to put strong boundaries on the exchange parameters of the system when the phase diagram of relevant parameters is considered.
From the susceptibility data presented above and as established by earlier work in the case of YbMgGaO 4 , YbZnGaO 4 can be characterized as easy-axis anisotropy while YbMgGaO 4 is easy-plane. Therefore, in the following we use two representative values of XXZ anisotropy for the model (1), the easy-plane case ∆ = 0.8 that is related to YbMgGaO 4 , and the easy-axis case ∆ = 1.1 as related to YbZnGaO 4 .
First, we explore the parameter region that would allow for a transition from the four-sublattice to the three-sublattice state at some finite value of the out-of-plane magnetic field H c using classical energy minimization. We fix the second-nearestneighbor coupling J 2 to the value 0.05J (red line in Figure 4) and the XXZ anisotropy in the model (1) to the two values discussed above. That leaves the bond-dependent anisotropies J ±± and J z± and the field as the parameters to scan through. We highlight our findings in Figure 5  . The 120 • phase is a three-sublattice state already at zero field and remains such for all the fields, while most of the stripe phase regions remain four-sublattice throughout the entire field range.
Our key finding is illustrated by the gradient-color regions interpolating between the "only-three-" and "only-foursublattice" regions in Figure 5 (a) and (b). They demonstrate that already at the level of our classical energy analysis, the four-to-three sublattice transition indeed takes place at some value of H c < H s in a surprisingly extended region that emanates from the 120 • part of the phase diagram into the stripeyz phase and extends up to J z± ∼ J in both cases, with the intensity emphasizing how far this transition is from the saturation field H s .
It should be noted that for S = 1/2, quantum effects in zero field are known to broaden the stability region of the 120 • phase beyond the boundaries of the classical model consideration [19]. Therefore, one may also expect the region of the fourto-three sublattice transition to extend beyond the classical predictions in this work.
While one can expect that quantum effects will further stabilize and extend the field-induced three-sublattice states, we note that experimentally the "4-to-3" (or M-to-K) transition in both YbMgGaO 4 and YbZnGaO 4 occurs at a rather low field, H c < 0.5H s , which provides further restrictions on the possible parameter ranges. At the level of our approximations, this constraint puts YbMgGaO 4 and YbZnGaO 4 in a close proximity to the 120 • phase boundary, giving the upper bound J z± /J 0.4.
As one can see in Figure 5 (a), there is a narrow region inside the stripe-x phase where the transition of interest also occurs, but it involves another transition back to the four-sublattice state at higher fields, so we do not consider it as a relevant region to the materials in question here.

G. Field-induced states
To provide further insights into the effects of the field, we explore the new phases it can produce. In some cases, the field-induced transformation of the spin structures involves simple canting toward the field until a full saturation is reached at some H s . However, in frustrated spin systems, or systems with anisotropic interactions, the field-evolution is more complicated. The case of the triangular-lattice Heisenberg antiferromagnet is paradigmatic in this respect [36], showcasing a well-known and much-studied sequence of transitions from "Y" to "up-up-down" (UUD) plateau and "V" states in its field-evolution toward saturation. The XXZ extension of the same model also includes non-coplanar "umbrella" and coplanar "fan" states [46]. Next-nearest neighbor interactions and anisotropies also introduce a wider variety of four-sublattice field-induced states [47,48].
However, the field-evolution of the phases of the anisotropicexchange model (1), which combines frustration from the bonddependent terms with that of the triangular-lattice geometry, remains largely unexplored. In this work, we offer some essential understanding and make significant steps of such an exploration.
We use the same representative values of J 2 = 0.05J and of the easy-axis and easy-plane XXZ anisotropies ∆ = 1.1 and ∆ = 0.8 as in Figs. 5 (a) and (b) to provide an insight into the rich phase diagram of the model (1) in a field. In Figure 5  In case of the easy-axis XXZ anisotropy in Figure 5 (c), at the lower values of J z± one can observe the expected canonical sequence of the Y-UUD-V phase transitions of the threesublattice states in the triangular-lattice XXZ model [36,46]. Their corresponding spin structures are sketched in the figure.
The field-induced behavior of the stripe states also includes multiple transitions. At the lowest field, the collinear stripe-yz spin configuration, in which spins are tilted off the basal plane of the lattice, is deformed into a non-coplanar four-sublattice state with all four spins on the elementary plaquette having different tilt angles. There is a broad crossover from this state to a structure with three spins forming an "umbrella" and the fourth strictly antiparallel to the magnetic field, see the sketches of the spin order in Figure 5 (c). This latter state is stable in a wide field region. As the field increases further, at not too small J z± there is a spin-flop-like transition to a similar state, umbrella with a fourth spin parallel to the field. For the yet larger values of J z± , transition to the saturation occurs directly from this "umbrella+up" state via a first-order transition. The main feature in both Figs. 5 (c) and (d), which is also our key finding, is that the region of stability of the three-sublattice states, related to the experimentally observed intensity at the K point, expands at the larger values of the magnetic field. Therefore, there is a region of the model parameters where an evolution from the four-sublattice to the saturated state necessarily proceeds via a high-field three-sublattice state.
For the easy-axis case of Figure 5 (c), this high-field state is a coplanar "V" state. For the easy-plane XXZ anisotropy case of Figure 5 (d), the four-sublattice phases and all the discussed trends are the same, while the three-sublattice region, classically, is a single noncoplanar "umbrella" state, which is a canted 120 • structure. We note that in the quantum case and for not too small XXZ anisotropy ∆, this umbrella state is replaced by the same sequence of Y-UUD-V phases as in Figure 5 (c) [46]. In contrast with the large-J z± first-order transition from the four-sublattice state to saturation, the transition from both V and umbrella states to saturation is second-order. We note that the discussed transitions agree with the prior classical Monte-Carlo simulations [31] conducted in the context of the parameter search for YbMgGaO 4 for a narrower range of parameters.

H. Further evolution of the phases
Here we present further details on the field-evolution of various phases of the model (1) for different choices of J 2 and J ±± . Our Figs. 6 and 7 show the intensity plots of the magnetic susceptibility in the J z± -H plane for the same choices of the XXZ anisotropy ∆ = 1.1 and ∆ = 0.8 as in Figs. 5 (c) and (d), respectively. Each row of the graphs corresponds to a different constant-J ±± cut of Figs. 5 (a) and (b) and each column corresponds to a different constant-J 2 cut of the 3D phase diagram in Figure 4.
While the variations of the J ±± term simply give an elaborate dissection of what happens underneath the projected view of the intensity plots of Figure 5 (a) and (b) for J ±± = ±0.05J cuts, it is clear that the next-nearest-neighbor interaction J 2 works strongly against the field-induced three-sublattice states. This is in accord with Figure 4, which shows that at J 2 ≈ 0.125J the 120 • state is completely eliminated from the zero-field phase diagram. The sets of parameters where no three-sublattice state is observed at any field are highlighted with a red frame in Figure 7. There are other changes that include shrinking or expanding of the phases already described and the appearance of a new canted stripe state for large J 2 and negative J ±± that occur at J z± 0.2J and larger fields for ∆ = 1.1 and for all fields for ∆ = 0.8, see upper right panels of Figs. 6 and 7.
The main takeaway from Figs. 6 and 7 is that the region of the four-to three-sublattice transition is limited by the extent of J ±± already shown in Figs. 5 (a) and (b) and more strongly by the next-nearest-neighbor interaction J 2 . Therefore, the experimental observations that the critical field of such a transition H c H s in both YMGO and YZGO put strong bounds on the anisotropic-exchange parameters. This analysis also limits next-nearest neighbor interactions to J 2 0.1J as only four-sublattice states survive for larger J 2 .

I. S (Q, ω) for select parameters
In order to relate the above theoretical analysis to our experimental results for YbZnGaO 4 and YbMgGaO 4 , we used SpinW [49] to calculate linear spin-wave theory (LSWT) S (Q, ω) for representative sets of parameters from the identified regions of the phase diagram that are shown by dots in Figs. 5 (a) and (b). Since the phase diagrams are for the parameters in units of the overall scale J, the latter is set to match LSWT dispersion in high fields. Comparison of the data to calculations from SpinW is inherently challenging owing to the implicit assumption in the SpinW calculations of the absence of disorder (which is necessary to reproduce the broadening in energy and Q observed in experiment). Nonetheless, for the parameters indicated in Figure 5 (a) we find very good qualitative agreement between calculations and data (see Figure 3) in the absence of disorder. The SpinW optimization algorithm optmagsteep was used with iterative manual adjustment of the spin state until the stable spin state (free of imaginary modes) was found [50].
In particular, comparing the diffuse scattering from Figure  3 (a) 0 T and the resolution-convoluted dispersion calculated from SpinW in Figure 3 (f), we note the enhanced intensity in the vicinity of M and K clearly present in both. Furthermore, the greater intensity residing in a range of energy from near the elastic line to about 0.2 meV in the calculation is clearly present, albeit undoubtedly broadened by disorder, in the ex- perimental data. As discussed in previous works [32,33,35], disorder "smears" the intensity in Q and E, which would also account for the intensity visible at higher energies in the INS data.
In Figure 3 (g) we show the calculated dispersion using SpinW at 2.5 T, where the spins reside in the 'V' state. The dispersion has again been convoluted with the resolution of the instrument for best comparison to the experimental data. Note that the intensity at the Γ point (though suppressed as discussed below) is at approximately 0.5 meV, as expected from comparison to the 2 and 3 T data. Furthermore, the modes above the Brillouin zone edge (which we emphasize follow from the superposition of 3 domains in the clean limit) offer a hint as to the the origins of the the broad continuum apparent in the data. As the field is raised, the contributions from modes above the zone edge in the SpinW calculation are reduced (compare the intensity from modes above M and K at 0 and 2 T). When the system is fully saturated, and the 'V'-state gives way to the field-induced ferromagnetic state, the modes above the zone edge are completely absent (see the curve calculated from the analytic expression following from linear spin-wave theory superimposed on the 8 T INS data in Figure 3 (h)). This evolution of intensity is qualitatively consistent with that observed in the data, where broadening due to disorder smears out the modes. That is, the modes residing above the zone edge contribute to the continuum at low to intermediate field, but are reduced in intensity as the system enters the (subject-to-disorder) 'V'-state, and are replaced by the broadened, otherwise singular mode when the system is fully polarized.
We also point out that an apparent visibility of the magnon modes in LSWT S (Q, ω) in the vicinity of the Γ point at low fields, Figure 3(f),(g), and lack thereof in the experimental data, Figure 3(a-c), can be explained by a significant interaction of the single-magnon branch with the two-magnon continuum. While the quantitative calculation of such effects in  Figure 6 for ∆ = 0.8. Phase diagrams without three-sublattice states are marked with the red frame. the anisotropic-exchange models is quite involved [51,52] and is beyond the scope of the present work, we, nevertheless, provide an intuitive insight into it by showing the bottom of the two-magnon continuum for three magnetic domains in Figure 3(f) and (g). As shown by the dashed lines in Figure 3 for both 0 T and 2.5 T, the bottom of the continuum has lower energy than one magnon modes at the large portion of the Brillouin zone, including the Γ point. For S = 1/2 systems, such an overlap can lead to strong decays and near-complete disappearance of the well-defined magnon excitations in the corresponding range of momenta together with the other phenomena such as strong renormalization of the single-magnon branches [53,54]. These effects require significant coupling between one-and two-magnon sectors -an inevitable consequence of the anisotropic exchange terms in (1), which are allowed by the presence of strong spin-orbit coupling in the rare-earth-based and some transition-metal compounds [29,51,55].
We have also calculated the 0 T S (Q, ω) for the parame-ters identified as appropriate for the easy-plane character of YbMgGaO 4 , identified by a dot in Figure 5 (b), (see Supplementary Figure 10), which is in excellent qualitative agreement with Figure 2 (a) of [33]. We further calculated the S (Q, ω) for the field-induced polarized state, where the overall scaling of the model (given by J) was chosen for best agreement. Within the uncertainty of the broad scattering observed in experiment, this calculation also yields a very good qualitative agreement (Supplementary Figure 11).

III. DISCUSSION
The hunt for the QSL state has yielded numerous studies of interesting new compounds and underlying physical phenomena, even though the ultimate goal may remain elusive. The search for a QSL state in the context of YbZnGaO 4 and YbMgGaO 4 materials seems likely nearing its end, although one cannot claim yet with absolute certainty that the issue is completely settled [56]. In the present work, we have established a potentially close description of YbZnGaO 4 and YbMgGaO 4 in a hypothetical clean limit, allowing for the implications of the disordered, real-world materials to follow. While possibly denied the coveted QSL state, a deeper understanding of the rich and diverse phase diagram of the triangular-lattice antiferromagnets is still a highly satisfying reward that will likely serve as a road map for future studies, in particular for designing new materials or advancing these studies to explore this rich phase diagram.
The present work seeks to advance our understanding of the phase diagram describing YbMgGaO 4 , YbZnGaO 4 , and many related systems. We note that our specific conclusions about these materials are a part of a greater effort to narrow down possible magnetic exchange parameters in these systems [30,31,57,58], and provide a guideline on how to address the issue of chemical disorder in real-world materials. In a very recent work [31] by some of the same authors, the phase crossovers similar to those shown in Figure 1 and Figure 2 were used to constrain parameters of YbMgGaO 4 , with the focus on reproducing them in the disorder-free limit, but without a broader map of the phase diagram beyond the observations offered by an optimization algorithm. Indeed, the methods used in that prior work suggest a general approach to find exchange parameters for disordered systems. By contrast, the goal of this work is to provide an expanded context of the phase diagram in which the observed crossovers can occur in principle, offering a common description for the phenomena observed in both systems. As such, the present work will serve as an invaluable map in the search of QSL and other intriguing states and phenomena in frustrated triangular-lattice compounds, and will importantly benefit the materials design efforts in this field.
With this last point in mind, recent experiments in the ytterbium-based chalcogenides also show promising QSL features and numerous phase transitions induced by an applied magnetic field [2,[7][8][9][10][11]. The precise nature of these transitions is not yet known and needs to be analyzed in a manner and within the framework suggested in the present study. These recent studies also highlight potential further insights that can bring deeper understanding of YbZnGaO 4 and YbMgGaO 4 by experiments with the field along different directions.
Altogether, we have demonstrated the power of the experimental and theoretical insights in identifying relevant parameter spaces of YbMgGaO 4 and YbZnGaO 4 and, potentially, other related materials. We have investigated their fieldinduced phases and have shown that despite their disorderinduced pseudo-SL ground states, we can significantly narrow the allowed regions of their phase diagram that are compatible with the phenomenologies of these materials. Similar consideration can be applied to the other materials such as chalcogenides. More experimental and theoretical investigations, such as targeted materials design, neutron scattering for the in-plane field directions augmented by analytical and numerical studies for this setting, and utilizing alternative tuning parameters such as external pressure, all could be useful for exploring the diverse phase diagram of this group of compounds and shedding further light on their rich physics.

A. Synthesis
Samples were synthesized from finely mixed Yb 2 O 3 (99.9%), ZnO(99.9%), and Ga 2 O 3 (99.999%) powders at 1350 • C. High quality single-crystals (such as pictured in Supplementary Figure 2) were grown using the optical floating zone technique. A typical growth was conducted in 1 MPa O 2 atmosphere with a speed ranging from 4 to 10 mm/hour. Single crystal quality was confirmed via laue x-ray diffraction (Supplementary Figure 3) while powder xray diffraction (PXRD) was used to confirm the correct phase at every step of synthesis (see Supplementary Figure 1).

B. High Resolution Magnetization Measurements
High-resolution measurements of magnetization were achieved with the complimentary tunnel diode oscillator (TDO) technique [31]. 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. A nearly single-crystal sample with dimensions of 2 mm in length and 1 mm in diameter was wound inside a detection coil, with the c axis of the sample aligned with the coil axis. The sample and coil constitute the inductor of the LC circuit. With the application of field, the sample magnetization change induces a change in the inductance, and thus shifts the resonance frequency. This technique enables highly sensitive detection of changes of magnetic moments 10 −15 Am 2 . [59].
Magnetization was also directly measured via an in-house Cryogenic S700X SQUID magnetometer in temperatures down to 300 mK using a 3He probe. A 1.81 mg sample was mounted on a silver straw with vacuum grease in H c and H ⊥ c orientations.
ac-susceptibility measurements were carried out on a long, approximately rectangular sample with field perpendicular to the sample c axis in a dilution refrigerator with a base temperature of 20 mK.

Diffuse magnetic scattering at CORELLI
Diffuse neutron scattering data were collected at the CORELLI spectrometer at Spallation Neutron Source, Oak Ridge National Laboratory [60]. This instrument is a quasi-Laue TOF instrument equipped with a 2D detector, with a -20 • to +150 • in-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 constrained the out-of-plane coverage to ± 8 • . A 0.8 g 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. Neutronabsorbing Cd was used to shield the sample holder to reduce the background scattering. Experiments were conducted with applied fields at the base temperature of 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 [61].
To account for the temperature factor in our background subtraction in total-scattering mode, we compared the ratio of integrated intensities of a rectangular volume of reciprocal space at 5 T for both temperatures. The region was bounded by −0.45 < h < −0.35, 0.5 < k < 0.6, and −2 < l < 2. We used this region, away from the zone edge, and our 5 T data to ensure the comparison was unaffected by diffuse magnetic scattering. This integration approximates the ratio of Bose population factors -our scaling factor was 1.01 ± 0.005. To improve statistics, we used symmetry operations. All analysis and visualization were performed using Mantid and Python.

Inelastic neutron scattering at the CNCS, DCS, and SPINS
We conducted inelastic neutron scattering experiments at the Cold Neutron Chopper Spectrometer [38] (CNCS) at Oak Ridge National Laboratory, and the Disk-Chopper Spectrometer [39] (DCS) and Spin Polarized Inelastic Neutron Spectrometer (SPINS) at the National Institute of Standards and Technology. The same 1.4 g single crystal sample was aligned to use the [h, k, 0] scattering plane and with the field parallel to the sample c axis (along the [001] direction). All measurements were carried out in dilution refrigerators, with sample mounted on copper sample mounts such as pictured in Supplementary  Figure 2.
The base temperatures for CNCS, DCS, and SPINS measurements were 50 mK, 70 mK, and 60 mK, respectively. For CNCS and DCS, the sample was rotated through approximately 180 degrees. The incident energy used at CNCS, DCS, and SPINS was 3.9 meV, 3.55 meV, and 3.7 meV, respectively.
Analysis of CNCS data was carried out in large part using the HORACE software package [45]. Analysis of DCS data was carried out in large part using the DAVE software package [62].

Polarized neutron scattering at BT7
We measured using polarized neutrons at the BT7 [40] beamline at the National Instistute of Standards and Technology. We measured with a fixed final energy of 14.7 meV at about 0.5 meV such that the FWHM of the collimated beam was 1.08 meV and encompassed the energy range of the continuum as pictured for low fields in Figure 3. Samples were again aligned to use the [h,k,0] scattering plane, and measurements were conducted using a vertical guide field (parallel to the sample c-axis). A 3He cryostat was used. 0 T measurements were performed in the absence of a magnet, and then a 7 T magnet was added to perform measurements at 2 T. Flipping ratios ranged from 19 to 33 for the 0 T measurements and from 17 to 28 for the 2 T measurements. Polarization correction was performed using pbcor software.

Theory
For the LSWT S (Q, ω) we used SpinW calculations [49]. The global phase diagram was obtained using classical energies for the single-Q states, see Refs. [19,20], and for the fieldinduced phases, classical energy minimization of the three-and four-sublattice structures was used.
V. DATA AVAILABILITY All relevant data are available from the authors upon reasonable request. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan(http://energy.gov/downloads/doepublic-access-plan).  collected along the Brillouin zone shows generally greater intensity for 2 T (corresponding to the enhanced projection out-of-plane for higher fields) b) As in a) but for the cut extending from M to Γ. At low fields, most scattering is along the zone edge. c) As in a) but showing the spin-flip channel for 0 T and 2 T, showing greater intensity for in-plane components, consistent with a). d) As in b) but showing the spin-flip channel for 0 and 2 T. See methods for details. Uncertainties represent one standard deviation.
FIG. 9. Intensity versus energy at the M point using 0 T -8 T subtraction (both base temperature). Blue curve shows estimated scattering from the elastic line based on the 0.141 meV instrument resolution at the elastic line, which is approximately 16% of the total scattering. Over-subtraction at 1 meV is due to the single-magnon dispersion present in the 8 T data. Uncertainties represent one standard deviation.