Magnetic Anisotropy in Spin-3/2 with Heavy Ligand in Honeycomb Mott Insulators: Application to CrI3

Ferromagnetism in the two-dimensional CrI3 has generated a lot of excitement, and it was recently proposed that the spin-orbit coupling in Iodine may generate bond-dependent spin interactions leading to magnetic anisotropy. Here we derive a microscopic spin model of S=3/2 on transition metals surrounded by heavy ligands in honeycomb Mott insulators using a strong-coupling perturbation theory. For ideal octahedra we find the Heisenberg and Kitaev interactions, which favor the magnetic moment along the cubic axis via quantum fluctuations. When a slight distortion of the octahedra is present together with the spin-orbit coupling, three additional terms, the off-diagonal symmetric interactions Γ and Γ′, and on-site Ising interactions arise. They result in the magnetic anisotropy that pins the moment perpendicular to the honeycomb plane as observed in a single layer CrI3, revealing the strong impact of octahedra enviornment on the spin model. A gap in the spin-wave spectrum and comparison to the spin-orbit coupled Jeff= 1/2 and S=1 models are also presented.

Ferromagnetism in the two-dimensional CrI3 has generated a lot of excitement, and it was recently proposed that the spin-orbit coupling in Iodine may generate bond-dependent spin interactions leading to magnetic anisotropy. Here we derive a microscopic spin model of S=3/2 on transition metals surrounded by heavy ligands in honeycomb Mott insulators using a strong-coupling perturbation theory. For ideal octahedra we find the Heisenberg and Kitaev interactions, which favor the magnetic moment along the cubic axis via quantum fluctuations. When a slight distortion of the octahedra is present together with the spin-orbit coupling, three additional terms, the off-diagonal symmetric interactions Γ and Γ , and on-site Ising interactions arise. They result in the magnetic anisotropy that pins the moment perpendicular to the honeycomb plane as observed in a single layer CrI3, revealing the strong impact of octahedra enviornment on the spin model. A gap in the spin-wave spectrum and comparison to the spin-orbit coupled J eff = 1/2 and S=1 models are also presented.

