We show that for a d-dimensional model in which a quench with a rate
\tau^{-1} takes the system across a d-m dimensional critical surface, the
defect density scales as n \sim 1/\tau^{m\nu/(z\nu +1)}, where \nu and z are
the correlation length and dynamical critical exponents characterizing the
critical surface. We explicitly demonstrate that the Kitaev model provides an
example of such a scaling with d=2 and m=\nu=z=1. We also provide the first
example of an exact calculation of some multispin correlation functions for a
two-dimensional model which can be used to determine the correlation between
the defects. We suggest possible experiments to test our theory.Comment: 4 pages including 4 figures; generalized the discussion of the defect
  density scaling to the case of arbitrary critical exponents, and added some
  references; this version will appear in Physical Review Letter

Diptiman Sen

K. Sengupta

L. Landau

S. Sachdev

S. Suzuki

Shreyoshi Mondal

English

arXiv

We show that for a d-dimensional model in which a quench with a rate $\tau^{-1}$ takes the system across a (d - m)-dimensional critical surface, the defect density scales as $n \sim 1/ \tau^{m\nu/(z\nu+1)}$, where  and z are the correlation length and dynamical critical exponents characterizing the critical surface. We explicitly demonstrate that the Kitaev model provides an example of such a scaling with d = 2 and m = \nu = z = 1. We also provide the first example of an exact calculation of some multispin correlation functions for a two-dimensional model that can be used to determine the correlation between the defects. We suggest possible experiments to test our theory

Sengupta, K

Sen, Diptiman

Mondal, Shreyoshi

Open Access Repository of IISc Research Publications

Exact Results for Quench Dynamics and Defect Production in a Two-Dimensional Model

Crossref

Name not available

We show that for a d-dimensional model in which a quench with a rate τ−1 takes the system across a (d-m)-dimensional critical surface, the defect density scales as n~1/τmν/(zν+1), where ν and z are the correlation length and dynamical critical exponents characterizing the critical surface. We explicitly demonstrate that the Kitaev model provides an example of such a scaling with d=2 and m=ν=z=1. We also provide the first example of an exact calculation of some multispin correlation functions for a two-dimensional model that can be used to determine the correlation between the defects. We suggest possible experiments to test our theory

Sengupta, K.

Publications of the IAS Fellows

ar
X
iv
:0
71
0.
17
12
v2
  [
co
nd
-m
at.
sta
t-m
ec
h]
  1
5 J
an
 20
