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

B. Schmidt  J. Sichelschmidt  K. M. Ranjith  Max-Planck-Institut für chemische
Physik fester Stoffe, Dresden, Germany  Th. Doert  Fakultät für Chemie und
Lebensmittelchemie, Technische Universität Dresden, Germany  M. Baenitz  Max-
Planck-Institut für chemische Physik fester Stoffe, Dresden, Germany

(Typeset )

###### Abstract

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 Yb3+ 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 Yb3+ delafossites. The
delafossites in general emerged from the end of the 19th century onwards
111The original name was given by the French mineralogist Charles Friedel
(1832–1899) to the natural mineral CuFeO${}_{\text{2}}$ in honor of Gabriel
Delafosse (1796–1878) [83, 84] in 1873., 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 $\alpha$-NaFeO2
delafossite structure [2, 3] with space group $R\bar{3}m$.

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 $R$X6 octahedra. In this respect they can serve as model systems for
quantum magnetism in a perfect 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 Ti3+ has an effective
spin $s=1/2$ doublet ground state. Unfortunately it turns out that compounds
based on Ti such as NaTiO2 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 Ce3+ or Yb3+ 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.

Table 1: Yb delafossites. Lengths of the long ($\lambda$) and short ($\sigma$) edges of the triangles forming the distorted YbCh${}_{\text{6}}$ octahedra for some MYbCh${}_{\text{2}}$ compounds with M being a metal and Ch being a chalcogenide. Next column: tilting angle $\alpha=\cos^{-1}\left(\lambda/(\sqrt{3}\sigma)\right)$ of the octahedral rectangles with respect to the triangular plane. Last column: Those compounds marked with an asterisk have no magnetic long-range order down to the respective lowest investigation temperature. For all others, this has not (yet) been investigated. Last row: undistorted octahedron. compound | $\lambda/\textup{\AA}$ | $\sigma/\textup{\AA}$ | $\alpha$ | reference | space group | remark
---|---|---|---|---|---|---
LiYbS${}_{\text{2}}$ | $3.80800$ | $3.24599$ | $47.35^{\circ}$ | [12, 13] | $R\bar{3}m$ | *
NaYbO${}_{\text{2}}$ | $3.3510$ | $2.87472$ | $47.70^{\circ}$ | [14, 15] | $R\bar{3}m$ | *
NaYbS${}_{\text{2}}$ | $3.90400$ | $3.40822$ | $48.60^{\circ}$ | [12, 16] | $R\bar{3}m$ | *
NaYbSe${}_{\text{2}}$ | $4.05680$ | $3.92235$ | $53.33^{\circ}$ | [17, 18] | $R\bar{3}m$ | *
KYbO${}_{\text{2}}$ | $3.40010$ | $3.01711$ | $49.41^{\circ}$ | [19] | $R\bar{3}m$ | *
KYbS${}_{\text{2}}$ | $3.96800$ | $3.62569$ | $50.81^{\circ}$ | [12, 20] | $R\bar{3}m$ | *
KYbSe${}_{\text{2}}$ | $4.11100$ | $5.06809$ | $62.07^{\circ}$ | [17] | $R\bar{3}m$ |
RbYbO${}_{\text{2}}$ | $3.41000$ | $3.16297$ | $51.51^{\circ}$ | [21] | $R\bar{3}m$ |
RbYbS${}_{\text{2}}$ | $3.99100$ | $3.71481$ | $51.66^{\circ}$ | [22] | $R\bar{3}m$ |
CsYbS${}_{\text{2}}$ | $4.02200$ | $3.70764$ | $51.22^{\circ}$ | [23] | $R\bar{3}m$ |
CsYbSe${}_{\text{2}}$ | $4.15390$ | $3.88335$ | $51.86^{\circ}$ | [24, 25] | $P6_{3}/mmc$ | *
TlYbS${}_{\text{2}}$ | $3.93500$ | $3.62771$ | $51.23^{\circ}$ | [26] | $R\bar{3}m$ | *
| $3.9454$ | $3.726$ | $52.31^{\circ}$ | [27] | $P6_{3}/mmc$ | *
AgYbO${}_{\text{2}}$ | $3.44040$ | $2.99493$ | $48.45^{\circ}$ | [28, 29] | $R\bar{3}m$ | *, linear Ag coordination
CuYbSe${}_{\text{2}}$ | $4.01670$ | $3.93641$ | $53.91^{\circ}$ | [30] | $P\bar{3}m1$ | tetrahedral Cu coordination
undistorted | $1$ | $1$ | $54.74^{\circ}$ |  | (cubic) |

Our research of the available literature returned fourteen Yb delafossite
systems, see Table 1. We expect that the number of compounds will increase
over time due to the huge interest among the quantum magnetism community.
Starting from NaYbS2, we have established the series of NaYbCh2 delafossites
as potential quantum spin liquids and will discuss these systems in particular
for a comparative analysis, since we consider them to be prototypical [31, 15,
18].

The most remarkable property of these materials is the absence of magnetic
order down to lowest reached temperatures $T=50\,\rm 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
long-range 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 spin-orbit coupling, leads to unusually large spin
and exchange anisotropies. In detail, the triangular crystal field splits the
spin-orbit entangled Yb3+ 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
Yb3+ 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

###  II.1 Crystallography

Table 1 summarizes the known Yb delafossites. Characteristic for the crystal
structure of these is a layered composition of sheets of tilted
YbCh${}_{\text{6}}$ octahedra alternating with filler planes comprised of
alkaline/transition/boron group metal ions.

Figure 1: Illustration of the crystal structure of NaYbS${}_{\text{2}}$. The
blue polygons represent the distorted YbS${}_{\text{6}}$ octahedra with Yb3+
in the center of their basal planes, the yellow polygons the distorted NaS6
octahedra with the Na+ ions. The S2- ions are printed in yellow. (a) Side view
perpendicular to the c direction of the unit cell (thin black lines)
containing three layers of YbS${}_{\text{2}}$ planes. (b) View from top
parallel to the c direction onto the middle layer. The light blue triangles
perpendicular to c are equilateral (edge length $\lambda$) and form triangular
lattice planes.

Fig. 1 illustrates the crystal structure of NaYbS${}_{\text{2}}$ as an
example. The sulfur octahedra (blue) are tilted such that the Yb3+ 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 $\lambda$ and six small isosceles triangles with two edges of
length $\sigma$ and one edge of length $\lambda$. The large triangles are
those oriented perpendicular to the c direction, forming a perfect triangular
lattice. For an ideal octahedron we would have $\sigma=\lambda$, the tilting
angle of the octahedral axis with respect to the triangular plane would be
$\alpha=\cos^{-1}\left(1/\sqrt{3}\right)\approx 54.74^{\circ}$. In contrast
all delafossites have $\sigma\neq\lambda$, the tilting angle of the octahedral
axis then is given by $\alpha=\cos^{-1}\left(\lambda/(\sqrt{3}\sigma)\right)$,
see Table 1 for numbers.

The delafossite structure can have two polytypes according to the orientation
of the planar layer stacking. The space group of the rhombohedral 3R type
delafossites is $R\bar{3}m$ wheras the hexagonal $2H$ types have a space group
of $P6_{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
$R\bar{3}m$ space group (Table 1). Assigning typical oxidation states, we have
Yb3+ 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).

###  II.2 Yb${}^{\text{3+}}$ in a trigonal crystal field

Figure 2: Schematic energy levels of the Yb3+ ion (not true to scale). From
left to right: free ion, with spin-orbit coupling, in a cubic crystal field,
in a trigonal crystal field. The thicknesses of the horizontal lines are
proportional to the degeneracies of the respective energy levels, see text.

Fig. 2 schematically shows the energy levels of a single Yb3+ ion. The
fourteen 4f13 states of Yb3+ with $\ell=3$, $s=1/2$ are split by the spin-
orbit coupling into a $j=\ell+s$ octet and a $j=\ell-s$ sextet. In a perfect
octahedral (cubic) environment with ideal tilting angle $\alpha$, the $j=7/2$
states are further split into two doublets $\Gamma_{6}$ and $\Gamma_{7}$ and a
$\Gamma_{8}$ quartet, the $j=5/2$ states into a $\Gamma_{7}$ doublet and a
$\Gamma_{8}$ quartet. Distorting the octahedron along one of its trigonal axes
lowers the CEF environment to trigonal, and the local site symmetry of the
Yb3+ ions is $C_{3\text{v}}$ with the threefold axis parallel to the c
direction. This transforms the formerly cubic states like
$\Gamma_{6}\to\Gamma_{6}^{\text{T}}$, $\Gamma_{7}\to\Gamma_{6}^{\text{T}}$,
and
$\Gamma_{8}\to\Gamma_{4}^{\text{T}}+\Gamma_{5}^{\text{T}}+\Gamma_{6}^{\text{T}}$.
We note that although the two representations $\Gamma_{4}^{\text{T}}$ and
$\Gamma_{5}^{\text{T}}$ 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 Yb3+ ion at an arbitrary lattice site $i$ is then
given by

$\displaystyle{\cal H}_{\text{CEF}}(i)$ $\displaystyle=B_{2}^{0}O_{2}^{0}({\bf
J}_{i})+B_{4}^{0}O_{4}^{0}({\bf J}_{i})+B_{4}^{3}O_{4}^{3}({\bf J}_{i})$
$\displaystyle+B_{6}^{0}O_{6}^{0}({\bf J}_{i})+B_{6}^{3}O_{6}^{3}({\bf
J}_{i})+B_{6}^{6}O_{6}^{6}({\bf J}_{i})$ (1)

