Yb delafossites: unique exchange frustration of 4f spin 1/2 moments on a perfect triangular lattice

While the Heisenberg model for magnetic Mott insulators on planar lattice structures is comparatively well understood in the case of transition metal ions, the intrinsic spin-orbit entanglement of 4f magnetic ions on such lattices shows fascinating new physics largely due to corresponding strong anisotropies both in their single-ion and their exchange properties. We show here that the Yb delafossites, containing perfect magnetic Yb$^{3+}$ triangular lattice planes with pseudospin $s=1/2$ at low temperatures, are an ideal platform to study these new phenomena. Competing frustrated interactions may lead to an absence of magnetic order associated to a gapless spin liquid ground state with a huge linear specific heat exceeding that of many heavy fermions, whereas the application of a magnetic field induces anisotropic magnetic order with successive transitions into different long ranged ordered structures. In this comparative study, we discuss our experimental findings in terms of a unified crystal-field and exchange model. We combine electron paramagnetic resonance (EPR) experiments and results from neutron scattering with measurements of the magnetic susceptibility, isothermal magnetization up to full polarization, and specific heat to determine the relevant model parameters. The impact of the crystal field is discussed as well as the symmetry-compatible form of the exchange tensor, and we give explicit expressions for the anisotropic g factor, the temperature dependence of the susceptibility, the exchange-narrowed EPR linewidth and the saturation field.

While the Heisenberg model for magnetic Mott insulators on planar lattice structures is comparatively well understood in the case of transition metal ions, the intrinsic spin-orbit entanglement of 4f magnetic ions on such lattices shows fascinating new physics largely due to corresponding strong anisotropies both in their single-ion and their exchange properties. We show here that the Yb delafossites, containing perfect magnetic Yb 3+ triangular lattice planes with pseudospin s = 1/2 at low temperatures, are an ideal platform to study these new phenomena. Competing frustrated interactions may lead to an absence of magnetic order associated to a gapless spin liquid ground state with a huge linear specific heat exceeding that of many heavy fermions, whereas the application of a magnetic field induces anisotropic magnetic order with successive transitions into different long ranged ordered structures. In this comparative study, we discuss our experimental findings in terms of a unified crystal-field and exchange model. We combine electron paramagnetic resonance (EPR) experiments and results from neutron scattering with measurements of the magnetic susceptibility, isothermal magnetization up to full polarization, and specific heat to determine the relevant model parameters. The impact of the crystal field is discussed as well as the symmetry-compatible form of the exchange tensor, and we give explicit expressions for the anisotropic g factor, the temperature dependence of the susceptibility, the exchange-narrowed EPR linewidth and the saturation field.

I. INTRODUCTION
We would like to discuss our findings and review further available results about an interesting class of compounds -the Ytterbium Yb 3+ delafossites. The delafossites in general emerged from the end of the 19th century onwards [1], the name is structurally motivated and as such has little relationship to the properties of the individual compounds falling into this class -there are insulators, metals, superconductors, semimetals, quasi two dimensional highly conductive materials, magnetic materials, and more. For a large variety of ground states one may think of it is highly probable that we can find representatives in the delafossite class of minerals.
The delafossites form as A 1+ R 3+ X 2− 2 , where A is an alkaline metal (Li, Na, K, Rb, Cs) or a monovalent transition metal ion (Pd, Pt or Cu, Ag), R is a trivalent transition metal or rare earth ion which might be magnetic (like Ti, V, Cr, Fe, Ce or Yb) or nonmagnetic (Al, Ga, In, Tl or Co, Rh), and X stands for a chalcogen which is either oxygen, sulfur, or selenium. Most of the compounds form in the rhombohedral α-NaFeO 2 delafossite structure [2,3] with space group R3m.
Until now, studies on 4f delafossites have been very rare. The reason is that the 4f ions are relatively large and are difficult to incorporate into oxygen-based delafossite structures which are widely investigated. For delafossites with sulfur or selenium, however, this is quite possible due to the larger size of voids formed by the the chalcogen and even the growth of sizable single crystals is possible [4]. An essential characteristic of the delafossite structure is the presence of triangular planes composed of edge-sharing RX 6 octahedra. In this respect they can serve as model systems for quantum magnetism in a per-fect planar triangular lattice. Due to the ideal triangular structure, geometric frustration counteracts or even suppresses magnetic order at low temperatures, eventually supporting spin liquid behavior. Alternatively, at zero temperature, a magnetic order, the so-called 120-degree order, is also predicted [5][6][7][8].
Among the magnetic trivalent transition metal ions only Ti 3+ has an effective spin s = 1/2 doublet ground state. Unfortunately it turns out that compounds based on Ti such as NaTiO 2 show structural instabilities which result in phase transitions and symmetry reductions introducing additional complexity in the interpretation of the results obtained [9][10][11]. Here the importance of the 4f ions comes into play: Among them, the Kramers ions with an odd number of electrons or holes in the 4f shell like Ce 3+ or Yb 3+ have a pronounced doublet ground state due to a low-symmetry crystal electric field (CEF) and can be described with a pseudospin s = 1/2. In this respect we underscore that Yb delafossites are ideal model systems for the study of spin 1/2 triangular lattices.
Our research of the available literature returned fourteen Yb delafossite systems, see Table I. We expect that the number of compounds will increase over time due to the huge interest among the quantum magnetism community. Starting from NaYbS 2 , we have established the series of NaYbCh 2 delafossites as potential quantum spin liquids and will discuss these systems in particular for a comparative analysis, since we consider them to be prototypical [15,18,31].
The most remarkable property of these materials is the absence of magnetic order down to lowest reached temperatures T = 50 mK, suggesting that we might have an experimental realization of the theoretically predicted spin-liquid type ground state [32,33]. Another striking feature is that upon the application of a magnetic field the nonmagnetic ground state transforms into a longrange ordered antiferromagnetic state. Therefore with the Yb delafossites we are dealing with systems in the vicinity of magnetic order which might be tagged as critical spin liquids. This is the crucial difference to the known putative spin liquid candidates like the triangular lattice organic salts, the kagome type herbertsmithite or the recently discovered hyperkagome iridates which are all far away from magnetic order [34][35][36].
In the following sections, we will try to reconcile our theoretical considerations based on a crystal-field plus nearest-neighbor exchange model with the experimental results. In summary the Yb delafossite compounds are interesting unique systems with an ideal triangular lattice structure which, together with the strong spinorbit coupling, leads to unusually large spin and exchange anisotropies. In detail, the triangular crystal field splits the spin-orbit entangled Yb 3+ states into a series of Kramers doublets, the lowest of which in turn results in a complex correlated ground state with a pseudospin s = 1/2. As a consequence also the magnetic exchange between the Yb 3+ ions mediated via the orbitals of the surrounding p states of the chalcogen ions becomes complex and bond dependent, similar to the iridate compounds with honeycomb structure.

