Extracting quantum magnetic model parameters in a triangular rare-earth magnet: systematic characterization of crystal electric field effects from magnetic susceptibility and resonant torsion magnetometry

Christopher A. Pocs, Peter E. Siegfried, Jie Xing, Athena S. Sefat, Michael Hermele, 3 B. Normand, 5 and Minhyea Lee Department of Physics, University of Colorado, Boulder, Colorado 80309, USA Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Center for Theory of Quantum Matter, University of Colorado, Boulder, Colorado 80309, USA Paul Scherrer Institute, CH-5232 Villigen-PSI, Switzerland Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland (Dated: October 8, 2021)

A primary goal at the interface of theoretical and experimental quantum magnetism is the investigation of exotic spin states, mostly notably quantum spin liquids (QSLs), that realize phenomena including quasiparticle fractionalization, long-ranged entanglement, and topological order. Magnetic rare-earth ions go beyond the straightforward paradigm of geometrical frustration in Heisenberg antiferromagnets by introducing competing energy scales, and in particular their strong spin-orbit coupling creates multiple split crystal electric-field (CEF) levels, leading to anisotropic effective spin models with intrinsic frustration. While rare-earth delafossites have a triangular-lattice geometry, and thus have gained recent attention as candidates for hosting spin-1/2 QSL physics, the reliable extraction of effective spin models from the initial many-parameter CEF spectrum is a hard problem. Using the example of CsYbSe2, we demonstrate the unambiguous extraction of the Stevens operators dictating the full CEF spectrum of Yb 3+ by translating these into parameters with a direct physical interpretation. Specifically, we combine low-field susceptibility measurements with resonant torsion magnetometry (RTM) experiments in fields up to 60 T to determine a sufficiently large number of physical parameters -effective Zeeman splittings, anisotropic van Vleck coefficients, and magnetotropic coefficients -that the set of Stevens-operator coefficients is unique. Our crucial identification of the strong corrections to the Zeeman splitting of Kramers doublets as van Vleck coefficients has direct consequences for the interpretation of all anisotropic magnetic susceptibility measurements. Our results allow us to determine the nature and validity of an effective spin-1/2 model for CsYbSe2, to provide input for theoretical studies of such models on the triangular lattice, and to provide additional materials insight into routes for achieving magnetic frustration and candidate QSL systems in rare-earth compounds.

