Honeycomb rare-earth magnets with anisotropic exchange interactions

Zhu-Xi Luo and Gang Chen3,4∗ Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106 Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84102 Department of Physics and HKU-UCAS Joint Institute for Theoretical and Computational Physics at Hong Kong, The University of Hong Kong, Hong Kong, China and State Key Laboratory of Surface Physics and Department of Physics, Fudan University, Shanghai 200433, China (Dated: May 14, 2020)


Introduction
Spin liquid candidates are often being searched among geometrically frustrated systems, such as triangular [1], kagomé [2] or pyrochlore [3] lattices. This is quite reasonable as the geometrical frustration could lead to a large number of degenerate or nearly-degenerate classical ground states for commonly studied Heisenberg models and thus enhance quantum fluctuations when the quantum effects are included. However, the destabilization of simple magnetically ordered states and driving a disordered one can happen even on unfrustrated lattices, by exploiting the power of anisotropic interactions [4]; the Kitaev honeycomb model [5] is a representative example of the latter. Besides being an academic interest, anisotropic spin interactions are also inevitable in realistic magnetic materials, especially those with heavy atoms. A large number of spin liquid candidates are known experimentally to possess a significant spin-orbit coupling, leading to rather anisotropic spin interactions [6][7][8][9][10][11]. Beyond the current interest in the spin liquid physics, understanding the relationship between the magnetic properties and the anisotropic spin interactions is a frontier topic in the field of quantum magnetism.
The most commonly studied anisotropic magnets on unfrustrated lattices are the 4d/5d magnets [11,12] that include the honeycomb iridates and RuCl 3 . Because of the possible relevance to Kitaev physics, these materials were referred as Kitaev materials. Due to the spatial extension of the 4d/5d electron wavefunctions, the exchange interactions between the local moments are usually beyond the nearest neighbors. Moreover, the iridates often suffer from a strong neutron absorption such that the data-rich neutron scattering measurement can be difficult. In comparison, the rare-earth family has the advantages of much stronger spinorbit couplings and much more localized 4 f orbitals [13][14][15], and the exchange interactions often restrict to first neighbors. This makes the understanding of the modeling Hamiltonian more accessible. In addition, the rare-earth magnets do not have the neutron absorption issue that prevails in iridates. Furthermore, their smaller energy scales allow for the possibility to quantitatively understand their Hamiltonian through the external magnetic fields. However, rare-earth materials have only been well-investigated on frustrated lattices [9,16].
In this paper, we will study the rare-earth magnets on the unfrustrated honeycomb lattice, and pursue an understanding of the experimental consequence of the spin-orbital entanglement on the honeycomb structure. We start by exploring the thermodynamic properties of a generic model with the nearest neighbor interactions. It is well-known that the anisotropic exchange couplings could appear in the temperature dependence of the thermodynamic quantities such as the specific heat, spin susceptibility [17] and magnetotropic coefficients [18,19]. Especially for the spin susceptibility and magnetic torque, magnetic fields along different directions induce magnetization of different magnitudes, leading to the anisotropic spin susceptibility [20][21][22][23][24][25] and the angular dependence of the magnetic torque [26][27][28][29], and providing a natural detection of the intrinsic spin anisotropy in the system. To go beyond the thermody-namic properties, we further consider the electron spin resonance (ESR) measurement [30,31] of the system. The ESR measurement turns out to be a very sensitive probe of the magnetic anisotropy and is especially useful for the study of the strong spin-orbit-coupled quantum materials, and we compute the ESR linewidth to reveal the intrinsic spin anisotropy of the spin interactions.
Due to the small energy scale of the interaction between the rare-earth local moment, it is ready to apply a small magnetic field in the laboratory to change the magnetic state into a fully polarized one. For such a simple product state, the magnetic excitation can be readily worked out from the linear spin wave theory. We further consider the spin wave spectrum and explore the possibility of topological magnons [32][33][34][35][36][37][38][39]. We find the magnon spectrum supports nontrivial topological band structure. This feature can be manifested in thermal Hall transport measurements.
The remaining parts of the paper are organized as follows. In Sec. 2, we introduce the nearest-neighbor spin Hamiltonian, followed by the high-temperature analysis of heat capacity, spin susceptibilities and magnetic torque coefficient in Sec 3. Then we consider the ESR and calculate the influence of anisotropy on the ESR linewidth in Sec. 4. Next the linear spin wave theory of the system is exploited under strong external fields in Sec. 5 and the aspect of topological magnons is discussed. Finally in Sec. 6 we comment on possible materials YbCl 3 and TlYbS 2 .