II. ONE YTTERBIUM ION
A. Crystallography Table I summarizes the known Yb delafossites. Characteristic for the crystal structure of these is a layered composition of sheets of tilted YbCh 6 octahedra alternating with «filler» planes comprised of alkaline/transition/boron group metal ions. Fig. 1 illustrates the crystal structure of NaYbS 2 as an example. The sulfur octahedra (blue) are tilted such that the Yb 3+ ions inside form perfect triangular lattice planes perpendicular to the crystallographic c direction. Ideally, an octahedron has four equivalent threefold axes perpendicular to the eight pairwise parallel equilateral triangles forming the surface of it. However, the octahedra of all Yb delafossite compounds are distorted in the same manner: one threefold axis is shortened such that each former octahedron is comprised of two «large» parallel equilateral triangles with edge length λ and six «small» isosceles triangles with two edges of length σ and one edge of length λ. The «large» triangles are those oriented perpendicular to the c direction, forming a perfect triangular lattice. For an ideal octahedron we would have σ = λ, the tilting angle of the octahedral axis with respect to the triangular plane would be α = cos −1 1/ √ 3 ≈ 54.74 • . In contrast all delafossites have σ = λ, the tilting angle of the octahedral axis then is given by α = cos −1 λ/( √ 3σ) , see Table I for numbers. The delafossite structure can have two polytypes ac- cording to the orientation of the planar layer stacking. The space group of the rhombohedral 3R type delafossites is R3m wheras the hexagonal 2H types have a space group of P 6 3 /mmc. The difference between the two polymorphs is the stacking of the planar layers in c direction. Most of the Yb-delafossites belong to the R3m space group (Table I). Assigning typical oxidation states, we have Yb 3+ ions with one hole in the 4f shell and A + «filler» ions. The latter mostly are alkaline metals, only one Yb delafossite exists with a metal from the Boron group (Tl) and two Yb delafossites have transition metal filler sheets from the Copper group (Ag and Cu).
B. Yb 3+ in a trigonal crystal field Fig. 2 schematically shows the energy levels of a single Yb 3+ ion. The fourteen 4f 13 states of Yb 3+ with = 3, s = 1/2 are split by the spin-orbit coupling into a j = +s octet and a j = − s sextet. In a perfect octahedral (cubic) environment with ideal tilting angle α, the j = 7/2 states are further split into two doublets Γ 6 and Γ 7 and a Γ 8 quartet, the j = 5/2 states into a Γ 7 doublet and a Γ 8 quartet. Distorting the octahedron along one of its trigonal axes lowers the CEF environment to trigonal, and the local site symmetry of the Yb 3+ ions is C 3v with the threefold axis parallel to the c direction. This transforms the formerly cubic states like Γ 6 → Γ T 6 , Γ 7 → Γ T 6 , and Γ 8 → Γ T 4 + Γ T 5 + Γ T 6 . We note that although the two representations Γ T 4 and Γ T 5 are one-dimensional, due to Kramer's theorem they are complex conjugates and correspond to the same energy. Apart from a constant the Hamiltonian for a single Yb 3+ ion at an arbitrary lattice site i is then given by where B m n are crystal-field parameters and O m n (J) are Stevens operators, being polynomials of the components of the total-momentum operator [37,38]. They are reproduced in Appendix A.
To gain insight into the structure of the wavefunctions and energy levels, let's for a moment assume the Yb 3+ ion resides in an ideal octahedron. The local symmetry of the Yb 3+ ion then is cubic with O h symmetry, and additional relationships between the B m n crystal-field parameters apply. Choosing the trigonal axis introduced above as the quantization axis of the hypothetic ideal delafossite, Eq. (1) reduces to [37,38] An explicit expression for the matrix of this Hamiltonian for j = 7/2 in the |j, m basis is given in Appendix B.
Only two independent crystal-field parameters remain. The corresponding wavefunctions of H (cubic,3) CEF which are grouped in Kramers pairs consisting of time reversed states are given by operator, giving T j, ± 1 2 = (−) j±1/2 j, ∓ 1 2 with T 2 = −1 for half-integer total momentum j.
The fourth doublet is the pure state 7 2 , ± 3 2 , adiabatically (as function of the trigonal distortion) connected to the pure state in the Γ 8 quartet.
One doublet defined in Eqs. (4) will be the ground state. In principle 7 2 , ± 3 2 could be the ground-state doublet as well. However, we will see at the end of this section why in our case it is not. With the help of Appendix A, the matrix elements of the total momentum operator J within the ground-state doublet can be read off: All other matrix elements vanish. The total momentum operator J transforms like J = −T JT −1 under time reversal T . With T = U K, U = exp(iπJ y ) in our (standard) representation and K the complex conjugation, this requires the matrix elements of J z and J x to be real, those of J y to be purely imaginary, which is equivalent to φ α = φ γ +2πn, n ∈ Z. This indeed allows us to introduce a pseudospin S for the ground-state doublet by mapping The Zeeman splitting for a single Yb 3+ ion at site i is then obtained from the Hamiltonian where µ 0 is the magnetic permeability constant and µ B the Bohr magneton.
The above equations give an argument why the pure state 7 2 , ± 3 2 cannot be the ground state here: We would have transverse matrix elements 7 2 , ± 3 2 |J x,y | 7 2 , ∓ 3 2 = 0, i. e. g ⊥ ≡ 0, and no coupling to a magnetic field applied perpendicular to the threefold axis enforced by trigonal symmetry alone. As we will see, this is in contradiction to experiment. To the best of our knowledge, no Yb compound is known to have the 7 2 , ± 3 2 doublet as ground state.

