6 form where J ′ is the derivative of Heisenberg coupling J$_{cluster}$ between two spins ℓ and m with respect to their distance r$_{ℓm}$ , J ′ = d J$_{cluster}$/ d r$_{ℓm}$ ; Q E 1 , 2 are the generalized coor- dinates of these two modes; and the functions f E 1 , 2 are where λ (in fact λ/k ) is the expansion parameter. Consider the perturbation H$_{perturbation}$ = Q E 1 j · Q E 1 $_{k}$, which means a coupling between the Q E 1 distortion modes of the two tetrahedra. Integrate out the optical phonons, at lowest non-trivial order, it produces a term (3 J $^{$^{′}$2}$λ ) / (4 k $^{2}$) τ x j · τ x $_{k}$. This can be seen by minimizing separately the two cluster Hamiltonians with respect to Q E 1 , which gives Q E 1 $_{= (}$√$_{3}$$_{J}$ $^{′}$) / (2 k ) τ $^{x}$, then plug this into the perturbation term. Thus we have produced the J$_{x}$ term in the Kitaev model with J$_{x}$ = − (3 J $^{$^{′}$2}$λ ) / (4 k $^{2}$). Similarly the perturbation H$_{perturbation}$ = Q E 2 j · Q E 2 k will generate (3 J $^{$^{′}$2}$λ ) / (4 k $^{2}$) τ y j · τ y k at lowest non-trivial order. So we can make J$_{y}$ = − (3 J $^{$^{′}$2}$λ ) / (4 k $^{2}$). The τ z j · τ z k coupling is more difficult to get. We treat it as − τ x j τ y j · τ x $_{k}$τ y $_{k}$. By the above reasoning, we need an anharmonic coupling H$_{perturbation}$ = Q E 1 $_{j}$Q E 2 j · Q E 1 $_{k}$Q E 2 $_{k}$. It will produce at lowest non-trivial or- der (9 J $^{$^{′}$4}$λ ) / (16 k $^{4}$) τ x j τ y j · τ x $_{k}$τ y $_{k}$. Thus we have J$_{z}$ = (9 J $^{$^{′}$4}$λ ) / (16 k $^{4}$). Finally we have made up a spin-lattice model H$_{SL}$ , which involves only S$_{ℓ}$ · S$_{m}$ interaction for physical spins, J $^{′}$( Q E 1 f E 1 + Q E 2 f E 2 ) $_{(}$√$_{3}$$_{/}$$_{2)}$$_{J}$ $^{′}$( − Q E 1 τ x + Q E 2 τ $^{y}$) (12) According to TABLE I we have f E 1 = − $_{(}$√$_{3}$$_{/}$$_{2)}$$_{τ}$ x and f E 2 $_{= (}$√$_{3}$$_{/}$$_{2)}$$_{τ}$ $^{y}$. Then the coupling becomes f E 2 = (1 / 2)( S$_{2}$ · S$_{4}$ + S$_{1}$ · S$_{3}$ − S$_{1}$ · S$_{4}$ − S$_{2}$ · S$_{3}$ ) , f E 1 = √ 1 / 12( S$_{1}$ · S$_{4}$ + S$_{2}$ · S$_{3}$ + S$_{2}$ · S$_{4}$ + S$_{1}$ · S$_{3}$ − 2 S$_{1}$ · S$_{2}$ − 2 S$_{3}$ · S$_{4}$ ) . H$_{cluster}$ j, SL = H$_{cluster}$ j + k 2 ( Q E 1 $_{j}$) 2 + k 2 ( Q E 2 $_{j}$) 2 − √$_{3}$ 2 J $^{′}$( Q E 1 $_{j}$τ x j − Q E 2 $_{j}$τ y j ) , (13) H$_{SL}$ = ∑ cluster H$_{cluster}$$_{,}$ SL + ∑ x $^{−}$links <jk> λ$_{x}$ Q E 1 j · Q E 1 k + ∑ y $^{−}$links <jk> λ$_{y}$ Q E 2 j · Q E 2 k + ∑ z $^{−}$links <jk> λ$_{z}$ Q E 1 $_{j}$Q E 2 j · Q E 1 $_{k}$Q E 2 k H$_{jk,}$$_{SL}$ = H$_{cluster}$ j, SL + H$_{cluster}$ k, SL + λ H$_{perturbation}$ [ Q E 1 $_{j}$, Q E 2 $_{j}$, Q E 1 $_{k}$, Q E 2 $_{k}$] H$_{SL}$ = ∑ cluster j [ ( J$_{cluster}$/ 2)( S$_{j}$$_{1}$ + S$_{j}$$_{2}$ + S$_{j}$$_{3}$ + S$_{j}$$_{4}$ ) 2 + k 2 ( Q E 1 $_{j}$) 2 + k 2 ( Q E 2 $_{j}$) 2 + J ′ ( Q E 1 j S$_{j}$$_{1}$ · S$_{j}$$_{4}$ + S$_{j}$$_{2}$ · S$_{j}$$_{3}$ + S$_{j}$$_{2}$ · S$_{j}$$_{4}$ + S$_{j}$$_{1}$ · S$_{j}$$_{3}$ − 2 S$_{j}$$_{1}$ · S$_{j}$$_{2}$ − 2 S$_{j}$$_{3}$ · S$_{j}$$_{4}$ √$_{12}$ + Q E 2 j S$_{j}$$_{2}$ · S$_{j}$$_{4}$ + S$_{j}$$_{1}$ · S$_{j}$$_{3}$ − S$_{j}$$_{1}$ · S$_{j}$$_{4}$ − S$_{j}$$_{2}$ · S$_{j}$$_{3}$ 2 )] − ∑ x $^{−}$links <jk> 4 J$_{x}$k 2 3 J $^{′}$2 Q E 1 j · Q E 1 k − ∑ y $^{−}$links <jk> 4 J$_{y}$k 2 3 J $^{′}$2 Q E 2 j · Q E 2 k + ∑ z $^{−}$links <jk> 16 J$_{z}$k 4 9 J $^{′}$4 Q E 1 $_{j}$Q E 2 j · Q E 1 $_{k}$Q E 2 k (14) The spin-lattice(SL) Hamiltonian on a single cluster j is [equation (1.8) in Ref. $^{35}$], where k > 0 is the elastic constant for these phonon modes, J ′ is the spin-lattice coupling constant, Q E 1 j and Q E 2 j are the generalized coordinates of the Q E 1 and Q E 2 distortion modes of cluster j , H$_{cluster}$ j is (2). As al- ready noted in Ref. $^{35}$, this model does not really break the pseudo-spin rotation symmetry of a single cluster. Now we put two clusters j and k together, and in- clude a perturbation λ H$_{perturbation}$ to the optical phonon Hamiltonian, where Q E 1 j is the generalized coordinate for the Q E 1 mode on cluster j , and Q E 1 $_{k}$, Q E 2 $_{j}$, Q E 2 k are similarly defined; λ$_{x,y}$ = − (4 J$_{x,y}$k $^{2}$) / (3 J $^{$^{′}$2}$) and λ$_{z}$ = (16 J$_{z}$k $^{4}$) / (9 J $^{$^{′}$4}$); the single cluster spin-lattice Hamiltonian H$_{cluster}$$_{,}$ SL is (13). Collect the results above we have the spin-lattice Hamiltonian H$_{SL}$ explicitly written as, The single cluster spin-lattice Hamiltonian [first three lines in (14)] is quite natural. However we need some harmonic(on x - and y -links of honeycomb lattice) and an- harmonic coupling (on z -links) between optical phonon