08
Exact results for quench dynamics and defect production in a two-dimensional model
K. Sengupta1, Diptiman Sen2 and Shreyoshi Mondal1
1 TCMP division, Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700 064, India
2 Center for High Energy Physics, Indian Institute of Science, Bangalore, 560 012, India
(Dated: February 2, 2008)
We show that for a d-dimensional model in which a quench with a rate τ−1 takes the system
across a d−m dimensional critical surface, the defect density scales as n ∼ 1/τmν/(zν+1), where ν
and z are the correlation length and dynamical critical exponents characterizing the critical surface.
We explicitly demonstrate that the Kitaev model provides an example of such a scaling with d = 2
and m = ν = z = 1. We also provide the first example of an exact calculation of some multispin
correlation functions for a two-dimensional model which can be used to determine the correlation
between the defects. We suggest possible experiments to test our theory.
PACS numbers: 73.43.Nq, 05.70.Jk, 64.60.Ht, 75.10.Jm
Quantum phase transitions have been studied exten-
sively for several years [1]. Such a transition is accom-
panied by diverging length and time scales [1] leading to
the absence of adiabaticity close to the quantum critical
point. Thus the system fails to follow its instantaneous
ground state when some parameter in its Hamiltonian is
varied in time at a finite rate 1/τ which takes the sys-
tem across the critical point. Since the final state of the
system does not conform to the ground state of its final
Hamiltonian, defects are produced [2, 3]. The defect den-
sity n depends on the quench time τ as n ∼ 1/τdν/(νz+1),
where ν and z are the correlation length and dynamical
critical exponents at the critical point [4]. A theoretical
study of such a quench dynamics requires a knowledge of
the excited states of the system. Such studies have there-
fore been mostly restricted to phase transitions in exactly
solvable models in one or infinite dimensions [5, 6, 7, 8].
Experimental studies of defect production due to quench-
ing of the magnetic field in a two-dimensional (2D) spin-1
Bose condensate have been undertaken [9]. However, ex-
act studies of quench dynamics have not been carried out
so far for 2D spin models.
In this Letter, we carry out such a study for the 2D
Kitaev model [10]. We show that when the quench takes
the system across a critical (gapless) line, the density of
defects scales as 1/
√
τ . The Kitaev model has d = 2 and
ν = z = 1; hence, the above scaling is in contrast to
the n ∼ 1/τ behavior expected when the system passes
through a critical point [4]. In this context we provide
a general discussion of the scaling of the defect density
for a general d-dimensional model with arbitrary ν and
z; we show that when the quench takes such a system
through a d − m dimensional gapless (critical) surface,
the defect density scales as n ∼ 1/τmν/(zν+1). This re-
sult is a generalization of the result of Ref. [4] and hence
constitutes a significant extension of our current under-
standing of defect production due to a quench. We also
compute exactly some multispin correlation functions for
our model, and use them to study the variation of defect
correlations with the quench rate and the model param-
eters. Such an exact analysis of defect correlations has
not been carried out so far for 2D systems.
The Kitaev model is a spin-1/2 model on a 2D honey-
comb lattice with the Hamiltonian [10]
H =
∑
j+l=even
(J1σ
x
j,lσ
x
j+1,l + J2σ
y
j−1,lσ
y
j,l + J3σ
z
j,lσ
z
j,l+1),
(1)
where j and l denote the column and row indices of
the lattice. This model, sketched in Fig. 1, is known
to have several interesting features [12, 13, 14, 15]. It
is a rare example of a 2D model which can be exactly
solved [10, 12, 14]. It supports a gapless phase for
|J1 − J2| ≤ J3 ≤ J1 + J2 [10] which is possibly con-
nected to a spin liquid state and demonstrates fermion
fractionalization at all energy scales [13]. In certain pa-
rameter regimes, the ground state exhibits topological
order and the low-energy excitations carry Abelian and
non-Abelian fractional statistics; these excitations can
be viewed as robust qubits in a quantum computer [11].
There have been proposals for experimentally realizing
this model in systems of ultracold atoms and molecules
trapped in optical lattices [16]; such systems are known
to provide easy access to the study of non-equilibrium
dynamics of the underlying model. However, in spite of
several studies of the phases and low-lying excitations of
the Kitaev model, its non-equilibrium dynamics has not
been studied so far. We will study what happens in this
model when J3 is varied from −∞ to ∞ at a rate 1/τ ,
keeping J1 and J2 fixed.
One of the main properties of the Kitaev model which
makes it theoretically attractive is that, even in 2D, it
can be mapped onto a non-interacting fermionic model
by a suitable Jordan-Wigner transformation [12, 14],
HF = i
∑
~n
[J1b~na~n− ~M1 + J2b~na~n+ ~M2 + J3D~nb~na~n], (2)
where a~n and b~n are Majorana fermions sitting at the
top and bottom sites respectively of a bond labeled ~n,
~n =
√
3iˆ n1 + (
√
3
2 iˆ+
3
2 jˆ) n2 denote the midpoints of the
2FIG. 1: Schematic representation of the Kitaev model on
a honeycomb lattice showings the bonds J1, J2 and J3.
Schematic pictures of the ground states, which correspond
to pairs of spins on vertical bonds locked parallel (antiparal-
lel) to each other in the limit of large negative (positive) J3,
are shown at one bond on the left (right) edge respectively.
~M1 and ~M2 are spanning vectors of the lattice, and a and b
represent inequivalent sites.
vertical bonds shown in Fig. 1, and n1, n2 run over all
integers. The vectors ~n form a triangular lattice. The
x, y coordinates of the triangular lattice sites are given
by x =
√
3(n1 + n2/2) and y = 3n2/2. [We will refer
to the sites of the honeycomb lattice either as (j, l) as
in Eq. (1), or as (a, ~n) and (b, ~n) as in Eq. (2); j + l
is even (odd) for a (b) sites respectively.] The vectors
~M1 =
√
3
2 iˆ+
3
2 jˆ and
~M2 =
√
3
2 iˆ− 32 jˆ are spanning vectors
for the reciprocal lattice. The operator D~n can take the
values ±1 independently for each ~n and commutes with
HF , so that the states can be labeled by the values of D~n
on each bond; the ground state corresponds to D~n = 1
on all bonds [10, 12, 13, 14] irrespective of the sign of J3
due to a special symmetry of the model [10]. Since D~n is
a constant of motion, the dynamics of the model starting
from the ground state never takes the system outside the
manifold of states with D~n = 1.
For D~n = 1, Eq. (2) can be diagonalized as HF =∑
~k ψ
†
~k
H~kψ~k, where ψ
†
~k
= (a†~k, b
†
~k
) are Fourier trans-
forms of a~n and b~n, the sum over ~k extends over half
the Brillouin zone (BZ) of the triangular lattice formed
by the vectors ~n, and H~k can be expressed in terms of
the Pauli matrices σi (where σ3 is diagonal) as H~k =
2[J1 sin(~k · ~M1)−J2 sin(~k · ~M2)]σ1+2[J3+J1 cos(~k · ~M1)+
J2 cos(~k · ~M2)]σ2. The spectrum consists of two bands
with energies E±~k = ±E~k, where
E~k = 2[{J1 sin(~k · ~M1)− J2 sin(~k · ~M2)}2
+{J3 + J1 cos(~k · ~M1) + J2 cos(~k · ~M2)}2]1/2(3)
For |J1−J2| ≤ J3 ≤ J1+J2, the bands touch each other,
and the energy gap ∆~k = E
+
~k
− E−~k vanishes for special
FIG. 2: Plot of defect density n versus Jτ and α =
tan−1(J2/J1). The density of defects is maximum at J1 = J2.
values of ~k leading to a gapless phase [10, 12, 14, 15].
We will now quench J3(t) = Jt/τ from −∞ to ∞ at a
fixed rate 1/τ , keeping J , J1 and J2 fixed at some positive
values. The ground states of HF corresponding to J3 →
−∞(∞), schematically shown in Fig. 1, are gapped and
have σzj,lσ
z
j,l+1 = 1(−1) for all lattice sites (j, l) of type b.
To study the time evolution of the system, we note that
after an unitary transformation U = exp(−iσ1π/4), we
obtain HF =
∑
~k ψ
′†
~k
H ′~kψ
′
~k
, where H ′~k = UH~kU
† is given
by H ′~k = 2[J1 sin(
~k · ~M1) − J2 sin(~k · ~M2)]σ1 + 2[J3(t) +
J1 cos(~k · ~M1)+ J2 cos(~k · ~M2)]σ3. Hence the off-diagonal
elements of H ′~k remain time independent, and the quench
dynamics reduces to a Landau-Zener problem for each
~k. The defect density can then be computed following a
standard prescription [17]: n = (1/A)
∫
~k
d2~k p~k, where
p~k = exp[−2πτ{J1 sin(~k · ~M1)− J2 sin(~k · ~M2)}2/J ] (4)
is the probability of defect production for the state la-
beled by momentum ~k, and A = 4π2/(3
√
3) denotes the
area of half the BZ over which the integration is carried
out. A plot of n as a function of the quench time Jτ
and an angle α is shown in Fig. 2; here we have taken
J1[2] = J cosα[sinα]. We note that the density of defects
produced is maximum when α = π/4 (J1 = J2). This oc-
curs because the length of the gapless line through which
the system passes during the quench is maximum for
J1 = J2. Hence the system remains in the non-adiabatic
state for the maximum time during the quench, leading
to the maximum density of defects. Note that unlike
some other models [5], the defects here do not corre-
spond to topological defects since the dynamics always
keeps D~n = 1 on all bonds.
3For sufficiently slow quench 2πJτ ≫ 1, p~k is expo-
nentially small for all values of ~k except near the line
J1 sin(~k · ~M1) = J2 sin(~k · ~M2); the contribution to the
momentum integral in the expression for n comes from
this region. Note that the line p~k = 1 precisely corre-
sponds to the zeros of the energy gap ∆~k as J3 is varied
for fixed J1, J2. By expanding p~k about this line, we see
that for a very slow quench, the defect density scales as
n ∼ 1/√τ . This demonstrates that the scaling of n with
τ crucially depends on the dimensionality of the critical
surface since for a quench which takes the system through
a critical point instead of a critical line, the defect density
of the Kitaev model, which has d = 2 and ν = z = 1 [18],
is expected to scale as 1/τ [4]. This observation leads to
the following general conclusion.
Consider a d-dimensional model with ν = z =
1 which is described by a Hamiltonian Hd =∑
~k ψ
†
~k
(
σ3ǫ(~k)t/τ +∆(~k)σ+ +∆∗(~k)σ−
)
ψ~k, where σ
±
= (σ1 ± iσ2)/2. Suppose that a quench takes the system
through a critical surface of d −m dimensions. The de-
fect density for a sufficiently slow quench is given by [17]
n = (1/Ad)
∫
BZ
ddke−πτf(~k) ≃ (1/Ad)
∫
BZ
ddk exp[ −
πτ
∑m
α,β=1 gαβkαkβ ] ∼ 1/τm/2, where Ad is the area of
the d-dimensional BZ, f(~k) = |∆(~k)|2/|ǫ(~k)| vanishes
on the d − m dimensional critical surface, α, β denote
one of the m directions orthogonal to that surface, and
gαβ = [∂
2f(~k)/∂kα∂kβ ]~k∈critical surface. Note that this
result depends only on the property that f(~k) vanishes
on a d − m dimensional surface, and not on the pre-
cise form of f(~k). For general values of ν and z, we
note that a Landau-Zener type of scaling argument yields
∆ ∼ 1/τzν/(zν+1), where ∆ is the energy gap [4]. When
one crosses a d−m dimensional critical surface during the
quench, the available phase space Ω for defect production
scales as Ω ∼ km ∼ ∆m/z ∼ 1/τmν/(zν+1); this leads to
n ∼ 1/τmν/(zν+1). For a quench through a critical point
where m = d, we retrieve the results of Ref. [4].
Next we study the correlations between the defects
produced during the quench. To this end, we define
the operators O~r = ib~na~n+~r. In the spin language,
O~0 = σ
z
j,lσ
z
j,l+1. For ~r 6= ~0, O~r can be written as a prod-
uct of spin operators going from a b site at ~n to an a site
at ~n+~r: the product begins with a σx or σy at (j, l) and
ends with a σx or σy at (j′, l′) with a string of σz’s in
between, where the forms of the initial and final σ ma-
trices depend on the positions of j + l and j′ + l′. Note
that for J3 → −∞(∞), where the z component of each
spin is locked with that of its vertically nearest neighbor,
〈ψ−∞(∞)|O~r|ψ−∞(∞)〉 = ±δ~r,~0. For the Kitaev model,
it is known that the spin correlations between sites ly-
ing on different bonds vanish, e.g., 〈σza,~nσzb,~n+~r〉 = 0 for
~r 6= ~0 [13]. Therefore 〈O~r〉 are the only non-vanishing
two-point correlators of the model [14]. A non-zero value
of 〈O~r〉 for ~r 6= ~0 in the final state provides a signature of
FIG. 3: Plot of 〈O~r〉, sans the δ-function peak at the origin,
as a function of ~r for four values of J2/J , for Jτ = 5.
the defects. In particular, a plot of 〈O~r〉 versus ~r gives an
estimate of the correlations between the defects. [Since
O2~r = 1, all the moments of O~r can be found trivially:
〈On~r 〉 = 〈O~r〉 if n is odd and = 1 if n is even.]
After the quench, the system, for each momentum ~k,
is described by a combination of ψ−∞~k with probability
p~k and ψ∞~k with probability 1 − p~k, where ψ±∞~k are
the eigenstates of H~k for J3 → ±∞. Hence 〈O~r〉 can be
computed in a straightforward manner:
〈O~r〉 = − δ~r,~0 +
2
A
∫
d2~k p~k cos(
~k · ~r), (5)
where the integral runs over half the BZ with area A.
For large values of τ , the dominant contribution comes
from the region near the line J1 sin(~k · ~M1) = J2 sin(~k ·
~M2) where p~k = 1. Introducing the variables k‖ and k⊥
which vary along and perpendicular to this line (along
the directions nˆ‖ and nˆ⊥ respectively), we see that the
integrand in Eq. (5) takes the form exp[−a(~k0)τk2⊥ ±
i(~k0 + k‖nˆ‖ + k⊥nˆ⊥) · ~r], where a(~k0) is a number de-
pending on ~k0. The evaluation of the integral over k⊥
gives a factor of exp
[−(~r · nˆ⊥)2/(4aτ)] /√aτ . We thus
see that the magnitude of the defect correlations goes
as 1/
√
τ , while the spatial extent of the correlations
goes as
√
τ . This is confirmed by the following relation:∑
~r ~r
2〈O~r〉 = −2(∇2~kp~k)~k=~0 = 24πτ(J21 + J22 + J1J2)/J .
To study the correlations between defects, we evaluate
Eq. (5) numerically; the ~r dependence of 〈O~r〉 is shown
in Fig. 3 for several values of J2/J1, where J1 = J and
Jτ = 5. Here we have omitted the δ-function peak at
~r = 0 in Eq. 5 so as to make the correlation at ~r 6= ~0
clearly visible. The plot of 〈O~r〉 reflects the defect corre-
lations. To understand the variation of these correlations
with the ratio J2/J1, we note that for large Jτ , the max-
imum contribution to 〈O~r〉 comes from around the wave
4FIG. 4: Plot of 〈O~r〉 at the points (1, 0) on the n1 axis (black
solid line), (0, 2) on the n2 axis (blue dotted line), and (2,−2)
along the −45◦ line in the n1 − n2 plane (red dashed line) as
a function of α = tan−1(J2/J1), for J
2 = 1 and Jτ = 5.
vectors ~k0 for which p(~k0) = 1. For J2 ≫ (≪)J1, this im-
plies sin[~k · ~M2( ~M1)] = 0 which yields ~k0 ∼
√
3ˆi± jˆ. The
maximum contribution to 〈O~r〉 comes from cos(~k0·~r) = 1,
i.e., ~k0 · ~r = 0. Thus for J2 ≫ (≪)J1, 〈O~r〉 is expected
to be maximal along the lines y = −(+)√3x, namely,
n1 = −n2 (n1 = 0) in the n1−n2 plane. This expectation
is confirmed in Fig. 3 where 〈O~r〉 can be seen to be max-
imal along n1 = −n2 (n1 = 0) for J2 = 4(0.2)J1. Such a
strong anisotropy can be understood by noting that the
Kitaev model reduces to a one-dimensional model when
J2 ≫ (≪)J1. For intermediate values of J1/J2, a grad-
ual evolution of the defect correlations can be seen in
Fig. 3. We note that if the Kitaev model can be real-
ized using ultracold atoms in an optical lattice [16], such
an evolution of the defect correlations with J1/J2 can,
in principle, be experimentally detected by spatial noise
correlation measurements as pointed out in Ref. [19].
We can obtain a different view of the spatial anisotropy
of the defect correlations by studying 〈O~r〉 as a function
of α = tan−1(J2/J1). As α changes, the ratio J2/J1
varies from 0 to∞ while fixing J21 +J22 = J2 = 1. A plot
of 〈O~r〉 at three representative points (n1, n2) = (1, 0) (on
the n1 axis), (0, 2) (on the n2 axis), and (2,−2) (along the
−45◦ line in the n1−n2 plane) as a function of α, shown
in Fig. 4, reveals the nature of the defect correlations.
We see that as J2/J1 = tanα is varied from 0 to ∞, the
magnitude of the correlation at the point (1, 0) on the
n1 axis increases till it reaches a maximum at J1 = J2
(α = π/4), and then decays to 0 as α approaches π/2. For
the points (0, 2) on the n2 axis and (2,−2) along the line
with slope−45◦, the correlation becomes maximum when
J2 ≪ J1 (α = 0) and J2 ≫ J1 (α = π/2) respectively,
as expected from Fig. 3. We conclude that the spatial
anisotropy of the correlations between the defects 〈O~r〉
depends crucially on the ratio J2/J1.
To conclude, we have shown that the density of defects
produced during a quench of the Kitaev model through
a critical line scales with the quench time as 1/
√
τ , in-
stead of the 1/τ behavior expected for a quench through
a critical point. We have provided a general result for
the defect density which reproduces the result of Ref. [4]
as a special case. We have also discussed the variation
of the defect correlations with the model parameters and
pointed out the possibility of detection of these varia-
tions in experiments. These results significantly improve
our general understanding of the scaling of the density of
defects and their correlations in 2D systems.
We thank A. Dutta and A. Polkovnikov for stimulating
discussions.
[1] S. Sachdev, Quantum Phase Transitions (Cambridge
University Press, Cambridge, England, 1999).
[2] T. W. B. Kibble, J. Phys. A 9, 1387 (1976); W. H. Zurek,
Nature (London) 317, 505 (1985).
[3] B. Damski, Phys. Rev. Lett. 95, 035701 (2005).
[4] A. Polkovnikov, Phys. Rev. B 72, 161201(R) (2005); A.
Polkovnikov and V. Gritsev, arXiv:0706.0212.
[5] K. Sengupta, S. Powell, and S. Sachdev Phys. Rev. A
69, 053616 (2004); P. Calabrese and J. Cardy, J. Stat.
Mech: Theory Expt P04010 (2005), and Phys. Rev. Lett.
96, 136801 (2006); J. Dziarmaga, Phys. Rev. Lett. 95,
245701 (2005), and Phys. Rev. B 74, 064416 (2006).
[6] A. Das et al., Phys. Rev. B 74, 144423 (2006).
[7] R. W. Cherng and L. S. Levitov, Phys. Rev. A 73, 043614
(2006); V. Mukherjee et al., Phys. Rev. B 76, 174303
(2007).
[8] B. Damski and W. H. Zurek, Phys. Rev. A 73, 063405
(2006); F. M. Cucchietti et al., Phys. Rev. A 75, 023603
(2007); T. Caneva, R. Fazio, and G. E. Santoro, Phys.
Rev. B 76, 144427 (2007).
[9] L. E. Sadler et al., Nature (London) 443, 312 (2006).
[10] A. Kitaev, Ann. Phys. 321, 2 (2006).
[11] A. Kitaev, Ann. Phys. 303, 2 (2003).
[12] X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett.
98, 087204 (2007).
[13] G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev.
Lett. 98, 247201 (2007).
[14] H.-D. Chen and Z. Nussinov, arXiv:cond-mat/0703633.
[15] D.-H. Lee, G.-M. Zhang, and T. Xiang, arXiv:07053499;
K. P. Schmidt, S. Dusuel, and J. Vidal, arXiv:07093017.
[16] L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev.
Lett. 91, 090402 (2003); A. Micheli, G. K. Brennen, and
P. Zoller, Nature Physics 2, 341 (2006).
[17] See for example, L. Landau and E. M. Lifshitz, Quan-
tum Mechanics: Non-relativistic Theory, 2nd Ed. (Perg-
amon Press, Oxford, 1965); S. Suzuki and M. Okada in
Quantum Annealing and Related Optimization Methods,
Eds. by A. Das and B. K. Chakrabarti (Springer-Verlag,
Berlin, 2005).
[18] Note that as we approach the gapless phase, ∆E ∼ |J3 −
J3c| and ∆E ∼ |k|; hence ν = z = 1 for the model [4].
[19] E. Altman, E. Demler, and M. D. Lukin, Phys. Rev. A
570, 013603 (2004).


