Strategy to Extract Kitaev Interaction using Symmetry in Honeycomb Mott Insulators

The Kitaev spin liquid, a ground state of the bond-dependent Kitaev model in a honeycomb lattice has been a centre of attraction, since a microscopic theory to realize such an interaction in solid-state materials was discovered. A challenge in real materials though is the presence of the Heisenberg and another bond-dependent Gamma interactions detrimental to the Kitaev spin liquid, and there have been many debates on their relative strengths. Here we offer a strategy to extract the Kitaev interaction out of a full microscopic model by utilizing the symmetries of the Hamiltonian. Two tilted magnetic field directions related by a two-fold rotational symmetry generate distinct spin excitations originating from a specific combination of the Kitaev and Gamma interactions. Together with the in- and out-of-plane magnetic anisotropy, one can determine the Kitaev and Gamma interactions separately. Dynamic spin structure factors are presented to motivate future experiments. The proposed setups will advance the search for Kitaev materials.


INTRODUCTION
An electron's orbital motion in an atom generates a magnetic field which influences its spin moment, known as spin-orbit coupling. When the coupling is strong in heavy atoms, the effective Hamiltonian is described by the spin-orbit-entangled pseudospin wave-function and the interactions among magnetic ions are highly anisotropic different from the standard Heisenberg interaction [1][2][3][4][5][6] . A fascinating example is the Kitaev model with a bond-dependent interaction in a two-dimensional honeycomb lattice, whose ground state is a quantum spin liquid (QSL) with Majorana fermions and Z 2 vortex excitations 7 . There have been extensive studies on the model because in the Kitaev QSL non-Abelian excitations emerge under a magnetic field, and their braidings provide topological computation. Since a microscopic mechanism to generate such an interaction was uncovered 8 , intense efforts toward finding QSLs including a variety of candidate materials from spin S =1/2 9-18 to higher-spin S systems have been made [19][20][21][22] . Despite such efforts, a confirmed Kitaev QSL is still missing.
One challenge in finding the Kitaev QSL in magnetic materials is the presence of other spin interactions which may generate magnetic orderings or other disordered phases [23][24][25][26][27][28][29] . A generic nearest neighbour (n.n.) model in an ideal honeycomb was derived which revealed the isotropic Heisenberg interaction and another bond-dependent interaction named the Gamma (Γ) 25 . Furthermore, there exist further neighbour interactions such as second and third n.n.
Heisenberg interactions, which makes it difficult to single out the Kitaev interaction itself.
There have been many debates on the relative strengths, especially between the dominant Kitaev and Gamma interactions in Kitaev candidate materials 13,18,28,30 , and an experimental guide on how to extract the Kitaev interaction out of a full Hamiltonian is highly desirable.
In this work, we present a symmetry-based experimental strategy to determine the Kitaev interaction. Our proposal is based on the π-rotation around the a-axis perpendicular to one of the bonds in the honeycomb plane, denoted by C 2a symmetry that is broken by a specific combination of the Kitaev and Γ interactions. This broken C 2a can be easily detected with the help of a magnetic field applied within the a−c plane where the c-axis is perpendicular to the honeycomb plane; spin excitations under the two field angles of θ and −θ, measured away from the honeycomb plane as shown in Fig. 1(a), are distinct due to the combination of the Kitaev and Gamma interactions. The two field angles are related by the π-rotation around a-axis, i.e. C 2a operation. Such differences are based on the symmetry and signal the relative strengths of these interactions. A magnetic ordering that further enhances the broken C 2a symmetry does not alter the asymmetry, but quantifying the interaction strengths requires the size of the magnetic ordering. For this reason, a polarized state in the high-field region would be ideal for our purpose.
To determine each of the interactions, one needs to use the conventional in-vs. outof-plane anisotropy in spin excitations. We note that the Gamma interaction affects the conventional anisotropy, but the Kitaev does not when the field is large enough to compensate the order by disorder effect 31 . Thus subtracting the Gamma contribution deduced from the conventional anisotropy allows us to estimate the Kitaev interaction from the measured spin excitations under the field angles of θ and −θ. Both the conventional anisotropy and the π-rotation-related spin excitations can be measured by angle-dependent ferromagnetic resonance (FMR) or inelastic neutron scattering (INS) techniques while sweeping the magnetic field directions in the a − c plane containing the C 2a rotation axis.
Below we present the microscopic model and main results based on the π-rotation symmetry around a-axis. To demonstrate our theory, we also show the FMR and dynamical spin structure factors (DSSF) obtained by exact diagonalization (ED). We analyze the different spin excitations under the two field angles at finite momenta using the linear spin wave theory (LSWT), which further confirms our results based on the symmetry argument. Our results will guide a future search of Kitaev materials.