I. INTRODUCTION
Transition metal trihalides (TMT) are layered materials composed of transition metals (M) and halides (X) of the group 9 in a 1:3 ratio. They have a honeycomb layered structure, and depending on the filling of d-orbitals in the transition metals, some are semiconductors and some are metals. [1] Among them, RuCl 3 , VI 3 and CrI 3 are Mott insulators. Magnetic orderings in these systems further establish the importance of electronic correlations and call for microscopic understanding of spin models. For exampe, based on a strong-coupling perturbation theory of the generic spin model [2][3][4], it was shown that α-RuCl 3 described by the effective spin J eff = 1/2 has dominant bond-dependent Kitaev and off-diagonal symmetric Γ interactions. [5,6] RuCl 3 has become an emergent candidate of spin-1/2 Kitaev spin liquid [7]. Intense research activities on various properties of RuCl 3 have been carried out [6,[8][9][10][11][12][13][14] and recently a magnetic-field induced spin liquid was suggested. [15][16][17][18][19][20][21] In parallel theoretical interest in the ground state of higher-spin Kitaev models was initiated by classical model studies. [22,23] The classical Kitaev model has a macroscopic degeneracy named a classical spin liquid [22], but the higher-spin quantum Kitaev model is not exactly solvable, and the ground state is currently unknown. Various numerical studies such as exact diagonalization on S=1 suggested that the ground state is possibly a spin liquid with gapless excitations. [24] These studies were mainly theoretical interests, until a microscopic derivation of the S=1 Kitaev model in multi-orbital systems was found. [25] A heavy ligand spin-orbit coupling (SOC) and strong Hund's coupling in e g -orbitals is a way to generate S=1 bond-dependent Kitaev interaction. The magnetic field effects on the S=1 Kitaev model have also been investigated. [26][27][28] The bond-dependent interactions have recently been adopted into TMT systems, because the nearest neighbour (n.n.) Heisenberg J, Kitaev K and Γ interactions * hykee@physics.utoronto.ca are allowed based on the symmetry of the lattice. [4,29,30] In particular, ferromagnetism in a single-layer CrI 3 has generated excitement in recent years. [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] CrI 3 is a ferromagnetic (FM) insulator with Tc ∼ 61K for bulk samples. [46][47][48] A single layer CrI 3 was successfully synthesized, which showed an FM ordering at Tc ∼ 45K. [49] The two-dimensional FM Heisenberg model is insufficient to explain finite T c , i.e., the Mermin-Wagner theorem [50], and several theoretical models were proposed to explain the magnetic anisotropy. They include the XXZ model [51,52], Ising anisotropy and Kitaev [53], and large Kitaev and small symmetric off-diagonal Γ interactions [54].
While the J, K, and Γ interactions are allowed by the symmetry, and found to be significant in the earlier derivations for lower-spins [3,4,25], they may not be applicable to S=3/2. Thus, a microscopic derivation of S=3/2 model is necessary to find the sources of the magnetic anisotropy. Here we derive a n.n. spin model for S=3/2 with three electrons in t 2g orbitals at transition metal sites and strong SOC in p-orbitals of ligands. We take into account strong electron-electron interaction in multi-orbital systems including Hund's coupling and effects of octahedra distortions present in R3 rhombohedral lattice. Contributions from e g orbitals are also included. The minimal n.n. model includes J, K, Γ, another symmetric off-diagonal Γ [55], and on-site Ising denoted as A c along c-axis.
The rest of the paper is organized as follows. In Sec.II, the on-site Kanamori interaction is reviewed, and tight binding Hamiltonian is derived. In Sec.III, we derive the spin model for ideal octahedra environment. In Sec.IV we analyze the effects of spin interactions on the FM moment direction. In Sec.V, we study the spin model for the distorted octahedra present in R3 rhombohedral lattice. This includes the distortion-induced hopping matrix elements and three additional spin interactions generated via combined effects of SOC and distortion. In Sec.VI we discuss the origin of the spin gap, finite T c , and spin wave spectrum within JKΓΓ A c model including the second n.n. Dzyaloshinskii-Moriya (DM) interaction. Finally in Sec.VII, we summarize our results and compare with J eff =1/2 and S=1 spin models. The detailed calculations are presented in the Appendix. Edge-shared octahedra honeycomb structure unit cell a, b, in global coordinates xyz. Transition metal sites M in gray and non-magnetic ligands X in purple. The n.n. bonds X, Y, Z are related by C3 symmetry. The sites M1, M2, X1, and X2 are involved in the second order strong-coupling expansion on the Z bond. Indirect hopping integrals t0, t1 and t2 are shown, and other direct hoppings can be found in the Appendix.

II. KANAMORI INTERACTION AND TIGHT BINDING HAMILTONIAN
The honeycomb network is made of metal (M) d-orbital sites with half filled t 2g orbitals and the octahedra cages of non-magnetic ligand (X) sites with fully occupied porbitals. The full Hamiltonian is composed of the on-site Kanamori interacton and tight binding Hamiltonian between two sites.
The on-site Hamiltonian of the M sites is described by the Kanamori interaction [56] as well as crystal field spitting: where the density operator n ασ is given by c † ασ c ασ , and c † ασ is the creation operator with α orbital and spin σ. U and U are the intra-orbital and inter-orbital density-density interaction respectively, and J H is the Hund's coupling for the spin-exchange and pair-hopping terms. ∆ c is a crystal field splitting on the M sites, originated from the surrounding octahedra, leading to the splitting of the d-orbitals into t 2g and e g orbitals. In a d 3 system one has half-filled t 2g orbitals, where the Hund's coupling J H selects for the S=3/2 configuration as the ground state. A table of the excited state energy spectra is show in the Appendix 1. The energies of the exited state are larger than the hopping in-tegrals, which allows us to treat the tight-binding hopping parameters as a pertubation.
In the edge shared octahedra structure, each bond between n.n. M sites involves two adjacent ligands, as shown in Fig.1. A tight binding Hamiltonian between two transition metal sites M 1 and M 2 including the two adjacent ligands X 1 and X 2 is given below.
where 0 nxn refers a nxn null matrix. The basis is chosen as are three porbitals at ligand site X m . Each block of indirect hopping between M i and X m is denoted by T M i Xm and the direct hopping between M sites by T M1M2 .

III. IDEAL HONEYCOMB STRUCTURE
To understand the microscopic origin of the spin model, we start with the ideal honeycomb network surrounded by perfect edge shared octahedra. We focus on the Z bond of the honeycomb Fig.1, as the other two bonds are related by C 3 symmetry.

A. Superexchange path: indirect hopping
Indirect hoppings indetween the M and X sites are the largest hopping parameters. The indirect hopping matrix between M 1 and X 1 sites in Fig.1 is: where t 1 , t 2 are related to the σ Slater-Koster integral t 1 = √ 3t 2 = t pdσ √ 3/2 < 0 and t 0 is the π Slater-Koster integral t 0 = t pdπ > 0. The non-zero hoppings of Eq.(3) are depicted in Fig.1. The other indirect bonds T M2X2 , T M2X1 , T M1X2 are obtained from T M1X2 by consecutively applying the symmetry operation of the octahedron C 4 (0, 0, 1).
To account for the indirect d to p hoppings we first integrate out the p orbitals through a pertubative procedure truncated at second order, leading to an effective d to d hopping model: where (a, m) represent a sum over all single hole states a of all sites X m . The hole states are SOC states, thus creating two energy costs ∆E a = ∆−λ p /2 or ∆+λ p , where ∆ = d − p is the atomic energy difference between M and X sites, and λ p is the SOC in p-orbitals. Carrying out the procedure for t 2g − p − t 2g paths results in the following hopping matrix in block form: where σ i with i = x, y, z are the Pauli matrices carrying the spin degrees of freedom and σ 0 is the 2 × 2 identity matrix. The nature of the intermediate hole states, as well as the form of our indirect hoppings T M i Xm , have selected only certain σ matrices in the effective hopping matrix. For example in the limit λ p → 0, only the d yz − d xz term t eff σ 0 → t 2 0 /∆ is present which contributes to the direct hopping channel. The new terms present for non-zero λ p are the spin-flip (SF) terms between d yz /d xz and d xy . Such terms would normally not appear in the second order perturbation process, as they involve a d yz /d xz − p z hopping followed by a p x /p y − d xy hopping, which will only occur if p z is entangled with p x /p y . The SOC among the p x , p y , p z generates such entanglement. Furthermore, when SOC is the dominant energy scale of the hole states, the wavefunctions are inevitably mixtures of p-orbitals and their spin, leading to σ i dependence proportional to the r ratio of the SOC and atomic energy difference.
We now perform the strong coupling perturbation processes where the on-site Hamiltonian given by Eq.(1) is the dominant term, and the effective hopping matrix Eq.(4) is considered as a perturbation. It is straightforward to check that the symmetry of the edge-shared octahedra crystal allows Heisenberg J, Kitaev K, and symmetric off-diagonal Γ interactions [4,29,30]. Truncating at second order in perturbation we arrive at the following Heisenberg-Kitaev (JK) spin model for the ideal honeycomb where γ = x, y, z bond, and J 0 and K 0 have contributions from t 2g and e g . The superexchange process involving only t 2g orbitals will be captured by the effective t 2g − t 2g hoppings of Eq.(5) leading to The spin-dependent hoppings have generated a S z i S z j Kitaev interaction. This can be rudimentary understood by a kinetic argument. To simplify the argument lets focus on one spin 1/2 electron hopping through only a SF hopping. Imagining two sites starting in (↑, ↑) state, the SF hopping can lower the energy by the process: , at a energy cost of −(rt eff ) 2 /U , while imagining the two sites starting in (↑, ↓) the SF process is forbidden from Pauli exclusion principle. Thus the path of starting and ending at (↑, ↑) lowers the energy, and their is a FM S z i S z j preference in the system. Carrying out the details for the S=3/2 leads to the exact forms of Eq. (7).
Among the symmetry allowed terms, the symmet- does not occur in the second order perturbation results. This operator would connect (↑, ↑) to (↓, ↓), which is definitely possible from a SF process, however, there is a subtle cancellation. Having a single Pauli matrix in the SF term creates such cancellation within the spin block, resulting in a null Γ term. Thus it becomes finite only when higher order perturbation terms are included, or when the octahedra are no longer ideal, as we will show in the later Sec.V. Before we present the direct hopping contribution, we show the contribution from e g orbitals. Given that the p orbitals hybridization with e g is larger than t 2g , this contribution is essential. The final form of T eff M1M2 (e g ⊗t 2g ) hopping includes the spin-dependent terms in the effective d yz , d xz to e g blocks, as well as spin-independent −2t eff (t 2 /t 0 )σ 0 hopping between d xy and d 3z 2 −r 2 . The detailed matrix can be found in Eq.(A.3) in the Appendix, and carrying out the strong coupling expansion with this hopping leads to the additional contribution to the Heisenberg and Kitaev interactions Similar to the t 2g − t 2g case, the spin-dependent hoppings contribute to the Kitaev interaction, while the spin-independent hopping, −2t eff (t 2 /t 0 )σ 0 generates a FM Heisenberg term. The FM Heisenberg interaction is origniated from the competition of two exited states separated by Hund's coupling, i.e, e g paths consistent with the earlier findings obtained by the first principle calculations. [45,57] B. Direct hopping Let us consider the contribution of direct hoppings to the spin model. The direct hopping, T M1M2 between M 1 and M 2 is given by The decomposition of the Slater-Koster is explained in the Appendix. Since there is no spin-dependent hopping terms, this does not generate the Kitaev interaction, but changes the Heisenberg interaction. Combining both indirect and direct hopping contributions, the two exchange interactions for the ideal octahedra environment Eq.(6) are found to be Note that the Kitaev interaction has the prefactor of r 2 = (2λ p /(2∆+λ p )) 2 . This leads to a smaller Kitaev compared to Heisenberg interaction, unless the contribution from e g path reduces the overall strength of the Heisenberg interaction.
Below we first show that the magnetic moment of JK model is pinned along the cubic axis, different from the experimental moment direction ofĉ = [111]/ √ 3, i.e, c-axis in CrI 3 . This leads to our motivation to study the effects of distortions of the octahedra.

IV. MAGNETIC ANISOTROPY
We have shown how the ideal honeycomb structure leads to Heisenberg-Kitaev (JK) model up to second order in perturbation. The magnetic moment of the ferromagnetic state obtained with the FM JK model is pinned along the cubic axis such as [100] and C 3 equivalent directions via quantum fluctuations [58,59]. We show three examples in Fig.2 to illustrate the moment pinning direction.
Following the method in Ref. [59], we perform an exact diagonalization calculation on a eight site honeycomb (details of the cluster in the Appendix) to get the ground state wavefunction |GS , and then compute the probability distribution P = | Ψ F M (θ, φ)| GS | 2 , where |Ψ F M (θ, φ) is a ferromagnetic ansatz with moment direction pointing at (θ, φ) on the sphere. Results are show in Fig.2 for different values of JKΓ. Independent of ratio of J and K, the moment is along the cubic axis as shown in the panel (a) and (b). On the other hand, when a small FM Γ is introduced, the moment is along theĉ direction as shown in the panel (c). This effect can be anticipated from the classical analysis. The classical JKΓ model in the FM state with moment S 0 = (S x 0 , S y 0 , S z 0 ) has an energy density per unit cell ] ⇒ S 0 =ĉ leading to a moment pined on the [1,1,1]. We choose Γ/K 1 to show that a tiny FM Γ anisotropy results in theĉ-axis moment shown in Fig.2 (c).
The JK model does not have a magnetic anisotropy, and thus anisotropic interactions beyond J and K are necessary to understand the magnetic anisotropy in CrI 3 . The Γ interaction allowed by the symmetry can be finite if higher order perturbation terms are included. However, one may ask if there are other interactions allowed by a slight distortion of the lattice within the second order perturbation theory without invoking higher order terms. Indeed TMT materials do not have ideal octahedra, but have either rhombohedral R3 or monoclinic C2/m structures, and their magnetism strongly depends on structural differences and number of layers. [60][61][62][63] Below we study the effects of distorted octahedra, which induce additional hopping integrals which were forbidden without the distortion, and investigate if the new spin interactions occur linear in the distortion-induced hoppings. If so, they are not insignificant, and play an important role in finite anisotropy responsible for a finite T c , and a finite spin gap at Γ-point in the inelastic neutron scattering measurement. [64,65] V. EFFECTS OF DISTORTION: DISTORTED OCTAHDERA CrI 3 goes through the structural transition from C2/m to R3 structure at low temperature. [48] In the rhomhedral structure, there are two types of X ligand position deviations from the ideal octahedra structure. As shown in Fig.3(a), a single octahedron can be viewed as two shaded triangles. One distortion is the staggered rotations of the two triangles denoted by δx blue arrows, with displacements of X sites perpendicular to theĉ direction. The other distortion is the compression of the distance between these two triangles along the c-axis denoted by δx orange arrows, with displacements along theĉ direction. Here the dimensionless parameters δx and δx are in unit of the distance between the n.n. M sites d M . Analytic formulas of the new positions of X sites under staggered rotations and compression are found in Appendix 3. In the R3 spacegroup there are other type of distortion, namely the M 1 −X 2,4,6 bond length can be different from the M 1 −X 1,3,5 bond length. This type of distortion is generally exceedingly small compared to δx and δx and we neglect it in the following analysis. Fig.3(b) shows a top view of the the honeycomb unit cell with two such distorted octahedra forming the Z bond, and the staggered rotation of the right octahedron is the mirror image of the left octahedron. To obtain the hopping parameters induced by the distortions one makes use of the Slater-Koster rotated bond formulas [66], with details of the procedure described in the Appendix, leading to distorted hopping matrices T MX .
The distortion-induced hopping matrices T MX have all elements non-zero as a result of lowering the local symmetry of the octahedron from O h to D 3 . We denote the new allowed hopping under distortion as δt i . Starting with the distortion-induced hoppings, we follow the procedure described in Section III. Using the distortion-induced T MX matrices, we derive the effective T eff M1M2 (t 2g ⊗ t 2g ) and T eff M1M2 (e g ⊗ t 2g ) (details in Appendix 4) which are then treated perturbativly in the strong coupling pertubation. The minimal spin model for the n.n. is finally given by where α, β, (γ) refers to the γ bond taking α and β spin components. [4,55] If we included only T eff M1M2 (t 2g ⊗ t 2g ) hopping contributions, the strengths of the three new terms are: where δt α and δt b,c are some combinations of different δt i . A full form of the spin model including JKΓΓ A c to leading order in δt i as well as e g contributions are listed in Table IV in the Appendix. Note that both Heisenberg and Kitaev interactions are renormalized by the distortion, but Heisenberg has a linear term in δt, while Kitaev does not. Given that there are five spin interactions only within the nearnest neighbor model, we do not intend to pin down quantative values of these exchange terms. However, it would be useful to present an example, and discuss how they are affected by U , J H , ∆, ∆ c , λ p , t pdσ , t pdπ , and δt i . For CrI 3 crystal structure in Ref. [48] it is straightforward to extract the values of (δx, δx ) (0.0297, 0.0137). Setting δx = 2δx , we compute JKΓΓ A c as a function of δx. An example case is shown in Fig.4 with U , U , J H , ∆ c , λ p , ∆, t pdπ , t pdσ , t ddσ , t ddπ and t ddδ set to values mentioned in the figure caption. The result is obtained wihtout linearlizing δx. Note that when δx = 0, Γ Γ and A c are absent, and they grow linearly in δx for a small δx. The new hopping paths under δx have generated feromagnetic Γ, Γ and A c which would lead to the moment pined alonĝ c direction. The FM Heisenberg term is the largest term while the Kitaev is AFM and remains small even at δx = 0. While the values of spin exchange terms are strongly influenced by the ratios of J H /U , ∆ c /U , or t pdσ /t pdπ , there are general trends: (1) The FM Heisenberg is originated from e g contribution, and thus it becomes further negative when ∆ c is smaller or Hund's coupling is larger. This affects the Kitaev interaction, which then becomes more positive. (2) A smaller λ p makes K, Γ, Γ and A c smaller. (3) When the ratio of t pdσ /t pdπ becomes smaller, the e g contributions are suppressed, leading to enhancement of the Heisenberg interaction towards AFM, while the Kitaev interaction towards the FM sign.