Model
We begin with the following microscopic spin model, that is the most general nearest neighbor Hamiltonian on a honeycomb lattice with the (usual) Kramers doublet effective spin-1/2 local moments [6,10,15,40], with γ i j taking e 2iπ/3 , e −2iπ/3 , and 1 on the bonds along a 1 , a 2 , a 3 directions respectively, as shown in Fig. 1. The first two terms give the usual XXZ Hamiltonian, while the latter two terms are bond-dependent and constitute the spin-orbit interaction. The spin components are defined in the global coordinate system in Fig. 1. This is possible because the system is planar and has an unique rotational axis. This differs from the rare-earth pyrochlore materials where the spins are often defined in the local coordinate system for each sublattice. This model applies to the rare-earth local moment such as the Yb 3+ ion. For non-Kramers doublet like Pr 3+ or Tb 3+ ion, the J ±z term is not allowed by symmetry, and the model becomes further simplified. In fact, a non-Kramers doublet based rare-earth honeycomb magnet arises from the triangular lattice magnet TbInO 3 after 1/3 of the Tb 3+ ions becomes inactive magnetically [41]. For the rare-earth local moments, the 4 f electrons are much localized, and most often, one only needs to consider the nearest-neighbor interactions, and occasionally, one would like to include the further neighbor dipole-dipole interactions. In contrast, for the 4d/5d systems, one may need to worry about further neighbor exchange interactions because of the large spatial extension of the electron wavefunctions. An alternative and often used parametrization of the Hamiltonian is that of the J-K-Γ -Γ model [42]: where α, β, γ take values in {x , y , z }. In the latter coordinate system, our unit vectors of Fig. 1 can be expressed byx = (−1, −1, 2)/ 6,ŷ = (1, −1, 0)/ 2 andẑ = (1, 1, 1)/ 3. The spin components in the above equation are

Thermodynamics
The highly anisotropic nature of the exchange interaction first impacts the thermodynamic properties of the system. Here we explicitly calculate the specific heat and the magnetic susceptibilities of the system from the generic exchange Hamiltonian, with details presented in Appendix B. Using the high-temperature series expansion [43,44], we find the heat capacity  to be where k B is the Boltzmann constant, and Due to the spin-orbit entanglement, the coupling of the local moment to the external magnetic field is also anisotropic. The Landé factors g's are different for the in-plane and out-plane magnetic fields, and the Zeeman coupling is given as h x , h y are the in-plane components of the external magnetic field, and h the z-component. µ 0 is the vacuum permeability and µ B the Bohr magneton. Again using high-temperature series expansion, we compute the parallel and perpendicular spin susceptibilities up to O(T −3 ) In the SU(2)-symmetric point, J zz = 2J ± , J ±± = J ±z = 0, the two expressions coincide. For the rare-earth local moments with non-Kramers doublets, g ⊥ = 0 so χ ⊥ = 0. In Fig. 2, we plot the magnetic susceptibilities and show the deviation from the simple Curie-Weiss law due to the high order anisotropic terms. In addition to the simple thermodynamics such as C v and χ, the magnetic torque measurement is proved to be quite useful in revealing the magnetic anisotropy. Intrinsically, this is because the induced magnetization is generically not parallel to the magnetic field. Thus, when the sample has an anisotropic magnetization, the system would experience a torque τ = M × H = −∂ F /∂ θ in an external magnetic field. The magnetotropic coefficient k = ∂ 2 F /∂ θ 2 , defined as the second derivative of the free energy to the angle θ between the sample and the applied magnetic field, can be introduced to quantify such anisotropy. It can be  (14)) on J ±± /J zz and J ±z /J zz , with fixed values of J ± /J zz . Left: J ± /J zz = 1, middle: In the first two panels, the linewidths are minimal at the "most isotropic" points J ±± = J ±z = 0. In the right panel, we depict the case J ± /J zz < 0, where the in-plane exchange is ferromagnetic. The linewidths ∆H in the three plots are in the units µ B g(θ )/ 2π. directly measured using the resonant torsion magnetometry [18]. Under the high temperature expansion, we find the magnetotropic coefficient k is given as where we have defined h 2 = h 2 x +h 2 y +h 2 z . Note that there are two different sources of a nontrivial torque magnetometry: the anisotropy of the g-tensor and the anisotropy of the exchange. In the limit of g ⊥ = g , there is still a nonzero contribution to k, detailed in Appendix B. The coefficient will only vanish if the Heisenberg limit is further taken: J zz = 2J ± , J ±± = J ±z = 0.