RESULT
Model -The generic spin exchange Hamiltonian among magnetic sites with strong spinorbit coupling for the ideal edge sharing octahedra environment in the octahedral x − y − z axes shown in Fig. 1(a) contains the Kitaev (K), Gamma (Γ), and Heisenberg (J) where S = 1 2 σ with ≡ 1 and σ is Pauli matrix, ij denotes the nearest neighbor (n.n.) magnetic sites, and αβ(γ) denotes the γ bond taking the α and β spin components (α, β, γ ∈ {x,y,z}). The x-, y-, and z-bonds are shown in red, blue, and green colours, respectively in Fig. 1(a). Further neighbour interactions and trigonal-distortion allowed interactions, and their effects will be discussed later.
To analyze the symmetry of the Hamiltonian, we rewrite the model in the a − b − c axes 32-34 : where φ γ = 0, 2π 3 , and 4π 3 for γ = z-, x-, and y-bond respectively, and the exchange interactions are given by The Hamiltonian H is invariant under π-rotation around the b-axis denoted by C 2b and 2π 3rotation around the c-axis by C 3c in addition to the inversion and time-reversal symmetry.
Our proposed experimental design is based on the observation that the H is not invariant under π-rotation about the a-axis C 2a due to the presence of only J ac , i.e., if J ac = 0, C 2a is also a symmetry of H. Since the C 2a is broken by J ac , if there is a way to detect the broken C 2a , that will signal the strength of J ac . We note that the magnetic field sweeping from the c-axis to a-axis within the a − c-plane does the job. The fields with angles of θ (blue line) and −θ (red line) for 0 < θ < π 2 shown in Fig. 1(b) and (c) are related by C 2a rotation, and thus measuring the spin excitation difference between these two field directions will detect the strength of J ac .
To prove our symmetry argument, we consider a full model with a magnetic field. Under a magnetic field, the total Hamiltonian including the Zeeman term is given by where the external field h has the polar angle θ measured away from the a − b honeycomb plane and the azimuthal angle φ from the a-axis as shown in Fig. 1 We focus on the lowest energy excitation n = 1 which gives a dominant resonance at low temperatures, and drop the n in ω n from now on for simplicity, even though our proposal works for all n. We define the excitation anisotropy between the magnetic field with angles of θ and −θ as δω K (θ) ≡ ω(θ) − ω(−θ) for 0 < θ < π 2 , and the conventional anisotropy between in-and out-of-plane fields as δω A ≡ ω(θ = 0) − ω(θ = π 2 ). Below we first show how δω K arises from J ac under the field in the a − c plane based on the symmetry.
Symmetry Analysis -To understand the origin of a finite δω K for φ = 0 under the magnetic field sweep, we first begin with a special case when φ = π 2 , i.e, when the external Direction of the external magnetic field h in abc axes where θ is measured from the a−b plane, and φ is from the a-axis. The blue arrow M represents the magnetic moment direction with the angle in the a − c plane is the difference in the spin excitation energies ω between two field directions: ω(θ) (blue) and ω(−θ) (red). C 2b maps ω(θ) to ω(π +θ), so δω K (π −θ) = −δω K (θ).
field is in the b − c plane. This is a special case where δω K = 0 for the following reason.
The Zeeman terms due to the field with the angle θ and with −θ are related by a π rotation of the field about theb axis, denoted by The same can be achieved by a π-rotation of the lattice, which also indicates H is invariant under C 2b . While H B breaks the C 2b symmetry of H, the total Hamiltonian H + H B (θ) and H + H B (−θ) are related by C 2b and therefore, share the same eigenenergies, i.e., δω K = 0. The difference due to the field is simply removed by a π rotation of the eigenstates about theb axis. The magnetic field sweeping from θ to −θ in the other planes equivalent to b − c plane by C 3c symmetry also gives δω K = 0. Now let us consider when the magnetic field sweeps in the a − c plane. Similarly, the magnetic field directions θ and −θ are related by Considering a π rotation of the lattice about theâ axis, we find J XY , J Z , J ab , terms are invariant under C 2a , while the J ac terms transform as magnetic field angle. Here, for simplicity, we calculate the excitation energy probed by the RF field (details can be found in the Methods) with a set magnetic field strength for spin 1 2 using ED on a C 3 -symmetric 24-site cluster.
We set our units the magnetic field h = 1 and g = µ B ≡ 1, leading to the excitation energy of a free spin, ω 0 = gµ B h = 1, so the excitation energies calculated are normalized by ω 0 . A few sets of different interaction parameters (in units of ω 0 ) are investigated. Figure   2(a) shows the J = −1 and K = Γ = 0.5 case with no δω K (θ) between −π/2 < θ < 0 (red line) and 0 < θ < π/2 (blue line), since J ac = 0. The conventional anisotropy δω A is finite, because the Γ interaction generates a strong anisotropy between the plane θ = 0 and the c-axis θ = π/2, i.e., J XY = J Z due to a finite Γ contribution. The black line is for only J = −1 showing a uniform FMR independent of angles which serves as a reference. Kitaev candidate materials 28 . Clearly, δω K (θ) is significant due to a finite J ac , and δω A is also large due to a finite Γ. While a magnetic field of strength h = 1 is used to polarize the ground state where the finite-size effect is small as shown in Supplementary Note 2, our symmetry argument works for any finite field. However, we note that the finite-size effect of ED is minimal when the ground state is polarized.  in the a − c plane. The square boxes denote the excitation energies obtained by the ED, and the colour bars indicate the intensity of DSSF α S αα (q, ω) (details can be found in the Methods). The structure factor is convolved with a Gaussian of finite width to emulate finite experimental resolution. We observe a clear difference between the two field directions, δω K at every momentum points. In particular, δω K is the largest at M 2 -point, while it is tiny at the K 1 -point. Note that M 1 and M 3 are related by the C 2b and inversion.
To gain more insights of δω K (θ) at finite momenta obtained by ED, we also perform LSWT calculations with the magnetization making an angle θ M as indicated in Fig. 1(b). with the ED results in Fig. 3(a). The mismatch between LSWT and ED is visible at every momentum, which implies the significant effects of nonlinear terms 47 .
However, when the field increases, the difference should decrease, since the magnetic polarization increases at a higher field. In Fig. 3(b), we show both ED and LSWT with h = 8 and θ M ∼ 25.8 • , where the two results match well as expected, and the nonlinear terms become less significant. In particular, the anisotropy δω K at the K-point at the high field limit given by the leading terms in 1/h, is simplified as where θ M (θ) → θ when h → ∞. This shows that both J ac and J ab should be finite for a finite δω K at the K-point, which explains no splitting of δω K at the K-point in Fig. 3(b), as our choice of parameters gives J ab = 0, i.e, Γ = −K/2. On the other hand, at the M 2 -point, there is no simple expression, but the leading terms of δω K (θ) in δθ a/c around the a-and c-axis (δθ a = 0 − θ and δθ c = θ − π/2) are given by where A and C are functions of other interactions given in Supplementary Note 3. Clearly, δω K (θ) appears as odd powers of J ac and δθ a/c , consistent with the symmetry analysis presented above.
So far, we have focused on the ideal octahedra environment. However, trigonal distortion is often present, albeit small, which introduces extra exchange interactions. Below we discuss other contributions to δω A complicating the isolation of K from J ac and our resolution of such complication in order to estimate the Kitaev interaction out of a full Hamiltonian.
Effects of trigonal distortion and further neighbour interactions -In principle, there are other small but finite interactions; few examples in δH include where Γ is introduced when a trigonal distortion is present 26 ; J 2 and J 3 are the second and third n.n. Heisenberg interactions respectively. It is natural to expect that they are smaller than the n.n. Kitaev, Gamma, and Heisenberg interactions 18,27,28 . Several types of interlayer exchange interactions are present, but they are even smaller than the terms considered in Eq. 12 18 .
Let's investigate how they affect the above analysis done for the ideal n.n. Hamiltonian.
First of all, the isotropic interactions such as further neighbour J 2 , J 3 , and the interlayer Heisenberg do not make any change to our proposal, since they do not contribute to δω A nor δω K . On the other hand, the Γ modifies the exchange parameters as follows: The conventional anisotropy δω A is now due to Γ + 2Γ obtained from J Z − J XY . Thus to single out the Kitaev interaction, one has to find both Γ and Γ , as J ac is a combination of K, Γ and Γ . Once the trigonal distortion is present, the g-factor also becomes anisotropic, i.e., the in-plane g a is different from the c-axis g c , which affects δω A .
However, the g-factor anisotropy does not affect the δω K , since the field angles of θ and −θ involve the same strength of in-and out-of-plane field components, i.e, h(θ) = h aâ + h cĉ and h(−θ) = −h aâ + h cĉ . Thus we wish to extract the information of K and Γ − Γ from δω K , as it is free from the g-factor anisotropy.
We note that δω K at the K-point, Eq. (10) offers both J ac and J ab from the first term independent of the field and the next term proportional to 1/h eff (h eff = h g 2 a cos 2 θ + g 2 c sin 2 θ).
Once J ac and J ab are deduced, K and Γ − Γ can be estimated from Eq. (13). The measurements of δω K at the K-point with a large magnetic field then determine K and Γ − Γ separately. Further neighbor Heisenberg interactions, J 2 and J 3 do not modify Eq. (10) in the high-field limit, so they do not affect our procedure.