VI. SPIN GAP AND FINITE TRANSITION TEMPERATURE
The preceding results suggest that it is likely that CrI 3 has dominant FM Heisenberg and smaller Kitaev interaction while the magnetic anistoropy is originated from the distortion of the octahedra together with the SOC of heavy ligands. The magnetic ordering pattern changes from bulk to films, suggesting a strong coupling between magnetism and crystal structures, which further implies the importance of the distorted octahedra in the presence of SOC. Our microscopic model has some overlap with the previous studies. The Kitaev and on-site Ising anistoropy in addition to the Hisenberg interaction were found in Ref. [53] using the density functinal theory. If Γ = Γ , the JΓΓ model maps to the XXZ model [51,52] in the a − b − c crystallographic coordinate where J ab = J − Γ and J c = 3Γ. Since the exchange parameters strongly depend on J H /U , ∆ c /U , and tight binding parameters, it is useful to obtain experimental inputs to determine some parameters. The inelastic neutron scattering experiments and magneto-Raman spectroscopy have reported a spin gap of approximately 0.36 meV at the Brillouin zone (BZ) center Γ-point. [64,68] This is also consistent with a small anistoropy found in the ferromagnetic resonance experiment [54], which is about 0.07 meV leading to the spin gap of 0.3 meV.
Based on our spin wave analysis using the JKΓΓ A c model including the second n.n. DM interaction, the spin wave dispersion ω k is expressed as ω 0 + ρ k 2 around the Γ-point. Here ω 0 and ρ are the spin gap and stiffness, respectively, and they are given by The details of the spin wave spectrum is presented in Appendix 7.
The spin gap is rather small, as expected because it is originated from a combination of slightly distorted octahedra and SOC, i.e., Γ, Γ , and A c . While tiny, it is essential for a finite T c in a single layer CrI 3 . In the FM ordered phase, at low temperatures (T ) the magnons are excited and their number is given by Without the spin gap ω 0 , N s diverges in two-dimension at any temperature except T = 0, i.e., the celebrated Mermin-Wanger theorem. Thus one can understand the essential role of ω 0 which cuts the divergence, and allows the FM ordering at finite temperatures, as long as ω 0 (T ) remains finite for T ≤ T c . While quantifying the transition temperature requires further analysis including magnon-magnon interactions at finite temperatures, the temperature dependence of ρ(T ) and ω 0 (T ) from the inelastic neutron scattering measurement [64] indicates the crucial role of ω 0 (T ) which vanishes at T c . Another important parameter is the Kitaev interaction K, which leads to a gap at the the BZ corner Kpoint known as Dirac gap [54], reported in the neutron scattering. [64,65] However, the second n.n. DM term also generates the Dirac gap. [64,65] Note that neither contribute to the spin gap at Γ-point. We would like to point out that Γ and Γ also play a part in the Dirac gap as shown in Eq.(A.14). Further analysis on the individual role of K and DM interactions remain to be resolved in future studies.