Electron spin resonance
In the thermodynamic properties, the leading contributions come from the J zz and J ± terms, while the J ±± and J ±z terms are subleading. Arising from spin-orbital entanglement and completely breaking the U(1) rotational symmetry, these terms play important roles in the potential quantum spin liquid behavior. To resolve them, we now turn to the electron spin resonance.
Electron spin resonance measures the absorption of electromagnetic radiation by a sample subjected to an external static magnetic field. For a SU(2) invariant system, the absorption is completely sharp, i.e. described by a delta function located exactly at the Zeeman energy [30]. Therefore, the broadening of the resonance spectrum has to arise from the magnetic anisotropy. To understand the contribution of the anisotropy of the nearest-neighbor spin interaction to the ESR linewidth, we decompose the Hamiltonian Eq. (1) into the isotropic Heisenberg part and the anisotropic exchange part where the Heisenberg coupling J = (J zz + 4J ± )/3, and the anisotropic part Here Γ i j is a traceless and symmetric exchange coupling matrix, satisfying Under the Zeeman term of Eq. (8), the ESR can be computed using the Kubo-Tomita formalism [45], yielding a Lorentzian-shaped spectrum. The corresponding linewidth at high temperatures is given by the second and the fourth moments of the ESR line-shape function [46][47][48] where θ is again the angle between the external field and the sample, and M 2 and M 4 are the second and the fourth moments, respectively, and M ± ≡ i S ± i . The expectation "〈· · · 〉"in the above equations is taken with respect to high temperatures. Specifically, we find that The high-temperature ESR linewidths (14), expressed in terms of the different exchanges (17), can be compared with future experiments on the rare-earth based honeycomb magnets in order to extract the anisotropic exchanges. The J zz and J ± exchanges are easier to indicate from analyzing the experimental data of susceptibilities, as they have lower-order effects. Using their extracted values, one can infer the corresponding J ±± and J ±z from the ESR linewidth ∆H. In Fig. 3, we depict the three-dimensional plots that explicitly demonstrate the dependence of the ESR linewidth on the anisotropic couplings J z± and J ±± for three different choices of J ± .