I. INTRODUCTION
The quantum spin liquid (QSL) is a nonmagnetic many-body ground state in which the spin correlations have long-ranged quantum entanglement [1]. Despite many theoretical studies of these enigmatic phases, which have served to drive detailed investigations of a wide range of magnetic compounds, neither a universally agreed QSL phase in a real material nor an unambiguous set of experimental QSL criteria has yet emerged. The strong quantum fluctuations responsible for producing exotic spin states are a consequence of generalized magnetic frustration, which leads to a highly degenerate manifold of competing states. Early examples of material realizations of candidate models for hosting QSL states were based on geometrical frustration in structures with triangular motifs, including kagome [2], pyrochlore [3,4], and triangular lattices [5][6][7][8], mostly of real S = 1/2 spins.
Magnetic insulators with strong spin-orbit coupling are now widely recognized as a platform for extending very significantly the nature of frustration and the variety of quantum many-body phases (including QSLs) and phenomena that can be realized. When compared to spin-1/2 magnetic insulators based on 3d ions, these sys-tems tend to exhibit complex microscopic physics even at the single-ion level. Magnetic materials containing 5d and 4d transition-metal ions possess interactions that are anisotropic in both spin space and real space, leading to complex phenomenology in pyrochlore systems [3,9,10] and in proximate Kitaev materials [11][12][13][14][15]. Compounds based on 4f rare-earth ions that combine the geometric frustration of pyrochlore and triangular lattices with strong spin-orbit coupling have also provided fertile ground for quantum magnetism research [16][17][18][19]. However, in these materials a detailed understanding of the microscopic physics is a prerequisite for developing effective spin models that serve as a basis for theories of many-body phenomena.
In insulating 4f materials, the key to a microscopic description of the magnetic interactions is a reliable determination of the single-ion crystal electric-field (CEF) Hamiltonian (H CEF ) [3,4]. In many rare-earth compounds, CEF effects split the degeneracy of the ground multiplet of the total angular momentum, J. For halfodd-integer J, the result is a set of Kramers doublets. At sufficiently low temperatures, a restriction of the dynamics to the lowest such doublet may be invoked to justify using an effective pseudospin-1/2 model [16,[20][21][22]. From an experimental standpoint this allows a significant diversification of both materials and magnetic phenomena, while on the theoretical side a minimal model may be relatively simple, or even a realization of one of the paradigm models in frustrated magnetism.
The leading system-specific parameters in H CEF are commonly determined by inelastic neutron scattering (INS), which provides direct information about the spectrum of CEF multiplets. However, there are in general more symmetry-allowed parameters in H CEF than there are gaps between CEF multiplets, resulting in an underdetermined fit of the CEF parameters and hence to uncertainties in the appropriate pseudospin-1/2 model that can be qualitative rather than merely quantitative corrections. Methods allowing the unambiguous determination of H CEF with higher reliability are thus very desirable. In this article, we present a high-fidelity determination of the CEF parameters of a selected 4f system, obtained by using a combination of low-field magnetic susceptibility and high-field (up to 60 T) resonant torsion magnetom-etry (RTM) measurements.
Yb-based triangular-lattice compounds offer an excellent combination of geometric frustration, quasi-twodimensional nature, half-odd-integer J, and a strongly split CEF spectrum that suggests the validity of pseudospin-1/2 models of quantum magnetism. In particular, the family of AYbX 2 delafossites (with A = Na, Cs and X = O, S, Se) has attracted intensive interest [21][22][23][24][25][26][27][28][29], and a very recent summary was compiled in Ref. [30]. Unlike the material YbMgGaO 4 [31][32][33][34][35][36], they are free from potential site disorder [37,38] and no members of the family have been found to exhibit magnetic ordering at temperatures down to tens of mK [30], while several have been reported to exhibit continua of low-energy magnetic excitations [21,22,39]. Their magnetic response is highly anisotropic between the in-and out-of-plane directions, providing a valuable opportunity to use the magnetic field to vary the free-energy landscape. Thus the trianguar-lattice delafossites present an ideal test case for unambiguous fitting of the fieldinduced CEF spectra and subsequent establishment of the effective spin-1/2 states in rare-earth compounds.
As an example material we focus on The complete CEF spectrum is specified by the coefficients of six Stevens operators, and thus the challenge is to measure enough independent physical quantities beyond the low-energy limit. From the susceptibility we extract not only the linear Zeeman coefficients for both of the primary, high-symmetry field directions but also the van Vleck (VV) coefficients for the ground doublet, which are the second-order corrections (i.e. quadratic in the applied field). While these VV coefficients are often used to describe the magnetic response of rare-earth compounds in terms of an additive contribution to the susceptibility [25,40], in fact they are embedded in a nontrivial way in the full expression for this quantity. We will show how the multiple roles of the ground-state VV coefficients and their significance as stand-alone physical quantities, with direct implications for the high-field magnetic anisotropy, have yet to be appreciated in connection with CEF fitting. They facilitate the bridge to the field range covered by our RTM measurements, from which we extract the full field-and temperature-dependence of the magnetotropic coefficients for the same two highsymmetry field directions. These round out a complete set of eight observables, allowing us to determine a unique set of microscopic CEF parameters and thus the full CEF spectrum. As Fig. 1(c) makes clear, our results reveal an intricate energy landscape and level-repulsion behavior in the CEF spectrum of CsYbSe 2 up to high magnetic field values and for both field directions.
The physical content of our analysis is to interpret the essential role of the VV coefficients in determining magnetic properties, even at low fields and temperatures, where they are often neglected in pseudospin-1/2 models. Quantifying the VV corrections allows us to demarcate the field-temperature range over which a pseudospin-1/2 approximation is appropriate in CsYbSe 2 , and to describe the high-field limit accurately. The resulting full characterization of the CEF spectra is essential for investigating their experimental consequences, for example when the application of intense magnetic fields causes multiple real or avoided level-crossings in narrowly spaced CEF spectra. It is also a prerequisite to study further mechanisms leading to different forms of magnetic frustration, in which additional degrees of freedom, such as phonons, hybridize with the electronic spectrum to cause profound effects to appear in low-temperature thermodynamics and transport properties [41,42]. Finally, while we have focused on one example material, our analysis can be applied widely to localized 4f -electron systems.
The structure of this paper is as follows. In Sec. II we introduce CsYbSe 2 , our experimental methods, the complete Stevens operator formalism for CEF levels and some approximate treatments. In Sec. III we show all the results of our susceptibility measurements and their analysis for the two primary field directions. Section IV presents our RTM measurements and extraction of the magnetotropic coefficients, allowing us to determine the full set of Stevens operators. In Sec. V we discuss the physical interpretation of the results and their consequences for experimental analysis, effective pseudospin-1/2 models in frustrated quantum magnetism, and materials selection for candidate QSL systems. Section VI contains a brief summary of our contribution and four Appendices provide additional information concerning details of the analysis.

II. MATERIALS AND METHODS
While most of the AYbX 2 family crystallizes in the R3bm space group [30], the layer stacking sequence of CsYbSe 2 gives a P 6 3 /mmc structure [43,44], in which the triangular-lattice planes are constructed from edgesharing YbSe 6 octahedra, as illustrated in Figs. 1(a-b). The ratio of inter-and intra-layer distances separating the Yb 3+ ions results in a quasi-2D magnetic system [43]. Within each plane, the Se atoms mediate identical AFM superexchange interactions between all nearest-neighbor Yb 3+ ion pairs, resulting in highly frustrated triangular magnetism. In fact the system does not order at zero field for temperatures down to 0.4 K, although an applied in-plane field induces an up-up-down ordering, as demonstrated by the observation of a 1/3 plateau in the magnetization, M (H), at this temperature [43].
Magnetization and susceptibility measurements were performed on single-crystalline samples using a Quantum Design MPMS. All susceptibilities we report are obtained from χ ab,c = dM ab,c dH , where the indices denote measurements performed with the field oriented in the triangularlattice plane (H ab) or perpendicular to it (H c). Beyond the 7 T upper limit of our magnetization measurements, we employ resonant torsion magnetometry (RTM) to probe the nature of the magnetic anisotropy. The sample is mounted on the tip of a vibrating cantilever and the measured shifts in the resonant frequency (f 0 ≈ 40 kHz) reflect changes in the magnetic rigidity caused by changes in the direction of the applied magnetic field, which are quantified by a tensor of magnetotropic coefficients, k( H). We focus on k(H, T ) at the high-symmetry angles θ = π/2 (to measure k ab ) and 0 (k c ), where θ is the polar angle of the applied field measured from the crystalline c-axis. The measurement configuration is summarized in Sec. IV and described in detail elsewhere [45,46]. All RTM data were taken using the capacitive magnet of the NHMFL pulsed-field facility at Los Alamos National Laboratory.

A. CEF Analysis and Model Hamiltonian
Yb 3+ ions subject to the CEF interactions of their surrounding anion charge distribution have a total J = 7/2 ground-state multiplet of allowed electronic states. Unlike the triangular-lattice material YbMgGaO 4 , the Yb 3+ ions in CsYbSe 2 are largely free from any site disorder associated with the non-magnetic ions. To describe CsYbSe 2 , we parametrize the single-ion CEF interaction as a linear combination of the six symmetry-allowed Stevens operators,Ô m n , for the D 3d site symmetry of the YbSe 6 octahedral environment [21,[23][24][25][26][28][29][30]  (1) This Hamiltonian splits the J = 7/2 multiplet into four Kramers doublets, |n ± with n = 0, 1, 2, 3, whose energies we use to define the separations from the groundstate doublet (n = 0) as ∆ 10 , ∆ 20 , and ∆ 30 [ Fig. 1(c)]. The corresponding wave functions are obtained by diagonalizing Eq. (1). A CEF spectrum for the zero-field limit can be identified using spectroscopic probes, particularly INS and Raman spectroscopy, but to date little information is available with which to investigate the high-field reorganization of the energy spectrum.
The symmetries of the triangular lattice of Yb 3+ ions permit a nearest-neighbor superexchange interaction with XXZ spin symmetry [47]. Additional symmetryallowed and bond-dependent anisotropic pseudo-dipolar exchange terms [48] are found to give vanishing contributions in a standard Weiss mean-field approximation. Thus we restrict our analysis to the minimal XXZ spin model,Ĥ XXZ , describing nearest-neighbor interactions between adjacent in-plane J = 7/2 moments in terms of two interactions, J ⊥ and J z , which are the respective couplings of spin components transverse and parallel tô z, i.e. sites andĴ i,γ , with γ = x, y, z, labels the components of spin J = 7/2 moments on site i.
To account for Zeeman coupling to the external field we add the termĤ where the Landé g-factor, g J = 8/7, is used for Yb 3+ moments. We note that the quantization axis defining theẑ direction of the chosen basis of spin operators in H CEF andĤ XXZ is identically the crystallographic c-axis, whence the component H z refers to a field (H c) applied along the c-axis in experiment. Similarly, H ⊥ = H ab refers to a field component perpendicular to the c axis, which lies precisely in the hexagonally symmetric abplane. In none of our experiments (magnetization, susceptibility, RTM) did we find a discernible difference in the response for different in-plane field directions, and hence we do not distinguish between these. The two separate contributions to the spin response under an external magnetic field, arising from the single-ion anisotropy and the superexchange anisotropy (both of which have their origin in the CEF spectrum), are then captured by the Hamiltonian

B. Weiss Mean-Field Approximation
We treat the physics of the system at finite temperature within a self-consistent Weiss mean-field approximation for the Yb 3+ spins, which we will find captures the bulk magnetic behavior of CsYbSe 2 exceptionally well at all measurement temperatures (T ≥ 2 K). In the Weiss mean-field treatment, Eq. (4) reduces to a system of N decoupled Yb 3+ ions each subject to the Hamiltonian where q = 6 is the nearest-neighbor coordination number on the triangular lattice and the mean-field expectation values of the spin operators are determined self-consistently as Because we find that the bulk magnetic susceptibility of CsYbSe 2 is almost entirely uniaxial, we make the additional simplifying assumption that the external magnetic field lies in the xzplane, such that H = H ⊥x + H zẑ .
The precise determination of the mean-field Hamiltonian given in Eq. as noted above measuring the three CEF level-splittings, ∆ 10 , ∆ 20 , and ∆ 30 [ Fig. 1(c)], is not sufficient to determine six unknowns. Electron spin resonance (ESR) [27] probes a much lower frequency range, making it the method of choice for determining effective g-factor values in the ground doublet, while the T -dependence of the line width can be used to estimate ∆ 10 [23,25,30]; in very well characterized systems, the line width may also be used to discuss the effective exchange anistropy between J ⊥ and J z [27].
Thus the parameter space for fitting the measured CEF spectra usually remains highly degenerate, even with accurate spectral measurements in an applied field and well-constrained g-tensors, and hence the uniqueness of a fitted set of Stevens-operator coefficients cannot be guaranteed [28,29,37]. The method we apply to solve this problem has two key components. First we apply a detailed consideration of the second-order corrections in the field-dependence of the CEF spectrum, encoded as the VV coefficients we extract from the low-field susceptibility. Then we leverage extensive RTM data providing systematic temperature-and field-dependent information about the magnetotropic coefficients for the two highsymmetry field directions to fix a unique set {B n m }.

III. ANISOTROPIC MAGNETIC SUSCEPTIBILITIES
A. Experiment: deviations from Curie-Weiss In Fig. 2 we show the temperature-dependences of the two magnetic susceptibilities, χ ab (T ) and χ c (T ), measured for CsYbSe 2 in the presence of an external magnetic field µ 0 H = 1 T applied respectively within the ab plane and along the c axis. There is no indication of long-range ordering down to T = 2 K. The crossing of the two curves around T = 35 K reflects a crossover from easy-plane behavior (∆χ = χ ab − χ c > 0) at low temperatures to easy-axis anisotropy (∆χ < 0) at high T .
Both susceptibilities increase rapidly as the temperature is lowered below 50 K, and the corresponding inverse quantities shown in the inset of Fig. 2 make clear the susceptibility anisotropy, with χ −1 c (T ) exhibiting a significantly sharper downward trend as T decreases. While a qualitative inspection suggests that a Curie-Weiss form might capture χ −1 ab (T ) for T < 50K, this is clearly impossible for χ −1 c (T ). This type of low-T downturn in χ −1 (T ) is observed quite commonly in similar rare-earth magnetic materials [16,20,23,25,40,49], but to date lacks a detailed analysis. Next we show that this behavior can be explained entirely by analyzing the field-induced evolution of the CEF energy levels, where the leading deviations from a linear (effective Zeeman) form are contained in two strongly anisotropic ground-state VV coefficients.

B. Van Vleck coefficients and low-field limit
The splitting of the CEF levels in a finite applied magnetic field can be understood systematically by perturbation theory. We express the H-linear and -quadratic corrections to the energy eigenvalues ofĤ CEF due toĤ Z in the form where E 0 n refers to the four CEF energy levels at zero field (and setting the ground-state energy, E 0 0 , to zero gives E 0 n = ∆ n0 for n = 1, 2, and 3). The second term describes the H-linear Zeeman splitting with a generalized g-tensor for all levels (n = 0, 1, 2, 3) defined as g n ⊥(z) = 2g J n ± |Ĵ x(z) |n ∓(±) , where the subscripts ⊥, z denote a field applied respectively within the ab-plane or along the c-axis. The conventional effective spin-1/2 g-tensor components are g 0 ⊥ and g 0 z , to which we refer henceforth, without the superscript 0, as the g-factors of the system.
We define all of the phenomena obtained at second order in the perturbative effect ofĤ Z on the zero-field eigenstates ofĤ CEF as the VV contribution. To describe the full field-and temperature-dependence of the susceptibility we define the VV coefficients, α n ⊥ and α n z , for the nth CEF level by where Although it is often stated that one may define a "VV susceptibility," χ VV , that is a small, additive, temperature-independent, and paramagnetic contribution to the total susceptibility, this is rarely an accurate approximation. By inspecting the form of the VV coefficients, one observes that the n = 0 terms should be negative, giving the expected type of 2nd-order correction to the ground state. The physics content of the anisotropic VV coefficients can be read from Fig. 1(c), where the field-dependence of the CEF levels is quite strongly nonlinear due to level-repulsion effects between adjacent doublets. In CsYbSe 2 this repulsion, which is equivalent to a negative curvature of the lower (and positive of the upper) branch in each case, is strongest between the n = 0 and 1 doublets for H c.
We obtain analytical formulas for the low-field magnetic susceptibilities of the system by using the Nparticle partition function calculated with Eq. (5). The full expression for χ ab(c) (T ), which includes the contributions of all four CEF levels, is presented in Eq. (A1) of App. A. Here we focus on the regime of temperatures sufficiently small that only the lowest-lying Kramers doublet need be considered, i.e. k B T ≪ ∆ 10 , and write the inverse susceptibility as We have defined the quantities Θ CW ⊥(z) = qJ ⊥(z) /k B in order to obtain an adapted Curie-Weiss (CW) form for the two applied-field directions and in the second line we have used the definition of the g-factors (above) to make this form more transparent. We note that the regime of validity of Eq. (8) is also that in which the system can be approximated by an effective pseudospin-1/2 description.
It is clear from Eq. (8) that the contrast to a conventional CW form, χ −1 ∝ T + Θ, is the additional Tlinear term in the denominator, whose prefactor is the corresponding VV coefficient. An explicit comparison is presented in App. B. The key advantage of our formulation is to observe that the same (VV) coefficients describing the H 2 correction that gives the leading nonlinear contribution to the CEF levels at low temperatures and finite fields [Eq. (6)] are those describing the nonlinear, "beyond-CW" form of the susceptibility at zero field and finite temperatures [Eq. (8)]. Thus one may TABLE I. Comparison of the lowest zero-field CEF level-splitting, ∆10, g-tensor components, and the VV coefficients for n = 0 in a number of Yb-based triangular-lattice compounds. For CsYbSe2 we show two sets of g-factors and VV coefficients, one obtained by a full calculation [Eq. (7)] using the B n m coefficients and wave functions, i.e. assisted by fitting to the RTM data as shown in Sec. IV, and the other from fitting χ ab,c (T ) using Eq. (8). For the other compounds, ∆10 and the g-factors are quoted from the respective references and we calculated the VV coefficients using the reported B n m coefficients. The compound nominally most similar to CsYbSe2, NaYbSe2, displays slightly less anisotropy in its VV coefficients and this is consistent with the differing forms of χ ab (T ) and χc(T ) reported in Ref. [25]. For NaYbO2 we note that the two sets of Stevens coefficients (Fit 1 and Fit 2) deduced in Ref. [28] yield dramatically different values of α 0 ⊥ and α 0 z , which underlines the crucial role of the VV coefficients in a complete and consistent characterization of HCEF. YbMgGaO4 is found to be least anisotropic among these materials, and its small VV coefficients are consistent with the reported Curie-Weiss form of the susceptibility.
conclude that the latter effect is also related to levelrepulsion between doublets, the fact that the susceptibility is a second field derivative of the free energy meaning that second-order perturbative effects do not vanish in the limit H → 0.
Returning to Fig. 2, the gray dashed lines in the inset show fits to Eq. (8), made in the regime T < 45 K, for each field direction. From these fits we obtain two exchange energies, J ⊥ and J z in Eqs. (2) and (5), and two VV coefficients for the ground-state doublet. A complete fit is deferred to Sec. V. Starting with the VV coefficients, we find the values α 0 ⊥ = −(2.9 ± 0.1) × 10 −4 meV/T 2 from χ −1 ab and α 0 z = −(17.9±0.1)×10 −4 meV/T 2 from χ −1 c . We stress that α 0 z is six times larger than α 0 ⊥ , which for all fields beyond 5 T becomes clearly manifest as a much larger downward level-repulsion of the lowest Kramers doublet [ Fig. 1(c) and the quantitative analysis of Sec. V]. Turning to the magnetic interactions, we find J ⊥ = 0.54 ± 0.01 K from χ ab and J z = 0.61 ± 0.01 K from χ c , both encapsulated in Θ CW ⊥(z) . If one reduces the system to an effective pseudospin-1/2 model, the corresponding interaction terms are J ′ ⊥ = 5.12 K (≃ 0.44 meV) and J ′ z = 0.84 K (≃ 0.07 meV), as detailed in App. C. We discuss the physical implications of these interaction parameters in Sec. VB.
As noted above, downward curvature of the lowtemperature χ −1 has been observed in other rare-earth compounds and our fitting results demonstrate that this feature should be characterized by using the VV coefficients as a part of a full description of the anisotropic magnetism. The validity of Eq. (8) as a replacement for the CW form of the susceptibility is confirmed by capturing the different T -dependences correctly for the two field directions with two VV coefficients that are consistent with Eq. (6). Although the VV coefficients are not parameters appearing directly in the system Hamiltonian, they impose constraints that are essential for a unique determination of the full parameter set, a topic we discuss further in the next section.
In Table 1 we compare the VV coefficients and other physical characteristics for a variety of Yb delafossites. For all compounds listed other than CsYbSe 2 , the ∆ 10 and g-tensor parameters are taken from experimental data. All α 0 ⊥ and α 0 z values were calculated from Eq. (7) using the Stevens coefficients (B n m ) provided by each reference, where they were obtained from fits to the CEF spectrum obtained by INS. We stress again the fact that, in several of the studies cited, different sets of Stevens coefficients can provide equally good descriptions ("Fit 1" and "Fit 2") of the same INS data due to the underconstrained nature of the problem. The anisotropy of CsYbSe 2 , α 0 z /α 0 ⊥ ≈ 6, is strikingly higher than the values reported for other compounds in the same family. Although all of the NaYbX 2 materials seem to show considerable directional anisotropy [30], this may be less severe for NaYbO 2 , except that the two sets of proposed B n m parameters yield wildly different VV coefficients, indicative of an underlying ambiguity of the type we demonstrate how to resolve. In this regard the sole non-delafossite in Table 1, YbMgGaO 4 , is a nearly isotropic outlier.
Before proceeding, we reiterate two important attributes of the VV coefficients as a mean of characterizing the magnetic anisotropy of a CEF system. First, nonzero VV coefficients are immediately evident in χ(T ), as strong deviations from a CW form, and thus failure to account for them means that the interaction parameters, J ⊥,z , cannot be determined correctly. Second, it is evident from Eq. (8) that caution is required in applying an effective pseudospin-1/2 description, because even at low T , where only the lowest Kramers doublet is thermally populated, the repulsion from the higher CEF levels can-not be ignored. Hence the twin roles of the anisotropic VV coefficients in dictating H-and T -dependent physical properties generic to many 4f electronic systems (Table  1) must be taken into account to obtain a meaningful description of the spin physics.

A. Dependence on Field and Temperature
Having characterized the magnetic response at low fields using χ(T ), and thereby obtained four independent physical quantities to include in the fitting procedure, we turn for more information to the magnetropic coefficients. The RTM method allows these to be measured over wide ranges of both field and temperature, which we will show provides enough input for an unambiguous determination of all the remaining free parameters in Eq. (5). The magnetotropic coefficient is defined as k(H, ϑ) = ∂ 2 F (H, ϑ)/∂ϑ 2 , where F (H, ϑ) is the portion of the Helmholtz free energy depending on the magnitude and orientation of the magnetic field (ϑ is the angular direction ofĤ measured in the plane of vibration of the sample [45,46]). k quantifies the magnetic rigidity of a material, whose origin lies in the energy cost of rotating a sample with an anisotropic magnetic free energy in a finite field, and RTM measurements constitute a highly sensitive probe of this magnetoanisotropy.
The variation of the magnetotropic coefficients at finite field produces a shift in the resonant frequency of a system composed of a cantilever and an attached sample given by with K the intrinsic bending stiffness of the cantilever at zero field [45,46]. In the low-field limit, where χ is constant at constant temperature, both the torque and the magnetotropic coefficient are not only straightforward functions of the polar angle, θ, but are also quadratic in H. Thus the magnetotropic coefficient can be expressed in terms of the difference, ∆χ = χ ab − χ c , between the principal components of the susceptibility tensor in the plane of vibration, as k = ∆χH 2 cos 2θ. More detailed derivations of field-dependent expressions for the highsymmetry angles (k ab for θ = 0 and k c for θ = π/2) can be deduced by introducing the transverse susceptibility, as shown in App. D. Figures 3(a) and 3(b) show the frequency shifts, ∆f ab and ∆f c , which are directly proportional to the magnetotropic coefficients, k ab and k c [Eq. (9)], as functions of field in the H ab and H c geometries. We note that the vibration plane of the cantilever and the rotation plane of H are the same, meaning that in this experimental configuration the angle ϑ of Eq. (D2) is identically the polar angle, θ, relative to the crystallographic c-axis. At low fields, where χ ab and χ c are H-independent constants (i.e. the magnetizations m ab and m c are linear in H), both k ab and k c are indeed proportional both to H 2 and to ∆χ with opposing signs. As the temperature is increased, ∆χ(T ) in Fig. 2 changes its sign above 50 K, and this is reflected in the sign-changes of ∆f ab and ∆f c between T = 30 and 50 K.
As the field is increased, ∆f ab (H) deviates from an H 2 dependence and becomes nonmonotonic, with a local maximum at low temperatures [T ≤ 30 K in Fig. 3(a)]. The origin of this behavior lies primarily in m ab (H) saturating in this field range, which we confirm from our calculations in Sec. V. By contrast, the behavior of ∆f c (H) remains monotonic in the same T range, where m c (H) continues to increase with H. We comment that ∆f c (H) does exhibit a weak maximum in H at T = 70 K, whereas no such behavior is found in ∆f ab at this temperature. This feature is also captured qualitatively by calculating the eigenstates of the full mean-field Hamiltonian [Eq. (5)], as we show next.

B. Fitting k ab and kc
Taking the magnetic interaction parameters J ⊥ and J z determined from χ ab,c (T ) [Fig. 2], we fit the measured data for k ab (H) and k c (H) using the full meanfield Hamiltonian based on Eq. (5), in which the six coefficients of the Stevens operators are free parameters to be determined. However, obtaining the ground-state VV coefficients, α 0 ⊥ and α 0 z , from the fit to Eq. (8) reduces the number of fitting parameters to four. Given the wide ranges of T and H covered by the RTM data, the remaining unknowns can be determined with an unprecedentedly high level of confidence by using k ab and k c . Figure 4 displays the magnetotropic coefficients converted from the measured ∆f data in both geometries.  It is clear that the fits capture the full field-dependence of k ab (H) and k c (H) in an excellent manner, despite the minimal spin model and the mean-field approximation.
In particular, the low-field H 2 curvature is fully consistent with |∆χ| (Sec. IVA). At µ 0 H > 30 T, beyond the range of some of the data, the agreement is no longer quantitative for all temperatures simultaneously, but the model continues to capture the majority of the fielddependent behavior for both field directions. One possible source of these deviations would be high-field magnetostrictive effects, which can distort the lattice structure and thus modify both the magnetic interactions and possibly even the sizes of the CEF coefficients at sufficiently large applied fields [50]. We use the optimal fits to obtain the six Stevens-operator coefficients [Eq. (1)] shown in Table 2 and we discuss the physical implications of having this fully determined spin Hamiltonian [Eq. (5)] in the next section.

A. CEF spectrum and VV coefficients
Having determined a full set of eight fitting parameters from the RTM data, we first verify the self-consistency of this fit against the measured magnetic susceptibilities, which are shown as χ −1 ab (T ) and χ −1 c (T ) over the full temperature range in Fig. 5. The solid lines show the same quantities calculated from the mean-field Hamiltonian of Eq. (5). The agreement of the fit with the data is quantitatively excellent over the entire measured T range, including temperatures allowing significant population of the higher CEF levels. The dot-dashed lines were calculated using the pseudospin-1/2 approximation, i.e. the contributions from the ground-state doublet [Eq. (8)], which as expected captures the T -dependence below a specific energetic cut-off; by inspection we find this cutoff at T ≈ 60 K in CsYbSe 2 for both directions ofĤ.
As discussed in Sec. III, the difference in T -dependence for the two primary field directions is accounted for in the modelling procedure by the large discrepancy in the VV coefficients. While this anisotropy is also reflected in the very different prefactors of the H-quadratic part of the energy spectrum for the two different field geometries, we stress again the fact that it affects the susceptibility strongly even at zero field. As a consequence, the susceptibility fits in Fig. 2 provide an independent determination of α 0 ⊥ and α 0 z , as well as of the squared matrix elements | |0 ± |Ĵ x(z) |0 ∓(±) | 2 . These constitute additional constraints on the allowed Stevens-operator coefficients that are crucial for reducing the enormous degeneracy of the six-variable parameter space describing the CEF levels. It is well known that suitable sets of coefficients are often highly degenerate, in the sense of yielding many indistinguishable fits of INS data for the CEF spectrum and g-tensor values [18,28,37]. In Table 1, we compare the ground-state g-tensor values and VV coefficients obtained by fitting k ab,c (H) (top line, Fig. 4) and χ −1 (second line, Fig. 5), where we find agreement at the 10% level.
Next we use our fit coefficients to calculate the CEF energy levels, which are shown in Fig. 6 as a function of field at zero temperature. The black solid lines show the spectra calculated, forĤ oriented in the ab-plane [panel (a)] and along the c-axis [panel (b)], by the direct diagonalization ofĤ CEF +Ĥ Z . Focusing first on zero field, we obtain the level-splittings ∆ 10 = 13.6 meV, ∆ 20 = 29.6 meV, and ∆ 30 = 52.6 meV. Thus ∆ 10 in CsYbSe 2 is similar to that of NaYbSe 2 (≈ 17.7 meV [29]), but is significantly lower than the values reported for delafossites in which the Yb 3+ ion has an oxygen environment [28,37]. The g-factor values we obtain for the ground state, g ab = 3.52 and g c = 1.33, are comparable to those measured in all of the Yb-based delafossites [30], as summarized in Table 1.
Turning to finite fields, a comparison of the groundstate doublet (n = 0) in Figs. 6(a) and 6(b) shows the expected strong anisotropy, which lies predominantly in the fact that the H-quadratic component is much larger for H c (the aforementioned factor-6 difference between α 0 z and α 0 ⊥ ). Because the quadratic curvatures are a consequence of repulsion between adjacent CEF levels, the n = 1 doublet exhibits a clear opposing curvature for H c (dominating the behavior of both doublet components), whereas in the H ab configuration this effect is weaker than the higher-order corrections. The dotdashed lines in Fig. 6 illustrate the fidelity of a fit made only at the level of second-order energy corrections for each n [Eq. (6)], i.e. by using all the VV coefficients, which agrees well until µ 0 H > 50 T in bothĤ directions.
Moving up in the CEF spectrum, our results also capture the unique properties of the n = 2 state. This dipolar-octupolar doublet, |J, m J = ±3/2 [20,49,51,52], does not transform as a magnetic dipole because the threefold rotational symmetry excludes any mixing between states with m J = ±3/2 and those with other values of m J [53]. Instead this doublet combines the features of dipolar and octupolar moments [20,49]. In particular, only one component of the dipole moment (in this case oriented along the c-axis) appears in the vector of pseudospin-1/2 operators, which makes the Zeeman splitting of the n = 2 level appear very different between H ab and H c. In the latter case, the m J = 3/2 state remains an exact eigenstate up to arbitrarily large field, so the Zeeman splitting is exactly H-linear with no higher-order perturbative corrections [ Fig. 6(b)], i.e. α 2 z = 0. By contrast, for H ab the matrix elements 2 ± |J x |2 ∓ are zero and the H-linear Zeeman splitting vanishes (g 2 ⊥ = 0). Thus these doublet components remain nearly degenerate [ Fig. 6(a)] until the field is sufficiently strong that higher-order perturbative corrections become visible.
A further key experimental quantity we consider is the magnetization. The solid lines in Fig. 7 show the magnetization of the system at several different temperatures calculated for fields in the two primary directions using the full Hamiltonian matrix [Eq. (5)]. In the low-field range, m ab and m c reproduce exactly our experimental data measured up to µ 0 H = 7 T. At higher fields, m ab calculated from the full spectrum changes slope at µ 0 H ≈ 12 T (at T = 4 K), beyond which it continues to increase more slowly with increasing field. This hint of saturation behavior is consistent with the broad maximum in the RTM frequency shift, ∆f ab (H), at low temperatures (Fig. 3). No such behavior is found in m c , although the overall slope does decrease as H is raised beyond approximately 25 T, and nor is it immediately evident in the RTM measurements; however, it can be found within the detailed H-dependence of the free en-ergy, which we parameterize using the transverse susceptibilities, χ T ab (T ) and χ T c (T ), in App. D. Finally, the magnetization presents an excellent test case for the validity of a pseudospin-1/2 description of CsYbSe 2 , meaning the extent to which any physical property is explained by the contributions from the groundstate doublet (n = 0) alone. In Figs. 7(a) and 7(b) we show in addition the magnetizations obtained using a pseudospin-1/2 approximation with (dotted lines) and without (dashed lines) the VV correction. For reference, the efficacy of the VV corrections in reproducing the full field-induced evolution of the ground doublet can be gauged in Fig. 6. For m ab (H), a pseudospin-1/2 model with no VV correction captures the field-dependence only up to µ 0 H ≈ 10 T before saturating to the generic tanh function of a Zeeman doublet. By contrast, a model with VV correction follows the full solution very closely over the entire field range to 60 T. This contrasts with the situation for m c (H), where even the VV-corrected model shows a clear departure from the full solution that sets in around 15 T, beyond which the approximate treatment separates systematically, and it is clear that higher-order corrections to the doublet spectrum are required. For the model with no VV correction, m c (H) does not even reproduce the low-field regime. Recalling the very significant magnetic anisotropy (α 0 z /α 0 ⊥ ≈ 6), the effectiveness of the VV corrections for the two directions is no surprise. As a function of temperature, we observe that the VV-corrected magnetizations in both geometries show increasing discrepancies at the same field as T is increased to 70 K (Fig. 7), even though the nominal thermal cutoff imposed by the next CEF level is ∆ 10 /k B ≈ 160 K.
Thus our experiments and analysis demonstrate two key results. The first is that all of the nontrivial magnetic properties of a material are captured by a single concept, the VV coefficients. The second is that the contrast between conventional and nontrivial magnetic properties can be found in a single material due to anisotropic VV coefficients that differ strongly between the field geometries H ab and H c. We have shown that the origin of the VV coefficients lies in the large matrix elements governing the mixing and consequent repulsion of adjacent CEF levels. Although the ground-state VV coefficients, α 0 ⊥ and α 0 z , are zero-field quantities whose effects can be found in non-CW behavior of the susceptibility, their strongest impact emerges as H increases (Fig. 6). Our analysis also demonstrates that the success or failure of the popular pseudospin-1/2 approximation to the magnetic properties of a Kramers-doublet system depends directly on the magnitude of the VV coefficients.
Our conclusions rely completely on determining all of the CEF parameters with high accuracy, which allows a detailed characterization of all the magnetic properties throughout the (H, θ, T ) parameter space within a single and self-consistent model. The simplest indication for a significant contribution from the VV coefficients is a deviation of the susceptibility from a CW form, which is visible most clearly as a downward curvature in χ −1 (T ) as T → 0; this type of behavior is obvious in χ −1 c (T ) in Fig. 5, but not in χ −1 ab (T ). As noted above, a more detailed characterization of the VV coefficients requires experiments under significant applied fields. Quantitatively, the origin of the VV coefficients in second-order perturbation theory [Eq. (7)] gives them a systematic dependence on ∆ −1 n0 , as well as on the expectation values of J x(z) . Although the latter terms are primarily responsible for the strong directional anisotropy in the Yb-based delafossites whose reported ∆ 10 values and extracted VV coefficients are collected in Table 1, the related material YbMgGaO 4 presents an example where the large value of ∆ 10 suppresses the VV coefficients, leading to a predominantly CW-type behavior of χ −1 (T ) [34].
We stress again that our full microscopic model is directly applicable to the computation of all aspects of the magnetic response of a material. Here we have illustrated the situation for our own magnetic susceptibility, magnetization, and RTM data, and we await its extension to describe magnetic specific-heat, torque, ESR, and magnetocalorimetric measurements. We have also distilled the properties of the full model to the physically relevant VV coefficients, which can be understood as governing the nonlinear field-dependence of the CEF energy levels (Fig. 6) and thus affect all of the magnetic properties. This linking function allows the use of our analysis to resolve a number of contradictions that have emerged in recent studies of delafossites by different techniques.
One example is the report from single-crystal INS measurements on CsYbSe 2 of no CEF levels below 20 meV [43], in direct contradiction to the present conclusions (Fig. 6). We have calculated the scattering cross-section for the 0 → 1 transition by using the diagonalized CEF matrix in order to compare the intensity obtained in the crystalline geometry of the measurement with that expected for a polycrystalline sample. Indeed we find the former to be smaller than the latter by four orders of magnitude, indicating that the initial conclusion was the consequence of an unfortunate choice of geometry, and subsequent reports suggest that a CEF level has been found around 15 meV [54]. We note also that the two sets of Stevens-operator coefficients proposed [21] for NaYbO 2 based on INS measurements of the CEF spectrum have widely divergent VV coefficients, whereas similar fits to the INS spectra of NaYbSe 2 show much more agreement when refined with complementary Raman scattering data [25].

B. Magnetic interaction parameters and triangular-lattice spin models
We turn our discussion from the Stevens operators determining the CEF levels to the magnetic interaction parameters that govern the low-energy physics, and hence the extent to which the Kramers-doublet system can be used as an effective realization of any of the paradigm S = 1/2 models in quantum magnetism. As Sec. VA made clear, an adequate understanding of the CEF energies and their evolution in an applied field is a prerequisite in the search for phenomena including field-induced phase transitions and candidate QSL phases [1]. In particular, narrowly spaced CEF levels that undergo significant mutual repulsion at a specific field scale offer a complex and correlated energy landscape that could accommodate unconventional spin states.
Focusing on the triangular geometry, the nearestneighbor triangular-lattice Heisenberg model is the original [55] and still one of the deepest problems in frustrated quantum magnetism [56,57]. Although the ground state of this model has a modest amount of magnetic order, triangular lattices have attracted extensive interest from a number of angles over the decades. Not only does this geometry have multiple materials realizations, but each generation of materials has opened a new dimension in research [56]. Triangular organic compounds drove a discussion of proximity to the Mott transition [5][6][7][8], Cs 2 CuCl 4 and Cs 2 CuBr 4 [58][59][60][61][62] drove studies of the spatially anisotropic (J-J ′ ) triangular-lattice models [57,63], and the cobaltates Ba 3 CoNb 2 O 9 [64,65], Ba 3 CoSb 2 O 9 [66][67][68], and Ba 8 CoNb 6 O 24 [69] spurred the consideration of spin-anisotropic triangular lattices with XXZ symmetry. In recent years, Yb-based triangularlattice materials have sparked very strong interest in further spin anisotropies, in the form of J ++ and J +z terms [38,48,70], all of which widen considerably the scope for finding QSL states. On the theoretical side, it has been shown that second-neighbor Heisenberg interactions also drive a QSL state, whose gapped or gapless nature remains undetermined at present [71][72][73][74][75]. In the nearest-neighbor model, the systematic treatment of parton-based formulations has led to qualitative advances in calculating the dynamical excitation spectrum by both Schwinger-boson [76] and pseudofermion methods [73,77]. Long-standing questions about the thermodynamic properties may soon be answered by DMRG methods [78], despite the constraints of working on a rather narrow cylinder, and by tensor-network methods [79] despite the challenge posed by the high connectivity of the triangular lattice.
In Sec. II we restricted our considerations to a triangular-lattice model of XXZ form [Eq. (2)], meaning that we allowed only a minimal spin anisotropy of Ising or XY form. Working at the mean-field level [Eq. (5)], in Sec. IIIB we obtained the results J ⊥ = 0.54 ± 0.01 K and J z = 0.61 ± 0.01 K for the magnetic interactions in the system of J = 7/2 Yb 3+ ions. As shown in App. C, the interaction terms of the corresponding effective pseudospin-1/2 model are J ′ ⊥ = 5.12 K (≃ 0.44 meV) and J ′ z = 0.84 K (≃ 0.07 meV), meaning that CsYbSe 2 is very strongly in the XY limit (J ′ z /J ′ ⊥ ≃ 0. 16). Although this contrasts with most assumptions made at present, we note that XY character has also been deduced for each of the NaYbX 2 materials by using a pseudospin-1/2 treatment [30]. We expect further investigations of CsYbSe 2 to confirm this result. As the most straightforward indicator of XY or Ising physics, we suggest that the width in field of the regime over which the 1/3 plateau (i.e. the state of up-up-down spin order) in the magnetization is stabilized constitutes a quantity rather sensitive to the ratio J ′ z /J ′ ⊥ [47,80]. This plateau is evident in the data of Ref. [43], although a lower temperature would assist a quantitative analysis. However, the relevant regime of parameter space is yet to be investigated theoretically for a field applied in the ab plane of the system. Finally, while we cannot exclude terms of J ++ and J +z type from the spin Hamiltonian, the accuracy of our fit indicates that the effect of any missing terms is extremely small in CsYbSe 2 ; either they are genuinely small or they are relevant only at the lowest temperatures where a treatment beyond the present meanfield level would be required.
On the materials side, it is also fair to say that the continuing lack of experimental signatures for QSL ground states can be blamed on two primary issues, namely the scarcity of candidate materials and the paucity of measurable physical quantities offering unambiguous signatures or predictors for QSL properties (such as fractional quantum numbers and nonlocal entanglement). To address the first of these, our analysis provides a definitive guide to the field-induced physics of highly spin-anisotropic systems such as the Yb-delafossites, and hence to the regions of parameter space where competing energy scales establish an environment conducive to the occurrence of exotic spin states. Although we cannot solve the second issue, we can provide a comprehensive understanding of the single-ion energy spectrum that identifies the extent to which a given material replicates a target pseudospin-1/2 model, thereby streamlining the experimental search for QSL fingerprints by available experimental methods.

VI. SUMMARY
We have investigated the anisotropic magnetic response of an insulating 4f electronic system by measuring two key thermodynamic quantities, the magnetic susceptibility in the low-field limit and the magnetotropic coefficients over very wide field and temperature ranges (up to µ 0 H = 60 T and T = 70 K). We have shown that the anisotropies in both quantities can be formulated within a set of anisotropic van Vleck (VV) coefficients, which arise as the second-order perturbative corrections of the Zeeman interaction to the zero-field crystal electric field (CEF) spectrum. This leads to the essential finding that the VV coefficients constitute independent physical quantities that describe the crucial magnetic properties of 4f spin systems across the full range of applied fields and extant anisotropies. A proper account of the groundstate VV coefficients is indispensable for an accurate and unambiguous determination of the microscopic parameters governing the CEF Hamiltonian, a process for which otherwise few routes are known to date. The VV coef-ficients fulfill the vital function of unifying the low-field, low-temperature magnetic susceptibility with the highfield magnetotropic coefficients and CEF levels, and in this sense their role as stand-alone physical quantities allowing a full interpretation of magnetic anisotropies has not been appreciated before.
Our experimental results highlight the value of the resonant torsion magnetometry (RTM) method, which is accurate and profoundly powerful in terms of the parameter ranges it accesses. The magnetotropic coefficients we extract over these broad field and temperature ranges play the key role in obtaining a unique set of Stevens operators describing the microscopic CEF Hamiltonian with unprecedent fidelity. We reiterate that a fitting analysis must provide complete consistency from zero to high field and at all relevant temperatures, and our fits meet this challenge. With the full CEF spectrum in hand, we can examine the validity of different and popular approximations that have been applied to many materials. Specifically, we identify the limits of a CW fit to the temperature-dependence of the magnetic susceptibility and the boundaries of the effective pseudospin-1/2 description for systems with a ground-state Kramers doublet.
We have developed and applied our analysis for the material CsYbSe 2 , which is a member of a family of Ybdelafossites displaying triangular-lattice geometry. Because the CEF levels of the Yb 3+ ion are four Kramers doublets, these compounds are leading candidates in the search for quantum spin-liquid (QSL) behavior, and indeed the full CEF spectrum we obtain up to high fields [ Fig. 6] reveals an intricate and anisotropic energy landscape amenable to unconventional magnetism. This spectrum allows one to construct a maximally informed pseudospin-1/2 model for the low-energy physics of the system, and within a minimal XXZ spin Hamiltonian we conclude that CsYbSe 2 is a strongly XY triangular antiferromagnet. While we await further experimental confirmation of this result, we note again that our analysis is applicable to a wide range of 4f materials with complex CEF spectra and especially with ground-state doublets allowing an effective spin-1/2 description, which should expand significantly the scope of the search for QSL phases. Appendix A: Exact form of magnetic susceptibility From the N -particle partition function calculated with Eq. (5), we derive analytical formulas for the low-field magnetic susceptibilities of our model using the thermodynamic relation where β = 1/k B T , Θ CW ⊥(z) = qJ ⊥(z) /k B as in Eq. (8), and α n i are the VV coefficients defined in Eq. (7). The form of Eq. (A1) is exact within the mean-field approximation, and by considering only n = 0 (the ground-state doublet) it reduces to Eq. (8).
Appendix B: Comparison to the Curie-Weiss magnetic susceptibility of a spin-1/2 system To illustrate how the VV coefficients alter even the low-T limit of the full susceptibility expression, where only the lowest Kramers doublet is relevant, we repeat Eq. (8) in the form Only when the VV term is negligible does this expression reduce to the familiar CW form, with a CW temperature for the effective S = 1/2 model given by qJ ′ ⊥(z) /4k B , where the parameters J ′ ⊥ and J ′ z are defined in App. C below.
Appendix C: Effective spin-1/2 projection At the lowest temperatures, where the physics is dominated by the properties of the lowest Kramers doublet, {|0 ± }, it is helpful to formulate a pseudospin-1/2 model in terms of effective spin-1/2 operators with the action The equivalent pseudospin-1/2 model for a Hamiltonian H is obtained formally by projecting onto the lowest Kramers-doublet subspace using the projection operator from which it follows that the equivalent g-factors and exchange constants of the pseudospin-1/2 model are related to J = 7/2 operator matrix elements and exchange constants by the expressions Applying these relations using g J = 8/7 gives the gfactors g ⊥ = 3.77 and g z = 1.76 shown in the second line of Table 2, and using the results J ⊥ = 0.54 ± 0.01 K and J z = 0.61 ± 0.01 K obtained in Sec. IIIB for the J = 7/2 system gives the values J ′ ⊥ = 5.12 K (≃ 0.44 meV) and J ′ z = 0.84 K (≃ 0.07 meV) for the interaction parameters of the pseudospin-1/2 model.