VII. SUMMARY AND DISCUSSION
We have shown a microscopic derivation of the n.n. spin model for honeycomb Mott insulators with three electrons in t 2g orbitals at M sites surrounded by octahedral cages of ligands X with strong SOC. Using the standard strongcoupling perturbation theory, we found that there are only Heisenberg J and Kitaev K interaction for the ideal honeycomb lattice among the three symmetry allowed interactions (J, K, Γ), because Γ is zero up to the second order perturbation term. The exchange paths between t 2g and t 2g vs. t 2g and e g via ligands generate opposite signs for both J and K interactions. The Heisenberg interaction is of order t 2 eff /U , while the Kitaev is smaller by a factor of r 2 ∼ (λ p /∆) 2 . The FM Heisenberg interaction is originated from the e g paths, with the hopping integral between e g and p orbitals being larger compared to t 2g and p orbitals.
The FM Heisenberg and Kitaev interaction leads to a FM ordering, but the moment direction is pinned along the cubic x, y, or z axis, e.g., [100] ( C 3 equivalent direction) via quantum fluctuations. The c-axis, [111] moment pinning found in CrI 3 should thus originated from other interactions, which are also responsible for the spin gap at the Γ-point in the neutron scattering measurement [64]. Including the distorted octahedra present in the rhombohedral structure, three additional spin interactions are found.
Inspecting the linear order in the distortion-induced hopping paths within the second order perturbation theory, Γ, Γ and on-site Ising anisotropy A c contain terms linear in the distortion-induced hopping integrals. The Heisenberg interaction also contains such additional linear term, but the Kitaev does not, implying that it is possible to finetune a system closer to the Kitaev-dominant regime via octahedra distortions.
Comparison to J eff = 1/2, S=1, and 3/2 spin systems would be useful. The SOC is necessary to generate the bond-dependent interaction, as spin and orbital should be entangled to get such a directional dependent spin interaction. However, the presence of SOC is not enough to find an exotic phase like a spin liquid, because the dominant interaction is often the Heisenberg interaction. To compare different spin cases, a summary of the ideal honeycomb exchange interactions, J 0 , K 0 , and Γ 0 including the effective indirect hopping (t eff ) and direct hopping (t d ) integrals, only up to the second order perturbation terms, is shown in the following Table I. Focusing on the ideal octahedra and n.n. model via second order superexchange processes, J eff = 1/2 is unique because the Heisenberg term is absent. On the other hand, for the S=1 model with d 2 in e g , the heavy ligand SOC generates the Kitaev interaction, which has the same order of magnitude with J. In fact, K = −2J, if only e g paths are considered. [25] For S=3/2 case, we found that J is order 1 in units of roughly t 2 eff /U , while K is smaller by r 2 . Thus it is hard to compete with the Heisenberg interaction. We speculate that this is valid for spins equal or higher than 3/2. Unlike S=1/2 and S=1 cases, the Heisenberg interaction is dominant in S=3/2, but we cannot rule out a possibility of cancellation among different contributions to the Heisenberg interaction, which may let the Kitaev interaction overtake a major place. In particular, the Heisenberg interaction is more sensitive to the distortion-induced hopping integrals than the Kitaev term as shown in the table in the Appendix, manipulating ocathedra may be a way to tune the system to a desired Kitaev dominant regime.