DISCUSSION
We propose an experimental setup to single out the Kitaev interaction for honeycomb Mott insulators with edge-sharing octahedra. In an ideal octahedra cage, the symmetryallowed n.n. interactions contain the Kitaev, another bond-dependent Γ and Heisenberg interactions. We prove that the magnetic anisotropy related by the π-rotation around the a-axis denoted by δω K occurs only when a combination of K and Γ, i.e. K − Γ, is finite.
This can be measured from the spin excitation energy differences under the magnetic field of angle sweeping from above to below the honeycomb plane using the FMR or INS techniques.
Since the in-and out-of-plane magnetic anisotropy, δω A is determined solely by Γ, one can estimate Γ strength first from δω A and then extract the Kitaev interaction from δω K .
While the trigonal distortion introduces an additional interaction, the Kitaev interaction is unique as it is the only interaction that contributes to δω K without altering δω A . Our theory is applicable to all Kitaev candidate materials including an emerging candidate RuCl 3 .
In particular, since the two dominant interactions are ferromagnetic Kitaev and positive Γ interactions in RuCl 3 3,5,18,27 , leading to a large J ac and a small J ab , we predict that δω K independent of the g-factor anisotropy is significant except at the K-point. Supplementary Note 4 shows the FMR and INS of a set of parameters with a small negative Γ interaction to stabilize a zero-field zig-zag ground state as in RuCl 3 18,27,28 . Another relevant perturbation in some materials is the effect of monoclinic structure which loses the C 3c symmetry of R3, making the z-bond different from the x-and y-bonds. The current theory of finite δω K due to a finite J ac still works for C2/m structure. However, since the z-bond of J z ac (= K z /3−Γ z /3) is no longer the same as the x-and y-bonds of J x ac (= J y ac ) and C 2a symmetry relates between the x-and y-bonds, the anisotropy δω K at different momenta, detecting both J x ac = K x /3 − Γ x /3 and J z ac , is required to determine different x-and z-bond strengths. The symmetry-based theory presented here is also valid for higher spin models with the Kitaev interaction such as S = 3/2 CrI 3 including a nonzero single-ion anisotropy 19,21,22,48 which generates a further anisotropy in δω A but does not affect the δω K . The next near-est neighbor Dzyaloshinskii-Moriya interaction with the d-vector along the c-axis 49 is also invariant under the C 2a symmetry. Further studies for higher-spin models remain to be investigated to identify higher-spin Kitaev spin liquid. We would like to emphasize that the proposed set-up is suitable for other experimental techniques such as low-energy terahertz optical and nuclear magnetic resonance spectroscopies that probe spin excitations in addition to the angle-dependent FMR and INS spectroscopy shown in this work as examples.