where $B_{n}^{m}$ are crystal-field parameters and $O_{n}^{m}({\bf 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 Yb3+ ion resides in an ideal octahedron. The
local symmetry of the Yb3+ ion then is cubic with $O_{\text{h}}$ symmetry, and
additional relationships between the $B_{n}^{m}$ 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]

$\displaystyle{\cal H}_{\text{CEF}}^{(\text{cubic},3)}$
$\displaystyle=B_{4}^{(3)}\left(O_{4}^{0}-20\sqrt{2}O_{4}^{3}\right)$
$\displaystyle\phantom{=}+B_{6}^{(3)}\left(O_{6}^{0}+\frac{35\sqrt{2}}{4}O_{6}^{3}+\frac{77}{8}O_{6}^{6}\right).$
(2)

An explicit expression for the matrix of this Hamiltonian for $j=7/2$ in the
$\left|j,m\right\rangle$ basis is given in Appendix B.

Only two independent crystal-field parameters remain. The corresponding
wavefunctions of ${\cal H}_{\text{CEF}}^{(\text{cubic},3)}$ which are grouped
in Kramers pairs consisting of time reversed states are given by

$\displaystyle\Gamma_{6}$
$\displaystyle:\quad\mp\frac{1}{3}\sqrt{\frac{35}{6}}\left|\frac{7}{2},\pm\frac{5}{2}\right\rangle-\frac{1}{3}\sqrt{\frac{7}{3}}\left|\frac{7}{2},\mp\frac{1}{2}\right\rangle\pm\frac{1}{3}\sqrt{\frac{5}{6}}\left|\frac{7}{2},\mp\frac{7}{2}\right\rangle,$
$\displaystyle\quad\left\langle\Gamma_{6}\right|J_{z}\left|\Gamma_{6}\right\rangle=\mathop{\rm
diag}\left(\pm\frac{7}{6}\right),$ $\displaystyle\Gamma_{7}$
$\displaystyle:\quad\mp\frac{1}{3}\sqrt{\frac{7}{2}}\left|\frac{7}{2},\pm\frac{7}{2}\right\rangle+\frac{\sqrt{5}}{3}\left|\frac{7}{2},\pm\frac{1}{2}\right\rangle\pm\frac{1}{3\sqrt{2}}\left|\frac{7}{2},\mp\frac{5}{2}\right\rangle,$
$\displaystyle\quad\left\langle\Gamma_{7}\right|J_{z}\left|\Gamma_{7}\right\rangle=\mathop{\rm
diag}\left(\pm\frac{3}{2}\right),$ $\displaystyle\Gamma_{8}$
$\displaystyle:\quad\left\\{\begin{aligned}
&\frac{1}{3}\sqrt{\frac{14}{3}}\left|\frac{7}{2},\pm\frac{7}{2}\right\rangle\pm\frac{1}{3}\sqrt{\frac{5}{3}}\left|\frac{7}{2},\pm\frac{1}{2}\right\rangle+\frac{2}{3}\sqrt{\frac{2}{3}}\left|\frac{7}{2},\mp\frac{5}{2}\right\rangle\\\
&\left|\frac{7}{2},\pm\frac{3}{2}\right\rangle\end{aligned}\right.,$
$\displaystyle\quad\left\langle\Gamma_{8}\right|J_{z}\left|\Gamma_{8}\right\rangle=\mathop{\rm
diag}\left(\pm\frac{7}{6},\pm\frac{3}{2}\right).$ (3)

They are independent of the crystal-field parameters and determined by
symmetry only. The notation $\mathop{\rm diag}(A)$ means that the matrix
representation
$\left\langle\Gamma_{\alpha}\left|J_{z}\right|\Gamma_{\alpha}\right\rangle$ of
the $J_{z}$ operator has eigenvalues $A$ with wavefunctions as noted.

###  II.3 Ground state of the Yb${}^{\text{3+}}$ ion

Lowering the local symmetry of the crystal field from cubic to trigonal by
compressing the threefold axis parallel to the c direction splits the
$\Gamma_{8}$ quartet obtained so far into two Kramers doublets and all six
crystal field parameters in Eq. (1) become independent. Also the CEF
wavefunctions for the four doublets now depend on the crystal-field potential,
see Appendix B for a $j=7/2$ matrix representation of the corresponding
Hamiltonian. This Hamiltonian couples only $|j,m\rangle$ states with $\Delta
m=\pm 3$, such that three of the four resulting Kramers doublets may be
written as [39, 40]

$\displaystyle\left|\psi^{+}\right\rangle$ $\displaystyle=-\alpha{\rm e}^{{\rm
i}\phi_{\alpha}}\left|\frac{7}{2},\frac{7}{2}\right\rangle+\beta\left|\frac{7}{2},\frac{1}{2}\right\rangle+\gamma{\rm
e}^{-{\rm i}\phi_{\gamma}}\left|\frac{7}{2},-\frac{5}{2}\right\rangle,$
$\displaystyle\left|\psi^{-}\right\rangle$ $\displaystyle=\alpha{\rm e}^{-{\rm
i}\phi_{\alpha}}\left|\frac{7}{2},-\frac{7}{2}\right\rangle+\beta\left|\frac{7}{2},-\frac{1}{2}\right\rangle-\gamma{\rm
e}^{{\rm i}\phi_{\gamma}}\left|\frac{7}{2},\frac{5}{2}\right\rangle$ (4)

where $\alpha$, $\beta$, and $\gamma$ are real with
$\alpha^{2}+\beta^{2}+\gamma^{2}=1$. We note that time reversal $T$ doesn’t
change the sign of
$\left|\frac{7}{2},-\frac{1}{2}\right\rangle=T\left|\frac{7}{2},\frac{1}{2}\right\rangle$.
This is a general feature of the time-reversal operator, giving
$T\left|j,\pm\frac{1}{2}\right\rangle=(-)^{j\pm
1/2}\left|j,\mp\frac{1}{2}\right\rangle$ with $T^{2}=-1$ for half-integer
total momentum $j$.

The fourth doublet is the pure state
$\left|\frac{7}{2},\pm\frac{3}{2}\right\rangle$, adiabatically (as function of
the trigonal distortion) connected to the pure state in the $\Gamma_{8}$
quartet.

One doublet defined in Eqs. (4) will be the ground state. In principle
$\left|\frac{7}{2},\pm\frac{3}{2}\right\rangle$ 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 ${\bf J}$ within the ground-state doublet can be read off:

$\displaystyle\left\langle\psi^{\pm}\left|J_{z}\right|\psi^{\pm}\right\rangle$
$\displaystyle=\pm\frac{1}{2}\left(7\alpha^{2}+\beta^{2}-5\gamma^{2}\right),$
$\displaystyle\left\langle\psi^{-}\left|J_{x}\right|\psi^{+}\right\rangle$
$\displaystyle=\sqrt{7}\alpha\gamma{\rm e}^{{\rm
i}(\phi_{\alpha}-\phi_{\gamma})}+2\beta^{2}=\left\langle\psi^{+}\left|J_{x}\right|\psi^{-}\right\rangle^{*},$
$\displaystyle\left\langle\psi^{-}\left|J_{y}\right|\psi^{+}\right\rangle$
$\displaystyle={\rm i}\left(\sqrt{7}\alpha\gamma{\rm e}^{{\rm
i}(\phi_{\alpha}-\phi_{\gamma})}+2\beta^{2}\right)$
$\displaystyle=\left\langle\psi^{+}\left|J_{y}\right|\psi^{-}\right\rangle^{*}.$

All other matrix elements vanish. The total momentum operator ${\bf J}$
transforms like ${\bf J}=-T{\bf J}T^{-1}$ under time reversal $T$. With
$T=UK$, $U=\exp({\rm i}\pi 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 $\phi_{\alpha}=\phi_{\gamma}+2\pi n$, $n\in\mathbb{Z}$. This
indeed allows us to introduce a pseudospin ${\bf S}$ for the ground-state
doublet by mapping

$\displaystyle g_{j}J_{z}$ $\displaystyle\to g_{\parallel}S_{z},$
$\displaystyle g_{j}J_{x}$ $\displaystyle\to g_{\perp}S_{x},$ $\displaystyle
g_{j}J_{y}$ $\displaystyle\to g_{\perp}S_{y},$ $\displaystyle g_{\parallel}$
$\displaystyle=g_{j}\left(7\alpha^{2}+\beta^{2}-5\gamma^{2}\right),$
$\displaystyle g_{\perp}$
$\displaystyle=g_{j}\left(2\sqrt{7}\alpha\gamma+4\beta^{2}\right),$
$\displaystyle g_{j}$
$\displaystyle=1+\frac{j(j+1)+s(s+1)-\ell(\ell+1)}{2j(j+1)}$
$\displaystyle=\frac{8}{7}\approx 1.14\quad\mbox{(Land\'{e} factor)}.$

The Zeeman splitting for a single Yb3+ ion at site $i$ is then obtained from
the Hamiltonian

$\displaystyle{\cal
H}_{\text{Zeeman}}(i)=-\mu_{0}g_{j}\mu_{\text{B}}\sum_{\alpha}J_{i}^{\alpha}H_{\alpha}$
$\displaystyle\to-\mu_{0}\mu_{\text{B}}\left[g_{\parallel}S_{i}^{z}H_{z}+g_{\perp}\left(S_{i}^{x}H_{x}+S_{i}^{y}H_{y}\right)\right],$
(5)

where $\mu_{0}$ is the magnetic permeability constant and $\mu_{\text{B}}$ the
Bohr magneton.

The above equations give an argument why the pure state
$\left|\frac{7}{2},\pm\frac{3}{2}\right\rangle$ cannot be the ground state
here: We would have transverse matrix elements
$\left\langle\frac{7}{2},\pm\frac{3}{2}\left|J_{x,y}\right|\frac{7}{2},\mp\frac{3}{2}\right\rangle=0$,
i. e. $g_{\perp}\equiv 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
$\left|\frac{7}{2},\pm\frac{3}{2}\right\rangle$ doublet as ground state.

###  II.4 Electron paramagnetic resonance

We have done thorough electron paramagnetic resonance studies on all three
NaYbCh${}_{\text{2}}$ delafossites [15, 16, 18]. Table 2 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 NaYbS2 in Fig.3. In a certain temperature range, the
resonance intensity $I_{\text{EPR}}$ is antiproportional to the temperature
[29]. Regarding $I_{\text{EPR}}$ as a measure for the resonant susceptibility
$\chi_{\text{R}}$, we can correspondingly assign characteristic temperatures
$\Theta_{\parallel,\perp}$, see Table 2. We note that these temperatures are
not necessarily identical with the Curie-Weiss temperatures obtained from
susceptibility measurements (see below).

Figure 3:  EPR spectra (upper frames) and anisotropy of the EPR $g$ factor
(lower frame) of single crystalline NaYbS2 at $T=19\,\text{K}$ and the
microwave field $b_{\text{mw}}\bot\text{c}$ axis ($\nu=9.4\,\text{GHz}$). The
spectra taken at external fields $B\bot\text{c}$ axis and $B\|\text{c}$ axis
were fitted by a Lorentzian lineshape (dashed lines). The anisotropy of the
EPR $g$ factor can be described by
$g(\Theta)=\sqrt{g_{\|}^{2}\cos^{2}\Theta+g_{\bot}^{2}\sin^{2}\Theta}$ with
$g_{\bot}=3.19(5)$ and $g_{\|}=0.57(3)$. The sample was rotated around an axis
lying in the basal plane parallel to $b_{\text{mw}}$.

Towards high temperatures the temperature dependence of the resonance
linewidth $\Delta 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 $\Delta 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
$\Delta E$ and yields for $\Delta E\gg k_{\text{B}}T$ an approximate
temperature dependence $\propto\exp(-\Delta E/k_{\text{B}}T)$. Identifying
$\Delta E$ with the energy $\Delta E_{12}$ of the first excited crystal-field
split state of the Yb3+ ion the data analysis of $\Delta B(T)$ gives a rough
estimate (within $\pm 2.5\,\rm meV$) of $\Delta E_{12}$, values of which are
also shown in Table 2.

###  II.5 Inelastic neutron scattering

Table 2: Characteristic values for the NaYbCh${}_{\text{2}}$ delafossites,
determined by different experimental methods discussed in the text. For
NaYbO${}_{\text{2}}$, only results on powder samples are available which are
reproduced in the respective $\parallel$ row of each quantity. An asterisk in
one of the Ref. columns denotes this publication.
$\begin{array}[]{cc|c|c|c|c|c|c|l}\text{observable}&\text{quantity}&\text{NaYbO}_{\text{2}}&\text{Ref.}&\text{NaYbS}_{\text{2}}&\text{Ref.}&\text{NaYbSe}_{\text{2}}&\text{Ref.}&\text{remarks}\\\
\hline\cr\hline\cr\text{EPR}&\left|g_{\parallel}\right|&1.75(3)&\text{\cite[cite]{[\@@bibref{Number}{ranjith:19}{}{}]}}&0.57(3)&\text{\cite[cite]{[\@@bibref{Number}{sichelschmidt:19}{}{}]}}&1.01(1)&\text{\cite[cite]{[\@@bibref{Number}{sichelschmidt:20}{}{}]}}&T=20\,\text{K}\\\
&\left|g_{\perp}\right|&3.28(8)&&3.19(5)&&3.13(4)&&\\\
&\Theta_{\parallel}&-9\,\text{K}&&-15.2\,\text{K}&&-14.3\,\text{K}&&\\\
&\Theta_{\perp}&&&-14.8\,\text{K}&&-14.0\,\text{K}&&\\\ &\Delta
E_{12}&27\,\text{meV}&&17\,\text{meV}&&14\,\text{meV}&&\\\
\hline\cr\hline\cr\text{INS}&\Delta
E_{12}&34.8\,\text{meV}&\text{\cite[cite]{[\@@bibref{Number}{ding:19}{}{}]}}&12\,\text{meV}&\text{\cite[cite]{[\@@bibref{Number}{baenitz:18}{}{}]}}&15.8\,\text{meV}&\text{\cite[cite]{[\@@bibref{Number}{zhang:21}{}{}]}}&T=5\,\text{K}\\\
&\Delta E_{13}&58.5\,\text{meV}&&23\,\text{meV}&&24.3\,\text{meV}&&\\\ &\Delta
E_{14}&83.1\,\text{meV}&&39\,\text{meV}&&30.5\,\text{meV}&&\\\
\hline\cr\hline\cr\chi(T)&\Theta_{\parallel}&-6\,\text{K}&\text{\cite[cite]{[\@@bibref{Number}{ranjith:19}{}{}]}}&-1.8\,\text{K}&*&-3.5\,\text{K}&\text{\cite[cite]{[\@@bibref{Number}{ranjith:19a}{}{}]}}&\\\
&\Theta_{\perp}&&&-11.2\,\text{K}&&-7\,\text{K}&&\\\
&\mu_{\parallel}&2.6\,\mu_{\text{B}}&&1.2\,\mu_{\text{B}}&&1.1\,\mu_{\text{B}}&&\\\
&\mu_{\perp}&&&2.87\,\mu_{\text{B}}&&2.43\,\mu_{\text{B}}&&\\\
\hline\cr&\hat{\Theta}_{\parallel,\perp}&-100\,\text{K}&*&-66\,\text{K}&*&-66\,\text{K}&*&\\\
&\mu_{\text{eff}}&4.6\,\mu_{\text{B}}&&4.6\,\mu_{\text{B}}&&4.6\,\mu_{\text{B}}&&\\\
\hline\cr\hline\cr
M(H)&\mu_{0}H_{\text{sat}}^{\parallel}&12\,\text{T}&\text{\cite[cite]{[\@@bibref{Number}{ranjith:19}{}{}]}}&\--&&25\,\text{T}&\text{\cite[cite]{[\@@bibref{Number}{ranjith:19a}{}{}]}}&T=470\,\text{mK}\\\
&\mu_{0}H_{\text{sat}}^{\perp}&&&14.7\,\text{T}&\text{\cite[cite]{[\@@bibref{Number}{luther:19}{}{}]}}&12\,\text{T}&&\\\
&M_{\text{sat}}^{\parallel}&1.36\,\mu_{\text{B}}/\text{Yb}^{3+}&&0.3\,\mu_{\text{B}}/\text{Yb}^{3+}&\text{\cite[cite]{[\@@bibref{Number}{baenitz:18}{}{}]}}&0.49\,\mu_{\text{B}}/\text{Yb}^{3+}&&\\\
&M_{\text{sat}}^{\perp}&&&1.6\,\mu_{\text{B}}/\text{Yb}^{3+}&\text{\cite[cite]{[\@@bibref{Number}{luther:19}{}{}]}}&1.5\,\mu_{\text{B}}/\text{Yb}^{3+}&&\\\
\hline\cr\hline\cr\end{array}$

Table 2 contains results on inelastic neutron scattering as well. For NaYbO2
and NaYbSe2, three clear maxima in energy scans at temperature $T=5\,\text{K}$
have been observed [43, 44], associated with the energy differences $\Delta
E_{1j}$ between the ground state $\left|\psi_{1}^{\pm}\right\rangle$ and the
excited states $\left|\psi_{j}^{\pm}\right\rangle$, $j=2,3,4$. Values are
quoted in the table.

Figure 4: (a) Excitation spectrum of a NaYbS${}_{\text{2}}$ powder sample
measured at $T=5\,\rm K$ with an incident neutron energy $E_{\text{i}}=50\,\rm
meV$. (b) The spectrum integrated over a momentum transfer range of $2\dots
4\,\text{\AA}^{-1}$ at different temperatures. Dashed line: fit of the $5\,\rm
K$ data. Figure taken from Ref. [31].

We investigated the excited CEF doublets of NaYbS${}_{\text{2}}$ in a time-of-
flight neutron scattering experiment [31], see Fig. 4 for a summary of the
results. At $T=5\,\rm K$ we have found only two clear maxima in the energy-
dependent scattering intensity, pointing to excited doublets located $23\,\rm
meV$ and $39\,\rm meV$ above the ground state. A third maximum is missing.
However at higher temperatures $T\geq 50\,\rm K$ additional intensity appears
at about $11\,\rm meV$ and $27\,\rm 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
$\Delta E_{12}\approx 12\,\rm 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${}_{\text{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 minimum exchange model discussed
in the following.

##  III Many Ytterbium ions

###  III.1 High-temperature magnetic susceptibility

The dimensionless uniform magnetic susceptibility $\chi(T)$ of a crystal with
volume $V$ is given by the change of the magnetisation ${\bf M}$ with the
magnetic field ${\bf B}=\mu_{0}{\bf H}$ with components

$\chi_{\alpha}=\mu_{0}\frac{\partial M_{\alpha}}{\partial
B_{\alpha}}=\frac{\mu_{0}}{V}\frac{\partial\bar{\mu}_{\alpha}}{\partial
B_{\alpha}}=-\frac{\mu_{0}}{V}\frac{\partial^{2}F}{\partial
B_{\alpha}^{2}},\quad\alpha=\parallel,\perp$ (6)

where $\bar{\mu}_{\alpha}=-\partial F/\partial B_{\alpha}$ is the total
magnetic moment in spatial direction $\alpha$ either parallel or perpendicular
to the c axis and $F=(1/\beta)\log\cal Z$ the canonical free energy,
$1/\beta=k_{\text{B}}T$ the inverse temperature and $k_{\text{B}}$ the
Boltzmann constant. For the molar susceptibility Eq. (6) has to be multiplied
by $N_{\text{L}}/(\nu/V)$ where $N_{\text{L}}$ is Avogadro’s number and
$\nu/V$ the volume density of Yb3+ ions.

We link the free energy in the usual way with statistical mechanics through
the partition function ${\cal Z}=\exp(-\beta{\cal H})$ where

${\cal H}=\sum_{i=1}^{\nu}\left[{\cal H}_{\text{CEF}}(i)+{\cal
H}_{\text{Zeeman}}(i)+{\cal H}_{\text{exc}}(i)\right]$ (7)

with ${\cal H}_{\text{CEF}}$ given by Eq. (1), ${\cal H}_{\text{Zeeman}}$
given by Eq. (5), and the exchange Hamiltonian for an arbitrary but fixed site
$i$

${\cal H}_{\text{exc}}(i)=\frac{1}{2}\sum_{\langle
ij\rangle}\sum_{\alpha\beta}J_{i}^{\alpha}\hat{\cal
J}_{ij}^{\alpha\beta}J_{j}^{\beta}$ (8)

where the sum is taken over the $z=6$ bonds connecting sites $j$ and site $i$
with an exchange tensor $\hat{\cal J}$ having the respective components
$\hat{\cal J}_{ij}^{\alpha\beta}$ 222We use the hat-notation to indicate that
the corresponding symbol relates to the full angular momentum $J$ with
$j=7/2$.. (See below in Sec. III.2 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 $\beta\to 0$ we can expand $\chi_{\alpha}$ (Eq.
(6)) in powers of $\beta$. We obtain a Curie-Weiss law

$\chi_{\alpha}=\frac{\nu}{V}\mu_{0}g_{j}^{2}\mu_{\text{B}}^{2}\frac{j(j+1)}{3}\beta\left(1+\beta
k_{\text{B}}\hat{\Theta}_{\alpha}\right)+{\cal O}(\beta^{3}),$ (9)

with $j=7/2$ and the Curie-Weiss temperatures given by

$\displaystyle k_{\text{B}}\hat{\Theta}_{\parallel}$
$\displaystyle=-\frac{4}{5}\left(j-\frac{1}{2}\right)\left(j+\frac{3}{2}\right)B_{2}^{0}-\frac{j(j+1)}{3}z\hat{J}_{\parallel},$
$\displaystyle k_{\text{B}}\hat{\Theta}_{\perp}$
$\displaystyle=+\frac{2}{5}\left(j-\frac{1}{2}\right)\left(j+\frac{3}{2}\right)B_{2}^{0}-\frac{j(j+1)}{3}z\hat{J}_{\perp}.$
(10)

Here $z\hat{J}_{\parallel}$ and $z\hat{J}_{\perp}$ are the contributions of
the exchange tensor parallel and perpendicular to the c direction. Other
components of the exchange tensor do not appear in $\chi_{\alpha}$ up to order
$\beta^{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_{2}^{0}$
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: Temperature dependence of the magnetic susceptibility of NaYbSe2
obtained in a field of $1\,\rm T$ for $\chi_{\perp}(T)$ and of $3\,\rm T$ for
$\chi_{\parallel}(T)$.

Figure 5 shows the temperature dependence of the susceptibilities
$\chi_{\parallel,\perp}(T)$ for NaYbSe${}_{\text{2}}$ in the full temperature
range $T=0.5\ldots 400\,\text{K}$ accessible to us where we have applied the
magnetic field in directions parallel to the c axis (label $H\parallel c$) and
perpendicular to it. At temperatures $T\gtrsim 120\,\text{K}$, the inverse
$\chi_{\parallel,\perp}^{-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 $\chi_{\perp}(T)$ according to Eq. (9) for $T\geq
150\,\text{K}$. From this fit we obtain $\mu_{\text{eff}}\approx
4.6\,\mu_{\text{B}}$, reflecting the full effective moment
$\mu_{\text{eff}}=g_{j}\mu_{\text{B}}\sqrt{j(j+1)}\approx
4.54\,\mu_{\text{B}}$. We also obtain a Curie-Weiss temperature
$\hat{\Theta}_{\perp}\approx-66\,\text{K}$. In the temperature range where the
Curie-Weiss law is a good approximation, $\chi_{\alpha}(T)$ is essentially
isotropic with $\hat{\Theta}_{\parallel}\approx\hat{\Theta}_{\perp}$,
independent of the direction of $\bf H$. The same fits for NaYbO2 (powder) and
NaYbS2 yield the full moment and isotropic Curie-Weiss temperatures as well.
Results are noted in Table 2.

Below $T=80\,\rm K$ the $s=1/2$ pseudospin state emerges. After subtracting a
temperature independent Van-Vleck contribution (obtained from high field
magnetization measurements) the susceptibility $\chi(T)$ for $10\,{\rm K}\leq
T\leq 40\,\rm K$ can be fitted with a Curie-Weiss law (right-hand plot of Fig.
5) which yields a Curie-Weiss temperature $\Theta_{\perp}=-7\,\rm K$, and an
effective moment $\mu_{\perp}=2.43\mu_{\text{B}}$ and a Curie-Weiss
temperature $\Theta_{\parallel}=-3.5\,\rm K$ and a moment of
$\mu_{\parallel}=1.1\mu_{\text{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.

Figure 6: Temperature dependence of the magnetic specific heat of the NaYbCh2
delafossites [31, 15, 18].

The pronounced maximum in the susceptibility (Fig. 5) corresponds to the
maximum found in the temperature dependence of the magnetic specific heat
$c_{\text{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 NaYbCh2 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_{\text{m}}(T)=\gamma T$ down to
the lowest accessible temperatures with a large residual value $\gamma\approx
1\,{\rm J}/\rm molK^{2}$ is found. This linear temperature dependence is well
known for heavy-fermion systems, for example YbRh2Si2 [49] or Yb4As3 [50, 51],
yet typical for a gapless spin liquid with fermionic excitations [36].

###  III.2 Pseudospin exchange model

We restrict ourselves now to the CEF ground-state doublet of each Yb3+ 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
$\langle ij\rangle$. We choose our coordinate system such that one bond is
parallel to the $x$ direction and the $z$ axis is perpendicular to the
triangular-lattice planes. For this bond, we can write

${\cal J}_{ij}=\begin{pmatrix}J_{\perp}&0&0\\\ 0&J_{\perp}&0\\\
0&0&J_{\parallel}\end{pmatrix}+\begin{pmatrix}J_{\Delta}&0&0\\\
0&-J_{\Delta}&J_{yz}\\\ 0&J_{yz}&0\end{pmatrix},$ (11)

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 for $\tilde{\cal J}_{ij}=U^{-1}{\cal J}_{ij}U$, $U$
unitary, $\tilde{\cal J}_{ij}$ diagonal. This is a direct consequence of the
trigonal distortion (tilting angle
$\alpha\neq\cos^{-1}\left(1/\sqrt{3}\right)$ defined in Sec. II.1) of the
YbCh6 octahedra.

The full pseudospin Hamiltonian with this parametrization, expressed with
ladder operators instead of cartesian spin operators then reads

$\displaystyle{\cal H}$ $\displaystyle=\sum_{\left\langle
ij\right\rangle}\left\\{\frac{1}{2}J_{\perp}\left(S_{i}^{+}S_{j}^{-}+S_{i}^{-}S_{j}^{+}\right)+J_{\parallel}S_{i}^{z}S_{j}^{z}+\frac{1}{2}J_{\Delta}\left({\rm
e}^{{\rm i}\phi_{ij}}S_{i}^{-}S_{j}^{-}+{\rm e}^{-{\rm
i}\phi_{ij}}S_{i}^{+}S_{j}^{+}\right)\right.$
$\displaystyle\phantom{=\sum_{\left\langle
ij\right\rangle}}\left.+\frac{1}{2\rm i}J_{yz}\left[{\rm e}^{{\rm
i}\phi_{ij}}\left(S_{i}^{z}S_{j}^{+}+S_{i}^{+}S_{j}^{z}\right)-{\rm e}^{-{\rm
i}\phi_{ij}}\left(S_{i}^{z}S_{j}^{-}+S_{i}^{-}S_{j}^{z}\right)\right]\right\\}+{\cal
H}_{\text{Zeeman}}$ (12)

with ${\cal H}_{\text{Zeeman}}$ given by Eqs. (5) and the direction-dependent
phases are

$\phi_{ij}=\left\\{\begin{array}[]{rl}0,&{\bf R}_{i}-{\bf R}_{j}=(\pm
1,0,0)\\\ \frac{2\pi}{3},&{\bf R}_{i}-{\bf
R}_{j}=\pm\left(-\frac{1}{2},\frac{\sqrt{3}}{2}\right)\\\
-\frac{2\pi}{3},&{\bf R}_{i}-{\bf
R}_{j}=\pm\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right)\end{array}\right..$
(13)

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 $\Delta s_{z}=1$ or $\Delta s_{z}=2$.

The three-dimensional ground-state phase diagram of this Hamiltonian has been
addressed by several authors [53, 54, 55, 56, 32, 33, 57] using slightly
differing parameterizations (see Appendix C for examples). Zhu et al. [33]
have shown that indeed spin-liquid type regions in the ground-state 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 1 at least nine do not show
any signature of long-range magnetic order down to the lowest measurement
temperature, typically either $0.4\,\rm K$ or $2\,\rm 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 YbCh6 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 Yb3+ ions, each being member of a different
YbCh6 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
Yb3+ 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
NaYbS2 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_{\text{N}}$. According to Ref. [58],
frustration ratios $f=\left|\Theta_{\text{CW}}\right|/T_{\text{N}}$ [59, 60]
of the order of $10$ can be achieved already for an isotropic Heisenberg
exchange model on the triangular lattice with $J_{\perp}=J_{\parallel}$ and
$J_{\Delta}=J_{yz}=0$ and a small interplane exchange coupling
$J_{\text{inter}}$. With typical Curie-Weiss temperatures of a few Kelvin, it
might be that $T_{\text{N}}$ falls into the few $100\,\text{mK}$ range,
however this requires extremely small $J_{\text{inter}}={\cal
O}(10^{-4}J_{\parallel,\perp})$.

Third, we cannot exclude off-stoichiometric Yb3+ 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 Yb3+ ions in the
ab plane (see Table 2), the impact of the long-ranged dipole-dipole
interaction might be a further reason for a suppression of $T_{\text{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.

###  III.3 Electron paramagnetic resonance

The exchange-narrowed linewidth of the EPR resonance in general is given by

$\mu_{0}\Delta H(\theta)=\frac{\rm
const.}{\mu_{\text{B}}g(\theta)}\,M_{2}\sqrt{\frac{M_{2}}{M_{4}}},$ (14)

where ${\rm const.}=\pi/\sqrt{3}$ for a cutoff Lorentzian lineshape, and ${\rm
const.}=\sqrt{2\pi}$ for a Lorentzian $\times$ Gaussian lineshape [65].
$\theta$ 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]

$\displaystyle M_{2}$ $\displaystyle=\frac{\left\langle\left[{\cal
H},S_{\text{total}}^{+}\right]\left[S_{\text{total}}^{-},{\cal
H}\right]\right\rangle}{\left\langle
S_{\text{total}}^{+}S_{\text{total}}^{-}\right\rangle},$ (15) $\displaystyle
M_{4}$ $\displaystyle=\frac{\left\langle\left[{\cal H},\left[{\cal
H},S_{\text{total}}^{+}\right]\right]\left[\left[S_{\text{total}}^{-},{\cal
H}\right],{\cal H}\right]\right\rangle}{\left\langle
S_{\text{total}}^{+}S_{\text{total}}^{-}\right\rangle},$ (16) $\displaystyle
S_{\text{total}}^{\pm}$ $\displaystyle=\sum_{i}S_{i}^{\pm}.$ (17)

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

$\displaystyle M_{4}^{\text{Gauss}}$
$\displaystyle=3\left(M_{2}^{\text{Gauss}}\right)^{2},$ (18)
$\displaystyle\mu_{0}\Delta H_{\text{Gauss}}(\theta)$
$\displaystyle=\frac{\text{const.}}{\mu_{\text{B}}g(\theta)}\sqrt{\frac{M_{2}^{\text{Gauss}}}{3}}.$
(19)

Using the Hamiltonian (12), we obtain

$\displaystyle\lim_{T\to\infty}M_{2}^{\parallel}$
$\displaystyle=\frac{3}{4}\left[2\left(J_{\parallel}-J_{\perp}\right)^{2}+2J_{\Delta}^{2}+5J_{yz}^{2}\right],$
(20) $\displaystyle\lim_{T\to\infty}M_{2}^{\perp}$
$\displaystyle=\frac{3}{8}\left[2\left(J_{\parallel}-J_{\perp}\right)^{2}+10J_{\Delta}^{2}+7J_{yz}^{2}\right],$
(21)

where the symbols $\parallel$, $\perp$ denote the direction of the applied
magnetic field $\mu_{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}^{\parallel}=2M_{2}^{\perp}$ for a rotationally
invariant exchange ($J_{\Delta}=J_{xy}=0$) which can be understood in the
following way: Thermal spin fluctuations are energetically favorable
perpendicular to the field $\mu_{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

$\displaystyle\lim_{T\to\infty}M_{4}^{\parallel}$
$\displaystyle=3\left(J_{\perp}-J_{\parallel}\right)^{2}\left(7J_{\parallel}^{2}-6J_{\parallel}J_{\perp}+11J_{\perp}^{2}\right)$
$\displaystyle\phantom{=}{}+\frac{3}{2}\left[22J_{\Delta}^{4}+71J_{\Delta}^{2}J_{yz}^{2}+76J_{yz}^{4}\right.$
$\displaystyle\phantom{=\frac{3}{2}}+J_{\Delta}^{2}\left(44J_{\perp}^{2}-56J_{\perp}J_{\parallel}+24J_{\parallel}^{2}\right)$
$\displaystyle\phantom{=\frac{3}{2}}+\left.J_{yz}^{2}\left(71J_{\perp}^{2}-104J_{\perp}J_{\parallel}+63J_{\parallel}^{2}\right)\right],$
(22) $\displaystyle\lim_{T\to\infty}M_{4}^{\perp}$
$\displaystyle=\frac{3}{2}\left(J_{\perp}-J_{\parallel}\right)^{2}\left(7J_{\parallel}^{2}-6J_{\parallel}J_{\perp}+11J_{\perp}^{2}\right)$
$\displaystyle\phantom{=}{}+\frac{3}{4}\left[182J_{\Delta}^{4}+193J_{\Delta}^{2}J_{yz}^{2}+92J_{yz}^{4}\right.$
$\displaystyle\phantom{=\frac{3}{4}}+J_{\Delta}^{2}\left(60J_{\perp}^{2}-56J_{\perp}J_{\parallel}+56J_{\parallel}^{2}\right)$
$\displaystyle\phantom{=\frac{3}{4}}+\left.J_{yz}^{2}\left(81J_{\perp}^{2}-104J_{\perp}J_{\parallel}+65J_{\parallel}^{2}\right)\right].$
(23)

Similar to $M_{2}$, also $M_{4}$ vanishes for a fully isotropic exchange
Hamiltonian, and $M_{4}^{\parallel}=2M_{4}^{\perp}$ for a rotationally
invariant exchange.

We assume $J_{\Delta},J_{yz}\ll J_{\perp},J_{\parallel}$. Taking into account
only finite exchange constants $J_{\perp}$ and $J_{\parallel}$, we get for the
high-temperature EPR linewidth, Eq. (14)

$\displaystyle\lim_{T\to\infty}\mu_{0}\Delta H_{\parallel}$
$\displaystyle=\frac{\rm
const.}{\mu_{\text{B}}g_{\parallel}}\frac{3}{2\sqrt{2}}\frac{\left(J_{\perp}-J_{\parallel}\right)^{2}}{\sqrt{7J_{\parallel}^{2}-6J_{\parallel}J_{\perp}+11J_{\perp}^{2}}},$
(24) $\displaystyle\lim_{T\to\infty}\mu_{0}\Delta H_{\perp}$
$\displaystyle=\frac{\rm
const.}{\mu_{\text{B}}g_{\perp}}\frac{3}{4\sqrt{2}}\frac{\left(J_{\perp}-J_{\parallel}\right)^{2}}{\sqrt{7J_{\parallel}^{2}-6J_{\parallel}J_{\perp}+11J_{\perp}^{2}}}$
(25)

for field parallel and perpendicular to the c axis. An estimate of
$\lim_{T\to\infty}\Delta H_{\parallel,\perp}$ 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\to\infty}\Delta H_{\parallel}/\Delta
H_{\perp}=2g_{\perp}/g_{\parallel}$. This relation is roughly consistent with
the experimental EPR data. Estimating $\lim_{T\to\infty}\Delta H$ with the
smallest value the linewidth reaches in its temperature dependence we obtain
for NaYbCh2, $\text{Ch}=\text{O}$, S, Se: $\lim_{T\to\infty}\Delta
H_{\parallel}/\Delta H_{\perp}=1.7/10/4.6$. From the EPR data in Table 2, we
obtain $2g_{\perp}/g_{\parallel}=3.7/11.1/6.2$.

###  III.4 Magnetization and susceptibility

In order to learn more about the size of the exchange constants, we have
measured the temperature-dependent uniform magnetic susceptibility $\chi(T)$
of the pseudospins at sufficiently low temperatures $T<30\,\rm K$ and the
magnetization $M(H)$ in applied magnetic fields $B=\mu_{0}H$ up to $B=30\,\rm
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 $\beta\to 0$ and obtain

$\chi_{\alpha}=\frac{\nu}{V}\mu_{0}g_{\alpha}^{2}\mu_{\text{B}}^{2}\frac{s(s+1)}{3}\beta\left(1+\beta
k_{\text{B}}\Theta_{\alpha}\right)+{\cal O}\left(\beta^{3}\right)$ (26)

where the Curie-Weiss temperatures $\Theta_{\parallel,\perp}$ are given by

$\displaystyle k_{\text{B}}\Theta_{\alpha}$
$\displaystyle=-\frac{s(s+1)}{3}\sum_{n=1}^{6}{\cal J}_{i,i+n}^{\alpha},$
$\displaystyle k_{\text{B}}\Theta_{\parallel}$
$\displaystyle=-\frac{3}{2}J_{\parallel},\quad
k_{\text{B}}\Theta_{\perp}=-\frac{3}{2}J_{\perp}$ (27)

with $s=1/2$ for a field applied parallel and perpendicular to the
crystallographic c direction. Here ${\cal J}_{i,i+n}^{\alpha}$ denotes the
component of the exchange energy between site $i$ and its $n$th neighbor along
the field direction $\alpha$. We note that $\Theta_{\parallel,\perp}$ do not
depend on the direction-dependent terms in the Hamiltonian (12) – the exchange
constants $J_{\Delta}$ and $J_{yz}$ only appear in higher orders of the
expansion (compare also Ref. [53]).

Table 3: Characteristic Hamiltonian parameters for the NaYbCh${}_{\text{2}}$
delafossites, calculated from the experimental values in Table 2 as discussed
in the text.
$\begin{array}[]{cc|c|c|c}\text{observable}&\text{parameter}&\text{NaYbO}_{\text{2}}&\text{NaYbS}_{\text{2}}&\text{NaYbSe}_{\text{2}}\\\
\hline\cr\hline\cr\chi(T)&g_{\parallel}&3.0&1.4&1.3\\\ &g_{\perp}&&3.3&2.8\\\
&J_{\parallel}&0.34\,\text{meV}&0.10\,\text{meV}&0.20\,\text{meV}\\\
&J_{\parallel}/k_{\text{B}}&4\,\text{K}&1.2\,\text{K}&2.3\,\text{K}\\\
&J_{\perp}&&0.64\,\text{meV}&0.40\,\text{meV}\\\
&J_{\perp}/k_{\text{B}}&&7.5\,\text{K}&4.7\,\text{K}\\\ \hline\cr\hline\cr
M(H)&g_{\parallel}&2.7&0.6&1.0\\\ &g_{\perp}&&3.2&3.0\\\
&J_{\parallel}&0.42\,\text{meV}&\--&0.25\,\text{meV}\\\
&J_{\parallel}/k_{\text{B}}&4.9\,\text{K}&\--&2.9\,\text{K}\\\
&J_{\perp}&&0.61\,\text{meV}&0.46\,\text{meV}\\\
&J_{\perp}/k_{\text{B}}&&7.0\,\text{K}&5.4\,\text{K}\\\
\hline\cr\hline\cr\end{array}$

Table 2 holds our findings from susceptibility measurements [15, 18]. It
contains the measured Curie-Weiss temperatures $\Theta_{\parallel,\perp}$ and
the effective moments $\mu_{\parallel,\perp}$ obtained from a Curie-Weiss fit
to $\chi(T)$ at $T<30\,\text{K}$. For NaYbO${}_{\text{2}}$ only powder samples
were available, the corresponding averaged values are listed in the respective
$x_{\parallel}$ 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
$\mu_{\parallel,\perp}=g_{\parallel,\perp}\sqrt{s(s+1)}$, we can determine the
effective g factors. Results are listed in Table 3.

Table 2 also holds the values for the effective moment $\mu_{\text{eff}}$ and
the Curie-Weiss temperatures $\hat{\Theta}_{\parallel,\perp}$ for
$T>150\,\text{K}$, in the high-temperature limit the latter are given by Eqs.
(10). From the derivation of the ground state pseudospin in Sec II.3 we have a
relationship

$J_{\alpha}=\left(\frac{g_{\alpha}}{g_{j}}\right)^{2}\hat{J}_{\alpha},\quad\alpha=\parallel,\perp$
(28)

with the exchange constants introduced in Eq. (8) for the full angular
momentum. Together with the EPR g factors given in Table 2, this allows us to
roughly estimate the $B_{2}^{0}$ CEF parameter of NaYbS${}_{\text{2}}$ and
NaYbSe${}_{\text{2}}$ to $B_{2}^{0}\approx-1\ldots{-0.5}\,\text{meV}$.

###  III.5 Saturation field

Figure 7: Illustration of the spin configuration near saturation. Solid arrows
represent the sublattice moments, the thick arrows the respective magnetic
field. Left: field applied parallel to the c direction. The sublattice moments
form a cone around the field direction. Middle: field applied perpendicular to
the c direction. For $J_{\parallel}<J_{\perp}$, the sublattice moments lie in
the ab plane indicated by the gray rectangle. Right: field applied
perpendicular to the c direction, and $J_{\parallel}>J_{\perp}$. The
sublattice moments lie in a plane containing the c axis and the field axis,
indicated by the gray rectangle.

Table 2 also contains the values for the saturation fields
$H_{\text{sat}}^{\parallel,\perp}$ where the field dependent magnetization
$M(H)$ per Yb3+ ion reaches its saturation values $M_{\text{sat}}$. For
NaYbO${}_{\text{2}}$ again only powder samples were available. The saturation
field $H_{\text{sat}}^{\parallel}$ of NaYbS${}_{\text{2}}$ parallel to the c
axis was too high to be reached in our experiments.

$H_{\text{sat}}$ is defined as an instability of the fully polarized state
towards $\Delta 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={\cal E}/(\nu s)$ is a function of three pairs of polar and
azimuthal angles of the sublattice moments.

####  III.5.1 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_{\perp}/J_{\parallel}\lesssim
3/2$), $\pi$-coplanar ($3/2\lesssim J_{\perp}/J_{\parallel}\lesssim 2.2$) and
umbrella phase ($J_{\perp}/J_{\parallel}\gtrsim 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, $\theta_{i}=\xi$. Near saturation with $\xi\ll 1$, the classical
ground-state energy density is given by

$\displaystyle e_{\parallel}$ $\displaystyle\approx
3sJ_{\parallel}-h_{\parallel}-\xi^{2}\left(\frac{3}{2}sJ_{\perp}+3sJ_{\parallel}-\frac{1}{2}h_{\parallel}\right),$
(29) $\displaystyle h_{\parallel}$
$\displaystyle=g_{\parallel}\mu_{\text{B}}\mu_{0}H.$

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 ground-state energy per
spin for the fully polarized state: We gain energy $h_{\parallel}$, but there
is an energy loss $3sJ_{\parallel}$ because of the perfect misalignment (i. e.
ferromagnetic alignment) of the moments on all three sublattices. The term
$\propto\xi^{2}$ shows what happens when the umbrella opens a bit: We gain
energy $(3/2)sJ_{\perp}\xi^{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_{\parallel}\xi^{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_{\parallel}\xi^{2}$. With the saturation magnetization
$M_{\text{sat}}^{\parallel,\perp}=sg_{\parallel,\perp}\mu_{\text{B}}/\text{Yb}^{3+}$
(see Table 2 for the experimental values), we obtain

$\mu_{0}H_{\text{sat}}^{\parallel}=\frac{3s^{2}(2J_{\parallel}+J_{\perp})}{M_{\text{sat}}^{\parallel}}.$
(30)

####  III.5.2 Field in the ab plane

For $J_{\parallel}\neq J_{\perp}$, the magnetic phase diagram has three
intermediate phases between $H=0$ and saturation [70, 71, 72]. Like in the
case ${\bf H}\parallel\text{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
$\theta_{i}$ and $\phi_{i}$. Near saturation, $\theta_{i}\to\pi/2$,
$\phi_{i}\to\alpha$ where $\alpha$ is the angle of the applied field relative
to the a axis, and we obtain

$\mu_{0}H_{\text{sat}}^{\perp}=\frac{9s^{2}J_{\perp}}{M_{\text{sat}}^{\perp}}$
(31)

for the saturation field in the ab plane. This is the case for
$J_{\parallel}<J_{\perp}$. The independence of $H_{\text{sat}}^{\perp}$ from
$J_{\parallel}$ suggests that a planar spin configuration with all spins in
the ab plane is energetically favorable at and at least infinitesimally below
$H_{\text{sat}}$. At finite $H<H_{\text{sat}}$, a distorted cone might form.
Writing $\delta_{i}=\pi/2-\theta_{i}$, $\epsilon_{i}=\alpha-\phi_{i}$, the six
angles are a solution to Eq. (38) with $\delta_{i},\epsilon_{i}\neq 0$.

A geometric interpretation can be obtained by looking at two possible planar
configurations near $H_{\text{sat}}$, sketched in the middle and right
illustration of Fig. 7: the sublattice moments form a fan either in the ab
plane setting $\delta_{i}=0$ or in a plane perpendicular to it containing the
c and the field axis with $\epsilon_{i}=0$. For the former, the gradient near
$H_{\text{sat}}$ is given by

$\left(\frac{\partial
e_{\perp}}{\partial\epsilon_{i}}\right)_{\delta_{i}=0}\approx\left(\frac{h_{\perp}}{3}-2sJ_{\perp}\right)\begin{pmatrix}\epsilon_{1}-u\left(\epsilon_{2}+\epsilon_{3}\right)\\\
\epsilon_{2}-u\left(\epsilon_{3}+\epsilon_{1}\right)\\\
\epsilon_{3}-u\left(\epsilon_{1}+\epsilon_{2}\right)\end{pmatrix}$ (32)

with $u=sJ_{\perp}/\left[2sJ_{\perp}-(1/3)h_{\perp}\right]$. Minimizing this
gives

$\epsilon_{1}=\epsilon_{2}=\frac{u}{1-u}\epsilon_{3},\quad u=-1\mbox{ or
}\frac{1}{2}.$ (33)

From $u=-1$ we obtain $h_{\perp}=9sJ_{\perp}$ (first factor in Eq. (44)), for
$u=1/2$ we obtain the unphysical solution $h_{\perp}=0$. The ground-state
energy density is $e_{\perp}=3sJ_{\perp}-h_{\perp}$. Similar to the case ${\bf
H}\parallel\text{c}$, we lose energy $3sJ_{\perp}$ due to the ferromagnetic
alignment of all spins, and we gain energy $h_{\perp}$. Deviating slightly
from full polarization, but still with all sublattice moments in the ab plane
we gain energy $\Delta e_{\perp}\propto 3sJ_{\perp}|\epsilon_{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_{\text{sat}}$ for a spin configuration in the plane
containing the c axis and the field ($\epsilon_{i}=0$) is given by

$\left(\frac{\partial
e_{\perp}}{\partial\delta_{i}}\right)_{\epsilon_{i}=0}\approx\left(\frac{h_{\perp}}{3}-2sJ_{\perp}\right)\begin{pmatrix}\delta_{1}-u\left(\delta_{2}+\delta_{3}\right)\\\
\delta_{2}-u\left(\delta_{3}+\delta_{1}\right)\\\
\delta_{3}-u\left(\delta_{1}+\delta_{2}\right)\end{pmatrix}$ (34)

with $u=sJ_{\parallel}/\left[2sJ_{\perp}-(1/3)h_{\perp}\right]$. Minimizing
this gives

$\delta_{1}=\delta_{2}=\frac{u}{1-u}\delta_{3},\quad u=-1\mbox{ or
}\frac{1}{2}.$ (35)

From $u=-1$ we obtain $h_{\perp}=3s\left(J_{\parallel}+2J_{\perp}\right)$
(second factor in Eq. (44)). This is the case for $J_{\parallel}>J_{\perp}$.
Deviating from full polarization now results in an energy gain $\Delta
e_{\perp}\propto s(J_{\parallel}+2J_{\perp})|\delta_{i}|$, which is, for
$J_{\parallel}<J_{\perp}$, smaller than the energy gain when opening in the ab
plane. For $u=1/2$ we obtain
$h_{\perp}=6s\left(J_{\perp}-J_{\parallel}\right)$. This solution is
unphysical because this would include an energy loss $\propto 2sJ_{\parallel}$
when deviating from full polarization.

In the last six rows of Table 3 we note the values obtained from the equations
above. For NaYbS2, $J_{\parallel}$ cannot be determined because of the missing
value for $H_{\text{sat}}^{\parallel}$. If we take the exchange constants
derived from $I_{\text{EPR}}(T)$ and $\chi(T)$, we can roughly estimate
$56\,\text{T}\lesssim\mu_{0}H_{\text{sat}}^{\parallel}\lesssim 112\,\text{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 NaYbCh2 compounds contain nearly
perfect magnetic Yb3+ triangular lattice planes. The determination of the
crystal field levels by neutron scattering has shown that at temperatures
$T\lesssim 100\,\rm 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 Yb3+ 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 Yb3+ triangular lattice system YbGaMgO4 [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 Yb3+ ions also plays a role like for example in the YbAlO3 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_{\parallel}$
and $J_{\perp}$ in the NaYbCh2 materials are comparatively large for an Yb3+
triangular lattice. In YbMgGaO4, NaBaYb(BO3)2 and Rb3Yb(PO4)2 the low-
temperature Curie-Weiss temperatures and thus the exchange constants are more
than one order of magnitude smaller compared to the NaYbCh2 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_{\text{m}}(T)$. Here a Schottky peak develops in an applied
magnetic field due to the Zeeman splitting of the Yb3+ 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 NaYbCh2 systems the saturation fields exceed $10\,\rm T$.

Furthermore, we have found a pronounced and large linear temperature
dependence of the magnetic specific heat $c_{\text{m}}(T)=\gamma T$ with a
residual value $\gamma\approx 1\,{\rm J}/\rm molK^{2}$ for the NaYbCh2
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 ($\mu$SR), nuclear magnetic resonance (NMR) and finally inelastic
neutron scattering. For the other Yb3+ triangular lattice materials mentioned
above, such residual contributions are absent or negligibly small. The
presence of the large $\gamma$ 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 Yb3+ 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 NaYbCh2 delafossites have
their own fascinating physics which differs significantly from the previously
known planar Yb3+ 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 NaYbS2 with a ratio of the coupling coefficients of
$J_{\perp}/J_{\parallel}=6.25$ and a ratio of the EPR g factors of
$g_{\perp}/g_{\parallel}=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.

###### Acknowledgements.

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).

##  Appendix A Stevens operators

The Stevens operators in Eq. (1) have the explicit form [82]

$\displaystyle O_{2}^{0}$ $\displaystyle=3J_{z}^{2}-j(j+1),$ $\displaystyle
O_{4}^{0}$
$\displaystyle=35J_{z}^{4}-30j(j+1)J_{z}^{2}+25J_{z}^{2}-6j(j+1)+3\left[j(j+1)\right]^{2},$
$\displaystyle O_{4}^{3}$
$\displaystyle=\frac{1}{4}\left[J_{z}\left(J_{+}^{3}+J_{-}^{3}\right)+\left(J_{+}^{3}+J_{-}^{3}\right)J_{z}\right],$
$\displaystyle O_{6}^{0}$
$\displaystyle=231J_{z}^{6}-315j(j+1)J_{z}^{4}+735J_{z}^{4}+105\left[j(j+1)\right]^{2}J_{z}^{2}-525j(j+1)J_{z}^{2}+294J_{z}^{2}$
$\displaystyle\phantom{=}-5\left[j(j+1)\right]^{3}+40\left[j(j+1)\right]^{2}-60j(j+1),$
$\displaystyle O_{6}^{3}$
$\displaystyle=\frac{1}{4}\left\\{\left[11J_{z}^{3}-3j(j+1)J_{z}-59J_{z}\right]\left(J_{+}^{3}+J_{-}^{3}\right)+\left(J_{+}^{3}+J_{-}^{3}\right)\left[11J_{z}^{3}-3j(j+1)J_{z}-59J_{z}\right]\right\\},$
$\displaystyle O_{6}^{6}$
$\displaystyle=\frac{1}{2}\left(J_{+}^{6}+J_{-}^{6}\right),$

and we use the standard definitions

$\displaystyle J^{2}\left|j,m\right\rangle$
$\displaystyle=j(j+1)\left|j,m\right\rangle,$ $\displaystyle
J_{z}\left|j,m\right\rangle$ $\displaystyle=m\left|j,m\right\rangle,$
$\displaystyle J_{+}\left|j,m\right\rangle$
$\displaystyle=\sqrt{j(j+1)-m(m+1)}\left|j,m+1\right\rangle,$ $\displaystyle
J_{-}\left|j,m\right\rangle$
$\displaystyle=\sqrt{j(j+1)-m(m-1)}\left|j,m-1\right\rangle,$ $\displaystyle
J_{x}$ $\displaystyle=\frac{1}{2}\left(J_{+}+J_{-}\right),$ $\displaystyle
J_{y}$ $\displaystyle=\frac{1}{2\rm i}\left(J_{+}-J_{-}\right).$

##  Appendix B Crystal-field Hamiltonian matrices

In a cubic environment with a threefold quantization axis, the Hamiltonian
matrix derived from Eq. (2) for $j=7/2$ is given by

$\displaystyle H_{\text{CEF}}^{(3)}$
$\displaystyle=\left(\begin{smallmatrix}420\left(B_{4}^{(3)}+3B_{6}^{(3)}\right)&0&0&30\sqrt{70}\left(21B_{6}^{(3)}-4B_{4}^{(3)}\right)\\\
0&-60\left(13B_{4}^{(3)}+105B_{6}^{(3)}\right)&0&0\\\
0&0&-180\left(B_{4}^{(3)}-63B_{6}^{(3)}\right)&0\\\
30\sqrt{70}\left(21B_{6}^{(3)}-4B_{4}^{(3)}\right)&0&0&180\left(3B_{4}^{(3)}-35B_{6}^{(3)}\right)\\\
0&-15\sqrt{10}\left(16B_{4}^{(3)}+147B_{6}^{(3)}\right)&0&0\\\ 0&0&0&0\\\
3465\sqrt{7}B_{6}^{(3)}&0&0&15\sqrt{10}\left(16B_{4}^{(3)}+147B_{6}^{(3)}\right)\\\
0&3465\sqrt{7}B_{6}^{(3)}&0&0\end{smallmatrix}\right\rfloor$
$\displaystyle\phantom{=}{}\left\lceil\begin{smallmatrix}0&0&3465\sqrt{7}B_{6}^{(3)}&0\\\
-15\sqrt{10}\left(16B_{4}^{(3)}+147B_{6}^{(3)}\right)&0&0&3465\sqrt{7}B_{6}^{(3)}\\\
0&0&0&0\\\ 0&0&15\sqrt{10}\left(16B_{4}^{(3)}+147B_{6}^{(3)}\right)&0\\\
180\left(3B_{4}^{(3)}-35B_{6}^{(3)}\right)&0&0&30\sqrt{70}\left(4B_{4}^{(3)}-21B_{6}^{(3)}\right)\\\
0&-180\left(B_{4}^{(3)}-63B_{6}^{(3)}\right)&0&0\\\
0&0&-60\left(13B_{4}^{(3)}+105B_{6}^{(3)}\right)&0\\\
30\sqrt{70}\left(4B_{4}^{(3)}-21B_{6}^{(3)}\right)&0&0&420\left(B_{4}^{(3)}+3B_{6}^{(3)}\right)\end{smallmatrix}\right).$

Its eigenvalues are

$\displaystyle E_{\Gamma_{6}}$
$\displaystyle=-315\left(4B_{4}^{(3)}+45B_{6}^{(3)}\right),$ $\displaystyle
E_{\Gamma_{7}}$
$\displaystyle=\phantom{-}405\left(4B_{4}^{(3)}-21B_{6}^{(3)}\right),$
$\displaystyle E_{\Gamma_{8}}$
$\displaystyle=-180\left(B_{4}^{(3)}-63B_{6}^{(3)}\right),$

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

$\displaystyle H_{\text{CEF}}$
$\displaystyle=\left(\begin{smallmatrix}21\left(B_{2}^{0}+20B_{4}^{0}+60B_{6}^{0}\right)&0&0&\sqrt{35}\left(B_{4}^{3}+12B_{6}^{3}\right)\\\
0&3\left(B_{2}^{0}-260B_{4}^{0}-2100B_{6}^{0}\right)&0&0\\\
0&0&-9\left(B_{2}^{0}-30B_{4}^{0}+1260B_{6}^{0}\right)&0\\\
6\sqrt{35}\left(B_{4}^{3}+12B_{6}^{3}\right)&0&0&-15\left(B_{2}^{0}+36B_{4}^{0}-420B_{6}^{0}\right)\\\
0&12\sqrt{5}\left(B_{4}^{3}-21B_{6}^{3}\right)&0&0\\\ 0&0&0&0\\\
360\sqrt{7}B_{6}^{6}&0&0&12\sqrt{5}\left(21\sqrt{5}B_{6}^{3}-B_{4}^{3}\right)\\\
0&360\sqrt{7}B_{6}^{6}&0&0\\\ \end{smallmatrix}\right\rfloor$
$\displaystyle\phantom{=}{}\left\lceil\begin{smallmatrix}0&0&360\sqrt{7}B_{6}^{6}&0\\\
12\sqrt{5}\left(B_{4}^{3}-21B_{6}^{3}\right)&0&0&360\sqrt{7}B_{6}^{6}\\\
0&0&0&0\\\ 0&0&12\sqrt{5}\left(21\sqrt{5}B_{6}^{3}-B_{4}^{3}\right)&0\\\
-15\left(B_{2}^{0}+36B_{4}^{0}-420B_{6}^{0}\right)&0&0&-6\sqrt{35}\left(B_{4}^{3}+12B_{6}^{3}\right)\\\
0&-9\left(B_{2}^{0}-30B_{4}^{0}+1260B_{6}^{0}\right)&0&0\\\
0&0&3\left(B_{2}^{0}-260B_{4}^{0}-2100B_{6}^{0}\right)&0\\\
-6\sqrt{35}\left(B_{4}^{3}+12B_{6}^{3}\right)&0&0&21\left(B_{2}^{0}+20B_{4}^{0}+60B_{6}^{0}\right)\\\
\end{smallmatrix}\right).$

Needless to say that closed-form expressions for the eigenvalues, apart from
the pure $\left|7/2,\pm 3/2\right\rangle$ 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 $\delta B_{n}^{m}$ like

$\displaystyle E_{\Gamma_{6}}$ $\displaystyle\to
E_{\Gamma_{6}}-\frac{70}{3}\delta\left[14B_{4}^{0}-\sqrt{2}B_{4}^{3}+10\left(24B_{6}^{0}+\sqrt{2}B_{6}^{3}+2B_{6}^{6}\right)\right],$
$\displaystyle E_{\Gamma_{7}}$ $\displaystyle\to
E_{\Gamma_{7}}+10\delta\left[42B_{4}^{0}-3\sqrt{2}B_{4}^{3}-14\left(24B_{6}^{0}+\sqrt{2}B_{6}^{3}+2B_{6}^{6}\right)\right],$
$\displaystyle E_{\Gamma_{8}}$ $\displaystyle\to
E_{\Gamma_{8}}\left\\{\begin{aligned}
&-9\delta\left[B_{2}^{0}+20\left(B_{4}^{0}-63B_{6}^{0}\right)\right]\\\
&+9\delta
B_{2}^{0}+\frac{20}{3}\delta\left(13B_{4}^{0}+\sqrt{2}B_{4}^{3}-357B_{6}^{0}+56\sqrt{2}B_{6}^{3}+112B_{6}^{6}\right)\end{aligned}\right.,$

in particular splitting the $\Gamma_{8}$ quartet into two Kramers doublets.

##  Appendix C Hamiltonian parametrization

Table 4: Common exchange parametrizations for the triangular-lattice
Hamiltonian
$\begin{array}[]{l|llll}\mbox{Authors}&J_{\perp}&J_{z}&J_{\Delta}&J_{xy}\\\
\hline\cr\hline\cr\mbox{Li et
al.~{}\cite[cite]{[\@@bibref{Number}{li:15}{}{}]}}&2J_{\pm}&J_{zz}&2J_{\pm\pm}&J_{z\pm}\\\
\mbox{Rau and
Gingras~{}\cite[cite]{[\@@bibref{Number}{rau:18}{}{}]}}&-2J_{\pm}&J_{zz}&2J_{\pm\pm}&2J_{z\pm}\\\
\mbox{Zhu et al.~{}\cite[cite]{[\@@bibref{Number}{zhu:18}{}{}]}}&J&\Delta
J&4J_{\pm\pm}&2J_{z\pm}\end{array}$

Our definition of the exchange constants occurring in the Hamiltonian (12) is
made such that it reduces to the standard XXZ model for $J_{\Delta}=0$ and
$J_{yz}=0$ which in turn reduces to the standard isotropic Heisenberg model
for $J_{\perp}=J_{z}$. Different authors use slightly different
parametrizations, shown in Table 4 for easy comparison with literature.

##  Appendix D Moment calculation

In the limit $T\to\infty$, the correlation functions between spins on
different sites factorize, only on-site spin-flip correlations retain a finite
value. The denominators of Eqs. (15) and (16) therefore are given by

$\displaystyle\left\langle
S_{\text{total}}^{+}S_{\text{total}}^{-}\right\rangle$
$\displaystyle=\sum_{ij}\left\langle
S_{i}^{+}S_{j}^{-}\right\rangle\xrightarrow{T\to\infty}\sum_{i}\left\langle
S_{i}^{+}S_{i}^{-}\right\rangle$ $\displaystyle=\nu\left\langle
S_{1}^{+}S_{1}^{-}\right\rangle=\nu\left(\frac{1}{2}+\left\langle
S_{i}^{z}\right\rangle\right)\to\frac{\nu}{2}.$

With a general short-range exchange Hamiltonian we have to evaluate an
expression like

$\displaystyle\left\langle\left[{\cal
H},S_{\text{total}}^{+}\right]\left[S_{\text{total}}^{-},{\cal
H}\right]\right\rangle$ $\displaystyle=\left\langle\left(\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\sum_{k}\left[S_{i}^{\alpha}S_{j}^{\beta},S_{k}^{+}\right]\right)\right.$
$\displaystyle\phantom{=}{}\left.\cdot\left(\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\sum_{k}\left[S_{k}^{-},S_{i}^{\alpha}S_{j}^{\beta}\right]\right)\right\rangle$
(36)

where $J_{\left\langle ij\right\rangle}^{\alpha\beta}$ is the exchange matrix
along (not necessarily nearest-neighbor) bond $\left\langle ij\right\rangle$
with spin indices $\alpha\beta$. The inner sum over sites $k$ can be removed
applying Jacobi’s identity

$\left[S_{i}^{\alpha}S_{j}^{\beta},S_{k}^{+}\right]=\left[S_{i}^{\alpha},\left[S_{j}^{\beta},S_{k}^{+}\right]\right]-\left[S_{j}^{\beta},\left[S_{i}^{\alpha},S_{k}^{+}\right]\right]$

and the standard spin commutation relations. Together with translational
invariance for a Bravais lattice with $\nu$ sites and $z$ neighbors the first
term in Eq. (36) simplifies to

$\displaystyle\left[{\cal H},S_{\text{total}}^{+}\right]$
$\displaystyle=\frac{\nu}{2}\sum_{n=2}^{z+1}\sum_{\alpha\beta}J_{n}^{\alpha\beta}\left(S_{1}^{\alpha}\left[S_{n}^{\beta},S_{n}^{+}\right]+\left[S_{1}^{\alpha},S_{1}^{+}\right]S_{n}^{\beta}\right),$

accordingly for its complex conjugate. The generalization to non-Bravais
lattices, not needed in the present context, should be obvious. For the
$R\bar{3}m$ triangular-lattice exchange Hamiltonian with $z=6$ nearest
neighbors given by Eq. (12), we eventually obtain for $M_{2}$, Eq. (15) in the
high-temperature limit the expressions (20) and (21).

Evaluating the commutators in

$\displaystyle\left[{\cal H},\left[{\cal
H},S_{\text{total}}^{+}\right]\right]$ $\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}$
$\displaystyle\phantom{=}{}\left(S_{m}^{\gamma}\left[S_{n}^{\delta},S_{i}^{\alpha}\left[S_{j}^{\beta},S_{j}^{+}\right]\right]\right.$
$\displaystyle\phantom{=}\left.{}+\left[S_{m}^{\gamma},S_{i}^{\alpha}\left[S_{j}^{\beta},S_{j}^{+}\right]\right]S_{n}^{\delta}\right.$
$\displaystyle\phantom{=}\left.{}+S_{m}^{\gamma}\left[S_{n}^{\delta},\left[S_{i}^{\alpha},S_{i}^{+}\right]S_{j}^{\beta}\right]\right.$
$\displaystyle\phantom{=}\left.{}+\left[S_{m}^{\gamma},\left[S_{i}^{\alpha},S_{i}^{+}\right]S_{j}^{\beta}\right]S_{n}^{\delta}\right)$

and its complex conjugate we also use Jacobi’s identity, the spin commutation
relations and translational invariance. With

$\displaystyle T_{1}$ $\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}S_{m}^{\gamma}\delta_{ni}\left[S_{i}^{\delta},S_{i}^{\alpha}\right]\left[S_{j}^{\beta},S_{j}^{+}\right]$
$\displaystyle=\frac{\nu}{4}\sum_{\delta
m=2}^{z+1}\sum_{\gamma\delta}J_{\delta m}^{\gamma\delta}\sum_{\delta
j=2}^{z+1}\sum_{\alpha\beta}J_{\delta j}^{\alpha\beta}S_{\delta
m}^{\gamma}\left[S_{1}^{\delta},S_{1}^{\alpha}\right]\left[S_{\delta
j}^{\beta},S_{\delta j}^{+}\right],$ $\displaystyle T_{2}$
$\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}S_{m}^{\gamma}\delta_{nj}S_{i}^{\alpha}\left[S_{j}^{\delta},\left[S_{j}^{\beta},S_{j}^{+}\right]\right]$
$\displaystyle=\frac{\nu}{4}\sum_{\delta
m=2}^{z+1}\sum_{\gamma\delta}J_{\delta m}^{\gamma\delta}\sum_{\delta
i=2}^{z+1}\sum_{\alpha\beta}J_{\delta i}^{\alpha\beta}S_{\delta
m}^{\gamma}S_{\delta
i}^{\alpha}\left[S_{1}^{\delta},\left[S_{1}^{\beta},S_{1}^{+}\right]\right],$
$\displaystyle T_{3}$ $\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\delta_{mi}\left[S_{i}^{\gamma},S_{i}^{\alpha}\right]\left[S_{j}^{\beta},S_{j}^{+}\right]S_{n}^{\delta}$
$\displaystyle=\frac{\nu}{4}\sum_{\delta
n=2}^{z+1}\sum_{\gamma\delta}J_{\delta n}^{\gamma\delta}\sum_{\delta
j=2}^{z+1}\sum_{\alpha\beta}J_{\delta
j}^{\alpha\beta}\left[S_{1}^{\gamma},S_{1}^{\alpha}\right]\left[S_{\delta
j}^{\beta},S_{\delta j}^{+}\right]S_{\delta n}^{\delta},$ $\displaystyle
T_{4}$ $\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\delta_{mj}S_{i}^{\alpha}\left[S_{j}^{\gamma},\left[S_{j}^{\beta},S_{j}^{+}\right]\right]S_{n}^{\delta}$
$\displaystyle=\frac{\nu}{4}\sum_{\delta
n=2}^{z+1}\sum_{\gamma\delta}J_{\delta n}^{\gamma\delta}\sum_{\delta
i=2}^{z+1}\sum_{\alpha\beta}J_{\delta i}^{\alpha\beta}S_{\delta
i}^{\alpha}\left[S_{1}^{\gamma},\left[S_{1}^{\beta},S_{1}^{+}\right]\right]S_{\delta
n}^{\delta},$ $\displaystyle T_{5}$ $\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}S_{m}^{\gamma}\delta_{ni}\left[S_{i}^{\delta},\left[S_{i}^{\alpha},S_{i}^{+}\right]\right]S_{j}^{\beta}$
$\displaystyle=\left(T_{2}\right)_{i\leftrightarrow
j}^{\alpha\leftrightarrow\beta},$ $\displaystyle T_{6}$
$\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}S_{m}^{\gamma}\delta_{nj}\left[S_{i}^{\alpha},S_{i}^{+}\right]\left[S_{j}^{\delta},S_{j}^{\beta}\right]$
$\displaystyle=\left(T_{1}\right)_{i\leftrightarrow
j}^{\alpha\leftrightarrow\beta},$ $\displaystyle T_{7}$
$\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\delta_{mi}\left[S_{i}^{\gamma},\left[S_{i}^{\alpha},S_{i}^{+}\right]\right]S_{j}^{\beta}S_{n}^{\delta}$
$\displaystyle=\left(T_{4}\right)_{i\leftrightarrow
j}^{\alpha\leftrightarrow\beta},$ $\displaystyle T_{8}$
$\displaystyle=\sum_{\left\langle
mn\right\rangle}\sum_{\gamma\delta}J_{\left\langle
mn\right\rangle}^{\gamma\delta}\sum_{\left\langle
ij\right\rangle}\sum_{\alpha\beta}J_{\left\langle
ij\right\rangle}^{\alpha\beta}\delta_{mj}\left[S_{i}^{\alpha},S_{i}^{+}\right]\left[S_{j}^{\gamma},S_{j}^{\beta}\right]S_{n}^{\delta}$
$\displaystyle=\left(T_{3}\right)_{i\leftrightarrow
j}^{\alpha\leftrightarrow\beta}$

the commutator may be written as

$\left[{\cal H},\left[{\cal
H},S_{\text{total}}^{+}\right]\right]=\sum_{\ell=1}^{8}T_{\ell}=2\sum_{\ell=1}^{4}T_{\ell},$

where the last equality holds for inversion-symmetric exchange only (which is
the case here). Then

$\left\langle\left[{\cal H},\left[{\cal
H},S_{\text{total}}^{+}\right]\right]\left[\left[S_{\text{total}}^{-},{\cal
H}\right],{\cal
H}\right]\right\rangle=4\sum_{\ell=1}^{4}\sum_{\ell^{\prime}=1}^{4}\left\langle
T_{\ell}T_{\ell^{\prime}}^{\dagger}\right\rangle$

has to be evaluated. The final results for the expectation value for $M_{4}$
in the limit $T\to\infty$ are reproduced in Eqs. (22) and (23).

##  Appendix E Saturation field

To calculate the saturation field of the magnetization, we regard the
pseudospins of the triangular lattice as classical vectors living on a three-
sublattice structure. On each sublattice $i$, any two spins are aligned
parallel relative to each other. The direction of each spin defines the
$z^{\prime}$ axis of a local coordinate system, characterized by polar and
azimuthal angels $\Theta_{i}$ and $\phi_{i}$. The full Hamiltonian with this
parameterization

$\left(\begin{array}[]{c}S_{i}^{x}\\\ S_{i}^{y}\\\
S_{i}^{z}\end{array}\right)=\left(\begin{array}[]{ccc}\cos\phi_{i}&-\sin\phi_{i}&0\\\
\sin\phi_{i}&\cos\phi_{i}&0\\\
0&0&1\end{array}\right)\left(\begin{array}[]{ccc}\cos\theta_{i}&0&\sin\theta_{i}\\\
0&1&0\\\
-\sin\theta_{i}&0&\cos\theta_{i}\end{array}\right)\left(\begin{array}[]{c}S_{i}^{x^{\prime}}\\\
S_{i}^{y^{\prime}}\\\ S_{i}^{z^{\prime}}\end{array}\right)$

and ${\bf S}^{\prime}=(0,0,s)^{T}$ is given below.

###  E.1 Field parallel to the crystallographic c axis

With a magnetic field in $z$ direction, the classical energy density per spin
is

$\displaystyle e_{\parallel}=\frac{{\cal E}_{z}}{\nu s}$
$\displaystyle=sJ_{\perp}\left[\sin\theta_{1}\sin\theta_{2}\cos\left(\phi_{1}-\phi_{2}\right)+\sin\theta_{2}\sin\theta_{3}\cos\left(\phi_{2}-\phi_{3}\right)+\sin\theta_{3}\sin\theta_{1}\cos\left(\phi_{3}-\phi_{1}\right)\right]$
$\displaystyle+sJ_{\parallel}\left[\cos\theta_{1}\cos\theta_{2}+\cos\theta_{2}\cos\theta_{3}+\cos\theta_{3}\cos\theta_{1}\right]-\frac{1}{3}h_{\parallel}\left[\cos\theta_{1}+\cos\theta_{2}+\cos\theta_{3}\right],$
(37) $\displaystyle h_{\parallel}$
$\displaystyle=g_{\parallel}\mu_{\text{B}}\mu_{0}H,$

independent of $J_{\Delta}$ and $J_{yz}$: We assume a three-sublattice
structure and positive gyromagnetic ratios, the direction-dependent terms
$\propto J_{\Delta},J_{yz}$ cancel exactly upon summation over the three
sublattice pairs $(1,2)$, $(2,3)$, and $(3,1)$. According to Eq. (12),
$J_{\Delta}$ and $J_{yz}$ give the energy gain when coupling states differing
by $\Delta S_{z}=1$ or $\Delta S_{z}=2$. For the classical model, the total
spin $S_{z}^{\text{tot}}=\sum_{i=1}^{\nu}S_{i}^{z}$ is conserved, therefore
$J_{\Delta}$ 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.

The saturation field is given by a singularity of the Hessian matrix

$\displaystyle m_{\text{H}}(e)$
$\displaystyle=\frac{\partial^{2}e}{\partial\\{\delta_{i},\phi_{i}\\}\partial\\{\delta_{j},\phi_{j}\\}},$
$\displaystyle\det m_{\text{H}}(e)$ $\displaystyle=0.$ (38)

By symmetry, for the classical ground state of Eq. (37), polar angles are all
$\theta_{i}=x$, azimuthal angles are $\alpha$ and $\alpha\pm 2\pi/3$ such that
$|\phi_{i}-\phi_{j}|=2\pi/3$. The nonzero block of the Hessian matrix
$m_{\text{H}}$ reduces to a scalar,
$m_{\text{H}}(e_{\parallel})=\partial^{2}e_{\parallel}/\partial x^{2}$, and
Eq. (38) reduces to $(3/2)sJ_{\perp}+3sJ_{\parallel}-(1/2)h_{\parallel}=0$
which gives

$\mu_{0}H_{\text{sat}}^{\parallel}=\frac{3s}{g_{\parallel}\mu_{\text{B}}}\left(J_{\perp}+2J_{\parallel}\right).$
(39)

###  E.2 Field in the ab plane

Assume the magnetic field ${\bf H}$ lies in the ab plane at an angle $\alpha$
relative to the $x$ direction. The classical energy density per spin then has
the form

$\displaystyle\frac{{\cal E}_{\perp}}{\nu s}$
$\displaystyle=sJ_{\perp}\left[\sin\theta_{1}\sin\theta_{2}\cos\left(\phi_{1}-\phi_{2}\right)+\sin\theta_{2}\sin\theta_{3}\cos\left(\phi_{2}-\phi_{3}\right)+\sin\theta_{3}\sin\theta_{1}\cos\left(\phi_{3}-\phi_{1}\right)\right]$
$\displaystyle+sJ_{\parallel}\left[\cos\theta_{1}\cos\theta_{2}+\cos\theta_{2}\cos\theta_{3}+\cos\theta_{3}\cos\theta_{1}\right]$
$\displaystyle-\frac{1}{3}h_{\perp}\left[\sin\theta_{1}\cos\left(\alpha-\phi_{1}\right)+\sin\theta_{2}\cos\left(\alpha-\phi_{2}\right)+\sin\theta_{3}\cos\left(\alpha-\phi_{3}\right)\right],$
(40) $\displaystyle h_{\perp}$
$\displaystyle=g_{\perp}\mu_{\text{B}}\mu_{0}H,$ (41)

independent of $J_{\Delta}$ and $J_{yz}$ as well. We may assume an umbrella-
shape spin structure right below the saturation field, however for
$J_{\parallel}\neq J_{\perp}$ the umbrella cannot have rotational symmetry
around the axis set by the direction of the magnetic field, and according to
work done in the context of (Cs,Rb)CuCl3 on similar Heisenberg models
energetically favorable spin configurations at high fields are planar [73, 74,
75, 76, 77, 78]. Therefore we don’t make assumptions about relations between
the six angles $\\{\theta_{i},\phi_{i}\\}$. With
$\delta_{i}=\pi/2-\theta_{i}$, $\epsilon_{i}=\alpha-\phi_{i}$ we write

$\displaystyle e_{\perp}$ $\displaystyle=\frac{{\cal E}_{\perp}}{\nu
S}=e_{12}+e_{23}+e_{31},$ (42) $\displaystyle e_{ij}$
$\displaystyle=sJ_{\perp}\cos\delta_{i}\cos\delta_{j}\cos\left(\epsilon_{i}-\epsilon_{j}\right)+sJ_{\parallel}\sin\delta_{i}\sin\delta_{j}$
$\displaystyle\phantom{=}-\frac{1}{6}h_{\perp}\left(\cos\delta_{i}\cos\epsilon_{i}+\cos\delta_{j}\cos\epsilon_{j}\right)$
(43) $\displaystyle=e_{ji}.$

In the limit $\delta_{i},\epsilon_{i}\to 0$ (full polarization), Eq. (38) for
$m_{\text{H}}(e_{\perp})$ reads

$\displaystyle\left(h_{\perp}-9sJ_{\perp}\right)^{2}\left[h_{\perp}-3s\left(J_{\parallel}+2J_{\perp}\right)\right]^{2}$
$\displaystyle\times$
$\displaystyle\left[h_{\perp}-6s\left(J_{\perp}-J_{\parallel}\right)\right]h_{\perp}$
$\displaystyle=0.$ (44)

In our case, we have $J_{\perp}>J_{\parallel}>0$, and the first factor
determines the saturation field when lowering the field from values above
saturation. We obtain

$\mu_{0}H_{\text{sat}}^{\perp}=\frac{9sJ_{\perp}}{\mu_{\text{B}}g_{\perp}}.$
(45)

## References

  * Note [1] The original name was given by the French mineralogist Charles Friedel (1832–1899) to the natural mineral CuFeO${}_{\text{2}}$ in honor of Gabriel Delafosse (1796–1878) [83, 84] in 1873\.
  * Marquardt  _et al._ [2006] M. A. Marquardt, N. A. Ashmore, and D. P. Cann, Crystal chemistry and electrical properties of the delafossite structure, Thin Solid Films Proceedings of the Fourth International Symposium on Transparent Oxide Thin Films for Electronics and Optics (TOEO-4), 496, 146 (2006).
  * Mackenzie [2017] A. P. Mackenzie, The properties of ultrapure delafossite metals, Rep. Prog. Phys. 80, 032501 (2017).
  * Schleid and Lissner [1993] T. Schleid and F. Lissner, Single crystals of NaMS2 (M=Ho–Lu) from reactions of the lanthanides with sulfur in the presence of NaCl, Eur. J. Solid State Inorg. Chem. 30, 829 (1993).
  * Huse and Elser [1988] D. A. Huse and V. Elser, Simple Variational Wave Functions for Two-Dimensional Heisenberg Spin-1/2 Antiferromagnets, Phys. Rev. Lett. 60, 2531 (1988).
  * Mila [2000] F. Mila, Quantum spin liquids, Eur. J. Phys. 21, 499 (2000).
  * White and Chernyshev [2007] S. R. White and A. L. Chernyshev, Neél Order in Square and Triangular Lattice Heisenberg Models, Phys. Rev. Lett. 99, 127004 (2007).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Hirakawa  _et al._ [1985] K. Hirakawa, H. Kadowaki, and K. Ubukoshi, Experimental Studies of Triangular Lattice Antiferromagnets with S=1/2: NaTiO2 and LiNiO2, J. Phys. Soc. Jpn. 54, 3526 (1985).
  * Clarke  _et al._ [1998] S. J. Clarke, A. J. Fowkes, A. Harrison, R. M. Ibberson, and M. J. Rosseinsky, Synthesis, Structure, and Magnetic Properties of NaTiO2, Chem. Mater. 10, 372 (1998).
  * Ramirez [2001] A. P. Ramirez, Chapter 4 Geometrical frustration, in  _Handbook of Magnetic Materials_, Vol. 13 (Elsevier, 2001) pp. 423–520.
  * Cotter  _et al._ [1994] J. P. Cotter, J. C. Fitzmaurice, and I. P. Parkin, New routes to alkali-metal–rare-earth-metal sulfides, J. Mater. Chem. 4, 1603 (1994).
  * Ranjith  _et al._ [2020] K. M. Ranjith, P. Schlender, T. Doert, and M. Baenitz, Magnetic Properties of an Effective Spin-$\frac{1}{2}$ Triangular-Lattice Compound LiYbS2, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2019)_, JPS Conference Proceedings, Vol. 30 (Journal of the Physical Society of Japan, 2020).
  * Hashimoto  _et al._ [2003] Y. Hashimoto, M. Wakeshima, and Y. Hinatsu, Magnetic properties of ternary sodium oxides NaLnO2 (Ln=rare earths), Journal of Solid State Chemistry 176, 266 (2003).
  * Ranjith  _et al._ [2019a] K. M. Ranjith, D. Dmytriieva, S. Khim, J. Sichelschmidt, S. Luther, D. Ehlers, H. Yasuoka, J. Wosnitza, A. A. Tsirlin, H. Kühne, and M. Baenitz, Field-induced instability of the quantum spin liquid ground state in the ${{J}}_{\mathrm{eff}}=\frac{1}{2}$ triangular-lattice compound ${\mathrm{{NaYbO}}}_{2}$, Phys. Rev. B 99, 180401(R) (2019a).
  * Sichelschmidt  _et al._ [2019] J. Sichelschmidt, P. Schlender, B. Schmidt, M. Baenitz, and T. Doert, Electron spin resonance on the spin-1/2 triangular magnet NaYbS2, J. Phys.: Condens. Matter 31, 205601 (2019).
  * Gray  _et al._ [2014] A. K. Gray, B. R. Martin, and P. K. Dorhout, Crystal structure of potassium ytterbium(III) selenide, KYbSe2, Zeitschrift für Kristallographie - New Crystal Structures 218, 20 (2014).
  * Ranjith  _et al._ [2019b] K. M. Ranjith, S. Luther, T. Reimann, B. Schmidt, P. Schlender, J. Sichelschmidt, H. Yasuoka, A. M. Strydom, Y. Skourski, J. Wosnitza, H. Kühne, T. Doert, and M. Baenitz, Anisotropic field-induced ordering in the triangular-lattice quantum spin liquid NaYbSe2, Physical Review B 100, 224417 (2019b).
  * Dong  _et al._ [2008] B. Dong, Y. Doi, and Y. Hinatsu, Structure and magnetic properties of ternary potassium lanthanide oxides KLnO2 (Ln=Y, Nd, Sm–Lu), Journal of Alloys and Compounds 453, 282 (2008).
  * Iizuka  _et al._ [2020] R. Iizuka, S. Michimura, R. Numakura, Y. Uwatoko, and M. Kosaka, Single Crystal Growth and Physical Properties of Ytterbium Sulfide $\text{{KYbS}}_{\text{2}}$ with Triangular Lattice, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2019)_, JPS Conference Proceedings, Vol. 30 (Journal of the Physical Society of Japan, 2020).
  * Seeger and Hoppe [1969] K. Seeger and R. Hoppe, Ternäre Oxide der Seltenen Erden vom Typ RbSEO2, Zeitschrift für anorganische und allgemeine Chemie 365, 22 (1969).
  * Bronger  _et al._ [1996] W. Bronger, J. Eyck, K. Kruse, and D. Schmitz, ChemInform Abstract: Ternary Rubidium Rare-Earth Sulfides; Synthesis and Structure., ChemInform 27, 10.1002/chin.199630031 (1996).
  * Bronger  _et al._ [1993] W. Bronger, W. Brüggemann, M. von der Ahe, and D. Schmitz, Zur Synthese und Struktur ternärer Chalcogenide der Seltenen Erden AlnX2 mit A $\equiv$ Alkalimetall und X $\equiv$ Schwefel, Selen oder Tellur, Journal of Alloys and Compounds 200, 205 (1993).
  * Deng and Ibers [2005] B. Deng and J. A. Ibers, CsYbSe2, Acta Cryst E 61, i15 (2005).
  * Xing  _et al._ [2019] J. Xing, L. D. Sanjeewa, J. Kim, G. R. Stewart, A. Podlesnyak, and A. S. Sefat, Field-induced magnetic transition and spin fluctuations in the quantum spin-liquid candidate ${\mathrm{{CsYbSe}}}_{2}$, Phys. Rev. B 100, 220407(R) (2019).
  * Duczmal and Pawlak [1994] M. Duczmal and L. Pawlak, Magnetic properties of TlLnS2 compounds (Ln $\equiv$ Nd, Gd, Dy, Er and Yb), Journal of Alloys and Compounds 209, 271 (1994).
  * Ferreira  _et al._ [2020] T. Ferreira, J. Xing, L. D. Sanjeewa, and A. S. Sefat, Frustrated Magnetism in Triangular Lattice TlYbS2 Crystals Grown via Molten Flux, Front. Chem. 8, 10.3389/fchem.2020.00127 (2020).
  * Miyasaka  _et al._ [2009] N. Miyasaka, Y. Doi, and Y. Hinatsu, Synthesis and magnetic properties of ALnO2 (A=Cu or Ag; Ln=rare earths) with the delafossite structure, Journal of Solid State Chemistry 182, 2104 (2009).
  * Sichelschmidt  _et al._ [2020] J. Sichelschmidt, B. Schmidt, P. Schlender, S. Khim, T. Doert, and M. Baenitz, Effective Spin-1/2 Moments on a Yb3+ Triangular Lattice: An ESR Study, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2019)_, JPS Conference Proceedings, Vol. 30 (Journal of the Physical Society of Japan, 2020).
  * Daszkiewicz  _et al._ [2008] M. Daszkiewicz, L. D. Gulay, V. Y. Shemet, and A. Pietraszko, Comparative Investigation of the Crystal Structure of LnCuSe2 Compounds (Ln = Tb, Dy, Ho, Er, Tm, Yb and Lu), Zeitschrift für anorganische und allgemeine Chemie 634, 1201 (2008).
  * Baenitz  _et al._ [2018] M. Baenitz, P. Schlender, J. Sichelschmidt, Y. A. Onykiienko, Z. Zangeneh, K. M. Ranjith, R. Sarkar, L. Hozoi, H. C. Walker, J.-C. Orain, H. Yasuoka, J. van den Brink, H. H. Klauss, D. S. Inosov, and T. Doert, ${\mathrm{{NaYbS}}}_{2}$: A planar spin-$\frac{1}{2}$ triangular-lattice magnet and putative spin liquid, Phys. Rev. B 98, 220409(R) (2018).
  * Rau and Gingras [2018] J. G. Rau and M. J. P. Gingras, Frustration and anisotropic exchange in ytterbium magnets with edge-shared octahedra, Phys. Rev. B 98, 054408 (2018).
  * Zhu  _et al._ [2018] Z. Zhu, P. A. Maksimov, S. R. White, and A. L. Chernyshev, Topography of Spin Liquids on a Triangular Lattice, Phys. Rev. Lett. 120, 207203 (2018).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Savary and Balents [2017] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2017).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A Field Guide to Spin Liquids, Annu. Rev. Condens. Matter Phys. 10, 451 (2019).
  * Fulde [1979] P. Fulde, Crystal fields, in  _Handbook on the Physics and Chemistry of Rare Earths_, Alloys and Intermetallics, Vol. 2 (Elsevier, 1979) Chap. 17, pp. 295–386.
  * Hutchings [1964] M. T. Hutchings, Point-Charge Calculations of Energy Levels of Magnetic Ions in Crystalline Electric Fields, in  _Solid State Physics_, Vol. 16, edited by F. Seitz and D. Turnbull (Academic Press, 1964) pp. 227–273.
  * Sakai  _et al._ [2000] O. Sakai, M. Kohgi, H. Shiba, A. Ochiai, H. Aoki, K. Takegahara, and H. Harima, Local Symmetry and Crystalline Field Effects in Charge-Ordered Yb4As3, J. Phys. Soc. Jpn. 69, 3633 (2000).
  * Shiba  _et al._ [2000] H. Shiba, K. Ueda, and O. Sakai, Effective Hamiltonian for Charge-Ordered Yb4As3, J. Phys. Soc. Jpn. 69, 1493 (2000).
  * Abragam and Bleaney [1970] A. Abragam and B. Bleaney,  _Electron Paramagnetic Resonance of Transition Ions_ (Clarendon Press, Oxford, 1970).
  * Orbach [1961] R. Orbach, Spin-Lattice Relaxation in Rare-Earth Salts, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 264, 458 (1961).
  * Ding  _et al._ [2019] L. Ding, P. Manuel, S. Bachus, F. Grußler, P. Gegenwart, J. Singleton, R. D. Johnson, H. C. Walker, D. T. Adroja, A. D. Hillier, and A. A. Tsirlin, Gapless spin-liquid state in the structurally disorder-free triangular antiferromagnet ${\mathrm{{NaYbO}}}_{2}$, Phys. Rev. B 100, 144432 (2019).
  * Zhang  _et al._ [2021] Z. Zhang, X. Ma, J. Li, G. Wang, D. T. Adroja, T. P. Perring, W. Liu, F. Jin, J. Ji, Y. Wang, Y. Kamiya, X. Wang, J. Ma, and Q. Zhang, Crystalline electric field excitations in the quantum spin liquid candidate ${\mathrm{{NaYbSe}}}\\_{2}$, Phys. Rev. B 103, 035144 (2021).
  * Luther [2019] S. Luther (2019), private communication.
  * Note [2] We use the hat-notation to indicate that the corresponding symbol relates to the full angular momentum $J$ with $j=7/2$.
  * Jensen and Mackintosh [1991] J. Jensen and A. R. Mackintosh,  _Rare Earth Magnetism: Structures and Excitations_ (Clarendon Press, Oxford, 1991).
  * Schmidt and Thalmeier [2017a] B. Schmidt and P. Thalmeier, Frustrated two dimensional quantum magnets, Physics Reports 703, 1 (2017a).
  * Custers  _et al._ [2003] J. Custers, P. Gegenwart, H. Wilhelm, K. Neumaier, Y. Tokiwa, O. Trovarelli, C. Geibel, F. Steglich, C. Pépin, and P. Coleman, The break-up of heavy electrons at a quantum critical point, Nature 424, 524 (2003).
  * Fulde  _et al._ [1995] P. Fulde, B. Schmidt, and P. Thalmeier, Theoretical Model for the Semi-Metal Yb4As3, EPL 31, 323 (1995).
  * Schmidt  _et al._ [1996] B. Schmidt, P. Thalmeier, and P. Fulde, Excitations in spin chains and specific-heat anomalies in Yb44As3, EPL 35, 109 (1996).
  * Moriya [1960] T. Moriya, Anisotropic Superexchange Interaction and Weak Ferromagnetism, Phys. Rev. 120, 91 (1960).
  * Li  _et al._ [2015] Y. Li, G. Chen, W. Tong, L. Pi, J. Liu, Z. Yang, X. Wang, and Q. Zhang, Rare-Earth Triangular Lattice Spin Liquid: A Single-Crystal Study of ${\mathrm{{YbMgGaO}}}_{4}$, Phys. Rev. Lett. 115, 167203 (2015).
  * Li  _et al._ [2016] Y.-D. Li, X. Wang, and G. Chen, Anisotropic spin model of strong spin-orbit-coupled triangular antiferromagnets, Phys. Rev. B 94, 035107 (2016).
  * Liu  _et al._ [2016] C. Liu, R. Yu, and X. Wang, Semiclassical ground-state phase diagram and $\text{multi-}{Q}$ phase of a spin-orbit-coupled model on triangular lattice, Phys. Rev. B 94, 174424 (2016).
  * Luo  _et al._ [2017] Q. Luo, S. Hu, B. Xi, J. Zhao, and X. Wang, Ground-state phase diagram of an anisotropic spin-$\frac{1}{2}$ model on the triangular lattice, Phys. Rev. B 95, 165110 (2017).
  * Iaconis  _et al._ [2018] J. Iaconis, C.-X. Liu, G. Halász, and L. Balents, Spin Liquid versus Spin Orbit Coupling on the Triangular Lattice, SciPost Physics 4, 003 (2018).
  * Schmidt and Thalmeier [2017b] B. Schmidt and P. Thalmeier, Néel temperature and reentrant ${H}\text{$-$}{T}$ phase diagram of quasi-two-dimensional frustrated magnets, Phys. Rev. B 96, 214443 (2017b).
  * Obradors  _et al._ [1988] X. Obradors, A. Labarta, A. Isalgué, J. Tejada, J. Rodriguez, and M. Pernet, Magnetic frustration and lattice dimensionality in SrCr8Ga4O19, Solid State Commun. 65, 189 (1988).
  * Ramirez [1994] A. P. Ramirez, Strongly Geometrically Frustrated Magnets, Annu. Rev. Mater. Sci. 24, 453 (1994).
  * Dey  _et al._ [2020] S. Dey, E. C. Andrade, and M. Vojta, Destruction of long-range order in noncollinear two-dimensional antiferromagnets by random-bond disorder, Phys. Rev. B 101, 020411(R) (2020).
  * Wu  _et al._ [2019] L. S. Wu, S. E. Nikitin, M. Brando, L. Vasylechko, G. Ehlers, M. Frontzek, A. T. Savici, G. Sala, A. D. Christianson, M. D. Lumsden, and A. Podlesnyak, Antiferromagnetic ordering and dipolar interactions of ${\mathrm{{YbAlO}}}_{3}$, Phys. Rev. B 99, 195117 (2019).
  * Maksimov  _et al._ [2019] P. A. Maksimov, Z. Zhu, S. R. White, and A. L. Chernyshev, Anisotropic-Exchange Magnets on a Triangular Lattice: Spin Waves, Accidental Degeneracies, and Dual Spin Liquids, Phys. Rev. X 9, 021017 (2019).
  * Sedrakyan  _et al._ [2020] T. Sedrakyan, R. Moessner, and A. Kamenev, Helical spin liquid in a triangular XXZ magnet from Chern-Simons theory, Phys. Rev. B 102, 024430 (2020).
  * Castner and Seehra [1971] T. G. Castner and M. S. Seehra, Antisymmetric Exchange and Exchange-Narrowed Electron-Paramagnetic-Resonance Linewidths, Phys. Rev. B 4, 38 (1971).
  * Zorko  _et al._ [2008] A. Zorko, S. Nellutla, J. van Tol, L. C. Brunel, F. Bert, F. Duc, J.-C. Trombe, M. A. de Vries, A. Harrison, and P. Mendels, Dzyaloshinsky-Moriya Anisotropy in the Spin-1/2 Kagome Compound ${\mathrm{{ZnCu}}}_{3}(\mathrm{{OH}}{)}_{6}{\mathrm{{Cl}}}_{2}$, Phys. Rev. Lett. 101, 026405 (2008).
  * Yamamoto  _et al._ [2014a] D. Yamamoto, G. Marmorini, and I. Danshita, Quantum Phase Diagram of the Triangular-Lattice XXZ Model in a Magnetic Field, Phys. Rev. Lett. 112, 127203 (2014a).
  * Yamamoto  _et al._ [2014b] D. Yamamoto, G. Marmorini, and I. Danshita, Erratum: Quantum Phase Diagram of the Triangular-Lattice XXZ Model in a Magnetic Field [Phys. Rev. Lett. 112, 127203 (2014)], Phys. Rev. Lett. 112, 259901(E) (2014b).
  * Yamamoto  _et al._ [2017] D. Yamamoto, H. Ueda, I. Danshita, G. Marmorini, T. Momoi, and T. Shimokawa, Exact diagonalization and cluster mean-field study of triangular-lattice XXZ antiferromagnets near saturation, Phys. Rev. B 96, 014431 (2017).
  * Kawamura and Miyashita [1985] H. Kawamura and S. Miyashita, Phase Transition of the Heisenberg Antiferromagnet on the Triangular Lattice in a Magnetic Field, J. Phys. Soc. Jpn. 54, 4530 (1985).
  * Seabra  _et al._ [2011] L. Seabra, T. Momoi, P. Sindzingre, and N. Shannon, Phase diagram of the classical Heisenberg antiferromagnet on a triangular lattice in an applied magnetic field, Phys. Rev. B 84, 214418 (2011).
  * Wu  _et al._ [2020] M. Wu, D.-X. Yao, and H.-Q. Wu, Exact diagonalization study of the anisotropic Heisenberg model related to YbMgGaO4, arXiv:2008.08751 [cond-mat] (2020).
  * Kawamura [1984] H. Kawamura, Spin-Wave Analysis of the Antiferromagnetic Plane Rotator Model on the Triangular Lattice–Symmetry Breaking in a Magnetic Field, J. Phys. Soc. Jpn. 53, 2452 (1984).
  * Rastelli  _et al._ [1992] E. Rastelli, A. Tassi, A. Pimpinelli, and S. Sedazzari, Triangular planar antiferromagnet in an external magnetic field, Phys. Rev. B 45, 7936 (1992).
  * Jacobs  _et al._ [1993] A. E. Jacobs, T. Nikuni, and H. Shiba, Theory of Magnetic Structures of CsCuCl3 in Transverse Magnetic Field, J. Phys. Soc. Jpn. 62, 4066 (1993).
  * Ohyama and Jacobs [1995] T. Ohyama and A. E. Jacobs, Incommensurate state of ${\mathrm{{CsCuCl}}}_{3}$ in a transverse magnetic field, Phys. Rev. B 52, 4389 (1995).
  * Nikuni and Jacobs [1998] T. Nikuni and A. E. Jacobs, Quantum fluctuations in the incommensurate phase of ${\mathrm{{CsCuCl}}}_{3}$ in a transverse magnetic field, Phys. Rev. B 57, 5205 (1998).
  * Jacobs and Nikuni [2002] A. E. Jacobs and T. Nikuni, Magnetic structures of ${\mathrm{{RbCuCl}}}_{3}$ in a transverse field, Phys. Rev. B 65, 174405 (2002).
  * Li [2019] Y. Li, YbMgGaO4: A Triangular-Lattice Quantum Spin Liquid Candidate, Advanced Quantum Technologies 2, 1900089 (2019).
  * Guo  _et al._ [2019] S. Guo, A. Ghasemi, C. L. Broholm, and R. J. Cava, Magnetism on ideal triangular lattices in $\mathrm{{NaBaYb}}{(\mathrm{{B}}{\mathrm{{O}}}_{3})}_{2}$, Phys. Rev. Materials 3, 094404 (2019).
  * Guo  _et al._ [2020] S. Guo, R. Zhong, K. Górnicka, T. Klimczuk, and R. J. Cava, Crystal Growth, Structure, and Magnetism of the 2d Spin 1/2 Triangular Lattice Material Rb3Yb(PO4)2, Chem. Mater. 32, 10670 (2020).
  * Rotter [2017] M. Rotter, McPhase manual: Stevens Operators (2017).
  * Friedel [1873] C. Friedel, Sur une combinaison naturelle des oxydes de fer et de cuivre, et sur la reproduction de l’atacamite, Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences 77, 211 (1873).
  * Rogers [1913] A. F. Rogers, Delafossite, a cuprous metaferrite from Bisbee, Arizona, American Journal of Science s4-35, 290 (1913).