Spin
Heisenberg J0 Kitaev K0 symm. off-diagonal Γ0 In summary, in the ideal octahedra environment, we find that there are only two spin interactions, Heisenberg and Kitaev interactions. Kitaev interaction is generally weaker compared to the Heisenberg interaction in contrast to the lower spin models. Indirect hoppings among t 2g and t 2g vs. e g and t 2g have opposite contributions. A detailed balance between the two indirect and direct hopping contributions highly depends on the hopping integrals, Hund's coupling strengh, and crystal field spitting. Γ interaction is absent up to the second order due to a subtle cancellation, despite that it is allowed by the symmetry. We further show that octahedra distortion allows three additional interactions, two symmetric off-diagonal interactions Γ and Γ , and Ising anisotropy A c along the c-axis. They are all linearly proportional to a distortion-induced hopping integral. While they are much smaller than the Heisenberg interaction, they are essential for a spin gap in the FM phase of CrI 3 leading to a finite T c . Our study offers a microscopic route to the n.n. spin models, JKΓΓ A c interactions. Given that there are five exchange terms within the n.n. model, and second n.n. interactions including DM may be comparable to Γ, Γ and A c , further theoretical and experimental studies are required to determine the microscopic parameters of CrI 3 .  with the Anderse prediction t ddσ : t ddπ : t ddδ = −6 : 4 : −1. Similarly to the derivation of Eq.(5), the effective tight binding Hamiltonian for t 2g at site M 2 and e g at site M 1 is given by