Exact Diagonalization Simulations
where the Lehmann representation is used; |λ and |λ are the eigenstates with the thermal population factor p λ , and S α,β are the spin operators. In the low temperatures, we take |λ to be the ground state |0 and we are interested in the lowest energy excitation to |1 with a nonzero probability. For optical spectroscopies such as FMR, α = β = direction of the RF electromagnetic field and q = 0, so |0 and |1 belong to the same momentum sector. The structure factor simplifies to For INS, the finite q must match the difference in the momenta of |0 and |1 . For simplicity, we calculate the DSSF for α = β, α S αα (q, ω) = α 1|S α q |0 2 δ( ω + E 0 − E 1 ).
Linear Spin Wave Theory -The Hamiltonian in Eq. (2) is bosonized by the standard Holstein-Primakoff transformation 53 expanded to linear order in the spin S: where the quantization axis is parallel to the c-axis. The Fourier transforms are a j = 1 √ N k e ik·r j a k for sublattice A and b j = 1 √ N k e ik·(r j +δ) b k for sublattice B, where δ is the vector pointing to nearest neighbors. The resulting quadratic Hamiltonian has the form . Diagonalizing this BdG Hamiltonian following standard methods 54 gives two spin wave excitation branches.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code used to generate the data used in this study is available from the corresponding author upon reasonable request.