Quench dynamics and defect production in the Kitaev and extended Kitaev models
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Quench dynamics and defect production in the Kitaev and extended Kitaev models
Shreyoshi Mondal , Diptiman Sen and K. Sengupta TCMP division, Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700 064, India Center for High Energy Physics, Indian Institute of Science, Bangalore, 560 012, India (Dated: October 27, 2018)We study quench dynamics and defect production in the Kitaev and the extended Kitaev models.For the Kitaev model in one dimension, we show that in the limit of slow quench rate, the defectdensity n ∼ / √ τ where 1 /τ is the quench rate. We also compute the defect correlation function byproviding an exact calculation of all independent non-zero spin correlation functions of the model.In two dimensions, where the quench dynamics takes the system across a critical line, we elaborateon the results of earlier work [K. Sengupta, D. Sen and S. Mondal, Phys. Rev. Lett. , 077204(2008).] to discuss the unconventional scaling of the defect density with the quench rate. In thiscontext, we outline a general proof that for a d dimensional quantum model, where the quenchtakes the system through a d − m dimensional gapless (critical) surface characterized by correlationlength exponent ν and dynamical critical exponent z , the defect density n ∼ /τ mν/ ( zν +1) . We alsodiscuss the variation of the shape and the spatial extent of the defect correlation function with thechange of both the rate of quench and the model parameters and compute the entropy generatedduring such a quench process. Finally, we study the defect scaling law, entropy generation anddefect correlation function of the two-dimensional extended Kitaev model. PACS numbers: 73.43.Nq, 05.70.Jk, 64.60.Ht, 75.10.Jm
I. INTRODUCTION
Quantum phase transitions involve a fundamentalchange in the symmetry of the ground state of a quan-tum system. Such a transition usually takes place dueto the variation of some parameter λ in the Hamilto-nian of the system and is necessarily accompanied bydiverging length and time scales . A direct consequenceof such a diverging time scale is that a quantum systemfails to be in the adiabatic limit when it is sufficientlyclose to the quantum critical point. Thus a time evolu-tion of the parameter λ at a finite rate 1 /τ , which takessuch a system across a quantum critical point locatedat λ = λ c , leads to failure of the system to follow theinstantaneous ground state in a finite region around λ c .As a result, the state of the system after such a timeevolution does not conform to the ground state of its fi-nal Hamiltonian leading to the production of defects .It is well known that for a slow quench, the density ofthese defects n depends on the quench time τ accordingto n ∼ /τ dν/ ( νz +1) , where ν and z are the correlationlength and the dynamical critical exponents character-izing the critical point . A theoretical study of sucha quench dynamics requires a knowledge of the excitedstates of the system. As a result, early studies of thequench problem are mostly restricted to quantum phasetransitions in exactly solvable models such as the one-dimensional (1D) Ising model in a transverse field ,the infinite range ferromagnetic Ising model , the 1DXY model , and quantum spin chains . On theexperimental side, trapped ultracold atoms in optical lat-tices provide possibilities of realization of many of theabove-mentioned systems . Experimental studies of de-fect production due to quenching of the magnetic field ina spin-one Bose condensate has also been undertaken . Recently, Kitaev proposed a 2D spin-1/2 model on ahoneycomb lattice with a Hamiltonian H = X j + l =even ( J σ xj,l σ xj +1 ,l + J σ yj − ,l σ yj,l + J σ zj,l σ zj,l +1 ) , (1)where j and l denote the column and row indices ofthe honeycomb lattice. This model has several interest-ing features which led to a plethora of theoretical workson it . For example, it provides a rare examplewhere a 2D model can be exactly solved using a Jordan-Wigner transformation . Further, when J = 0,the model provides an example of an 1D spin model whichsupports a topological quantum phase transition with thecritical point at J = J . Moreover, in d = 2, the modelsupports a gapless phase for | J − J | ≤ J ≤ J + J which has a possible connection to a spin liquid stateand demonstrates fermion fractionalization at all energyscales . Finally, it has been shown in Ref. 18 that thepresence of magnetic field, which induces a gap in the2D gapless phase, leads to non-Abelian statistics of thelow-lying excitations of the model; these excitations canbe viewed as robust qubits in a quantum computer . Anextended version of this model has also been suggestedin Ref. 25 which has the Hamiltonian H = J " X j + l =odd σ yj,l σ zj +1 ,l σ xj +2 ,l + X j + l =even σ xj,l σ zj +1 ,l σ yj +2 ,l + H . (2)The quench dynamics of the 2D Kitaev model has beenstudied very recently in Ref. 26. It has been shown thatfor this model, quenching J takes the system through acritical line instead of critical point which leads to uncon-ventional variation of the defect density as a function ofthe quench rate. In this context, it has also been shownthat for a general d -dimensional model, where the quenchtake the system through a d − m dimensional hypersurfacecharacterized by the correlation length exponent ν anddynamical critical exponent z , the defect density obeys n d ∼ /τ mν/ ( zν +1) . The Kitaev model provides a con-crete example of such a quench for d = 2 and m = 1.The defect correlation function for such a quench hasalso been computed in Ref. 26.In this work, we extend and elaborate on the resultsof Ref. 26 and study the quench dynamics of the Kitaevmodel both in d = 1 and d = 2 and the extended Kitaevmodel in d = 2. The main results that we have obtainedare the following. First, we show that in 1D ( J = 0),where quenching J takes the system across the topologi-cal quantum critical point located at J = J , the densityof defects produced due to the quench scales as 1 / √ τ inthe limit of slow quench (large τ ). We also identify andcompute all independent non-zero spin-spin correlationfunctions and use them to elucidate the spatial extentof the defect correlation function. Second, we outline ageneral proof of the result reported in Ref. 26 that fora d dimensional quantum model, where the quench takethe system through a d − m dimensional hypersurfacecharacterized by the correlation length exponent ν anddynamical critical exponent z , the defect density obeys n d ∼ /τ mν/ ( zν +1) . Third, we elaborate on the variationof shape and size of the defect correlation function for the2D Kitaev model with the quench rate and the model pa-rameters. Fourth, we compute the entropy generated dueto such a quench and discuss its dependence on the modelparameters and the quench rate. Finally, we study thedefect scaling law, entropy generation and defect correla-tion function of the 2D extended Kitaev model describedby H .The organization of the paper is as follows. In Sec. II A,we analyze the quench dynamics of the Kitaev model in1D and obtain the quench rate dependence of the defectdensity. This is followed by Sec. II B, where we computethe 1D correlation functions and use them to discuss thenature of the defect correlation function. Next, in Sec.III A, we obtain the quench rate dependence of the defectdensity in 2D. The computation of the defect correlationfunction is detailed in Sec. III B and the entropy gener-ated during the quench process is computed in Sec. III C.This is followed by the study of quench dynamics of theextended Kitaev model in Sec. IV. Finally, we concludein Sec. V. FIG. 1: Schematic representation of the Kitaev model ona honeycomb lattice showings the bonds J , J and J .Schematic pictures of the ground states, which correspondto pairs of spins on vertical bonds locked parallel (antiparal-lel) to each other in the limit of large negative (positive) J ,are shown at one bond on the left (right) edge respectively. ~M and ~M are spanning vectors of the lattice, and a and b represent inequivalent sites. II. QUENCH IN 1DA. Defect density
For J = 0, the Kitaev model represents a spin-1/2model in 1D with the Hamiltonian H = X n (cid:0) J σ x n σ x n +1 + J σ y n − σ y n (cid:1) , (3)where n denotes site indices of a one dimensional chainwith N sites (we will assume N is a multiple of 4). Thelattice spacing a and the Planck constant ~ will be setequal to 1 in the rest of this work. The Hamiltonianin Eq. (3) can be exactly diagonalized using a standardJordan-Wigner transformation a n = n − Y j = −∞ σ zj σ y n ,b n = n Y j = −∞ σ zj σ x n +1 , (4)where a n and b n are independent Majorana fermions atsite n . They satisfy relations such as a † n = a n , { a m , a n } =2 δ m,n and { a m , b n } = 0. The label n for a n and b n goover N/ H can be written as H = i X n [ J b n a n + J b n a n +1 ]= 2 i π X k =0 [ b † k a k ( J + J e ik )+ a † k b k ( − J − J e − ik ) ] , (5)where the Majorana fermion creation and destruction op-erators a † k and a k are Fourier components of the a n ’s, a n = r N π X k =0 [ a k e ikn + a † k e − ikn ] . (6)The sum over k in Eq. (6) only goes over half the Bril-louin zone because a n describes a Majorana fermion; thenumber of modes lying in the range 0 ≤ k ≤ π is N/ k = 0 and π for which there is nodistinction between k and − k ; these two modes shouldhave a coefficient of p /N instead of p /N . However,we will ignore this correction here because we will be in-terested in the N → ∞ limit, and we will change from asum over k to an integral over k .] The operators a k and a † k satisfy the anticommutation relations { a k , a † k ′ } = δ kk ′ and { a k , a k ′ } = 0. One can now define a two componentfermionic creation operator ψ k = ( a k b k ), so that H can be written as H = π X k =0 ψ † k H k ψ k , where H k = 2 i (cid:18) − J − J e − ik J + J e ik (cid:19) . (7)From Eq. (7), we find that H can be diagonalized lead-ing to an energy spectrum consisting of two bands E ± k = ± q J + J + 2 J J cos k. (8)Note that the band gap vanishes at J = ± J for k = π and 0 respectively, where the bands touch each other. Itwas shown in Ref. 19 that this vanishing of the energygap signals a topological phase transition between thetwo phases of the model at J > J and J < J .To study the quench of the system across this criticalpoint, we will now consider what happens when we evolve J linearly in time at a rate 1 /τ from −∞ to ∞ , keeping J fixed at some positive value: we take J = J t/τ . Theground states of H in Eq. (7) have σ x n σ x n +1 = 1 and − t = −∞ and ∞ respectively for all values of n .In terms of the Hamiltonian in Eq. (7), the ground andexcited states for J → −∞ are respectively given by ψ k = 1 √ (cid:18) i (cid:19) and ψ k = 1 √ (cid:18) − i (cid:19) . (9)For J → ∞ , the ground and excited states are given by ψ k and ψ k respectively. FIG. 2: Defect density produced by quenching J in d = 1. By a change of basis, one can rewrite Eq. (7) in theform H = P k ψ ′ † k H ′ k ψ ′ k where H ′ k = 2 (cid:18) J + J cos k − J sin k − J sin k − J − J cos k (cid:19) . (10)Note that unlike Eq. (7), the off-diagonal elements of Eq.(10) do not change with time if J is held fixed. Asa result, the problem of quench dynamics is reduced tosolving a standard Landau-Zener problem for each mo-mentum k . The density of defect formation n can thusbe found to be n = Z π dkπ p k , where p k = e − πJ τ sin k (11)denotes the probability of the system to remain in theinitial ( J → −∞ ) ground state for momentum k . For J τ ≫
1, the contribution to n comes mainly from theregions near k = 0 and π where p k = 1. Thus one findsthat in the slow quench regime n ≃ / √ J τ . Such a1 / √ τ scaling of defect density conforms to the predictionof Ref. 4. For the present case, it is easy to see from Eq.(8), that the gap ∆( k ) = E + ( k ) − E − ( k ) vanishes linearlyat the critical point both with the quench parameter J and with momentum around k = 0 and π , so that zν = z = 1. Thus, n ∼ /τ dν/ ( zν +1) = 1 / √ τ in 1D.A plot of the defect density as a function of the quenchtime τ is shown in Fig. 2. The plot confirms the expectedresult, that the defect density is maximum for an infinitequench rate ( τ → n decreases quickly before settlingdown to a 1 / √ τ behavior for large τ .It is useful to note that the Hamiltonian H k in Eq. (7)can also be written, after a suitable change of basis, as H ′ k = 2 (cid:18) J − sin( k/ − iJ + cos( k/ iJ + cos( k/ − J − sin( k/ (cid:19) , (12)where J ± = J ± J . This form is useful if, for instance,one wants to study the effect of quenching J − from −∞ to ∞ keeping J + fixed. B. Correlation functions
Let us now consider how the system may be de-scribed at the final time t → ∞ when J = ∞ .In principle, the time evolution of the system is uni-tary, so that it will always be a pure state. However,for each momentum k , the wave function is given by √ − p k ψ k e − iE k t + √ p k ψ k e − iE k t , where E , k = ±∞ .As a result, the final density matrix of the system willhave off-diagonal terms involving ψ ∗ k ψ k and ψ ∗ k ψ k which vary extremely rapidly with time; their effects onphysical quantities will therefore average to zero. Hencethe final density matrix is effectively diagonal like thatof a mixed state , where the diagonal entries are time-independent as t → ∞ and are given by 1 − p k and p k .Such a density matrix is associated with an entropy whichwe will discuss in Sec. III C in the context of 2D Kitaevmodel.Using the above density matrix, we will now computethe correlation functions corresponding to the operators O r = ib n a n + r , where r is an integer. In terms of thespins, as can be seen from Eq. (4), the operator O r canbe written as O = σ x n σ x n +1 , O = σ y n +1 σ y n +2 ,O r = σ y n +1 n +2 r − Y j =2 n +2 σ zj σ y n +2 r for r ≥ , = σ y n +2 r n Y j =2 n +2 r +1 σ zj σ y n +1 for r ≤ − . (13)We will calculate the expectation values of these opera-tors shortly. In principle, one can also consider expecta-tion values of the operators ia n a n + r and ib n b n + r ; how-ever a direct calculation shows that these vanish if r = 0.Further, for the Kitaev model, it has been shown thatthe spin-spin correlations between sites lying on differentbonds vanish, i.e. , h σ x n σ x n + r i = 0 for | r | ≥ . There-fore h O r i are the only non-vanishing spin-correlators ofthe model .To compute h O r i we note that O r can be expressed interms of the fermion operators a k and b k . This will ingeneral involve summations over two different momenta k and k ′ . However, when h O r i is computed in a directproduct of states involving a k and b k , only terms in which k ′ = k will contribute. In the limit N → ∞ , the relevantpart of O r which contributes to the correlation functioncan be written as O r = − iN π X k =0 [ b † k a k e ikr − a † k b k e ikr ] . (14) < O r > r FIG. 3: Plot of correlation function h O r i as a function of r for J τ = 10 (red circles and red solid line) and J τ = 1(black squares and black dashed line). h O r i shows a dampedoscillatory behavior as a function of r . Using the wave functions given in Eq. (9), we find that h O r i = ± Z π dk cos( kr ) = ± δ r, , (15)where the + and − signs refer to the ground states of J = −∞ and ∞ respectively. This is expected since σ x n σ x n +1 = ± J = ∞ with probabilities 1 − p k and p k respectively,we find that h O r i = − δ r, + 2 π Z π dk p k cos( kr ) . (16)A plot of h O r i as a function of r is shown for represen-tative values of J τ = 1 ,
10 in Fig. 3. We find that h O r i shows a damped oscillatory behavior. Note that since h O r i = − δ r, for the ground state of H for J → ∞ , theplot of h O r i in the state of the system after the quenchprovides a direct measurement of the spatial extent ofthe correlation between the defects generated during thequench.For J τ ≫
1, the dominant contribution in the integralin Eq. (16) comes from the regions near k = 0 and π ascan be seen from the expression for p k in Eq. (11). Onecan combine these two regions, and write the expressionin (16) approximately as h O r i = − δ r, + 2 π Z ∞ dk e − πJ τk × [cos( kr ) + cos { ( π − k ) r } ]= − δ r, + 1 + ( − r π e − r / (8 πJ τ ) √ J τ . (17)Note that this vanishes if r is odd. For a given value of J τ , the expression in Eq. (17) decreases with increasing r , particularly for r > √ πJ τ . On the other hand,for a given large value of r , Eq. (17) has a maximumat τ = r / (4 πJ ). The fact that the crossover in bothcases occurs around r ∼ √ πJ τ signals the fact that theassociated length scale for the defect correlation functionis of order √ πJ τ . C. Sum rule
There is a sum rule that we can write down for h O r i .From Eq. (16), we see that O total ≡ ∞ X r = −∞ h O r i = − p , (18)where we have used the identity P r e ikr = 2 πδ ( k ) for − π < k < π . Going back to Eq. (10), we see that for k = 0, the Hamiltonians at different times commute witheach other irrespective of how J is varied in time from −∞ to ∞ . This means that if we start with the groundstate of J = −∞ , no transition will occur at any time,and we will have p = 1. Eq. (18) then implies that O total = 1. III. 2D KITAEV MODELA. Defect density
When J = 0, the Kitaev model with Hamiltoniangiven by Eq. (1) describes a spin model on a hexag-onal 2D lattice. Usually spin models are not exactlysolvable in two dimensions. One of the main proper-ties of the Kitaev model which makes it theoreticallyattractive is that, even in 2D, it can be mapped ontoa non-interacting fermionic model by a suitable Jordan-Wigner transformation . In terms of the Majo-rana fermions a jl and b jl one can write a jl = j − Y i = −∞ σ zil ! σ yjl for even j + l,b jl = j − Y i = −∞ σ zil ! σ xjl for odd j + l. (19) Such a transformation maps the spin Hamiltonian H inEq. (1) to a fermionic Hamiltonian given by H = i X ~n [ J b ~n a ~n − ~M + J b ~n a ~n + ~M + J D ~n b ~n a ~n ] , (20)where ~n = √ i n + ( √ ˆ i + ˆ j ) n denote the midpointsof the vertical bonds. Here n , n run over all integersso that the vectors ~n form a triangular lattice whose ver-tices lie at the centers of the vertical bonds of the un-derlying honeycomb lattice; the Majorana fermions a ~n and b ~n sit at the top and bottom sites respectively ofthe bond labeled ~n . The vectors ~M = √ ˆ i + ˆ j and ~M = √ ˆ i − ˆ j are spanning vectors for the reciprocallattice, and D ~n can take the values ± ~n . The crucial point that makes the solution of Ki-taev model feasible is that D ~n commutes with H , sothat all the eigenstates of H can be labeled by specificvalues of D ~n . It has been shown that for any value ofthe parameters J i , the ground state of the model alwayscorresponds to D ~n = 1 on all the bonds. Since D ~n is aconstant of motion, the dynamics of the model startingfrom any ground state never takes the system outside themanifold of states with D ~n = 1.For D ~n = 1, it is straightforward to diagonalize H in momentum space. We define Fourier transforms of theMajorana operators a ~n as a ~n = r N X ~k [ a ~k e i~k · ~n + a † ~k e − i~k · ~n ] (21)(and similarly for b ~n ), where N is the number of sites(assumed to be even, so that the number of unit cells N/ ~k extends over halfthe Brillouin zone of the 2D hexagonal lattice. We thenobtain H = P ~k ψ † ~k H ~k ψ ~k , where ψ † ~k = ( a † ~k , b † ~k ), and H ~k can be expressed in terms of Pauli matrices σ , , as H ~k = 2 [ J sin( ~k · ~M ) − J sin( ~k · ~M )] σ + 2 [ J + J cos( ~k · ~M ) + J cos( ~k · ~M )] σ . (22)The energy spectrum of H therefore consists of twobands with energies E ± ~k = ± J sin( ~k · ~M ) − J sin( ~k · ~M )) + ( J + J cos( ~k · ~M ) + J cos( ~k · ~M )) ] / . (23)We note for | J − J | ≤ J ≤ ( J + J ), these bands toucheach other so that the energy gap ∆ ~k = E + ~k − E − ~k vanishesfor special values of ~k leading to the gapless phase of themodel .We will now quench J ( t ) = Jt/τ at a fixed rate 1 /τ ,from −∞ to ∞ , keeping J , J and J fixed at some non-zero values; we have introduced the quantity J to fix thescale of energy. We note that the ground states of H corresponding to J → −∞ ( ∞ ) are gapped and have σ zj,l σ zj,l +1 = 1( −
1) for all lattice sites ( j, l ). To study thestate of the system after the quench, we first note thatafter an unitary transformation U = exp( − iσ π/ H = P ~k ψ ′ † ~k H ′ ~k ψ ′ ~k , where H ′ ~k = U H ~k U † isgiven by H ′ ~k = 2 [ J sin( ~k · ~M ) − J sin( ~k · ~M )] σ + 2 [ J ( t ) + J cos( ~k · ~M ) + J cos( ~k · ~M )] σ . (24)Hence the off-diagonal elements of H ′ ~k remain time inde-pendent and the problem of quench dynamics reduces toa Landau-Zener problem for each ~k . The defect densitycan then be computed following a standard prescription n = 1 A Z ~k d ~k p ~k ,p ~k = e − πτ [ J sin( ~k · ~M ) − J sin( ~k · ~M )] /J , (25)where A = 4 π / (3 √
3) denotes the area of half the Bril-louin zone over which the integration is carried out. Sincethe integrand in Eq. (25) is an even function of ~k , onecan extend the region of integration over the full Brillouinzone. This region can be chosen to be a rhombus withvertices lying at ( k x , k y ) = ( ± π/ √ ,
0) and (0 , ± π/ v , v ,each with a range 0 ≤ v , v ≤
1, one finds that k x = 2 π v + v − √ , k y = 2 π v − v . (26)Such a substitution covers the rhombus uniformly andfacilitates the numerical integration necessary for com-puting n .A plot of n as a function of the quench time Jτ and α =tan − ( J /J ) (we have taken J = J cos( α )[sin( α )]) isshown in Fig. 4. We note that the density of defects pro-duced is maximum when J = J . This is due to the factthat the length of the gapless line through which the sys-tem passes during the quench is maximum at this point.This allows the system to remain in the non-adiabaticstate for the maximum time during the quench, leadingto the maximum density of defects. For J /J > J /J ,the system does not pass through a gapless phase duringthe quench, and the defect production is exponentiallysuppressed.For sufficiently slow quench 2 πJτ ≫ p ~k is exponen-tially small for all values of ~k except in the region nearthe line J sin( ~k · ~M ) − J sin( ~k · ~M ) = 0 , (27)and the contribution to the momentum integral in (25)comes from values of ~k close to this line of zeroes. Wenote that the line of zeroes where p ~k = 1 precisely corre-sponds to the zeroes of the energy gap ∆ ~k as J is varied FIG. 4: Plot of defect density n as a function of the quenchtime Jτ and α = tan − ( J /J ). The density of defects ismaximum at J = J . for a fixed J and J . Thus the system becomes non-adiabatic when it passes through the intermediate gap-less phase in the interval | J − J | ≤ J ( t ) ≤ ( J + J ).It is then easy to see, by expanding p ~k about this linethat in the limit of slow quench, the defect density scalesas n ∼ / √ τ . We note that the scaling of the defectdensity with the quench rate in a quench where the sys-tem passes through a critical line in momentum space isdifferent from the situation where the quench takes thesystem through a critical point . In the latter case, forthe Kitaev model which has z = ν = 1, Ref. 4 predictsa defect density n ∼ /τ for d = 2. Thus the defectdensity crucially depends on the dimensionality of thecritical surface through which the system passes duringthe quench. This observation leads to a simple but gen-eral conclusion which we present below.Consider a d -dimensional model with z = ν = 1 de-scribed by a Hamiltonian H d = X ~k ψ † ~k (cid:18) ǫ ( ~k, t ) ∆( ~k )∆ ∗ ( ~k ) − ǫ ( ~k, t ) (cid:19) ψ ~k , (28)where ǫ ( ~k, t ) = ǫ ( ~k ) t/τ . Now let us assume that thequench takes such a system through a critical surface of d − m dimensions. The defect density for a sufficientlyslow quench can be expressed as n = 1 A Z BZ d d k p ( ~k ) , where p ( ~k ) = e − πτf ( ~k ) , ≃ A Z BZ d d k exp[ − τ X αβ =1 ,m g αβ k α k β ] ∼ /τ m/ , (29)where p ~k is the defect probability for momentum ~k , f ( ~k ) = | ∆( ~k ) | / | ǫ ( ~k ) | vanishes on the d − m dimen-sional critical surface, α, β denote one of the m or-thogonal directions to the critical surface and g αβ =( ∂ f ( ~k ) /∂k α ∂k β ) ~k ∈ critical surface . We note that this re-sult depends only on the property that f ( ~k ) has to vanishon a d − m dimensional surface, and not on the precisenature of f ( ~k ). For m = d , where the quench takes thesystem through a critical point, our results coincide withthat of Ref. 4.Finally we generalize our arguments for models wherethe d − m dimensional hypersurface is characterized bycorrelation length exponent ν and dynamical critical ex-ponent z . Let us assume that the system is describedby a Hamiltonian H [ λ ( t )] with quasi-energy eigenval-ues E ( ~k, t ) and that the time evolution of the param-eter λ ( t ) = λ ( t/τ ) takes the system through the criticalpoint λ = λ c at t = t . First we note that for large τ ,a non-vanishing probability of defect formation requiresthe non-adiabaticity condition | ∆( ~k ) | ∼ | ∂E ( ~k, t ) /∂t | [4]. Also, since ∂E ( ~k, t ) /∂t = ( ∂E ( ~k, t ) /∂λ ) τ − and nearthe critical point E ∼ λ zν , we get∆ ∼ τ − λ zν − (30)Further, as shown in Ref. 4, near any point on the criticalsurface, quite generally, one has ∆ ∼ | k | z , λ ∼ k /ν and k ∼ /τ ν/ ( zν +1) . Using these relations we find from Eq.30 that on any point near the gapless surface∆ ∼ /τ zν/ ( zν +1) (31)Next, let us consider the available phase space for for-mation of defects. When the quench takes the systemthrough a d − m dimensional hypersurface in momen-tum space, the available phase space is Ω ∼ k m ∼ ∆ m/z .Since this available phase space is directly proportionalto the defect density , we find, using Eq. (31), n ∼ Ω ∼ ∆ m/z ∼ /τ mν/ ( zν +1) (32)This generalizes the scaling law for defect density to ar-bitrary critical systems. Note that for z = ν = 1, werecover our earlier result n ∼ /τ m/ (Eq. (29)). For m = d , which represents quench through a critical point,we also recover the result of Ref. 4 ( n ∼ /τ dν/ ( zν +1) ) asa special case. B. Defect Correlation
The calculation of the correlation function can be ac-complished along similar lines as in 1D. First, we definethe operators O ~r = ib ~n a ~n + ~r . (33)In terms of the spin operators, we have O ~ = σ zj,l σ zj,l +1 .For ~r = ~ O ~r can be written as a product of spin op-erators going from a b site at ~n = ( j, l ) to an a site at ~n + ~r = ( j ′ , l ′ ): the product will begin with a σ x or σ y and end with a σ x or σ y with a string of σ z ’s in be-tween, where the choice of the initial and final σ matri-ces depends on whether the values of j + l and j ′ + l ′ areeven or odd. Note that O ~r for ~r = ~ h O ~r i versus ~r in the defect ground state provides an estimateof the shape and spatial extent of the defect correlationsproduced during the quench. Note that ( O ~r ) = 1,so that all the moments of O ~r can be found trivially: h ( O ~r ) n i = h O ~r i if n is odd and = 1 if n is even. O ~r can be written in terms of the Majorana fermionoperators a ~k and b ~k ; this again involves a sum over twodifferent momenta ~k and ~k ′ . However, the expectationvalue of O D~r in a direct product of states involving ~k only gets a contribution from terms in which ~k ′ = ~k . Itturns out that the relevant part of O D~r contributes tothe expectation values can be written as, O D~r = 4 iN X ~k [ b † ~k a ~k e i~k · ~r − a † ~k b ~k e − i~k · ~r ] . (34)The ground state and excited states for J = −∞ aregiven by ψ ~k and ψ ~k respectively, while the two statesare interchanged for J = ∞ . Using Eq. (9), we find that h O ~r i = ± N X ~k cos( ~k · ~r ) , (35)where the + and − signs refer to the ground states of J = −∞ and ∞ respectively. This confirms our earlierexpectation that in the ground states of J → −∞ ( ∞ ), h O ~r i = ± δ ~r,~ . Finally, in the state after quench, inwhich we have a mixture of the ground and excited statesof J = ∞ with probabilities 1 − p ~k and p ~k respectively,we find that h O ~r i = − δ ~r,~ + 2 A Z d ~k p ~k cos( ~k · ~r ) , (36)where the integral over momentum runs over half theBrillouin zone with area A . Note that the full Brillouinzone as well as p ~k remains invariant under a reflectionthrough the point ~k = ( π/ √ , k x → π/ √ − k x , k y → − k y . However, cos( ~k · ~r ) changes by a factor of( − n , if the components of ~r are given by x = √ n + n /
2) and y = 3 n /
2. Hence, h O ~r i = 0 for odd valuesof n .For large values of τ , substituting the expression inEq. (25) in the above integral, we find that the dominantcontribution comes from the region near the line given inEq. (27). Thus at every point ~k lying on that line, wecan introduce variables k k and k ⊥ which vary along theline and perpendicular to it along the directions ˆ n k andˆ n ⊥ respectively. Close to ~k , the integrand in Eq. (36)will take the form exp[ − aτ k ⊥ ± i ( ~k + k k ˆ n k + k ⊥ ˆ n ⊥ ) · ~r ], FIG. 5: Plot of O ~r sans the delta function peak at the originfor J = J = J and Jτ = 10 as a function of n and n (see text for details). The spatial anisotropy of the defectcorrelation function is clearly evident even for J = J . where a is a number of order 1 whose value dependson ~k . The integral over k ⊥ will give a factor ofexp (cid:2) − ( ~r · ˆ n ⊥ ) / (4 aτ ) (cid:3) / √ aτ . Thus we find that the den-sity of defects is of order 1 / √ τ in accordance with Eq.(29). This also leads us to expect that the spatial rangeof the defect correlation should go as √ τ .Next we consider the shape of the defect correlationfunction. For this purpose, we evaluate Eq. (36) numer-ically so as to obtain the ~r dependence of h O ~r i . Ingeneral, we expect the correlation will be anisotropic inspace if J /J ≫ ≪ J = 0 or J = 0 leads to the1D result derived in Sec. II B. A plot of the correlationfunction h O ~r i , without the delta function peak at ~r = 0,and as a function of n and n , where x = √ n + n / y = 3 n / h O ~r =0 i in or-der to make the correlations at ~r = ~ x direction, the correlations oscillate; the amplitude of os-cillations decays monotonically with x , in a qualitativelysimilar manner to the 1D correlation function O r shownin Fig. 2 for y = n = 0. The correlations decay in amonotonic way with y for x = n + n / θ = tan − ( − .
5) in the figure).Thus the correlations behave quite anisotropically evenfor J = J .We now aim at obtaining an understanding of the vari-ation of the spatial dependence of (cid:10) O r (cid:11) with the pa-rameters J and J . Such a variation can be analyticallyunderstood by noting that for Jτ ≫
1, the maximumcontribution to h O ~r i comes from around the wave vector ~k for which p ( ~k ) = 1. For J ≫ ( ≪ )1, this occurs whensin[ ~k · ~M ( ~M )] = 0 which yields ~k = π ( √ i ∓ ˆ j ) /
2. The
FIG. 6: Plot of ˙ O r ¸ sans the delta function peak at theorigin as a function of ~r for several representative values of J /J for J = J and Jτ = 5. The plot displays the change inthe shape of defect correlation function as a function of J /J (see text for details). maximum contribution to (cid:10) O ~r (cid:11) occurs where cos( ~k · ~n )is maximum, ie . , ~k · ~n = 0. Thus for J ≫ ( ≪ ) J , (cid:10) O ~r (cid:11) is expected to be maximal along the lines n + n =0( n = 0) in the n − n plane. This expectation is con-firmed as seen in Fig. 6 which shows (cid:10) O r (cid:11) for several rep-resentative values of J /J for a fixed J = J and Jτ = 5.We find that h O ~r i is maximal along n = 0( n + n = 0)line for J = 5(0 . J . This clearly shows that the de-fects produced in the quench will be highly anisotropicin this limits. For intermediate values of J and J , theanisotropy in (cid:10) O r (cid:11) can be similarly deduced by firstfinding ~k for which p ~k = 0 and then computing ~n forwhich ~k · ~n vanishes. The gradual evolution of the shapeof (cid:10) O r (cid:11) as we go from the limit J ≪ J to the limit J ≪ J can be seen from in Fig. 6.To obtain a more detailed picture of the spatialanisotropy of the defect correlations as a function of J /J we define a parameter α : J = J cos( α )[sin( α )].A variation of α therefore changes the ratio J /J from0 to ∞ while fixing J + J = J = 1. The plot of (cid:10) O ~r (cid:11) at points ( n , n ) = ( − ,
0) (on the x axis of the n − n plane), ( n , n ) = (2 , −
2) (along the − ◦ line inthe n − n plane) and ( n , n ) = (0 , −
2) (on the y axisof the n − n plane) as a function of α shown in Fig. 7clearly reveal the nature of the anisotropy of the correla-tion function. We find that as the ratio of J /J = tan( α )is varied from 0 to ∞ , the correlation on the representa-tive point (1 ,
0) along the x axis increases till it reachesthe point J = J ( α = π/
4) and then decays to 0 as α approaches π/
2. This signifies that the correlation alongthe x axis in the n − n plane becomes maximum for J = J . On the other hand, for the representative point(0 ,
2) on the y axis and 2 , − − ◦ , the correlation becomes maximum when J ≪ J FIG. 7: Plot of ˙ O r ¸ (sans the delta function peak) at rep-resentative points ( − ,
0) on the x axis (black solid line)(0 ,
2) on the y axis (blue dotted line) and (2 , −
2) along − ◦ in the n − n plane (red dashed line) as a functionof α = tan − ( J /J ) for fixed J = 1. ( α = 0) and J ≫ J ( α = π/
2) respectively, as expectedfrom Fig. 6. This lead us to conclude that the spatialanisotropy of the defect correlation function (cid:10) O ~r (cid:11) de-pends crucially on the ratio of J /J .Finally we note that we can obtain a measure of thespatial extent of the defect correlation function by calcu-lating h ~r i ≡ X ~r ~r h O ~r i . (37)To evaluate this, we first rewrite Eq. (36) as h O ~r i = − δ ~r,~ + 1 A Z d ~k p ~k e i~k · ~r , (38)where the integral now runs over the entire Brillouinzone. We then note that ~r e i~k · ~r = −∇ ~k e i~k · ~r , inte-grate by parts in Eq. (38) so as to make ∇ ~k act on p ~k , and use the identity P ~r e i~k · ~r = 2 Aδ ( ~k ), to obtain h ~r i = − ∇ ~k p ~k ) ~k = ~ = 24 πτ ( J + J + J J ) /J . Thisshows that the spatial extent of h O ~r i grows as √ τ forlarge τ . [Eq. (17) shows that we get the same behavior in1D.] Finally, we can get an idea of the spatial anisotropyof h O ~r i by computing h ~r i θ ≡ X ~r ( x cos θ + y sin θ ) h O ~r i , (39)where ~r = ( x, y ), and θ denotes a direction alongwhich the spatial extent is being calculated. Bywriting ( x cos θ + y sin θ ) e i~k · ~r = − (cos θ∂/∂k x + sin θ∂/∂k x ) e i~k · ~r , we can prove that h ~r i θ = 6 πτ [( J − J ) cos θ + √ J + J ) sin θ ] /J . We see that h ~r i θ hasa marked dependence on θ ; in fact, it vanishes in thedirection given by θ = tan − [( J − J ) / √ J + J )],and is maximum in the perpendicular direction. Thesestatements should be interpreted with some care; h ~r i θ may be small for some value of θ either due to a can-cellation between positive and negative correlations orbecause h O ~r i is small in that direction.We note that the sum rule discussed in Sec. II C isalso valid in 2D, and we get P ~r h O ~r i = − p ~ = 1regardless of how J is varied from −∞ to ∞ . C. Entropy
As discussed in Sec. II B, for each momentum ~k , thefinal density matrix is effectively diagonal, with entries1 − p ~k and p ~k . The density matrix of the entire systemtakes the product form ρ = N ρ ~k . The von Neumannentropy density corresponding to this state is given by s = − A Z d ~k [ (1 − p ~k ) ln(1 − p ~k ) + p ~k ln p ~k ] , (40)where the integral again goes half the Brillouin zone. Letus now consider the dependence of this quantity on thequenching time scale τ . If τ is very small, the systemstays in its initial state and p ~k will be close to 1 for allvalues of ~k ; for the same reason, h O ~ i will remain closeto 1. If τ is very large, the system makes a transition tothe final ground state for all momentum except near theline described in Eq. (27). Hence p ~k will be close to 0 forall ~k except near that line, and h O ~ i will be close to -1.In both these cases, the entropy density will be small.We therefore expect that there will be an intermediateregion of values of τ in which s will show a maximumand h O ~ i will show a crossover from − s and as a function of Jτ and α , shown in Fig. 8 confirmsthis expectation. We find that the entropy reaches amaximum for the intermediate value of Jτ where h O ~ i crosses over from − α . IV. EXTENDED KITAEV MODEL
The extended Kitaev model, described by H (Eq. (2)),can also be mapped, using the Majorana transformationgiven by Eq. (19), to a Fermionic Hamiltonian H ′ = iJ X ( j,l ) ∈ A [ a j,l a j +2 ,l − b j,l +1 b j +2 ,l +1 ] + H , (41)where H is given by Eq. (20). We note that in thismodel, just as for H , D n commutes with all the termsin the Hamiltonian and the ground state corresponds to D n = 1 for all links of the honeycomb lattice. Thus, in0 FIG. 8: Plot of the entropy density s as a function of Jτ and α = tan − ( J /J ). The entropy density peaks when h O ~ i crosses from − momentum space, H ′ reduces to a bilinear 2 by 2 matrixHamiltonian H ′ = P ~k ψ ( ~k ) † H ′ ( ~k ) ψ ( ~k ), where H ′ ( ~k ) = 2 ( [ J sin( ~k · ~M ) − J sin( ~k · ~M )] σ +[ J + J cos( ~k · ~M ) + J cos( ~k · ~M )] σ − J X k sin( √ k x ) σ ) . (42)This can be diagonalized to obtain the energy eigenvalues E ′ ± ~k = ± J sin ( √ k x ) + [ J + J cos( ~k · ~M )+ J cos( ~k · ~M )] + [ J sin( ~k · ~M ) − J sin( ~k · ~M )] ! / (43)Note that the presence of a non-zero J introduces a gapin the spectrum (except when √ k x = nπ ) for all valuesof J , J and J . Thus the quench of J ( J = J ( t/τ ))carries the system through a critical point at t = 0 pro-vided | J − J | ≤ J ≤ ( J + J ).The probability p ~k of defect formation in such aquench, where the system evolves according to Landau-Zenner dynamics, can be read off from Eqs. (42-43) as p ~k = e − πτ ( E ′ ± ~k ) | J / | J sin( √ k x ) | . (44)The density of defects is thus given by n = R d ~kp ~k /A ,where the integral is taken over half the Brillouin zone FIG. 9: Plot of the defect density as a function of η = J /J and Jτ for J = J = J . defined by the triangle with vertices lying at ( k x , k y ) =(2 π/ √ , , (0 , π/ , (0 , − π/
3) and A is the area of thisregion. A plot of the defect density as a function of thequench rate τ and η = J /J for J = J = J is shownin Fig. 9. Note that for large quench time τ , the maxi-mum contribution to the quench comes from around themomentum ~k = ( k x , k y ) for which E ′ ± ~k | J =0 vanishes.Around this point p ~k ∼ exp[ − πJτ P α,β = x,y f αβ ( ~k )( ~k − ~k ) α ( ~k − ~k ) β ] so that n ∼ /τ in accordance with theprediction of the general formula Eq. (32) for d = m = 2and ν = z = 1.Next, we look at the defect correlation functions forthe extended Kitaev model. To this end, we define theoperator O ext ~r = i ( a ~n a ~n + ~r − b ~n b ~n + ~r ) (45)and consider its expectation value for ~r = ~
0. Here ~r = ( √ n + √ n / , n /
2) (with integers n and n )specifies the sites of the honeycomb lattice. (For ~r = ~ O ext ~r vanishes since a ~n = b ~n = 1).For J → ∓∞ , the model reduces to a set of decoupledchains involving Majorana fermions on nearest neighborsites. For this model, it is known that for ~r = ~ (cid:10) O ext ~r (cid:11) = ∓ δ n , πn [( − n −
1] (46)in the ground states for J → ∓∞ respectively. Forgeneric values of J and for a mixed final state after the1quench characterized by a defect probability p ~k , we find (cid:10) O ext ~r (cid:11) = − N X ~k D a † ~k a ~k − b † ~k b ~k E sin( ~r · ~k )= δ n , πn [( − n − A Z d ~k sgn [sin( √ k x )] p ~k sin( ~r · ~k ) . (47)The sign of sin( √ k x ) appears in Eq. (47) becausefor J → ∞ , the ground state of Eq. (42) has D a † ~k a ~k − b † ~k b ~k E = ± √ k x ) > < J = J = J , J = ηJ and Jτ → ∞ . Note that oneneeds the condition 0 ≤ η ≤ h O ext ~r i in Eq. (47) comes from ~k = ~k = ((2 / √
3) cos − ( − η/ ,
0) where p ~k = ~k = 1. Thus for Jτ → ∞ one gets (cid:10) O ext ~r (cid:11) ≃ sin (cid:20) { n + n } cos − (cid:18) − η (cid:19)(cid:21) (48)where we have omitted the first term (proportional to δ n , ) in Eq. (47). Eq. (48) clearly brings out the de-pendence of the spatial anisotropy of the defect correla-tion function as a function of η . In particular, for η = 0, h O ext ~r i ∼ sin { ( n + n / π } , so that its sign alternates be-tween sites with odd and even values of n (if n is odd).Similarly, for η = 2, h O ext ~r i ∼ sin { (2 n + n ) π } ∼ h O ext ~r i for J = J = J , J = ηJ and Jτ = 3 as shown in Fig.10. We find that for η = 0 (top left plot of Fig. 10),it alternates between odd and even n sites, while for η close to 2 (bottom right plot in Fig. 10), the correlationfunction is much smaller than for η = 0.Finally, we compute the entropy generated due to sucha quench process given by Eq. (40) where p ~k is given byEq. (44). A plot of the entropy density as a function of Jτ and α = tan − ( J /J ) with J = J = J is shown inFig. 11. Once again we find, similar to that in the Kitaevmodel, that the entropy density peaks for intermediatevalue of τ . V. DISCUSSION
In conclusion, we have studied the quench dynamicsof the Kitaev model in 1D and 2D and the extended Ki-taev model in 2D. For the 1D Kitaev model and the 2Dextended Kitaev model, we have shown that the defect
FIG. 10: Plot of the defect correlation function (sans the firstterm with the delta function peak in Eq. (47)) with Jτ = 3and J = J = J for several representative values of η = J /J . See text for details.FIG. 11: Plot of the entropy density s as a function of quenchtime τ and η = J /J . density scales as 1 /τ d/ with the quench time τ , in accor-dance with the general results of Ref. 4. For the 2D Ki-atev model, where the quench takes the system througha gapless line, we found that the scaling of the defectdensity with τ changes due to the presence of a critical line instead of a critical point . In this context, we havepresented a general formula for the quench rate depen-dence of the defect density for a d dimensional systemwhen the quench takes such a system through a d − m dimensional critical surface. We have also computed the2defect correlation function for such quenches by an ex-act computation of all independent non-zero spin corre-lation functions in the defect ground state. In d = 2,we have found that the defect correlation function ex-hibit spatial anisotropy and studied the dependence ofthis anisotropy with the system parameter. Finally, wehave computed the entropy generated in such processesand have shown that the entropy peaks approximately atvalues of the quench rate for which the defect correlationfunction changes from − . If this can bedone, the evolution of the defect correlations with vari-ous parameters (such as J /J as shown in Fig. 7) can,in principle, be experimentally detected by spatial noisecorrelation measurements as pointed out in Ref. 32. Finally, we would like to note that the quench dynam-ics of the XXZ spin-1/2 chain has been recently studiedwith the Hamiltonian being varied along a line in pa-rameter space where the model is critical . In momen-tum space, the model only has a finite number of criticalpoints, but the system stays close to those critical pointsfor a long time. This is a different situation from the onethat we have analyzed in Sec. III where there is a line ofcritical points in momentum space; hence our results forthe scaling of the defect density are not applicable to thework in Ref. 33.We thank Amit Dutta and Anatoly Polkovnikov forstimulating discussions and several important sugges-tions. DS thanks DST, India for financial support underthe project SR/S2/CMP-27/2006. S. Sachdev,
Quantum Phase Transitions (Cambridge Uni-versity Press, Cambridge, 1999). T. W. B. Kibble, J. Phys. A , 1387 (1976); W. H. Zurek,Nature (London) , 505 (1985). B. Damski, Phys. Rev. Lett. , 035701 (2005). A. Polkovnikov, Phys. Rev. B , 161201(R) (2005). A. Polkovnikov and V. Gritsev, arXiv:0706.0212v2, to ap-pear in Nature Physics. The rate of defect production might change if the quenchtakes the system through some special multicritical point.See Ref. 12 for a discussion of this issue. K. Sengupta, S. Powell, and S. Sachdev Phys. Rev. A ,053616 (2004). J. Dziarmaga, Phys. Rev. Lett. , 245701 (2005), andPhys. Rev. B , 064416 (2006). P. Calabrese and J. Cardy, J. Stat. Mech: Theory ExptP04010 (2005), and Phys. Rev. Lett. , 136801 (2006). A. Das, K. Sengupta, D. Sen, and B. K. Chakrabarti, Phys.Rev. B , 144423 (2005). R. W. Cherng and L. Levitov, Phys. Rev. A , 043614(2006). V. Mukherjee, U. Divakaran, A. Dutta, and D. Sen, Phys.Rev. B , 174303 (2007). B. Damski and W. H. Zurek, Phys. Rev. A T. Caneva, R. Fazio, and G. E. Santoro, Phys. Rev. B ,144427 (2007). F. M. Cucchietti, B. Damski, J. Dziarmaga, and W. H.Zurek, Phys. Rev. A , 023603 (2007). For a review see I. Bloch, J. Dalibard, and W. Zwerger,arXiv:0704.3011v2, to appear in Rev. Mod. Phys. L. E. Sadler, J. M. Higbie, S. R. Leslie, M. Vengalat-tore, and D. M. Stamper-Kurn, Nature (London) , 312(2006). A. Kitaev, Ann. Phys. , 2 (2006). X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. , 087204 (2004). H.-D. Chen and Z. Nussinov, arXiv:cond-mat/0703633v5(unpublished). Z. Nussinov and G. Ortiz, arXiv:0709:2717v3 (unpub-lished). G. Baskaran, S. Mondal, and R. Shankar, Phys. Rev. Lett. , 247201 (2007). K. P. Schmidt, S. Dusuel, and J. Vidal, Phys. Rev. Lett. , 057208 (2008). A. Kitaev, Ann. Phys. , 2 (2003). D.-H. Lee, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. , 196805 (2007). K. Sengupta, D. Sen, and S. Mondal, Phys. Rev. Lett. ,077204 (2008). See for example, J. B. Kogut, Rev. Mod. Phys. , 659(1979). See for example, L. Landau and E. M. Lifshitz,
QuantumMechanics: Non-relativistic Theory , 2nd Ed. (PergamonPress, Oxford, 1965). See for example S. Suzuki and M. Okada in
Quantum An-nealing and Related Optimization Methods , Eds. by A. Dasand B. K. Chakrabarti (Springer-Verlag, Berlin, 2005). B. S. Shastry and D. Sen, Phys. Rev. B , 2988 (1997). L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. , 090402 (2003); A. Micheli, G. K. Brennen, and P.Zoller, Nature Physics , 341 (2006). E. Altman, E. Demler, and M. D. Lukin, Phys. Rev. A ,013603 (2004).33