Positions of X sites under distortions
In the R3 space group, octahedra made of two triangles as shown in Fig.3 are rotated aroundĉ-axis denoted by δx where X m lable sites as in Fig.3(a) andâ,b andĉ are unit vectors along the lattice vectors of Fig.3(b).

Distortion-induced hoppings
To get the distortion-induced hoppings one needs to find the directional cosines of the indirect d − p bonds using the ligand positions Eq.(A.4). With the metal site located at r M1 and ligand X m located at r m , one can then find the required directional cosines from r m − r M1 . They are the inputs in the Slater-Koster formula for rotated bonds Ref. [66]. Here we only show the leading term in each matrix element generated from the above procedure. Note that all elements become finite, and the distortion-induced hopping between t 2g and p oribtials are denoted by δt i and e g and p orbitals by δt i : where X 1 lable the site as in Fig.3(a). The distorted octahedron realizes the D 3 point group, which contains C 3 (1, 1, 1) and C 2 (−1, 1, 0) rotations. The hoppings T M1X3 and T M1X5 are recovered by applying C 3 (1, 1, 1) to T M1X1 . Further T M1X1 relates to T M1X2 by a C 2 (−1, 1, 0). Finally T M1X2 relates to T M1X4 and T M1X6 by C 3 (1, 1, 1). The direct hopping integrals denoted by T M1M2 is same as Eq. (9). Making use of Eq.(A.5) one can derive effective hoppings T eff M1M2 following the same method as in the ideal case Sec.III A. The distortion induced effective hoppings to leading order in δt and δt are listed in Table III.

On-site Ising spin interaction
Due to the distortion and SOC, the on-site anisotropic term is also generated via the hopping to anions which can hop back to create spin-dependent on-site terms denoted by T eff M1M1 = T eff M2M2 . Without the octahedra distortion, the effective hopping integrals between t 2g at M 1 are given by (A. 6) Similarly, effective hopping between t 2g and e g is found as TABLE III. Effective hoppings to leading order in distortion induced hoppings δt and δt of Eq.(A.5). Distortion induced hoppings are grouped in: δta = δt1 − δt4, δt b = 2δt1 + δt4, δtc = δt3 + δt5, δt d = δt1 − δt2, δte = δt1 + δt4, δt f = δt1 + δt2.