D. Electron paramagnetic resonance
We have done thorough electron paramagnetic resonance studies on all three NaYbCh 2 delafossites [15,16,18]. Table II contains our experimental findings on the g factors. Common remarkable feature is the strong anisotropy between in-plane and out-of-plane gyromagnetism. This is shown exemplarily for NaYbS 2 in Fig.3. In a certain temperature range, the resonance intensity I EPR is antiproportional to the temperature [29]. Regarding I EPR as a measure for the resonant susceptibility χ R , we can correspondingly assign characteristic temperatures Θ ,⊥ , see Table II. We note that these temperatures are not necessarily identical with the Curie-Weiss temperatures obtained from susceptibility measurements (see below).
Towards high temperatures the temperature dependence of the resonance linewidth ∆B(T ) can be understood with a spin-lattice relaxation through the modulation of the ligand field by the lattice vibrations. The basic mechanism behind this relaxation is the spin-orbit coupling which make the electron spins «feel» the ligand field modulation [41, p. 60 ff]. For the temperature dependence of this relaxation various processes are involved among which an exponential temperature dependence identifies a two-phonon process [42]. In this so-called Orbach process the thermal equilibrium of the Zeeman split ground doublet is achieved by a phonon absorption exciting the spin system to an upper state at energy ∆E and then a phonon emission back to the ground state, the absorption and emission energies differing by the Zeeman energy. This process is determined by the number of phonons at energy ∆E and yields for ∆E k B T an approximate temperature dependence ∝ exp(−∆E/k B T ). Identifying ∆E with the energy ∆E 12 of the first excited crystal-field split state of the Yb 3+ ion the data analysis of ∆B(T ) gives a rough estimate (within ±2.5 meV) of ∆E 12 , values of which are also shown in Table II. Table II contains results on inelastic neutron scattering as well. For NaYbO 2 and NaYbSe 2 , three clear maxima in energy scans at temperature T = 5 K have been observed [43,44], associated with the energy differences ∆E 1j between the ground state ψ ± 1 and the excited states ψ ± j , j = 2, 3, 4. Values are quoted in the table. We investigated the excited CEF doublets of NaYbS 2 in a time-of-flight neutron scattering experiment [31], see

E. Inelastic neutron scattering
FIG. 3. EPR spectra (upper frames) and anisotropy of the EPR g factor (lower frame) of single crystalline NaYbS2 at T = 19 K and the microwave field bmw⊥c axis (ν = 9.4 GHz). The spectra taken at external fields B⊥c axis and B c axis were fitted by a Lorentzian lineshape (dashed lines). The anisotropy of the EPR g factor can be described by g(Θ) = g 2 cos 2 Θ + g 2 ⊥ sin 2 Θ with g ⊥ = 3.19(5) and g = 0.57(3). The sample was rotated around an axis lying in the basal plane parallel to bmw. scattering intensity, pointing to excited doublets located 23 meV and 39 meV above the ground state. A third maximum is missing. However at higher temperatures T ≥ 50 K additional intensity appears at about 11 meV and 27 meV in the energy scan. If we attribute these features to excitations from the thermally populated first excited doublet to the two higher doublets already observed at low temperatures, we can assume the first excited doublet to be located at an excitation energy of ∆E 12 ≈ 12 meV above the ground state.
Two qualitative observations can be made here: The higher the g factor anisotropy, (i) the lower the CEF excitation energies, and (ii) the smaller the transition matrix element from the ground state to the first excited doublet is.
As we will see later from analyses of susceptibility data, the typical exchange energy scale of the NaYbCh 2 delafossites is of the order of a few Kelvin. With the CEF excitations being two orders of magnitude higher, this qualifies the pseudospin description of the ground-state doublet introduced in the prior section also for the min-  imum exchange model discussed in the following.