Strong field normal to the honeycomb plane
To further explore the effect of the anisotropic exchange interaction, we study the spin wave excitation with respect to the polarized states under the strong magnetic fields. This is clearly feasible in the current laboratory setting for the rare-earth magnets as the energy scales for them are usually rather small. For the 4d/5d magnets, there can be difficulty to achieve as the energy scale over there is much higher, typically by two orders. Our results here are relevant to the inelastic neutron scattering and thermal Hall transport measurements. We first consider the case of a strong magnetic field in the direction normal to the honeycomb plane such that the system is in the fully polarized paramagnetic phase and all the spins are aligned along the z direction. In this case, the magnon bands carry nontrivial Chern numbers for generic range of parameters, as found in reference [39] in the J −K−Γ −Γ presentation. We expand about this fully polarized state using the conventional Holstein-Primakoff transformations of the spin variables [49], which are Keeping only the bilinear terms of bosonic operators and taking the Fourier transformation, we arrive at and N being the number of sites in one sublattice. Here we have denoted k 1 = − 1 2 k x + 3 2 k y , k 2 = − 1 2 k x − 3 2 k y , and k 3 = k x that correspond to the y-, zand x-bonds, respectively. We further define f (k) = i e ik i , g 1 (k) = i e −ik i γ i , g 2 (k) = i e ik i γ i and u ≡ (g µ 0 µ B h − 3J zz /2), we then have for H k a block form where we have All the J ±z terms are not present. The spin wave dispersion relation for H k follows as where only the positive square root of ε(k) 2 is taken. Several simple limits of this expression can be checked: (1) in the Heisenberg limit J zz = 2J ± ≡ 2J, then ε(k) = ( when only J zz is finite, it reduces to the Ising case ε(k) = The two spin wave bands ω ± when a strong field in the z-direction is applied.
At high fields, the results can be simplified by the Schrieffer-Wolff transformation, with the commutator understood as and η is a diagonal matrix with entries (1, 1, −1, −1). Following the treatment of Ref. [39], we choose the transformation to be so that up to O(h −2 ), we have theH k to become A(k) →Ã(k), B(k) →B(k), At high fields, we can thus ignoreB(k) and focus on thẽ A(k) term. WritingÃ(k) = d 0 (k)1 + 1 2 d(k) · σ, with the three components being At each momentum k we have the eigenvalues The above spin wave bands Eq. (27) do not touch unless J ± = J ±± = 0, as we have depicted in Fig. 4. We further compute the Berry curvature as follows This is negative semi-definite in the Brillouin zone. The Chern numbers follow as This implies the presence of chiral magnon edge states and thermal Hall effect, resulting from the presence of magnon number non-conserving terms B(k) in the Hamiltonian [34,39]. The edge state for the open boundary condition is depicted in Fig. 5.

Strong field along the honeycomb plane
We now turn to a strong in-plane field, in the x-direction. This is relevant for the rare-earth local moments with the usual Kramers doublet, and does not apply to the non-Kramers doublet. The Holstein-Primakoff transformation for sublattice A is modified as , and that for sublattice B is obtained by substituting a by b. These are again bosonic operators satisfying

Keeping only the bilinear terms of bosonic operators and taking the Fourier transformation, we obtain
Define g 3 (k) = e ik 1 + e ik 2 − 2e ik 3 , g 4 (k) = e ik 1 − e ik 2 and recall f (k) = i e ik i . The H k is of the familiar form , but now with the A, B matrices given by Appealing again to the Schrieffer-Wolff transformation with then up to O(h −2 ⊥ ), we have the effectiveH k to be A(k) →Ã(k), B(k) →B(k). At high fields h ⊥ , we can ignoreB(k) and focus on theÃ(k) term. Rewritẽ A(k) = d 0 (k)1 + 1 2 d(k) · σ, with each component being We then arrive at the dispersions ω ± (k) = d 0 (k) ± 1 2 |d(k)|. The spectrum is plotted in Fig. 6. We find both bands have zero Chern numbers, and we have checked for many other parameter choices and also obtained trivial zero Chern number. Thus, the in-plane field magnon band structure is quite distinct from the topological magnon band structure for the normal-plane field case.

Discussion
We have studied the experimental consequences of the spin-orbital entanglement and the anisotropic spin exchange interactions in the honeycomb rare-earth magnets. These results can be directly compared with the experiments, thereby providing a useful guidance for the future study on candidate systems. One future direction would be to involve higher-lying crystal field states based on the information of specific materials.
One potential rare-earth candidate for the anisotropic honeycomb lattice model is YbCl 3 [50], which has a similar crystal structure to that of RuCl 3 . The Yb 3+ ions have nearly filled 4 f -orbitals, which, combined with the large crystal fields lead to Kramers doublet ground state manifold. This is modeled as an effective spin-1/2 local moment. Furthermore, its edge-shared octahedral structure gives simple exchange physics that is relatively well-understood according to a microscopic calculation in Ref. [15]. There is very limited information about this material in the literature apart from a very recent work [51]. Comparing their susceptibility data with our calculations, we are able to pin down the exchange parameters J zz ∼ 8K, J ± ∼ 6K. Another relevant material is TlYbS 2 with AB-stacking triangular structure, which is equivalent to the honeycomb lattice. Recent experiment [52] shows that it has no long-range magnetic order down to 0.4K, suggesting it to be a frustrated magnet and possible candidate for spin liquid.
In this paper, we have focused the analysis on the honeycomb lattice rare-earth magnets and its anisotropic interaction. It is noticed that, the generic model for the rare-earth honey-comb magnets contains a Kitaev interaction as one independent exchange interaction out of four. It is thus reasonable for us to consider the possibility of Kitaev materials among the honeycomb rare-earth magnets. In fact most rare-earth magnets have not been discussed along the line of Kitaev interactions, except the first few works [13][14][15]. In the previous work [13], we have illustrated this observation with the FCC rare-earth magnets. Since many non-honeycomb lattice iridates are claimed as Kitaev materials, it is thus reasonable to consider the rare-earth magnets with other crystal structures to be potential Kitaev materials beyond the previously proposed ones and the honeycomb one here [53,54]. The reason that these rare-earth magnets contain a Kitaev interaction is due to two facts. The first fact is the spin-orbital-entangled effective spin-1/2 local moment. The second fact is the three-fold rotation symmetry at the lattice site. This symmetry permutes the effective spin components and generates a Kitaev interaction. These two ingredients can be used as the recipe to search for other rare-earth Kitaev materials beyond the honeycomb one.
To summarize, we have focused on rare-earth honeycomb materials with nearest-neighbor interactions and computed the high-temperature thermodynamic properties, ESR linewidth, and spin-wave behaviors as the experimental consequences of the anisotropic spin interaction.

A Generic spin models and candidate states for higher spins
This spin model in Eq. (1) is designed for effective spin-1/2 local moments. It can be well extended to the high-spin local moments. For the honeycomb lattice with spin-1 local moments, the pairwise spin interaction is given as Because of the larger Hilbert space, a single-ion anisotropy (the D-term) is allowed and new states such as the quantum paramagnet can be favored here. Thus, the phase transitions between quantum paramagnet and other ordered phases can be interesting. Further neighbor exchange interaction, if included, could bring more frustration channel than the spin-orbit entanglement induced frustration. It is known that, simple J 1 -J 2 (first neighbor and second neighbor Heisenberg) model on honeycomb lattice could induce spiral spin liquids in two dimensions where the spiral degeneracy has a line degeneracy in the momentum space rather than the surface degeneracy. The presence of the anisotropic interaction in Eq. (36) would overcome the quantum/classical order by disorder effect and lift the degeneracy. In addition to the spin-1 local moments, the model in Eq. (36) also applies to the spin-3/2 systems. Since the honeycomb lattice contains three nearest neighbor bonds, one may consider the possibility of the AKLT states on the honeycomb lattice where the nearest neighbor bonds are covered with spin singles of the spin-1/2 states and three onsite spin-1/2 spins are combined back to a spin-3/2 local moments.
Here, we have only listed the pairwise spin interactions. Due to the spin-orbital entanglement and the spin-lattice coupling, the effective interaction for the spin-1 and spin-3/2 magnets can contain significant multipolar interactions. A simple example would be the biquadratic exchange −(S i · S j ) 2 that is induced effectively by the spin-lattice coupling. The presence of these multipolar interactions can significantly enhance quantum fluctuation by allowing the system to tunnel more effectively within the local spin Hilbert space and thus create more quantum states such as multipolar ordered phases and quantum spin liquids [55,56].

B Details of high temperature expansion
The high temperature expansion requires to take into account the commutation relations between different spin operators on a same site. To this end, we define a vertex function ν i (n x , n y , n z ) at each site i following Ref. [43], where n x , n y and n z have to be even integers for the function to be nonzero. Its explicit form can be calculated by introducing the generating function Expanding the exponential and using the definition of ν, we have ν(n x , n y , n z ) 2 n x +n y +n z ξ n x η n y ζ n z n x !n y !n z ! .
On the other hand, by diagonalizing the matrix of the exponential, we have ψ(ξ, η, ζ) = 2 cosh Expanding this and comparing with the previous equation, We note this function is symmetric under the permutation of n x , n y , n z . The heat capacity is related to the zero-field partition function in the following way where we have divided by the number of sites N to get the intensive quantity, and β = 1/k B T . Z 0 is given by where the 2 N factor results from the summation over all possible configurations. Susceptibility in direction a can be reduced to the following expectation values of two spin operators, Using the vertex function ν(n x , n y , n z ) defined above, we obtain for the parallel case, Here, the summation for {S i } is over the possible configurations of spins on all sites.
The expression above reduces to, in the limit g ⊥ = g , . (48)