Spin-1 Kitaev model in one dimension
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov Spin-1 Kitaev model in one dimension
Diptiman Sen , R. Shankar , Deepak Dhar and Kabir Ramola Center for High Energy Physics, Indian Institute of Science, Bangalore 560 012, India The Institute of Mathematical Sciences, CIT Campus, Chennai 600 113, India Department of Theoretical Physics, Tata Institute of Fundamental Research, Mumbai 400 005, India (Dated: October 15, 2018)We study a one-dimensional version of the Kitaev model on a ring of size N , in which there is aspin S > / J P n S xn S yn +1 . The cases where S is integer andhalf-odd-integer are qualitatively different. We show that there is a Z valued conserved quantity W n for each bond ( n, n + 1) of the system. For integer S , the Hilbert space can be decomposedinto 2 N sectors, of unequal sizes. The number of states in most of the sectors grows as d N , where d depends on the sector. The largest sector contains the ground state, and for this sector, for S = 1, d = ( √ /
2. We carry out exact diagonalization for small systems. The extrapolation of ourresults to large N indicates that the energy gap remains finite in this limit. In the ground statesector, the system can be mapped to a spin-1/2 model. We develop variational wave functions tostudy the lowest energy states in the ground state and other sectors. The first excited state ofthe system is the lowest energy state of a different sector and we estimate its excitation energy.We consider a more general Hamiltonian, adding a term λ P n W n , and show that this has gaplessexcitations in the range λ c ≤ λ ≤ λ c . We use the variational wave functions to study how theground state energy and the defect density vary near the two critical points λ c and λ c . PACS numbers: 75.10.Jm
I. INTRODUCTION
In recent years, there have been many studies of quan-tum spin systems which are characterized by a high de-gree of frustration and topological order. The word ‘frus-tration’ here refers to systems with competing interac-tions having a large number of states with energy nearthe minimum energy. Topological order implies the ex-istence of invariants which, for topological reasons, arerobust against a large class of perturbations. Such sys-tems are often associated with a novel structure of theground state and low-lying excitations, and are interest-ing from the point of view of possible applications inquantum computation [1–5]. A particularly interestingmodel in this context is the two-dimensional frustratedspin-1/2 model introduced by Kitaev [3]. This modelhas several fascinating properties which have been stud-ied in great detail [6–12]. For instance, the model and itsvariants constitute the only known class of spin modelsin two dimensions or more dimensions that is fully in-tegrable, being reducible to a system of non-interactingMajorana fermions. A similar model, called the compassmodel, although not exactly solvable, was introduced byKugel and Khomskii many years ago [13] to understandthe magnetic properties of transition metal oxides whichhave orbital degeneracies. Recently physical realizationsof the spin-1/2 Kitaev model have been proposed in op-tical lattice systems [14] and in quantum circuits [15].Variants of the model have also been studied in two di-mensions [16–24], three dimensions [25, 26] and also onquasi-one-dimensional lattices [27–29]. Finally, the spin- S Kitaev model has been studied in the large S limitusing spin wave theory [30], and the classical version ofthe Kitaev model has been studied at finite temperaturesusing analytical and Monte Carlo techniques [31]. Their results indicate that while the phenomenon of order-by-disorder [32–35] may occur in the quantum mechanicalKitaev model, it does not in the corresponding classicalmodel.For the Kitaev model with spin S > /
2, there is a Z invariant associated with each plaquette for arbitraryspin- S , which reduces to the conserved Z gauge fluxfor the spin-1/2 case [30]. However, the model does notseem to be fully integrable. While some differences inthe structure of the invariants between the models withhalf-odd-integer and integer spins have been pointed out[30], the issue of whether there are systematic differencesin the nature of the low-energy spectrum is also of in-terest. In the present paper, we approach this prob-lem by examining the spin-1 Kitaev model. The two-dimensional model appears difficult to analyze, but eventhe one-dimensional version of it has a lot of interestingstructure, as we proceed to show.The plan of this paper is as follows. In Sec. II, weconsider the spin- S Kitaev chain. In Sec. II A we showthat this model has local, mutually commuting conservedquantities W n , for integer S . The eigenvalues of W n are ±
1. For open boundary conditions, there are some ad-ditional conserved quantities at the ends of the system.The existence of these conserved quantities implies thatthe Hilbert space of a N -site system can be decomposedinto a sum of 2 N disjoint subspaces. The dimensions ofthese subspaces are not equal. In Sec. II B we developa formalism to compute the dimension of these sectors.For large N , the dimension varies as d N in most sectors,with the constant d depending on the sector. The sec-tors show complicated spatial structures, arising from thespatial structure of { W n } . We show this in Sec. II C, bycomputing the non-trivial spatial dependence of expec-tation values of spin operators in some sectors, averagedover all states in the sector. We then consider the spin-1model in Sec. III. In Sec. III A, we consider the groundstate and lowest excited state of the system. Exact diag-onalizations of small systems show that the ground statelies in a sector in which W n = +1 for all n . In this sec-tor, there is a gap between the ground state and the firstexcited state. The lowest excited state of the system isthe ground state of a different sector; and the energy-gap seems to approach a non-zero value in the limit ofthe system size going to infinity. In Sec. IV, we con-sider the sector containing the ground state, and showthat the Hamiltonian is equivalent to the Hamiltonian ofa deposition-evaporation process of a nearest-neighbor-exclusion lattice gas model, which can be written as of aspin-1 / λ P n W n ,and discuss its ground states as a function of λ . We showthat the ground state of this new Hamiltonian is gaplessfor a range of couplings λ c ≤ λ ≤ λ c , and gapped oth-erwise. We argue that for λ just above λ c , in the sectorcontaining the ground state, the density of negative W ’sis of order | λ − λ c ) | . For λ just below λ c , the densityof positive W ’s goes to zero as ( λ c − λ ) / . In the finalsection, we summarize our conclusions, and discuss therelationship of this model with the Fibonacci chain. II. ONE-DIMENSIONAL KITAEV MODEL
In this section, we will discuss a one-dimensional spin- S model which is obtained by considering a single row ofthe Kitaev model in two dimensions.Let us begin with the Kitaev model on the honeycomblattice. This is governed by the Hamiltonian H hexKit = J x X h ij i x S xi S xj + J y X h ij i y S yi S yj + J z X h ij i z S zi S zj , (1)where h ij i a denote the nearest-neighbor bonds in the a th direction. If we set J z = 0, we get a set of decoupledchains. We call this the Kitaev chain, and this is thetopic of this paper. The Hamiltonian is by H = X n (cid:0) J n − S x n − S x n + J n S y n S y n +1 (cid:1) . (2)In general, the couplings J m could be all different fromeach other. If some of the couplings are negative, we canchange the signs of those couplings by performing theunitary transformation S xm → − S xm , S ym → − S ym , and S zm → S zm (3) n n + 1 n + 2 W n FIG. 1: Picture of the Kitaev chain showing one of the con-served quantities W n . on appropriate sites. We consider the simpler case, whereall couplings have the same value, J m = J . Without anyloss of generality, we set J = 1. Finally, the Hamiltoniancan be unitarily transformed to a more convenient formby the following transformation on the even sites, S x n → S y n , S y n → S x n , and S z n → − S z n . (4)The Hamiltonian in Eq. (2) then takes the translationinvariant form H Kit = X n S xn S yn +1 . (5) A. Invariants
The Hamiltonian in Eq. (5) has the following localsymmetries for all S . Let us introduce the operators onsites Σ an = e iπS an , (6)and operators on bonds W n = Σ yn Σ xn +1 . (7)as shown in Fig. 1. We then find that[ W n , H ] = 0 . (8)The eigenvalues of Σ an are ± S and ± i forhalf-odd-integer S . Thus for any value of the spin S , theeigenvalues of W n are ± S . For integer values of S , all the matrices Σ an matrices commute with each other,whereas for half-odd-integer values, Σ an commutes withΣ bm for n = m but anticommutes with Σ bn for a = b . Con-sequently, for integer S , all the invariants W n commute,but for half-odd-integer S , W n anti-commutes with itsneighboring invariants, W n ± , and commutes with W m , m = n, n ±
1. We will now show that this implies that allthe eigenstates of the chain with half-odd-integer S are2 N/ fold degenerate.The invariants for half-odd-integer S can be combinedin the following way to form a set of mutually commutingangular momentum operators, one per every two bonds, µ zn = W n , µ xn = W n − Y m 1. We now add a site r + 1 to thechain, and also specify W r . Let the new set of { W } bedenoted by W ′ .Consider first the case W r = +1. Clearly, we can havetwo possibilities: Σ xr +1 = Σ yr = +1, or Σ xr +1 = Σ yr = − ν ( p, p ′ ) denote the number of states of a single sitewith Σ y = p , and Σ x = p ′ . Then, we clearly have therecursion equation Z r +1 ( y |W ′ ) = ν ( y, +1) Z r (+1 , W ) + ν ( y, − Z r ( − , W ) . (22)This equation can be written as a matrix equation (cid:20) Z r +1 (+1 |W ′ ) Z r +1 ( − |W ′ ) (cid:21) = T + (cid:20) Z r (+1 |W ) Z r ( − |W ) (cid:21) , (23)where T + is a 2 × T + = (cid:20) ν (+1 , +1) ν (+1 , − ν ( − , +1) ν ( − , − (cid:21) , (24)It then follows from Eqs. (18-20) that T + = 12 (cid:20) S − S + 1 S + 1 S + 1 (cid:21) for S odd , (25)= 12 (cid:20) S + 2 SS S (cid:21) for S even . (26)Similarly, when W r = − 1, the corresponding recursionequation is (cid:20) Z r +1 (+1 |W ′ ) Z r +1 ( − |W ′ ) (cid:21) = T − (cid:20) Z r (+1 |W ) Z r ( − |W ) (cid:21) , (27)where the matrix T − is given by T − = T + τ x , with τ x = (cid:20) (cid:21) . (28)It is then clear that for a given set of invariants W , thenumber of states can be written in terms of a product ofthe matrices T + and T − .For example, for an open chain of N sites, and W = { W N − , ...W , W , W } = { +1 , ... + 1 , − , − } , we have (cid:20) Z N (+1 | ... + −− ) Z N ( − | ... + −− ) (cid:21) = T + ... T + T − T − (cid:20) Z (+1 | φ ) Z ( − | φ ) (cid:21) , (29)where φ denotes the null string, and Z ( y | φ ) denotes thenumber of states of the spin at site 1 with Σ y = y . Thus Z (+1 | φ ) = S + 1, Z ( − | φ ) = S , when S is an eveninteger, and Z (+1 | φ ) = S , Z ( − | φ ) = S + 1 when S isan odd integer. The total number of states in this sectoris then given byΓ( W ) = Z N (+1 |W ) + Z N ( − |W ) . (30)For a closed chain, there is an additional invariant W N = y N x and the number of states in the sector be-comes Γ( W ) = Tr N Y n =1 T W n ! , (31)where T W n ≡ T ± for W n = ± Q Nn =1 is an orderedproduct of T ± matrices, from site 1 to N with the indexincreasing from right to left.We now calculate the dimensions of some sectors fora closed chain of length N . It is easy to get an explicitanswer for the two extreme limits when W n = ± n . In these cases, the number of states, Γ ± , isΓ ± = (cid:0) d ± (cid:1) N + (cid:0) d ± (cid:1) N , (32)where d ± ( S ) and d ± ( S ) are the larger and smaller eigen-values of T ± respectively. The eigenvalues can be com-puted to give, d +1(2) = 12 (cid:16) S ± p S + 2 S + 2 (cid:17) for S odd , (33)= 12 (cid:16) S + 1 ± p S + 1 (cid:17) for S even , (34) d − = 12 (cid:16) S + 1 ± p S − (cid:17) for S odd , (35)= 12 (cid:16) S ± p S + 2 S (cid:17) for S even . (36)For S = 1, d +1 is equal to the golden ratio, γ = (1 + √ / 2, and d +2 = − /γ . As N → ∞ , the dimension ofthe Hilbert space in the sector with all W n = 1 grows as γ N . On the other hand, d − = d − = 1. The dimension ofthe sector with all W n = − S = 1, the larger of the twoeigenvalues d ± is always greater than 1, and in the N →∞ limit, we have Γ ± ( S ) = (cid:0) d ± ( S ) (cid:1) N . (37) d ± ( S ) is referred to as the quantum dimension of thesector. As can be seen it is, in general, fractional for any S . In the limit S → ∞ , the quantum dimension tendsto S + 1 / C. Expectation values of the Σ operators indifferent sectors In this section we find the expectation values of theΣ an operators in various sectors. We will assume peri-odic boundary conditions. Our calculation will averageover all the states of a given sector considered with equalweight; this can be considered as a calculation in the limitthat the temperature T → ∞ , so that it does not dependon the Hamiltonian.We evaluate the expectation values of Σ an by insertingprojection operators at site n in the product of transfermatrices in Eq. (31). This yields the following expressionfor the expectation value of the Σ an operator in a generalsector with a W -configuration Wh Σ an +1 i W = Tr N Y j = n +1 T W j T aW n n − Y i =1 T W i ! / Γ( {W} )where T xW n = W n T W n τ z T yW n = τ z T W n T zW n = W n τ z T W n τ z (38)and τ z and τ x are the well-known Pauli matrices.We now compute the expectation values of Σ an in twosectors: the sector W with all W n = +1, and the sector W in which one of the W n = − W n =+1 (without loss of generality we pick W N = − h Σ an i W , ≡ h Σ an i , can be evaluated interms of the eigenvectors and eigenvalues of T + . The T + matrix is a linear combination of the Pauli matrices, τ z and τ x . Its eigenvectors are spinors polarized paralleland anti-parallel to a direction in the z x plane, formingan angle θ S with the z − axis, where θ S is defined by,cos θ S ≡ − p S + 1) for S odd , ≡ √ S for S even , (39)sin θ S ≡ S + 1 p S + 1) for S odd , ≡ S √ S for S even . (40)For the sector with all W n = +1 it is easy to see that h Σ xn i = h Σ yn i . For large N , we obtain h Σ x ( y ) n i = cos θ S , (41) h Σ zn i = cos θ S + d +2 d +1 sin θ S . (42)In the sector where W N = − + + + + − − − − + − − + + − + − n < ( S na ) > FIG. 2: Plot of h S xn i (dotted line), h S yn i (dashed line) and h S zn i (full line) as a function of n for S = 1, on a ring with16 sites in the sector + + + + − − − − + − − + + − + − , withperiodic boundary conditions (site 17 = site 1). to +1, we get, for large N , h Σ xn i = h Σ xn i − (cid:18) d +2 d +1 (cid:19) n − ! , (43) h Σ yn i = h Σ yn i − (cid:18) d +2 d +1 (cid:19) n ! , (44) h Σ zn i = h Σ zn i θ S d +1 cos θS + d +2 sin θ S (cid:18) d +2 d +1 (cid:19) n − ! . (45)Note that in the limit n → ±∞ , h Σ an i in Eqs. (43-45)approach the values given in Eqs. (41-42) exponentiallyquickly.While in general Σ an are complicated multi-spin opera-tors, for S = 1 we have Σ an = 1 − S an ) . Thus, for S = 1we are essentially computing the expectation values of( S an ) . To see what the spin textures are like in a typicalsector, we have plotted in Fig. 2 the expectation valuesof S xn , S yn and S zn for S = 1, as a function of the spatialcoordinate n for a ring of size 16, in the sector where thesequence of W ’s is + + + + − − − − + − − + + − + − .This sequence was chosen as it is a de Bruijn sequence[36] of length 16, in which each of the 16 possible bi-nary sequences of length 4 occur exactly once, taking theperiodic boundary conditions into account. III. S = 1 MODEL We will now focus on the Kitaev chain with spin-1’s ateach site. We will work with the natural spin-1 represen-tation in which ( S a ) bc = iǫ abc . (46) In this representation, the matrices Σ a are diagonal andare given by Σ x = − − , Σ y = − − , Σ z = − − . (47)We note that these matrices satisfy Σ x Σ y Σ z = I . Wedenote the basis vectors by | x i , | y i and | z i defined as | x i = , | y i = , | z i = . (48)We then see that the 9 possible states at sites ( n, n +1)are given by: | xy i , | xz i , | yx i , | zy i and | zz i with W n = 1 , (49)and | xx i , | yy i , | yz i and | zx i with W n = − . (50)From Eq. (46) we have, S x | x i = 0 , S y | x i = i | z i , S z | x i = − i | y i ,S x | y i = − i | z i , S y | y i = 0 , S z | y i = i | x i ,S x | z i = i | y i , S y | z i = − i | x i , S z | z i = 0 . (51)Eqs. (47) and (51) imply that ( S a ) = (1 − Σ a ) / W n = 1, wehave the following actions of the relevant term in theHamiltonian, S x S y | xy i = 0 ,S x S y | xz i = 0 ,S x S y | zy i = 0 ,S x S y | zz i = | yx i ,S x S y | yx i = | zz i . (52)For the 4 states in Eq. (50) satisfying W n = − 1, theactions of the relevant term in the Hamiltonian are givenby S x S y | xx i = 0 ,S x S y | yy i = 0 ,S x S y | yz i = −| zx i ,S x S y | zx i = −| yz i . (53)As mentioned earlier, for an open chain with site num-bers going from 1 to N , we find that S x and S yN commutewith H . We define the operators, τ ≡ iW S x , τ ≡ S x , (54) τ ≡ − S x W S x , τ ≡ (cid:0) − ( S x ) (cid:1) . (55)It can be verified that these operators obey a SU (2) × U (1) algebra. Exactly the same construction on the lastbond, with S x → S yN and W → W N , yields the samealgebra on that bond. A. Numerical studies We have carried out exact diagonalization studies ofsmall systems with periodic boundary conditions in or-der to find the energies of the ground state and the lowestexcited state of the spin-1 Kitaev chain. We find that theground state lies in the sector with all W n = 1 and haszero momentum (momentum is a good quantum numberin this sector since the values of the W n ’s are translationinvariant). The ground state energy per site as a func-tion of the system size N is presented in Table I. We seethat E /N shows odd-even oscillations as a function of N but seems to converge quite fast. The fast convergenceindicates that the ground state must have a fairly shortcorrelation length. The N -dependence of ¯ E N = E /N can be fitted to the form¯ E N = E ∞ + B ( − α ) N . (56)A simple plot of log | ¯ E N − E ∞ | versus N (Fig. 3), gives agood straight line for E ∞ = − . E ∞ . The corresponding valuesof B and α are 0 . 07 and 0 . 51. The estimated errors ofextrapolation are about 1 in the last significant digit. N E /N N E /N N . In the sector with all W n = 1, the first excited statehas momentum equal to π if N is even. We find thatthe gap separating it from the ground state is given by1 . N = 4 and 0 . N = 6. These valuesalso seem to be converging rapidly, and the large value isconsistent with a short correlation length. However, thisis not the lowest excited state of the system. Rather,we find that the state nearest in energy to the groundstate is the ground state of the sector with exactly one W n = − W n = 1. (We cannot usemomentum to classify the states in this sector since it isnot translation invariant). The energy gap ∆ E betweenthe lowest energy state in this sector and the ground stateof the sector with all W n = 1 is shown in Table II. We -7 -6 -5 -4 -3 -2 -1 1 2 4 6 8 10 12 14 16 18 20 | E / N - E ∞ | N 0.07 (0.51 N ) FIG. 3: Graph of | E /N − E ∞ | with N , where E ∞ = − . N . see that these also oscillate between even and odd valuesof N but seem to converge quite fast to a small but non-zero value. This is evidence that the spin-1 Kitaev chainhas a finite gap in the thermodynamic limit N → ∞ . N ∆ E N . IV. MAPPING THE SPIN-1 CHAIN TO ASPIN-1/2 CHAIN For a given value of the state of the spin at site n ,and a given value of W n , there are at most two choicesfor the spin state at site n + 1. Hence it is clear thatthe Hilbert space of a given sector can be mapped intothe Hilbert space of a spin-1/2 chain, with some statesexcluded which correspond to infinite energy. However,in general, the corresponding Hamiltonian would have arather complicated form, with long-ranged interactions.The mapping is easy to construct explicitly in the sectorwith all W n = +1, and the corresponding Hamiltonianhas only local interactions. This is what we now proceedto show.Consider the state zzzz · · · that belongs to the sectorwith all W n = +1. The only allowed process in this sectoris zz ⇋ yx [Eqs. (52)]. We may think of this process asa quantum dimer deposition-evaporation model. The z -spins are treated as empty sites; two empty sites can bechanged to being occupied by a dimer yx by a ‘deposition’process, and conversely, yx can ‘evaporate’ and become zz again. The dimers have a hard-core constraint, anda site cannot be shared by two dimers. The dimers areoriented: the ‘head’ x being to the right of the ‘tail’ y .This dimer deposition-evaporation model can alsobe described as a deposition-evaporation of a nearest-neighbor exclusion lattice gas. We just think of the headsas particles, and do not distinguish between the tails andempty sites, except for ensuring that we deposit a par-ticle at a site only if it is empty and both its nearestneighbors are also empty. Then this model is describedby the Hamiltonian H d = − X n (1 − σ zn − ) σ xn (1 − σ zn +1 ) . (57)We note that this model is different from the dimerdeposition-evaporation models studied earlier [37], inthat the two ends of the dimer are distinct, and there isno reconstitution. Also, this Hamiltonian does not havean interpretation as the evolution operator of a classicalMarkov process, as there are no diagonal terms corre-sponding to probability conservation.We have introduced a minus sign in the Hamiltonianfor later convenience. This does not change the eigen-value spectrum as the eigenvalues of H d occur in pairs ± e i . V. VARIATIONAL STUDY OF SECTOR WITHALL W n = 1 We will now use a variational approach to study theground state of the Hamiltonian H d with periodic bound-ary conditions. We use the z -basis, and denote the ↑ stateat the site i by an occupied site ( n i = 1), and the ↓ stateby an empty state ( n i = 0). Since two adjacent sites can-not be simultaneously occupied, the state space is thatof hard-core particles with nearest-neighbor exclusion ona line. A configuration C is specified by an N -bit binarystring 0010010101 · · · , which gives the values of all the N occupation numbers n i . We note that in the basis whereall the n i are diagonal, the Hamiltonian H d has all ma-trix elements non-positive. This implies that the (real)eigenvector corresponding to the lowest energy will haveall components of the same sign in this basis.For the ground state of H d , we consider a variationalwave function of the form | ψ i = X C p Prob( C ) | C i , (58)where Prob( C ) is chosen as the probability of the latticegas configuration C in some classical equilibrium ensem-ble corresponding to a suitably chosen lattice gas Hamil-tonian. Clearly, this trial vector is normalized, with h ψ | ψ i = 1 . (59)With this choice, Prob( C ) is also the probability of theconfiguration C in the quantum mechanical variationalstate | ψ i . The simplest choice of the lattice-gas Hamiltonian isthat of a classical lattice gas with nearest-neighbor ex-clusion, and a chemical potential µ , with a Hamiltoniangiven by H cl = + ∞ X i n i n i +1 − µ X i n i , (60)where we use the convention that 0 · ∞ = 0; hence thefirst term in Eq. (60) allows states with n i n i +1 = 0but disallows states with n i n i +1 = 1. Let us denote z = exp( βµ ). It is straightforward to determine variouscorrelation functions in the thermal equilibrium state cor-responding to H cl . The probability of a configuration C is given byProb( C ) = exp[ − βH cl ( C )] / Ω N ( z ) , (61)where Ω N ( z ) is the grand partition function for a ring of N sites.The grand partition function Ω N ( z ) can be determinedusing the standard transfer matrix technique. We findthe largest eigenvalue of the 2 × T given by T = (cid:20) z (cid:21) . (62)We now calculate h ψ | H d | ψ i . The matrix element of the i -th term is clearly zero, unless n i − = n i +1 = 0. Thenthe only non-zero matrix element is h H d i /N = − √ z Prob(000) = − √ z Prob(010) . (63)Here Prob(000) denotes the probability that randomlyselected three consecutive sites in the ring will be emptyin the classical ensemble, and similar definition forProb(010). This is easily calculated for the Hamiltonian H cl in the limit of large N . We getProb(010) = Prob(1) = ρ. (64)The largest eigenvalue Λ of T is given byΛ = (1 + √ z ) / , (65)and ρ is the density per site given by ρ = zd log(Λ) /dz .Extremizing h H d i with respect to z , we find that theminimizing value occurs for z = 0 . h H d i = − . E E ≤ − . . (66)This energy is somewhat higher than the energy obtainedin the previous section (see Fig. 3), indicating that thecorrelations in the classical Hamiltonian H cl do not ex-actly reproduce the correlations in the quantum groundstate of H d .We can make a better variational calculation by con-sidering a classical lattice gas with an additional next-nearest-neighbor interaction. The Hamiltonian of thislattice gas is H ′ cl = + ∞ X i n i n i +1 − K X i n i n i +2 − µ X i n i . (67)Let us denote z = exp( βµ ), and u = exp( βK ). In thiscase, the transfer matrix is a 3 × T = z zu . (68)The probability of the configuration C in the equilibriumensemble is given byProb( C ) = exp[ − βH ′ cl ( C )] / Ω N ( z, u ) , (69)where Ω N ( z, u ) is the grand partition function for a ringof N sites. We then get − h H d i /N = 2Prob(00000) √ z + 4Prob(10000) √ zu +2Prob(10001) √ zu , (70)Here Prob(00000) is the probability of finding a randomlyselected set of five consecutive sites all unoccupied inthe equilibrium ensemble corresponding to the Hamilto-nian H ′ cl . These probabilities are also easily calculated.Treating z and u as variational parameters, we find that h H d i is minimized for z = 0 . u = 1 . ρ = 0 . . . . E ≤ − . − . VI. STUDY OF GROUND STATES IN OTHERSECTORS We define a more general Hamiltonian H ( λ ) = H Kit + λ X n W n . (72)Since the W n ’s commute with H Kit , all the eigenvectorsof H Kit can be chosen to be simultaneous eigenvectorsof H ( λ ), for all λ . However, if we vary λ , we can getdifferent eigenvectors to have the lowest energy.Clearly, if λ is large and positive, the ground state willlie in the sector with all W n = − 1. Conversely, if λ islarge and negative, the ground state is the lowest energy eigenvector in the sector with all W n = +1. In boththese regions, the gap in the excitation spectrum is oforder | λ | . As we vary λ from −∞ to + ∞ , initially thegap decreases and becomes zero at some value λ c . Wethen expect a gap to open up again when λ is greaterthan a second critical point λ c > λ c . A. Sectors with most W n ’s positive Since the ground state for λ = 0 lies in the sector withall W n = +1, we have λ c > 0. In fact, if the lowestexcitation energy in the Hamiltonian H Kit is ∆ E , wehave λ c = ∆ E/ 2. At this point, the energy required tochange a single W n from +1 to − W n = 1.Without loss of generality, we may assume that in thissector, W N = − 1, and the rest of the W ’s are +1. Thebasis vectors in this sector are of type | xU i , or | V y i ,where U and V are all possible strings of length N − zzz . . . z of length N − 1, usingthe substitution rule zz → yx . Let | ψ i be the eigenvectorcorresponding to the lowest eigenvalue of H in this sec-tor. It is easy to verify that h xU | H | xU ′ i and h V y | H | V ′ y i are negative, for all U and U ′ , and V and V ′ . But, h V y | H | xU i are positive. This implies that h xU | ψ i and h xU ′ | ψ i have the same sign for all U and U ′ . Similarly h V y | ψ i and h V ′ y | ψ i have the same sign for all V and V ′ .This suggests a variational wave function of the form | ψ i = 1 √ "X U p Prob( U ) | xU i − X V p Prob( V ) | V y i . (73)Here Prob( U ) and Prob( V ) are arbitrary functions,satisfying the constraint X U Prob( U ) = X V Prob( V ) = 1 . (74)Each configuration U is in one-to-one correspondencewith the configurations of a nearest-neighbor-exclusionlattice gas on a linear chain of length ( N − C = { n i } . We put n i = 1 if andonly if there is a y in U in the position i + 1, otherwise n i = 0. Note that the last element of U cannot be a y .We specify C by a binary string of length ( N − C to U , we first add a single 0 to the binary stringof C at the right end, and then use the substitution rule10 → yx . The remaining zeros in C are replaced by z ’s.Similarly, we specify V also by a binary string of length( N − x → y, z → 0, and as the leftmostelement of the resulting string is always a zero, it may bedeleted.As in the previous calculation, we construct a classi-cal Hamiltonian to variationally estimate the parametersProb( C ). In this case, there is no translational symme-try, and in general, the lattice gas will have a non-trivialdensity profile. This is taken into account by making theactivities of the lattice gas in the classical Hamiltoniansite-dependent. We write H Ccl = + ∞ N − X i =1 n i n i +1 − N − X i =1 µ i n i . (75)The probability of each configuration C of the lattice gasis then given byProb(C) = exp[ − βH Ccl ( C )] / Ω N − ( { z i } ) , (76)with z i = exp( βµ i ), and Ω N − ( { z i } ) is the grand parti-tion function of the open chain of N − H is unchanged under thespace reflection i ↔ N + 1 − i , and at the same time ex-changing x and y . This can be built into our eigenvectorby assuming that if V are strings corresponding to latticegas configurations C , we setProb( V ) = Prob( U ) , (77)where U is the string corresponding the lattice gas con-figuration C T , the transpose of C .The rest of the calculation is done as before. By con-struction, we have h ψ var | ψ var i = 1 . (78)It is straightforward to express h ψ | H | ψ i in terms of themarginal probabilities of the different local configurationsof the lattice gas, remembering that there is no transla-tional invariance. For example, we get h ψ | S xN S y | ψ i = − Prob C ( n = 0) . (79)In the simplest case, we work with only two parame-ters, and set z = z ′ , and z i = z for i = 1. We wouldlike to estimate the difference of the ground state en-ergy in this sector and the ground state over all sectors.These energies are of order N , and to cancel the lead-ing linear N -dependence, we have to set z equal to theoptimal value z ∗ = 0 . N is large, so that onlythe term in the partition function corresponding to thelargest eigenvalue is kept. Extremizing over z ′ we obtain z ′ = 0 . h ψ var | H | ψ var i = − . N + 0 . . (80)This implies the following bound on the lowest eigenvaluein this sector E ′ ≤ N E + ∆ . (81)with ∆ = 0 . z ′ and z ′′ atthe two opposite ends. Extremizing with respect to theseparameters we find the energy gap to be 0 . z ′ =0 . z ′′ = 0 . z , z , z N − and z N − adjustable,and the rest of the z i ’s set equal to z ∗ . Table III showsthe improvement in the value of the energy gap with thenumber of parameters used.We thus obtain a variational estimate of the energygap of the first excited state from the ground state en-ergy. This matches quite well with the numerical esti-mates obtained in the previous section. Number of ∆parameters1 0.187512 0.164194 0.158456 0.156428 0.1557810 0.15556TABLE III: Estimate of the energy gap with the number ofparameters used. It is straightforward to extend this treatment to sectorswith two or more W n ’s negative. There is an energy ∆required to create a single negative W n . Thus λ c = ∆ / n defects, the energy would beminimized if the defects are equally spaced. Thus thedistance between the defects is N/n , and the energy costof creating n defects ∆ E ( n ) in H ( λ ), for small n , is wellapproximated by∆ E ( n ) ≈ − nλ + n ∆ + nA exp( − BN/n ) , (82)where A and B are some constants. This then impliesthat for λ = λ c + ǫ , the density of defects in the trueground state of H ( λ ) will vary as 1 / | log ǫ | . B. Sectors with most W n ’s negative We now discuss the behavior of the ground state energynear the critical point λ c . This depends on the behaviorof the ground state energy in sectors in which only a fewof the W n ’s are +1.For a ring of N sites, the sector with all W n = − xxxxx · · · and yyyyyy · · · . Thetwo are degenerate, with eigenvalue equal to − λN .Now consider the sector with only one W n = +1, say W = +1. Consider the state ψ = | zxxxx · · · i in thissector. From Eqs. (53), under H Kit , we have zx ⇋ yz ,0and this state can make a transition only to the state ψ = | yzxxx · · · i . And ψ can return to ψ or go to ψ = | yyzxx · · · i . Thus, the dynamics may be consideredas the dynamics of a particle z , which can hop to a nearestneighbor under the action of the Hamiltonian. Thereis a string of y ’s connecting the current position of theparticle to the leftmost allowed position which is n = 1.This string can become longer, or shorter, as the particlemoves, with no energy cost. When the z -spin is at the site N , it cannot move further to the right. The ground stateenergy E g − sector if this sector is seen to be the same asthat of a particle with nearest-neighbor hopping, confinedto move in the space 1 ≤ x ≤ N . It is thus given by E g − sector = − ( N − λ − J cos( πN + 1 ) . (83)Thus, we see that for large N , the state with all W n ’sequal to +1 is no longer the ground state for λ < J .We now consider a sector with exactly two of the W n ’sequal to +1, and the rest negative. Let us start withthe state | zxxx · · · zxxx · · · i , where the spins at two sites i = 1 and i = m + 1 ≤ N are in the state z (thesestates will be referred to as z -spins in the following). Thiscorresponds to W N = W m = +1. Then, under the actionof H ( λ ), this state mixes with other states where thepositions of the z -spins can change; the general state inthis sector may be labeled by the positions of the z -spins, r and r . We will write the vector as | r , r i , where1 ≤ r ≤ m < r ≤ N . Then, for 1 < r < m and m + 1 < r < N , we get H Kit | r , r i = −| r , r + 1 i − | r + 1 , r i−| r , r − i − | r − , r i . (84)If the first z -spin is at m and the second is not at m + 1, the first spin cannot move to m + 1, as that sitewould be in spin state y , and the state zy cannot change[Eqs. (52)]. Similarly, if r = N , and r = 1, then thesecond spin cannot move to the right. However, if thetwo z -spins are adjacent, then they can change to a state zz ⇋ yx [Eqs. (52)]. But from the state yx the state canonly return to zz .If we disallow the transitions to state yx , the z -spinsact as independent particles moving in two disjoint re-gions of space, 1 ≤ r ≤ m and m + 1 ≤ r ≤ N . Inthis case, the minimum energy of this system is just thesum of the energies of two particles. This energy is anupper bound on the true ground state energy of this sys-tem. Thus, we find that the ground state energy in thissector, E g − sector , has the upper bound E g − sector ≤ − J cos( πm + 1 ) − J cos( πN − m + 1 ) − λ ( N − . (85)Next, suppose that the state with the m -th site in the y -state and the ( m + 1)-th in the x -state is called thestate r = m + 1 , r = 1, and a similar definition for theother end. Then the range of r is at most m + 1, and the range of r is at most l − m + 1. By excluding somestates (here r = m + 1 , r = r ), the kinetic energy canonly increase, and hence we have E g − sector ≥ − J cos( πm + 2 ) − J cos( πN − m + 2 ) − λ ( N − . (86)For N, m ≫ 1, these bounds can be expanded in pow-ers of 1 /m , and have the same leading order correction.Also, the minimum energy corresponds to equally spaceddefects, with m = N/ W n ’s equal to +1. In case thelengths of the intervals between the positive W n ’s are m , m , m , · · · , m r , the bounds on the lowest energy inthis sector E gr − sector become − J r X i =1 cos( πm i + 2 ) − λ ( N − r ) ≤ E gr − sector ≤ − J r X i =1 cos( πm i + 1 ) − λ ( N − r ) . (87)Thus, we see that for λ > J , the ground state belongsto the sector with all W n ’s equal to − 1. If λ = J (1 − ǫ ),the ground state will be in the sector with n equispacedbonds with W n = +1, where the spacing ℓ between them ≈ N/n is given by ǫ − / . The minimum energy per siteof H ( λ ) for λ = J (1 − ǫ ) varies as ǫ / for small ǫ . Equiv-alently, if we restrict ourselves to sectors with only a frac-tion ǫ of W n ’s having the value +1, the minimum energyper site varies as − ǫ / . This is equivalent to the state-ment that for H Kit corresponding to λ = 0, in the sectorwith the fractional number of positive W n ’s being equalto ∆, the minimum energy per site varies as ∆ / . VII. DISCUSSION In this paper, we first analyzed the symmetries of aspin- S Kitaev chain. We found a Z invariant, W n , asso-ciated with every link ( n, n + 1), namely, N invariants forthe model defined on a ring with N sites. For integer S ,these invariants commute with each other and the Hamil-tonian. The Hilbert space can therefore be split into 2 N sectors, where the Hamiltonian is block diagonal. Forhalf-odd integer S , W n anti-commutes with W n ± andcommutes with the rest. We showed that this impliesthat all the eigenstates of the half-odd-integer spin mod-els are 2 N/ -fold degenerate, thus showing a qualitativedifference between the integer and half-odd-integer mod-els. We have developed a formalism to compute the di-mensions of the invariant sectors. We showed that the di-mension of most of the sectors can be calculated in termsof products of 2 × T + and T − . For S = 1 thequantum dimension of the sector with all W n = 1 is thegolden ratio, (1 + √ / 2. For S → ∞ , the quantum di-mension tends to S + 1 / W n = 1 and the W n = − λ , this would show gapless exci-tations in the range λ c ≤ λ ≤ λ c . We extended our vari-ational calculation to study how the ground state energyand the defect density would vary near the two criticalpoints λ c and λ c . At λ = λ c , Eq. (82) implies that theenergy of the lowest excited state in a system of length L goes as E ∼ exp( − BL ), corresponding to a state inwhich one W n = − W n = 1. By theusual scaling arguments, the gap to the first excited stategoes as 1 /L z , where z is the dynamical critical exponent.We therefore conclude that z = ∞ . At λ = λ c , the low-energy excitations form a low-density gas of hard-coreparticles. In one dimension, this can be mapped to asystem of non-interacting spinless fermions with a non-relativistic spectrum E ∼ k . Hence in a system of size L , the gap to the lowest energy states goes as 1 /L corre-sponding to k ∼ /L ; thus z = 2. It would be interestingto find the value of z in the critical region λ c < λ < λ c .Finally, we note that there is another interesting one-dimensional spin model called the golden or Fibonaccichain [38, 39], for which the number of states on a ringof size N is the same as that of the spin-1 Kitaev chainin the sector with all W n = 1. The Hamiltonian for thismodel is H GC = X i (( n i +1 + n i − − − n i − n i +1 ( γ − / σ xi + γ − n i + 1 + γ ) (cid:17) , (88)where n i = (1 − σ z ) / 2. It has been shown [38, 39] thatthis model is critical. Its long-range correlations are de- scribed by a SU (2) level 3 Wess-Zumino-Witten model,which is a conformally invariant field theory with cen-tral charge equal to 7 / 10. The Hamiltonian in Eq. (88)differs from the spin-1 Kitaev chain in the W n = 1 sec-tor by terms which are products of the n i operators. Wehave shown that the spin-1 Kitaev chain is gapped. Thusthese terms correspond to some relevant operators whichtake the golden chain Hamiltonian away from criticality.We can show that it is possible to add multi-spin termsto the minimal Kitaev chain which reduce to the extraterms in the W n = 1 sector. We need to add products ofthe n i operators to the minimal Kitaev chain to obtainthe golden chain in the sector with all W n = 1. The basisstates | ↑i , | ↓i that we use in Sec. IV are eigenstates ofthe n i operators with eigenvalues 1 and 0 respectively.The | ↑i state represents a state with the head, namely | x i . The | ↓i state represents either an empty site, | z i ,or a tail, | y i . It is clear from Eqs. (47) and (48) thatthe operator P x ≡ (1 + Σ x ) / | x i and 0 for | y i and | z i . Since all the Σ matrices commutefor integer S , they commute with the invariants and areblock diagonal within the invariant sectors. Thus, theHamiltonian, H KGC = γ − / H KC − X i (cid:0) − P xi +1 − P xi − + γ P xi − P xi +1 + γ − P xi − P xi P xi +1 (cid:1) , (89)when restricted to the W n = 1 sector, is exactly thegolden chain Hamiltonian discussed by Feiguin et al. andothers [38, 39]. We have thus constructed a realizationof the golden chain model as a spin-1 chain. Acknowledgments We thank G. Baskaran for interesting comments. DSthanks DST, India for financial support under ProjectNo. SR/S2/CMP-27/2006. DD thanks DST, Indiafor support through a J. C. Bose Fellowship underSR/S2/JCB-24/2006. [1] A. Kitaev, Ann. Phys. (N.Y.) , 2 (2003).[2] M. A. Levin and X.-G. Wen, Phys. Rev. B , 045110(2005).[3] A. Kitaev, Ann. Phys. (N.Y.) , 2 (2006).[4] S. Das Sarma, M. Freedman, and C. Nayak, Phys. Today (7), 32 (2006).[5] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S.Das Sarma, Rev. Mod. Phys. , 1083 (2008).[6] X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. , 087204 (2007).[7] G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev.Lett. , 247201 (2007).[8] D.-H. Lee, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. , 196805 (2007). [9] H.-D. Chen and Z. Nussinov, J. Phys. A , 075001(2008).[10] Z. Nussinov and G. Ortiz, Phys. Rev. B , 064302(2008).[11] K. P. Schmidt, S. Dusuel, and J. Vidal, Phys. Rev. Lett. , 057208 (2008).[12] K. Sengupta, D. Sen, and S. Mondal, Phys. Rev. Lett. , 077204 (2008); S. Mondal, D. Sen, and K. Sengupta,Phys. Rev. B , 045101 (2008).[13] K. I. Kugel and D. I. Khomskii, Sov. Phys. Usp. , 231(1982).[14] L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev.Lett. , 090402 (2003); A. Micheli, G. K. Brennen, andP. Zoller, Nat. Physics , 341 (2006). [15] J. Q. You, X.-F. Shi, X. Hu and F. Nori, Phys. Rev. B , 014505 (2010).[16] X.-G. Wen, Phys. Rev. D , 065003 (2003).[17] H. Yao and S. A. Kivelson, Phys. Rev. Lett. , 247203(2007).[18] S. Yang, D. L. Zhou, and C. P. Sun, Phys. Rev. B ,180404(R) (2007).[19] S. Dusuel, K. P. Schmidt, J. Vidal, and R. L. Zaffino,Phys. Rev. B , 125102 (2008).[20] G. Kells, N. Moran, and J. Vala, J. Stat. Mech.: TheoryExp. P03006 (2009).[21] H. Yao, S.-C. Zhang, and S. A. Kivelson, Phys. Rev. Lett. , 217202 (2009).[22] C. Wu, D. Arovas, and H.-H. Hung, Phys. Rev. B ,134427 (2009).[23] M. Kamfor, S. Dusuel, J. Vidal, and K. P. Schmidt, J.Stat. Mech.: Theory Exp. P08010 (2010).[24] X.-F. Shi, Y. Chen, and J. Q. You, Phys. Rev. B ,174412 (2010).[25] T. Si and Y. Yu, Nucl. Phys. B , 428 (2008).[26] S. Mandal and N. Surendran, Phys. Rev. B , 024426(2009).[27] V. Karimipour, Phys. Rev. B , 214435 (2009). [28] D. Sen and S. Vishveshwara, EPL , 66009 (2010).[29] A. Saket, S. R. Hassan, and R. Shankar, Phys. Rev. B , 174409 (2010).[30] G. Baskaran, D. Sen, and R. Shankar, Phys. Rev. B ,115116 (2008).[31] S. Chandra, K. Ramola, and D. Dhar, Phys. Rev. E ,031113 (2010).[32] J. Villain, R. Bidaux, J. P. Carton, and R. J. Conte, J.Phys. (Paris) , 1263 (1980).[33] R. Moessner and J. T. Chalker, Phys. Rev. B , 12049(1998).[34] C. L. Henley, Phys. Rev. B , 014424 (2005).[35] M. V. Gvozdikova and M. E. Zhitomirsky, JETP Lett. , 236 (2005).[36] See http://en.wikipedia.org/wiki/De Bruijn sequence[37] M. Barma, M. D. Grynberg, and R. B. Stinchcombe,Phys. Rev. Lett. , 1033 (1993).[38] A. Feiguin, S. Trebst, A. W. W. Ludwig, M. Troyer, A.Kitaev, Z. Wang, and M. H. Freedman, Phys. Rev. Lett. , 160409 (2007).[39] S. Trebst, M. Troyer, Z. Wang, and A. W. W. Ludwig,Prog. Theor. Phys. Supp.176