III. MANY YTTERBIUM IONS
A. High-temperature magnetic susceptibility The dimensionless uniform magnetic susceptibility χ(T ) of a crystal with volume V is given by the change of the magnetisation M with the magnetic field B = µ 0 H with components whereμ α = −∂F/∂B α is the total magnetic moment in spatial direction α either parallel or perpendicular to the c axis and F = (1/β) log Z the canonical free energy, 1/β = k B T the inverse temperature and k B the Boltzmann constant. For the molar susceptibility Eq. (6) has to be multiplied by N L /(ν/V ) where N L is Avogadro's number and ν/V the volume density of Yb 3+ ions. We link the free energy in the usual way with statistical mechanics through the partition function Z = exp(−βH) where with H CEF given by Eq. (1), H Zeeman given by Eq. (5), and the exchange Hamiltonian for an arbitrary but fixed site i where the sum is taken over the z = 6 bonds connecting sites j and site i with an exchange tensorĴ having the respective componentsĴ αβ ij [46]. for more details on the symmetry-allowed form of the exchange.) A factor (1/2) is included to compensate for double-counting the bonds when executing the sum over the lattice sites in Eq. (7). In the high-temperature limit β → 0 we can expand χ α (Eq. (6)) in powers of β. We obtain a Curie-Weiss law with j = 7/2 and the Curie-Weiss temperatures given by Here zĴ and zĴ ⊥ are the contributions of the exchange tensor parallel and perpendicular to the c direction.
Other components of the exchange tensor do not appear in χ α up to order β 2 . This result coincides with Ref. [47] where a four-parameter crystal field Hamiltonian has been treated. Remarkably, due to the orthogonality and tracelessness of the Stevens operators, only the B 0 2 CEF parameter enters the Curie-Weiss temperatures, independent of the form and symmetry of the crystal field otherwise as long as it contains a more-than-2-fold symmetry axis. Figure 5 shows the temperature dependence of the susceptibilities χ ,⊥ (T ) for NaYbSe 2 in the full temperature range T = 0.5 . . . 400 K accessible to us where we have applied the magnetic field in directions parallel to the c axis (label H c) and perpendicular to it. At temperatures T > ∼ 120 K, the inverse χ −1 ,⊥ (T ) (not shown) show a linear temperature dependence. The solid line in the left plot of Fig. 5 denotes a corresponding Curie-Weiss fit to χ ⊥ (T ) according to Eq. (9) for T ≥ 150 K. From this fit we obtain µ eff ≈ 4.6 µ B , reflecting the full effective moment µ eff = g j µ B j(j + 1) ≈ 4.54 µ B . We also obtain a Curie-Weiss temperatureΘ ⊥ ≈ −66 K. In the temperature range where the Curie-Weiss law is a good approximation, χ α (T ) is essentially isotropic withΘ ≈Θ ⊥ , independent of the direction of H. The same fits for NaYbO 2 (powder) and NaYbS 2 yield the full moment and isotropic Curie-Weiss temperatures as well. Results are noted in Table II. Below T = 80 K the s = 1/2 pseudospin state emerges. After subtracting a temperature independent Van-Vleck contribution (obtained from high field magnetization measurements) the susceptibility χ(T ) for 10 K ≤ T ≤ 40 K can be fitted with a Curie-Weiss law (right-hand plot of Fig. 5) which yields a Curie-Weiss temperature Θ ⊥ = −7 K, and an effective moment µ ⊥ = 2.43µ B and a Curie-Weiss temperature Θ = −3.5 K and a moment of µ = 1.1µ B for fields in the ab plane and in the c direction, respectively. These low-temperature effective moments are consistent with the measured g values from EPR within the s = 1/2 pseudospin model.
The pronounced maximum in the susceptibility (Fig. 5) corresponds to the maximum found in the temperature dependence of the magnetic specific heat c m (T ) (Fig. 6). Such a maximum is expected in the isotropic triangular lattice [48]. At the maximum, the thermal energy roughly corresponds to the exchange coupling energy of the spin system. In the susceptibility it is clearly visible that for fields in c direction the maximum is at lower temperatures. This is to be expected since the magnetic coupling is smaller in this direction. Fig. 6 shows the magnetic specific heat for all three NaYbCh 2 compounds. In contrast to the susceptibility, the maximum is broadened here. At temperatures right below the maximum, this specific heat decreases with T 2 as expected for two-dimensional magnon-like states, and then a linear dependence c m (T ) = γT down to the lowest accessible temperatures with a large residual value γ ≈ 1 J/molK 2 is found. This linear temperature dependence is well known for heavy-fermion systems, for example YbRh 2 Si 2 [49] or Yb 4 As 3 [50,51], yet typical for a gapless spin liquid with fermionic excitations [36].