Exact results for quench dynamics and defect production in a two-dimensional model

Indian Academy of Sciences

ar
X
iv
:0
71
0.
17
12
v2
  [
co
nd
-m
at
.s
ta
t-
m
ec
h]
  1
5 
Ja
n 
20
08
Exact results for quench dynamics and defect production in a two-dimensional model
K. Sengupta1, Diptiman Sen2 and Shreyoshi Mondal1
1 TCMP division, Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700 064, India
2 Center for High Energy Physics, Indian Institute of Science, Bangalore, 560 012, India
(Dated: February 2, 2008)
We show that for a d-dimensional model in which a quench with a rate τ−1 takes the system
across a d − m dimensional critical surface, the defect density scales as n ∼ 1/τmν/(zν+1), where ν
and z are the correlation length and dynamical critical exponents characterizing the critical surface.
We explicitly demonstrate that the Kitaev model provides an example of such a scaling with d = 2
and m = ν = z = 1. We also provide the first example of an exact calculation of some multispin
correlation functions for a two-dimensional model which can be used to determine the correlation
between the defects. We suggest possible experiments to test our theory.
PACS numbers: 73.43.Nq, 05.70.Jk, 64.60.Ht, 75.10.Jm
Quantum phase transitions have been studied exten-
sively for several years [1]. Such a transition is accom-
panied by diverging length and time scales [1] leading to
the absence of adiabaticity close to the quantum critical
point. Thus the system fails to follow its instantaneous
ground state when some parameter in its Hamiltonian is
varied in time at a finite rate 1/τ which takes the sys-
tem across the critical point. Since the final state of the
system does not conform to the ground state of its final
Hamiltonian, defects are produced [2, 3]. The defect den-
sity n depends on the quench time τ as n ∼ 1/τdν/(νz+1),
where ν and z are the correlation length and dynamical
critical exponents at the critical point [4]. A theoretical
study of such a quench dynamics requires a knowledge of
the excited states of the system. Such studies have there-
fore been mostly restricted to phase transitions in exactly
solvable models in one or infinite dimensions [5, 6, 7, 8].
Experimental studies of defect production due to quench-
ing of the magnetic field in a two-dimensional (2D) spin-1
Bose condensate have been undertaken [9]. However, ex-
act studies of quench dynamics have not been carried out
so far for 2D spin models.
In this Letter, we carry out such a study for the 2D
Kitaev model [10]. We show that when the quench takes
the system across a critical (gapless) line, the density of
defects scales as 1/
√
τ . The Kitaev model has d = 2 and
ν = z = 1; hence, the above scaling is in contrast to
the n ∼ 1/τ behavior expected when the system passes
through a critical point [4]. In this context we provide
a general discussion of the scaling of the defect density
for a general d-dimensional model with arbitrary ν and
z; we show that when the quench takes such a system
through a d − m dimensional gapless (critical) surface,
the defect density scales as n ∼ 1/τmν/(zν+1). This re-
sult is a generalization of the result of Ref. [4] and hence
constitutes a significant extension of our current under-
standing of defect production due to a quench. We also
compute exactly some multispin correlation functions for
our model, and use them to study the variation of defect
correlations with the quench rate and the model param-
eters. Such an exact analysis of defect correlations has
not been carried out so far for 2D systems.
The Kitaev model is a spin-1/2 model on a 2D honey-
comb lattice with the Hamiltonian [10]
H =
∑
j+l=even
(J1σ
x
j,lσ
x
j+1,l + J2σ
y
j−1,lσ
y
j,l + J3σ
z
j,lσ
z
j,l+1),
(1)
where j and l denote the column and row indices of
the lattice. This model, sketched in Fig. 1, is known
to have several interesting features [12, 13, 14, 15]. It
is a rare example of a 2D model which can be exactly
solved [10, 12, 14]. It supports a gapless phase for
|J1 − J2| ≤ J3 ≤ J1 + J2 [10] which is possibly con-
nected to a spin liquid state and demonstrates fermion
fractionalization at all energy scales [13]. In certain pa-
rameter regimes, the ground state exhibits topological
order and the low-energy excitations carry Abelian and
non-Abelian fractional statistics; these excitations can
be viewed as robust qubits in a quantum computer [11].
There have been proposals for experimentally realizing
this model in systems of ultracold atoms and molecules
trapped in optical lattices [16]; such systems are known
to provide easy access to the study of non-equilibrium
dynamics of the underlying model. However, in spite of
several studies of the phases and low-lying excitations of
the Kitaev model, its non-equilibrium dynamics has not
been studied so far. We will study what happens in this
model when J3 is varied from −∞ to ∞ at a rate 1/τ ,
keeping J1 and J2 fixed.
One of the main properties of the Kitaev model which
makes it theoretically attractive is that, even in 2D, it
can be mapped onto a non-interacting fermionic model
by a suitable Jordan-Wigner transformation [12, 14],
HF = i
∑
~n
[J1b~na~n− ~M1 + J2b~na~n+ ~M2 + J3D~nb~na~n], (2)
where a~n and b~n are Majorana fermions sitting at the
top and bottom sites respectively of a bond labeled ~n,
~n =
√
3î n1 + (
√
3
2 î+
3
2 ĵ) n2 denote the midpoints of the
2
FIG. 1: Schematic representation of the Kitaev model on
a honeycomb lattice showings the bonds J1, J2 and J3.
Schematic pictures of the ground states, which correspond
to pairs of spins on vertical bonds locked parallel (antiparal-
lel) to each other in the limit of large negative (positive) J3,
are shown at one bond on the left (right) edge respectively.
~M1 and ~M2 are spanning vectors of the lattice, and a and b
represent inequivalent sites.
vertical bonds shown in Fig. 1, and n1, n2 run over all
integers. The vectors ~n form a triangular lattice. The
x, y coordinates of the triangular lattice sites are given
by x =
√
3(n1 + n2/2) and y = 3n2/2. [We will refer
to the sites of the honeycomb lattice either as (j, l) as
in Eq. (1), or as (a, ~n) and (b, ~n) as in Eq. (2); j + l
is even (odd) for a (b) sites respectively.] The vectors
~M1 =
√
3
2 î+
3
2 ĵ and
~M2 =
√
3
2 î− 32 ĵ are spanning vectors
for the reciprocal lattice. The operator D~n can take the
values ±1 independently for each ~n and commutes with
HF , so that the states can be labeled by the values of D~n
on each bond; the ground state corresponds to D~n = 1
on all bonds [10, 12, 13, 14] irrespective of the sign of J3
due to a special symmetry of the model [10]. Since D~n is
a constant of motion, the dynamics of the model starting
from the ground state never takes the system outside the
manifold of states with D~n = 1.
For D~n = 1, Eq. (2) can be diagonalized as HF =
∑
~k ψ
†
~k
H~kψ~k, where ψ
†
~k
= (a†~k, b
†
~k
) are Fourier trans-
forms of a~n and b~n, the sum over ~k extends over half
the Brillouin zone (BZ) of the triangular lattice formed
by the vectors ~n, and H~k can be expressed in terms of
the Pauli matrices σi (where σ3 is diagonal) as H~k =
2[J1 sin(~k · ~M1)−J2 sin(~k · ~M2)]σ1 +2[J3+J1 cos(~k · ~M1)+
J2 cos(~k · ~M2)]σ2. The spectrum consists of two bands
with energies E±~k = ±E~k, where
E~k = 2[{J1 sin(~k · ~M1) − J2 sin(~k · ~M2)}
2
+{J3 + J1 cos(~k · ~M1) + J2 cos(~k · ~M2)}2]1/2(3)
For |J1−J2| ≤ J3 ≤ J1 +J2, the bands touch each other,
and the energy gap ∆~k = E
+
~k
− E−~k vanishes for special
FIG. 2: Plot of defect density n versus Jτ and α =
tan−1(J2/J1). The density of defects is maximum at J1 = J2.
values of ~k leading to a gapless phase [10, 12, 14, 15].
We will now quench J3(t) = Jt/τ from −∞ to ∞ at a
fixed rate 1/τ , keeping J , J1 and J2 fixed at some positive
values. The ground states of HF corresponding to J3 →
−∞(∞), schematically shown in Fig. 1, are gapped and
have σzj,lσ
z
j,l+1 = 1(−1) for all lattice sites (j, l) of type b.
To study the time evolution of the system, we note that
after an unitary transformation U = exp(−iσ1π/4), we
obtain HF =
∑
~k ψ
′†
~k
H ′~kψ
′
~k
, where H ′~k = UH~kU
† is given
by H ′~k = 2[J1 sin(
~k · ~M1) − J2 sin(~k · ~M2)]σ1 + 2[J3(t) +
J1 cos(~k · ~M1)+ J2 cos(~k · ~M2)]σ3. Hence the off-diagonal
elements of H ′~k remain time independent, and the quench
dynamics reduces to a Landau-Zener problem for each
~k. The defect density can then be computed following a
standard prescription [17]: n = (1/A)
∫
~k
d2~k p~k, where
p~k = exp[−2πτ{J1 sin(~k · ~M1)− J2 sin(~k · ~M2)}
2/J ] (4)
is the probability of defect production for the state la-
beled by momentum ~k, and A = 4π2/(3
√
3) denotes the
area of half the BZ over which the integration is carried
out. A plot of n as a function of the quench time Jτ
and an angle α is shown in Fig. 2; here we have taken
J1[2] = J cosα[sinα]. We note that the density of defects
produced is maximum when α = π/4 (J1 = J2). This oc-
curs because the length of the gapless line through which
the system passes during the quench is maximum for
J1 = J2. Hence the system remains in the non-adiabatic
state for the maximum time during the quench, leading
to the maximum density of defects. Note that unlike
some other models [5], the defects here do not corre-
spond to topological defects since the dynamics always
keeps D~n = 1 on all bonds.
3
For sufficiently slow quench 2πJτ ≫ 1, p~k is expo-
nentially small for all values of ~k except near the line
J1 sin(~k · ~M1) = J2 sin(~k · ~M2); the contribution to the
momentum integral in the expression for n comes from
this region. Note that the line p~k = 1 precisely corre-
sponds to the zeros of the energy gap ∆~k as J3 is varied
for fixed J1, J2. By expanding p~k about this line, we see
that for a very slow quench, the defect density scales as
n ∼ 1/√τ . This demonstrates that the scaling of n with
τ crucially depends on the dimensionality of the critical
surface since for a quench which takes the system through
a critical point instead of a critical line, the defect density
of the Kitaev model, which has d = 2 and ν = z = 1 [18],
is expected to scale as 1/τ [4]. This observation leads to
the following general conclusion.
Consider a d-dimensional model with ν = z =
1 which is described by a Hamiltonian Hd =
∑
~k ψ
†
~k
(
σ3ǫ(~k)t/τ + ∆(~k)σ+ + ∆∗(~k)σ−
)
ψ~k, where σ
±
= (σ1 ± iσ2)/2. Suppose that a quench takes the system
through a critical surface of d −m dimensions. The de-
fect density for a sufficiently slow quench is given by [17]
n = (1/Ad)
∫
BZ
ddke−πτf(
~k) ≃ (1/Ad)
∫
BZ
ddk exp[ −
πτ
∑m
α,β=1 gαβkαkβ ] ∼ 1/τm/2, where Ad is the area of
the d-dimensional BZ, f(~k) = |∆(~k)|2/|ǫ(~k)| vanishes
on the d − m dimensional critical surface, α, β denote
one of the m directions orthogonal to that surface, and
gαβ = [∂
2f(~k)/∂kα∂kβ ]~k∈critical surface. Note that this
result depends only on the property that f(~k) vanishes
on a d − m dimensional surface, and not on the pre-
cise form of f(~k). For general values of ν and z, we
note that a Landau-Zener type of scaling argument yields
∆ ∼ 1/τzν/(zν+1), where ∆ is the energy gap [4]. When
one crosses a d−m dimensional critical surface during the
quench, the available phase space Ω for defect production
scales as Ω ∼ km ∼ ∆m/z ∼ 1/τmν/(zν+1); this leads to
n ∼ 1/τmν/(zν+1). For a quench through a critical point
where m = d, we retrieve the results of Ref. [4].
Next we study the correlations between the defects
produced during the quench. To this end, we define
the operators O~r = ib~na~n+~r. In the spin language,
O~0 = σ
z
j,lσ
z
j,l+1. For ~r 6= ~0, O~r can be written as a prod-
uct of spin operators going from a b site at ~n to an a site
at ~n+~r: the product begins with a σx or σy at (j, l) and
ends with a σx or σy at (j′, l′) with a string of σz’s in
between, where the forms of the initial and final σ ma-
trices depend on the positions of j + l and j′ + l′. Note
that for J3 → −∞(∞), where the z component of each
spin is locked with that of its vertically nearest neighbor,
〈ψ−∞(∞)|O~r|ψ−∞(∞)〉 = ±δ~r,~0. For the Kitaev model,
it is known that the spin correlations between sites ly-
ing on different bonds vanish, e.g., 〈σza,~nσzb,~n+~r〉 = 0 for
~r 6= ~0 [13]. Therefore 〈O~r〉 are the only non-vanishing
two-point correlators of the model [14]. A non-zero value
of 〈O~r〉 for ~r 6= ~0 in the final state provides a signature of
FIG. 3: Plot of 〈O~r〉, sans the δ-function peak at the origin,
as a function of ~r for four values of J2/J , for Jτ = 5.
the defects. In particular, a plot of 〈O~r〉 versus ~r gives an
estimate of the correlations between the defects. [Since
O2~r = 1, all the moments of O~r can be found trivially:
〈On~r 〉 = 〈O~r〉 if n is odd and = 1 if n is even.]
After the quench, the system, for each momentum ~k,
is described by a combination of ψ−∞~k with probability
p~k and ψ∞~k with probability 1 − p~k, where ψ±∞~k are
the eigenstates of H~k for J3 → ±∞. Hence 〈O~r〉 can be
computed in a straightforward manner:
〈O~r〉 = − δ~r,~0 +
2
A
∫
d2~k p~k cos(
~k · ~r), (5)
where the integral runs over half the BZ with area A.
For large values of τ , the dominant contribution comes
from the region near the line J1 sin(~k · ~M1) = J2 sin(~k ·
~M2) where p~k = 1. Introducing the variables k‖ and k⊥
which vary along and perpendicular to this line (along
the directions n̂‖ and n̂⊥ respectively), we see that the
integrand in Eq. (5) takes the form exp[−a(~k0)τk2⊥ ±
i(~k0 + k‖n̂‖ + k⊥n̂⊥) · ~r], where a(~k0) is a number de-
pending on ~k0. The evaluation of the integral over k⊥
gives a factor of exp
[
−(~r · n̂⊥)2/(4aτ)
]
/
√
aτ . We thus
see that the magnitude of the defect correlations goes
as 1/
√
τ , while the spatial extent of the correlations
goes as
√
τ . This is confirmed by the following relation:
∑
~r ~r
2〈O~r〉 = −2(∇2~kp~k)~k=~0 = 24πτ(J
2
1 + J
2
2 + J1J2)/J .
To study the correlations between defects, we evaluate
Eq. (5) numerically; the ~r dependence of 〈O~r〉 is shown
in Fig. 3 for several values of J2/J1, where J1 = J and
Jτ = 5. Here we have omitted the δ-function peak at
~r = 0 in Eq. 5 so as to make the correlation at ~r 6= ~0
clearly visible. The plot of 〈O~r〉 reflects the defect corre-
lations. To understand the variation of these correlations
with the ratio J2/J1, we note that for large Jτ , the max-
imum contribution to 〈O~r〉 comes from around the wave
4
FIG. 4: Plot of 〈O~r〉 at the points (1, 0) on the n1 axis (black
solid line), (0, 2) on the n2 axis (blue dotted line), and (2,−2)
along the −45◦ line in the n1 − n2 plane (red dashed line) as
a function of α = tan−1(J2/J1), for J
2 = 1 and Jτ = 5.
vectors ~k0 for which p(~k0) = 1. For J2 ≫ (≪)J1, this im-
plies sin[~k · ~M2( ~M1)] = 0 which yields ~k0 ∼
√
3̂i± ĵ. The
maximum contribution to 〈O~r〉 comes from cos(~k0·~r) = 1,
i.e., ~k0 · ~r = 0. Thus for J2 ≫ (≪)J1, 〈O~r〉 is expected
to be maximal along the lines y = −(+)
√
3x, namely,
n1 = −n2 (n1 = 0) in the n1−n2 plane. This expectation
is confirmed in Fig. 3 where 〈O~r〉 can be seen to be max-
imal along n1 = −n2 (n1 = 0) for J2 = 4(0.2)J1. Such a
strong anisotropy can be understood by noting that the
Kitaev model reduces to a one-dimensional model when
J2 ≫ (≪)J1. For intermediate values of J1/J2, a grad-
ual evolution of the defect correlations can be seen in
Fig. 3. We note that if the Kitaev model can be real-
ized using ultracold atoms in an optical lattice [16], such
an evolution of the defect correlations with J1/J2 can,
in principle, be experimentally detected by spatial noise
correlation measurements as pointed out in Ref. [19].
We can obtain a different view of the spatial anisotropy
of the defect correlations by studying 〈O~r〉 as a function
of α = tan−1(J2/J1). As α changes, the ratio J2/J1
varies from 0 to ∞ while fixing J21 +J22 = J2 = 1. A plot
of 〈O~r〉 at three representative points (n1, n2) = (1, 0) (on
the n1 axis), (0, 2) (on the n2 axis), and (2,−2) (along the
−45◦ line in the n1 −n2 plane) as a function of α, shown
in Fig. 4, reveals the nature of the defect correlations.
We see that as J2/J1 = tanα is varied from 0 to ∞, the
magnitude of the correlation at the point (1, 0) on the
n1 axis increases till it reaches a maximum at J1 = J2
(α = π/4), and then decays to 0 as α approaches π/2. For
the points (0, 2) on the n2 axis and (2,−2) along the line
with slope −45◦, the correlation becomes maximum when
J2 ≪ J1 (α = 0) and J2 ≫ J1 (α = π/2) respectively,
as expected from Fig. 3. We conclude that the spatial
anisotropy of the correlations between the defects 〈O~r〉
depends crucially on the ratio J2/J1.
To conclude, we have shown that the density of defects
produced during a quench of the Kitaev model through
a critical line scales with the quench time as 1/
√
τ , in-
stead of the 1/τ behavior expected for a quench through
a critical point. We have provided a general result for
the defect density which reproduces the result of Ref. [4]
as a special case. We have also discussed the variation
of the defect correlations with the model parameters and
pointed out the possibility of detection of these varia-
tions in experiments. These results significantly improve
our general understanding of the scaling of the density of
defects and their correlations in 2D systems.
We thank A. Dutta and A. Polkovnikov for stimulating
discussions.
[1] S. Sachdev, Quantum Phase Transitions (Cambridge
University Press, Cambridge, England, 1999).
[2] T. W. B. Kibble, J. Phys. A 9, 1387 (1976); W. H. Zurek,
Nature (London) 317, 505 (1985).
[3] B. Damski, Phys. Rev. Lett. 95, 035701 (2005).
[4] A. Polkovnikov, Phys. Rev. B 72, 161201(R) (2005); A.
Polkovnikov and V. Gritsev, arXiv:0706.0212.
[5] K. Sengupta, S. Powell, and S. Sachdev Phys. Rev. A
69, 053616 (2004); P. Calabrese and J. Cardy, J. Stat.
Mech: Theory Expt P04010 (2005), and Phys. Rev. Lett.
96, 136801 (2006); J. Dziarmaga, Phys. Rev. Lett. 95,
245701 (2005), and Phys. Rev. B 74, 064416 (2006).
[6] A. Das et al., Phys. Rev. B 74, 144423 (2006).
[7] R. W. Cherng and L. S. Levitov, Phys. Rev. A 73, 043614
(2006); V. Mukherjee et al., Phys. Rev. B 76, 174303
(2007).
[8] B. Damski and W. H. Zurek, Phys. Rev. A 73, 063405
(2006); F. M. Cucchietti et al., Phys. Rev. A 75, 023603
(2007); T. Caneva, R. Fazio, and G. E. Santoro, Phys.
Rev. B 76, 144427 (2007).
[9] L. E. Sadler et al., Nature (London) 443, 312 (2006).
[10] A. Kitaev, Ann. Phys. 321, 2 (2006).
[11] A. Kitaev, Ann. Phys. 303, 2 (2003).
[12] X.-Y. Feng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett.
98, 087204 (2007).
[13] G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev.
Lett. 98, 247201 (2007).
[14] H.-D. Chen and Z. Nussinov, arXiv:cond-mat/0703633.
[15] D.-H. Lee, G.-M. Zhang, and T. Xiang, arXiv:07053499;
K. P. Schmidt, S. Dusuel, and J. Vidal, arXiv:07093017.
[16] L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev.
Lett. 91, 090402 (2003); A. Micheli, G. K. Brennen, and
P. Zoller, Nature Physics 2, 341 (2006).
[17] See for example, L. Landau and E. M. Lifshitz, Quan-
tum Mechanics: Non-relativistic Theory, 2nd Ed. (Perg-
amon Press, Oxford, 1965); S. Suzuki and M. Okada in
Quantum Annealing and Related Optimization Methods,
Eds. by A. Das and B. K. Chakrabarti (Springer-Verlag,
Berlin, 2005).
[18] Note that as we approach the gapless phase, ∆E ∼ |J3 −
J3c| and ∆E ∼ |k|; hence ν = z = 1 for the model [4].
[19] E. Altman, E. Demler, and M. D. Lukin, Phys. Rev. A
5
70, 013603 (2004).


Exact results for quench dynamics and defect production in a
  two-dimensional model

Exact results for quench dynamics and defect production in a two-dimensional model

Abstract

Similar works

Full text

Available Versions

Open Access Repository of IISc Research Publications

Crossref

Name not available

Publications of the IAS Fellows

Indian Academy of Sciences