Zero energy Majorana modes in spin ladders and a possible realization of the Kitaev model
aa r X i v : . [ c ond - m a t . s t r- e l ] A ug Zero energy Majorana modes in spin ladders and a possible realization of the Kitaevmodel.
A. A. Nersesyan , and A. M. Tsvelik The Abdus Salam International Center for Theoretical Physics, 34100, Trieste, Italy Center of Condensed Matter Physics, ITP, Ilia University, 0162, Tbilisi, Georgia Department of Condensed Matter Physics and Materials Science,Brookhaven National Laboratory, Upton, NY 11973-5000, USA
We show that in double-chain Mott insulators (ladders), disordered alternating ionic potentialsmay locally destroy coherence of magnetic excitations and lead to the appearance of spontaneouslydimerized islands inside the Haldane spin-liquid phase. We argue that a boundary between thedimerized and Haldane phases of a spin-1/2 ladder supports a localized zero-energy Majoranafermion mode. Based on these findings we suggest a realization of a generalized Kitaev modelwhere Majorana fermions can propagate in more than one dimension.
INTRODUCTION
A possible existence of Majorana zero modes in certaincondensed matter systems has been a subject of muchdiscussion in recent literature. The interest is caused bya strong nonlocality of Majorana zero modes which riseshopes of possible applications in such areas as quantumcomputation [1] and electron teleportation [2]. It hasbeen argued that stable Majorana zero energy modes(MZEM) may appear on surfaces [3] and in vortex coresof certain superconductors with a nontrivial topologicalinvariant (this includes spinless p x + ip y ones)[4], [5]. Asfar as magnetic systems are concerned, the existence ofMZEM has been suggested only for a model system (theKitaev model)[6] which material realization has not yetbeen achieved.This paper presents two main results. First, we showthat localized MZEM appear on a boundary betweenspontaneously dimerized (SD) and magnetically disor-dered phases in systems of two coupled spin-1/2 chains(spin ladders). As was demonstrated in [7],[8], the SDphase may exist in spin ladders with sufficiently strongfour-spin interaction. There are reasons to believe [10]that such interactions are present, for instance, in spinladder material CaCu O . Second, we suggest that a bal-ance between competing spin-disordered and SD phasescan be shifted by application of a staggered ionic po-tential which effectively enhances the four-spin interac-tion. This may be important for applications since real-ization of such potentials is quite feasible by chemistrymeans. Based on these findings we suggest a particularrealization of a generalized Kitaev model where MZEMare located on a manifold of arbitrary dimension. As inthe original Kitaev model, these dispersionless modes fa-cilitate propagation of Majorana fermions with nonzerodispersion so that such propagation indeed becomes be-comes akin to teleportation as described in [2].The idea that an alternating ionic potential can lead toSD was expressed in [11],[12] where the authors studieda transition from the Mott to a conventional ionic band insulator (IBI) in one dimension. It was found that sincethe Mott and IBI states differ under a local parity inver-sion, the transition driven by a change in the strength ofthe alternating ionic potential proceeds through an inter-mediate spin liquid phase characterized by a spontaneousdimerization. In the present paper we generalize theseideas for the case of (i) disordered ionic potentials, (ii)the systems consisting of two half filled chains coupledtogether. In the second part of the paper we construct ageneralized Kitaev model. STAGGERED IONIC POTENTIAL IN A SINGLEHUBBARD CHAIN
Let us briefly recall the single Hubbard chain case. Athalf filling such chain has a spectral gap in the chargesector and gapless excitations in the spin sector. At lowenergies the spin sector decouples from the charge oneand can be considered separetly. For reasons of conve-nience we will use the continuous approximation for theentire Hubbard model. This is legitimate if the Mott-Hubbard gap is much smaller than the bandwidth. Theeffective Euclidian Lagrangian of the Hubbard model inthis case is L = L c + L s (1) L c = 12 K h v − c ( ∂ τ Φ c ) + v c ( ∂ x Φ c ) i − G cos √ π Φ c (2) L s = 12 h v − s ( ∂ τ Φ s ) + v s ( ∂ x Φ s ) i (3)where Φ c,s are charge and spin bosonic fields, v c,s aretheir velocities, K is the Luttinger parameter and G isthe Umklapp coupling constant. All these parametersare related to the bare parameters of the Hubbard model.The cosine term in (2) is relevant for K < MH ) in the chargesector. The spin sector remains gapless and the correla-tion functions of the spin fields display power-law decayat T = 0. We omitted in (3) a marginally irrelevant termcontributed by backscattering of the electrons.Let us consider the influence of a staggered ionic po-tential. Using the well known textbook bosonization for-mulae [13],[14] for the staggered charge density we ob-tain the following contribution of such potential to theLagrangian density: V = Z d x V ( π, x ) ρ stag ( x ) (4) ∼ Z d x V ( π, x ) sin √ π Φ c cos √ π Φ s , where V ( π, x ) is an envelope made of the Fourier harmon-ics of the potential concentrated around wave vector π ,and ρ stag ( x ) is the continuum limit part of the staggereddensity ( − n C + nσ C nσ . Since the charge sector has a gapand the vacuum average of sin √ π Φ c is zero, the effectof the staggered potential shows up only in the secondorder of perturbation theory in V (2 k F = π ) [11]: δH = λ h ∂ x Φ sR ∂ x Φ sL − (2 πα ) − cos √ π Φ s i , (5)where Φ sR ; L are chiral components of the scalar field Φ s , λ ∼ V and α is a short-distance cutoff of the bosonictheory. More precisely, the coupling λ is expressed as aconvolution of the disorder average of the potentials withthe correlation function of the charge fields: λ ( x ) = Z d τ d y V ( π, x + y ) V ( π, y ) × (6) hh sin √ π Φ c ( τ, x + y ) sin √ π Φ c (0 , y ) ii ( τ + y /v s ) , and in principle is coordinate dependent. To provide theconvergence of (6) in time domain we need an exponentialdecay of the dynamical correlation function which occursonly if the charge gap is nonzero. To estimate the averagecoupling we may assume that h V ( π, x ) V ( π, i = V exp( −| x | /ξ ) . (7)In the most interesting case when the disorder correlationlength ξ ≪ v c / ∆ we have λ ∼ ξ h V i a ∆ . Interaction (5) ismarginally relevant and competes with the marginallyirrelevant interaction of the same form present in thespin sector of the original Hubbard model. Thereforesuch interaction will open a spin gap only when λ ex-ceeds some critical value related to the aforementionedmarginally irrelevant backscattering term. Refs. [11],[12]consider only a uniform potential, and one could get animpression that the onset of the gapped dimerized phasewould require doubling of the unit cell across the wholechain. However, from the above example it follows that,to achieve the desired effect, it is sufficient to break trans-lational symmetry only locally. DOUBLE CHAINS. THE CASE WHEN CHARGEGAPS >> THAN THE SPIN ONES
Now we consider a system of two coupled half-filledchains. For simplicity one may think about them as Hub-bard chains even though precise details of the charge sec-tor are not important. It was shown in [7],[8] that at halffilling, when the low-energy dynamics in the spin sectorof a single chain coincides with the one for the spin S=1/2Heisenberg model, the continuum limit description of theladder is given by the model of four Majorana fermionswith the Lagrangian density L = 12 χ aR ( ∂ τ − i v∂ x ) χ aR + 12 χ aL ( ∂ τ + i v∂ x ) χ aL + i m t X a =1 χ aR χ aL + i m s χ R χ L , (8)where the Majorana fields χ R , χ L in the path integralare real Grassmann numbers and parameters v, m t , m s depend on the bare interactions. This description workswell when the interchain exchange is much smaller thanthe exchange along the ladder. Then parameters m t , m s (masses of the triplet and the singlet Majorana fermions)are proportional to the linear combinations of the con-ventional J ⊥ and the four-spin J cycle exchange integrals[8],[9]: m t = AJ ⊥ + BJ cycle , m s = − AJ ⊥ + BJ cycle . (9)( A, B are numerical coefficients). The spectra of thetriplet and singlet Majorana fermions are ǫ t,s ( k ) = q ( vk ) + m t,s ; they do not depend on the signs of m t,s .However, since the spin operators are nonlocal in termsof the fermions, their correlation functions do depend onthese signs. In particular, if m t m s < | m t | . If,on the other hand, m t m s > O where the rung exchange is much smaller thanthe leg one the situation is different and the singlet mode(though appearing only in the magnetic continua) visiblyaffects the dynamical correlation functions [10].One may expect that by arranging a proper coordinatedependence of the interchain interactions one can create asituation when some or all of the Majorana masses changesign. This will lead to creation of Majorana zero modes[15]. In [16] it was shown that this does happen whenone puts a single vacancy on one of the chains; then allMajorana masses change sign. The result of such change,as it might be expected, is a creation of local spin 1/2localized at the vacancy. Here we would like to considera different possibility, namely, when only an odd numberof Majorana modes change sign. This happens either ona border between SD and magnetically disordered stateor, if the SU(2) symmetry in the spin sector is lost, bya change of sign in one the components of the Majoranatriplet. The latter can be achieved by applying an easyaxis anisotropy.To be definite we will consider the SD case. We suggestthat from the practical point of view a local change of m s can be achieved by applying staggered ionic potentialsacting on the staggered charge densities of both chains: δH = X i =1 , V i ( x )( − n C + σ,i ( n ) C σ,i ( n ) (10)where V i ( x ) is a slow function of x = na . In order toachieve the desired effect that only m s changes sign weneed V V <
0. In the opposite case V V > m t . For simplic-ity we consider only the former possibiity. To performactual calculations we adopt the following hierarchy ofenergy scales 4 t k ≫ V, ∆ MH ≫ t ⊥ , where ∆ MH is theMott-Hubbard gap for the individual chains, V is the in-terchain density-density interaction and t ⊥ is the inter-chain tunneling. With such assumptions the high energydegrees of freedom are the charge modes which can bedescribed by the following model: L = 12 K X i =1 , h v − c ( ∂ τ Φ ci ) + v c ( ∂ x Φ ci ) i + (11)+ V ∂ x Φ c ∂ x Φ c − G (cid:16) cos √ π Φ c + cos √ π Φ c (cid:17) , where fields Φ c , Φ c are the charge fields of individualchains, V being proportional to the interchain density-density forward scattering. Passing in (10) to the con-tinuum limit and bosonizing the resulting expression weobtain: δH ∼ V ( x ) sin( √ π Φ c )( ǫ + + ǫ − )+ V ( x ) sin( √ π Φ c )( ǫ + − ǫ − ) , (12)where ǫ ± are the symmetric and antisymmetric combi-nations of the chain spin-dimerization operators [17]. Inthe Majorana representation (8) these operators are ex-pressed as linear combinations of products of four orderor disorder parameters of the corresponding Ising models(each noninteracting Majorana fermion is equivalent tothe quantum Ising model; see [7] and also [13] for details): ǫ + ∼ µ µ µ µ , ǫ − ∼ σ σ σ σ . (13) At small distances the following fusion rules take place: ǫ ± ( R + ρ / ǫ ± ( R − ρ /
2) = (14)14 (cid:18) αρ (cid:19) Y i =1 [1 ± i πρ κ i ( R ) + · · · ]= const ±
14 i πα X i =1 κ i − π αρ X i =1 κ i ! + · · · , where κ i = χ iR χ iL . Using fusion rule (14) and keepingonly the most relevant terms we obtain in the secondorder in the potentials the following shift of all masses: m s,t = m (0) s,t + C ( x ) , (15) C ( x ) ∼ Z d y d τ V ( x + y ) V ( y ) ×hh sin √ π Φ c ( τ, x + y ) sin √ π Φ c (0 , y ) ii . Naturally, C vanishes in the absence of interchain in-teractions. Since in the spin-liquid state the massesof triplet and singlet Majorana fermions have differentsigns, such shift may lead to a change in sign of one ofthe masses. In the uniform case this would push the sys-tem to a spontaneously dimerized phase via QuantumCritical Point [8]. If the ionic potential is coordinatedependent the sign change will occur locally. For theHubbard model where the interchain coupling is antifer-romagnetic only the singlet mass changes sign. Then onphase boundaries between SD and disordered phases asingle Majorana mode has zero energy states, as it occurswith fermionic states located on the boundary of quasi-1D p-wave superconductor [3]. Another direct analogy isa boundary between ordered and disordered states in thequantum Ising model.A purely one-dimensional case (single ladder) will dis-play disorder effects studied in [18] and [16] where thecase of random potential was considered. It was shownthat when such boundaries are randomly distributed witha final concentration n i they give a singular contrubutionto the specific heat: C ( T ) ∼ n i ln − ( T /T ) . (16)We remark that the models considered in [18],[16] dealtwith the situation when the number of Majorana zeromodes on each domain wall was four. In that case, as wehave already mentioned, each domain wall carries localspin 1/2. HOW TO CONSTRUCT THE KITAEV MODEL
A much more interesting situation emerges when thenumber of modes is odd like in the case we consider. Letus consider the above ladder model with large areas of thedimerized phase separated by large areas of the Haldane (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)
FIG. 1: Construction of the Kitaev model. Double rails rep-resent spin ladders, the ellipses denote interactions (18) facil-itated by the magnetically active media of the substrate de-scribed in the main text. Areas with different filling representSD and Haldane phases. MZEMs are located at boundariesbetween the phases. phase, such that the wave functions of the zero modesfrom different domain walls do not overlap. Now let usput these ladders on top of a superparamagnet whosemagnetic susceptibility has a peak at ( π,
0) wave vector.So the strongest magnetic fluctuations in such a materialare smooth along the ladder direction and oscillate inthe direction perpendicular to the ladders. Then theystrongly couple to the K a = i( χ R χ aR + χ L χ aL ) , ( a = 1 , , , (17)the operator representing the difference between magne-tization densities of two legs of each ladder [7]. Integrat-ing over the bulk magnetic fluctuations we obtain theeffective interaction between Majorana modes from dif-ferent ladders: V int = 12 K a ( y, x ) J ab ( x − x ′ , r ) K b ( y + r, x ′ ) , (18) J ab ( x − x ′ ; r ) = hh N a ( y, x ) N b ( y + r, x ′ ) ii , where y is a location of the ladder in the transverse direc-tion, x, x ′ are the coordinates along the ladder and N a is the a -th component of the staggered magnetization ofthe paramagnet. Since the paramagnet is disordered, thisinteraction is short ranged. To get the Kitaev model wehave to arrange the positions of domain walls in such apattern that a wall on one ladder interacts only with oneother wall. Such pattern does not need to be regularthough disorder will lead to localization of the propagat-ing Majorana modes.As an example we may consider a pattern where eachladder contains alternating SD and Haldane phases of al-ternating lengths a and b . The ladders are shifted byamount c < a in the y direction (see Fig.1). All lengthsare assumed to be much larger than the domain wallwidth so that the wave functions of MZEMs on differentdomain walls do not overlap. Since the correlation func-tion J ab is short range, we can neglect all interactionsbetween domain walls except between the ones locatedon neighboring ladders. Hence we get H = (19) X y,a =1 , , Z d x h i v − χ aR ∂ x χ aR + χ aL ∂ x χ aL ) + i m t χ aR χ aL i y + X x c ,y J ab ( y, y + 1) × ( χ L χ aL + χ R χ aR ) y,x c ( χ L χ bL + χ R χ bR ) y +1 ,x c , where x c coordinates mark positions of domain walls. Foreach wall one has to specify whether it holds a zero modewith right or left chirality and replace χ operator withits zero mode component: (cid:18) χ R ( x ) χ L ( x ) (cid:19) y = γ ( x c , y ) N e ± R xxc m s ( x ′ ) d x ′ /v (cid:18) ∓ (cid:19) + ..., (20)where the sign depends on the sign of δm s = m s ( ∞ ) − m s ( −∞ ), γ ( x c , y ) is a zero energy Majorana fermionmode satisfying the Clifford algebra anticommutation re-lations, N is the normalization factor, the dots stand forcontributions of higher energy modes.Model (19) is a particular type of the Kitaev modelfermionized by means of the Jordan-Wigner transforma-tion, as in [19]. The latter method, when the original spinmodel is treated as a model of coupled chains and Jordan-Wigner transformation is done for each chain, representscertain advantages since it allows to avoid a cumbersomegauge fixing used in the original paper [1]. The most im-portant feature of the arrangement depicted on Fig. 1is that each MZEM in Hamiltonian (19) interacts onlywith one nearest neighbor exactly as in the original Ki-taev model. In this way the products of neighboring zeroMajorana modes γ ( x, y ) γ ( x, y + 1) are integrals of mo-tion and can be replaced by constants. Thus the fourfermion term in (19) becomes effectively a hopping termfor Majorana fermions χ a which now can propagate inthe transverse direction. The latter fermions have spec-tral gaps, but such sector also exists in the Kitaev model.To make sure that MZEMs are protected, we considerthe operators which couple these modes to local pertur-bations. The dangerous ones are those which can coupleMZEMs located at different positions. Neither K op-erator considered above nor the staggered componentsof magnetization can do this since they contain massivemodes and their correlations falls off exponentially. Thisleaves the staggered energy density operators ǫ ± whichare products of four order ( disorder) Ising model oper-ators (13). The antisymmetric staggered energy density ǫ − contains operators σ a (a=1,2,3) which have no groundstate average and therefore vanishes. Operator ǫ + has anonzero ground state average affected by MZEMs andtherefore may couple to it. The coupling originates fromthe change in µ ; the operators from the triplet sectordo not experience any change at the boundary. UsingEqs.(83,84) from [16] we get ( x < h ǫ + ( x ) i = ( a m t m s ) / exp {− f [ m s ( x − x c )] } , (21) f ( x ) = 12 θ ( x ) x + 18 h K ( | x | ) + K − ( | x | ) i , - - - mx e FIG. 2: ǫ + ( x ) /ǫ ( −∞ ) as a function of 2 m s ( x − x c ). This function is depicted on Fig.2. Here m s standsfor the asymptotic value of the singlet gap far from theboundary. As we see, at the boundary where the MZEMis located the symmetric staggered energy density van-ishes: h ǫ + ( x ) i ∼ | x − x c | / . This fact renders the inter-action with MZEM numerically weak. CONCLUSIONS
We have demonstrated that zero energy Majoranamodes (MZEMs) robust against external perturbationsexist on a boundary between the Haldane and SD phasesof spin S=1/2 ladders. We have also suggested a practicalway to manufacture systems with such phase boundaries.Namely, the effective interactions in spin ladders canbe manipulated by purely chemical means through stag-gered ionic potentials. There are two ways MZEMs fromdifferent domain walls interact between each other: bydirect overlap and via interaction of a particular Fouriercomponent of the magnetization (17). The former inter-action leads to formation of an impurity band inside ofthe spin gap with associated singularity in the specificheat (16). However, by maintaining a proper distancebetween the phase boundaries one can make this overlapsmall. For this case we have suggested a particular ar-rangement when the interactions between MZEMs fromdifferent ladders create a network similar to the one ex-isting in the Kitaev model in its gapped phase.We are grateful to B. L. Altshuler and V. Gritsev for valuable discussions. A. M. T. is grateful to Abdus SalamICTP for kind hospitality and also acknowledges a sup-port from the US DOE, Basic Energy Sciences, MaterialSciences and Engineering Division. A. N. N. is supportedby grants GNSF-ST09/4-447 and IZ74Z0-128058/1. [1] A. Y. Kitaev, Ann. Phys. , 2 (2003).[2] L. Fu, Phys. Rev. Lett , 056402 (2010).[3] K. Sengupta, I. Zutic, H.-J. Kwon, V. M. Yakovenko, andS. Das Sarma, Phys. Rev. B , 144531 (2001).[4] D. A. Ivanov, Phys. Rev. Lett. , 268 (2001).[5] R. Roy, Phys. Rev. Lett. , 185401 (2010).[6] A. Saket, S. R. Hazan, R. Shakar, Phys. Rev. B ,174409 (2010).[7] D. G. Shelton, A. A. Nersesyan and A. M. Tsvelik, Phys.Rev. B , 8521 (1996).[8] A. A. Nersesyan and A. M. Tsvelik, Phys. Rev. Lett. ,3939 (1997).[9] V. Gritsev, B. Normand, and D. Baeriswyl Phys. Rev. B69, 094431 (2004).[10] B. Lake, A. M. Tsvelik, S. Notbohm, et al . NaturePhysics , 50 (2010).[11] M. Fabrizio, A. O. Gogolin and A. A. Nersesyan, Phys.Rev. Lett. , 2014 (1999).[12] M. Fabrizio, A. O. Gogolin and A. A. Nersesyan,Nucl.Phys. B , 647, (2000).[13] A. O. Gogolin, A. A. Nersesyan and A. M. Tsvelik,”Bosonization and Strongly Correlated Systems”; Cam-bridge University Press, Cambridge, 1999; A. M. Tsvelik,”Quantum Field Theory in Condensed Matter Physics”,2nd edition, Cambridge University Press, Cambridge,2003.[14] T. Giamarchi ”Quantum Physics in One Dimension”,Oxford University Press, Oxford, 2004.[15] R. Jackiw and C. Rebbi, Phys. Rev. D , 3398 (1976).[16] A. O. Gogolin and A. A. Nersesyan, A. M. Tsvelik andLu Yu, Nucl. Phys. B , 705 (1999).[17] For the Heisenberg chains such operators are continuumlimits of the staggered spin energy densities ǫ + ± ǫ − ∼ ( − n ( S n S n +1 ) , .[18] D. G. Shelton and A. M. Tsvelik, Phys Rev. B , 14 242(1998).[19] X.-Y. Feng, G.-M. Zhang and T. Xiang, Phys. Rev. Lett.98