B. Pseudospin exchange model
We restrict ourselves now to the CEF ground-state doublet of each Yb 3+ ion and use the pseudospin description introduced above. The simplest model Hamiltonian includes nearest-neighbor exchange on the triangular lattice (six neighbors) only. Inspecting one Yb-Yb bond, we see that it contains a twofold rotation axis, a mirror plane in the middle of the bond, and a center of inversion in the middle of the bond. (For an ideal delafossite with undistorted octahedra, we have an additional mirror plane containing the basal plane of an octahedron and a further mirror plane perpendicular to it.) According to Moriya's rules [52], antisymmetric exchange must vanish due to the presence of the inversion center. This leaves us with four independent components of the exchange matrix (three for the ideal case) on any bond ij . We choose our coordinate system such that one bond is par-allel to the x direction and the z axis is perpendicular to the triangular-lattice planes. For this bond, we can write whereby we split the exchange matrix into a rotationally invariant part (rotations around the z axis) plus a traceless directional-dependent part. We note that for finite J yz the cartesian coordinate axes y and z are not the main axes of the exchange tensor, rather all three main axis components are different forJ ij = U −1 J ij U , U unitary,J ij diagonal. This is a direct consequence of the trigonal distortion (tilting angle α = cos −1 1/ √ 3 defined in Sec. II A) of the YbCh 6 octahedra.
The full pseudospin Hamiltonian with this parametrization, expressed with ladder operators instead of cartesian spin operators then reads with H Zeeman given by Eqs. (5) and the directiondependent phases are We note that the direction-dependent terms in Eq. (12) in contrast to the rotationally invariant terms contain at least one spin flip coupling states with ∆s z = 1 or ∆s z = 2. The three-dimensional ground-state phase diagram of this Hamiltonian has been addressed by several authors [32,33,[53][54][55][56][57] using slightly differing parameterizations (see Appendix C for examples). Zhu et al. [33] have shown that indeed spin-liquid type regions in the groundstate phase diagram of the model (12) without magnetic order may exist. However these nonmagnetic regions are comparatively small and require a special choice of the exchange constants. Much more common are magnetically ordered states like stripe phases and the 120-degree pattern known for the isotropic case. Nonetheless of the 14 compounds listed in Table I at least nine do not show any signature of long-range magnetic order down to the lowest measurement temperature, typically either 0.4 K or 2 K. Given the narrowness of the nonmagnetic regions of the model (12), it would be surprising if all compounds can be described with exchange parameters leading to ground states in those regions.
We are faced with a number of possible issues: First, crystallographic peculiarities. A symmetry-allowed buckling of the YbCh 6 octahedral planes possibly introduces additional exchange frustration not contained in the model above. Another problem are stacking faults: One unit cell contains three crystallographically equivalent Yb 3+ ions, each being member of a different YbCh 6 distorted-octahedra plane. Adjacent planes are stacked in an A-B-C like fashion where the projections along the c direction of the positions of the Yb 3+ ions of the «next» layer fall in the middle points of the triangular lattice formed by the «current» layer. This stacking might be distorted. Compatible with these two effects is a sample dependence of the EPR data for NaYbS 2 we have observed: The measurements were made for single crystals from two different batches [16]. While the resonance field is essentially identical in both measurements, a pronounced difference has been observed in the width of the EPR resonance. For the smaller crystal, a description of the latter using two Lorentzian lines has been necessary, indicating that roughly half of all spin probes has a larger linewidth than the other half. However, further investigations have to be undertaken to clarify this.
Second, a trivial reason for no magnetic order would be a frustration-induced extremely low ordering temperature T N . According to Ref. [58], frustration ratios f = |Θ CW | /T N [59,60] of the order of 10 can be achieved already for an isotropic Heisenberg exchange model on the triangular lattice with J ⊥ = J and J ∆ = J yz = 0 and a small interplane exchange coupling J inter . With typical Curie-Weiss temperatures of a few Kelvin, it might be that T N falls into the few 100 mK range, however this requires extremely small J inter = O(10 −4 J ,⊥ ).
Third, we cannot exclude off-stoichiometric Yb 3+ ions, and also Na vacancies might be present. Both effects introduce an unknown amount of disorder in the exchange constants, suppressing the magnetic order [61].
Fourth, the perfect threefold symmetry of the magnetic sublattice, including the inversion center in the middle of an Yb-Yb bond, might be distorted, introducing changes in bond angles and additional nonzero elements in the exchange matrix. Synchrotron data taken at low temperatures would clarify that.
Fifth, related to the comparatively large local moment of the Yb 3+ ions in the ab plane (see Table II), the impact of the long-ranged dipole-dipole interaction might be a further reason for a suppression of T N below our accessible temperature range [62].
Finally it might well be that a pure nearest-neighbor exchange model is not sufficient. Additional competing exchange between further neighbors might as well lead to a suppression of a magnetically ordered ground state [63,64]. Nevertheless we continue using this model for reasons which will become clear later.

C. Electron paramagnetic resonance
The exchange-narrowed linewidth of the EPR resonance in general is given by where const. = π/ √ 3 for a cutoff Lorentzian lineshape, and const. = √ 2π for a Lorentzian × Gaussian lineshape [65]. θ is the angle of the applied field H relative to a given crystallographic direction, for example the c axis. M 2 and M 4 denote the second and fourth moment of the EPR lineshape function, respectively, given by [66] Here the z axis not necessarily corresponds to the crystallographic threefold c axis but rather is defined by the direction of the applied field. Appendix D contains more details on the calculation of M 2,4 . For a Gaussian lineshape, all odd moments vanish, the higher even moments all factorize into powers of the second moment, and we have Using the Hamiltonian (12), we obtain lim T →∞ where the symbols , ⊥ denote the direction of the applied magnetic field µ 0 H relative to the crystallographic c direction. We note that for a fully isotropic exchange Hamiltonian the expressions above vanish identically. In this case, a finite linewidth is due to dipole-dipole interaction which by its very nature is anisotropic. Furthermore we read off M 2 = 2M ⊥ 2 for a rotationally invariant exchange (J ∆ = J xy = 0) which can be understood in the following way: Thermal spin fluctuations are energetically favorable perpendicular to the field µ 0 H, maximizing the energy gain due to the Zeeman energy by leaving the component of the total moment (anti-) aligned to the field unchanged. Thermal fluctuations out of the crystallographic ab plane are suppressed because they would break the threefold rotational symmetry around the c axis. For a field parallel to c, we therefore have two possible fluctuation directions, for a field in the ab plane only one fluctuation direction remains.
More algebra has to be done to calculate the fourth moment. We eventually obtain lim T →∞ Similar to M 2 , also M 4 vanishes for a fully isotropic exchange Hamiltonian, and M 4 = 2M ⊥ 4 for a rotationally invariant exchange.
We assume J ∆ , J yz J ⊥ , J . Taking into account only finite exchange constants J ⊥ and J , we get for the high-temperature EPR linewidth, Eq. (14) lim T →∞ for field parallel and perpendicular to the c axis. An estimate of lim T →∞ ∆H ,⊥ from our linewidth data is difficult to obtain, because in this limit, phonon-dominated relaxation mechanisms like the Orbach process discussed above might dominate. In particular we have lim T →∞ ∆H /∆H ⊥ = 2g ⊥ /g . This relation is roughly consistent with the experimental EPR data. Estimating lim T →∞ ∆H with the smallest value the linewidth reaches in its temperature dependence we obtain for NaYbCh 2 , Ch = O, S, Se: lim T →∞ ∆H /∆H ⊥ = 1.7/10/4.6. From the EPR data in Table II, we obtain 2g ⊥ /g = 3.7/11.1/6.2.

D. Magnetization and susceptibility
In order to learn more about the size of the exchange constants, we have measured the temperature-dependent uniform magnetic susceptibility χ(T ) of the pseudospins at sufficiently low temperatures T < 30 K and the magnetization M (H) in applied magnetic fields B = µ 0 H up to B = 30 T. The dimensionless magnetic susceptibility is given by Eq. (6). Here we expand the thermal traces for the pseudospin Hamiltonian Eq. (12) in the low-temperature limit β → 0 and obtain where the Curie-Weiss temperatures Θ ,⊥ are given by with s = 1/2 for a field applied parallel and perpendicular to the crystallographic c direction. Here J α i,i+n denotes the component of the exchange energy between site i and its nth neighbor along the field direction α. We note that Θ ,⊥ do not depend on the direction-dependent terms in the Hamiltonian (12) -the exchange constants J ∆ and J yz only appear in higher orders of the expansion (compare also Ref. [53]). Table II holds our findings from susceptibility measurements [15,18]. It contains the measured Curie-Weiss temperatures Θ ,⊥ and the effective moments µ ,⊥ obtained from a Curie-Weiss fit to χ(T ) at T < 30 K. For NaYbO 2 only powder samples were available, the corresponding averaged values are listed in the respective x rows. We can calculate the exchange constants for the parameters of the pseudospin Hamiltonian in Eq. (12) from the Curie-Weiss temperatures given by Eqs. (27), and with the effective moments given by µ ,⊥ = g ,⊥ s(s + 1), we can determine the effective g factors. Results are listed in Table III.  Table II also holds the values for the effective moment µ eff and the Curie-Weiss temperaturesΘ ,⊥ for T > 150 K, in the high-temperature limit the latter are given by Eqs. (10). From the derivation of the ground state pseudospin in Sec II C we have a relationship with the exchange constants introduced in Eq. (8) for the full angular momentum. Together with the EPR g factors given in Table II, this allows us to roughly estimate the B 0 2 CEF parameter of NaYbS 2 and NaYbSe 2 to B 0 2 ≈ −1 . . . −0.5 meV.  The saturation field H sat of NaYbS 2 parallel to the c axis was too high to be reached in our experiments.

E. Saturation field
H sat is defined as an instability of the fully polarized state towards ∆m s = 1 spin flips (magnons). It can be calculated within a classical approximation which is described in Appendix E. We parameterize the spins as moment vectors on three interpenetrating sublattices where on each sublattice, the moments are aligned pairwise parallel, and the classical energy density e = E/(νs) is a function of three pairs of polar and azimuthal angles of the sublattice moments.

Field parallel to the c direction
A magnetic field applied parallel to the crystallographic c direction preserves the rotational symmetry of the Hamiltonian (12). For s = 1/2, three possible configurations of the sublattice moments near saturation have been found [67][68][69], called 0-coplanar (0 < J ⊥ /J < ∼ 3/2), π-coplanar (3/2 < ∼ J ⊥ /J < ∼ 2.2) and umbrella phase (J ⊥ /J > ∼ 2.2) which we assume here. The left sketch in Fig. 7 illustrates this particular configuration. At any finite field value, the sublattice moments arrange on a cone such that the projections onto the ab plane form a 120-degree structure. Polar angles are all equal, θ i = ξ. Near saturation with ξ 1, the classical ground-state energy density is given by This demonstrates the competition between the best possible antiferromagnetic alignment of the spins with respect to each other and the alignment parallel to the field direction. The first two terms yield the groundstate energy per spin for the fully polarized state: We gain energy h , but there is an energy loss 3sJ because of the perfect misalignment (i. e. ferromagnetic alignment) of the moments on all three sublattices. The term ∝ ξ 2 shows what happens when the umbrella opens a bit: We gain energy (3/2)sJ ⊥ ξ 2 from the (small) 120-degree projections of the spins onto the ab plane perpendicular to the field. This energy gain is half of what we would obtain if it were possible to align the projections in a Néel like manner, impossible on a triangular lattice. We also gain energy 3sJ ξ 2 due to the (small) reduction of the pairwise ferromagnetic spin alignment in the direction parallel to the field on the three sublattices. For the same reason, we lose energy (1/2)h ξ 2 . With the saturation magnetization M ,⊥ sat = sg ,⊥ µ B /Yb 3+ (see Table II for the experimental values), we obtain

Field in the ab plane
For J = J ⊥ , the magnetic phase diagram has three intermediate phases between H = 0 and saturation [70][71][72]. Like in the case H c, we assume an umbrella-shape spin structure right below the saturation field, however the umbrella has no rotational symmetry around the axis set by the direction of the magnetic field and might degenerate to a planar configuration [73][74][75][76][77][78]. In Appendix E, we parameterize a spin on sublattice i with its polar and azimuthal angles θ i and φ i . Near saturation, θ i → π/2, φ i → α where α is the angle of the applied field relative to the a axis, and we obtain for the saturation field in the ab plane. This is the case for J < J ⊥ . The independence of H ⊥ sat from J suggests that a planar spin configuration with all spins in the ab plane is energetically favorable at and at least infinitesimally below H sat . At finite H < H sat , a distorted cone might form. Writing δ i = π/2 − θ i , i = α − φ i , the six angles are a solution to Eq. (E2) with δ i , i = 0.
A geometric interpretation can be obtained by looking at two possible planar configurations near H sat , sketched in the middle and right illustration of Fig. 7: the sublattice moments form a fan either in the ab plane setting δ i = 0 or in a plane perpendicular to it containing the c and the field axis with i = 0. For the former, the gradient near H sat is given by Minimizing this gives From u = −1 we obtain h ⊥ = 9sJ ⊥ (first factor in Eq. (E8)), for u = 1/2 we obtain the unphysical solution h ⊥ = 0. The ground-state energy density is e ⊥ = 3sJ ⊥ − h ⊥ . Similar to the case H c, we lose energy 3sJ ⊥ due to the ferromagnetic alignment of all spins, and we gain energy h ⊥ . Deviating slightly from full polarization, but still with all sublattice moments in the ab plane we gain energy ∆e ⊥ ∝ 3sJ ⊥ | i | due to both the small antiferromagnetic component in the spin alignment and the reduction of the moment parallel to the magnetic field. The gradient near H sat for a spin configuration in the plane containing the c axis and the field ( i = 0) is given by Minimizing this gives From u = −1 we obtain h ⊥ = 3s J + 2J ⊥ (second factor in Eq. (E8)). This is the case for J > J ⊥ . Deviating from full polarization now results in an energy gain ∆e ⊥ ∝ s(J + 2J ⊥ )|δ i |, which is, for J < J ⊥ , smaller than the energy gain when opening in the ab plane. For u = 1/2 we obtain h ⊥ = 6s J ⊥ − J . This solution is unphysical because this would include an energy loss ∝ 2sJ when deviating from full polarization. In the last six rows of Table III we note the values obtained from the equations above. For NaYbS 2 , J cannot be determined because of the missing value for H sat . If we take the exchange constants derived from I EPR (T ) and χ(T ), we can roughly estimate 56 T < ∼ µ 0 H sat < ∼ 112 T which is indeed a larger range than experimentally accessible for us.

IV. SUMMARY AND CONCLUDING REMARKS
In summary, we have shown that the series of NaYbCh 2 compounds contain nearly perfect magnetic Yb 3+ triangular lattice planes. The determination of the crystal field levels by neutron scattering has shown that at temperatures T < ∼ 100 K an application of a pseudospin model is justified. A consistent treatment in the framework of the anisotropic Heisenberg triangular lattice model with next nearest neighbor (nn) coupling of Yb 3+ pseudospins of 1/2 gives good agreement with some experimental results. In particular, the estimation of the main exchange energies (rotationally invariant elements of the exchange matrix) and the derived saturation fields are in good agreement. Nevertheless, the question to the origin of the absence of magnetic order and the emerging spin liquid state is not completely answered. In principle, off diagonal contributions in the exchange matrix can be responsible for it. However, it is not possible to estimate these off diagonal contributions on the basis of the available data. This is not unusual and is also true for the other prominent Yb 3+ triangular lattice system YbGaMgO 4 [53,79].
Furthermore, the next-nearest neighbor exchange (nnn) may also play a role. It is already known from the Heisenberg lattice without spin-orbit coupling that the additional frustration introduced with an antiferromagnetic coupling can lead to a quantum spin liquid phase where in most cases already a nnn coupling which is one order of magnitude smaller than the nearest-neighbor coupling is sufficient [64]. However, we cannot exclude that the dipolar interaction between the Yb 3+ ions also plays a role like for example in the YbAlO 3 quantum magnet [53]. This long-ranged interaction as an additional competing effect could be another reason for the strong frustration and the suppression of the order.
It must also be noted that the effective exchange constants J and J ⊥ in the NaYbCh 2 materials are comparatively large for an Yb 3+ triangular lattice. In YbMgGaO 4 , NaBaYb(BO 3 ) 2 and Rb 3 Yb(PO 4 ) 2 the lowtemperature Curie-Weiss temperatures and thus the exchange constants are more than one order of magnitude smaller compared to the NaYbCh 2 system [53,[79][80][81]. In this sense, these other systems with small exchange couplings can be regarded as nearly-single-ion systems (without significant exchange coupling between the ions). This is most evident in the magnetic part of the specific heat c m (T ). Here a Schottky peak develops in an applied magnetic field due to the Zeeman splitting of the Yb 3+ CEF ground-state doublet. This peak shifts with the field to higher temperatures. Also, the saturation field in the magnetization M (H) is in the range of a few Tesla, while in the NaYbCh 2 systems the saturation fields exceed 10 T. Furthermore, we have found a pronounced and large linear temperature dependence of the magnetic specific heat c m (T ) = γT with a residual value γ ≈ 1 J/molK 2 for the NaYbCh 2 systems. This is a crucial feature of the (gapless) quantum spin liquid ground state and is consistent with the residual fluctuations detected by muon spin relaxation (µSR), nuclear magnetic resonance (NMR) and finally inelastic neutron scattering. For the other Yb 3+ triangular lattice materials mentioned above, such residual contributions are absent or negligibly small. The presence of the large γ term is a clear evidence for a gapless spin liquid ground state. In analogy to correlated 4f heavy fermion systems with enhanced renormalized electronic density of states at the Fermi level we have an enhanced (renormalized) density of magnetic «fermion-like» states due to fluctuations associated with the generic spin liquid ground state.
Another difference to other Yb 3+ triangular lattices is the occurrence of field induced order. This shows that competing interactions are responsible for the spin liquid state and place the systems in the vicinity of a critical point. As already mentioned, this might originate from small but finite symmetry-compatible off-diagonal components in the Yb-Yb exchange matrix and/or a possible next-nearest neighbor interaction.
Taken together, we have successfully shown that the NaYbCh 2 delafossites have their own fascinating physics which differs significantly from the previously known pla-nar Yb 3+ triangular lattices. There is both exchange anisotropy and spin anisotropy in the systems which is an an essential ingredient and enhancement for frustration and absence of magnetic order. The anisotropy is particularly pronounced in NaYbS 2 with a ratio of the coupling coefficients of J ⊥ /J = 6.25 and a ratio of the EPR g factors of g ⊥ /g = 5.6. We are waiting to see further intriguing developments in the field of the 4f delafossites. For example, the series could be extended to include compounds with the chalcogen Te. These new compounds might exhibit a smaller band gap or even be metallic due to extended Te 5p orbitals. This would establish a bridge between the quantum spin liquid in the Mott insulator and a Fermi liquid in a correlated (semi) metal.

ACKNOWLEDGMENTS
We are grateful to Sasha Chernychev, Jeff Rau, Oliver Stockert, and Peter Thalmeier for elucidating discussions and valuable comments. This work has been partially funded by the German Research Foundation (DFG) within the Collaborative Research Center SFB 1143 (Projects B03, C02,C03, and A05). , the corresponding wavefunctions are given by Eqs. (3).
In a trigonal CEF, all six crystal-field parameters introduced in Eq. (1) are independent, and we obtain for the j = 7/2 matrix representation of the Hamiltonian Needless to say that closed-form expressions for the eigenvalues, apart from the pure |7/2, ±3/2 doublet, are lengthy and not very insightful.
If we regard the trigonal distortion of the ideal octahedron as small, we obtain corrections to the cubic eigenvalues to first order in the CEF parameters δB m n like , dependent terms ∝ J ∆ , J yz cancel exactly upon summation over the three sublattice pairs (1, 2), (2, 3), and (3,1). According to Eq. (12), J ∆ and J yz give the energy gain when coupling states differing by ∆S z = 1 or ∆S z = 2. For the classical model, the total spin S tot z = ν i=1 S z i is conserved, therefore J ∆ and J yz cannot contribute to the ground-state energy. The energy is symmetric with respect to the exchange of any two sublattice labels i = 1, 2, 3, so the saturation field (as any other ground-state property) cannot depend on it.
In the limit δ i , i → 0 (full polarization), Eq. (E2) for m H (e ⊥ ) reads In our case, we have J ⊥ > J > 0, and the first factor determines the saturation field when lowering the field from values above saturation. We obtain