Article pubs.acs.org/JPCA
Lattice Defects in the Kitaev Honeycomb Model John Brennan† and Jiří Vala*,†,‡ †
Department of Mathematical Physics, Maynooth University, Maynooth, County Kildare, Ireland School of Theoretical Physics, Dublin Institute for Advanced Studies, 10 Burlington Road, Dublin, Ireland
‡
ABSTRACT: The Kitaev honeycomb lattice system is an important model of topological materials whose phase diagram exhibits both abelian and non-abelian topological phases. The latter, a so-called Ising phase, is related to topological superconductors. Its quasiparticle excitations, which are formed by Majorana fermions attached to vortices, show non-abelian fractional statistics and are known as Ising anyons. We investigate dislocation defects in the Ising phase of the Kitaev honeycomb model. After introducing them to the system, we accordingly generalize our solution of this model to the situation with the defects. The important part of this effort is developing an appropriate Jordan−Wigner fermionization procedure. It is expected that the presence of defects manifests itself by the formation of fermionic zero-energy modes around the defect end points. We numerically confirm this expectation and further investigate properties of these modes. The computational potential of our technique is demonstrated for both diagonalization and dynamical simulations. The latter focuses on the process of fusion of the vortex zero-energy modes with the Majorana fermions attached to the defect. This process simulates fusion of non-abelian Ising anyons.
■
Delgado12 and Beigi and co-workers.13 Defects of the honeycomb model have also been studied by Petrova et al.14 who found that some defects carry Majorana fermions. However, this work has focused on the abelian phase of the model. The presented study of defects in the non-abelian phase of the honeycomb system requires a nontrivial generalization of the exact solution of the model without dislocations.9,15 The solution, which is based on Jordan−Wigner fermionization, usefully factorizes into the Bogoliubov−de Gennes description of superconductors and the stabilizer description of the topological +(2 ) vacuum, thus revealing the underlying toric code structure of the model. Consequently, this structure allows us to introduce the defects, studied by Kitaev and Kong within the toric code,10 to all phases of the honeycomb model, including the non-abelian phase in particular. We have implemented a methodology to calculate the energy eigenstates of the model with the defects numerically as their presence breaks the model’s translational invariance. We confirm the expectation that the defects affect the fermionic spectrum of the model the same way a pair of vortices do. Specifically, in the non-abelian phase, the dislocation defects are accompanied by zero-energy modes that are similar to the Majorana zero modes bound to vortices in the model.16 The non-abelian statistical properties of these composite particles have been verified by
INTRODUCTION The Kitaev honeycomb lattice model is an exactly solvable spin− lattice system that exhibits both abelian and non-abelian topological phases. It was first defined by Kitaev,1 who solved the system by a reduction to free fermions in a static 2 gauge field. The abelian phase is precisely the doubled-2 , or +(2 ), topological phase, which is equivalent to the toric code, another spin−lattice model with topological order defined by Kitaev2 in the framework of stabilizer quantum error correction.3 The nonabelian topological phase is of the Ising type and supports nonabelian Ising anyons formed as a Majorana fermion mode bound to a vortex. Majorana fermions, particles, or modes whose creation and annihilation operators satisfy γ = γ†, appear in condensed matter systems with particle−hole symmetry, like superconductors, as the zero-energy modes as implied by γ(E) = γ†(−E).4 Physical realizations of the Kitaev model were considered using trapped atoms in optical lattices5 and polar molecules.6 The model is also relevant to solid state materials like p-wave superconductors.7,8 A review of the topic of topological phases in the relevant class of lattice models can be found for example in Watts et al.9 We are interested in defects of the honeycomb lattice specifically in the non-abelian phase of the model. The defect we have chosen to investigate is a nontrivial dislocation of the lattice which, when tuned to the abelian phase, is equivalent to the toric code defect string, or domain wall, presented recently by Kitaev and Kong.10 Their work stems from the earlier study of boundaries in toric code models by Bravyi and Kitaev.11 Generalizations of the toric code model with boundaries and domain walls have also been investigated by Bombin and Martin© XXXX American Chemical Society
Special Issue: Ronnie Kosloff Festschrift Received: January 6, 2016 Revised: February 17, 2016
A
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A Berry phase calculations17,18 to be consistent with the nonabelian fractional statistics of the Ising anyons. The exact solution of the model also enables us to simulate its dynamics.19 We demonstrate these capabilities by simulating the system with four zero-energy modes, two attached to the end points of a defect and two attached to a pair of vortices. During the dynamical process, one of the vortices is translated across the lattice from its original position to one of the two end points. At the end of the process, both zero modes are annihilated as the zero-energy modes and their wave functions combine and “melt” into the bulk fermionic spectrum. This process is related to the fusion of the Ising anyons.20
Vp = σ2yσ1zσ6x + σ3zσ2xσ1y + σ4xσ3yσ2z + σ5yσ4zσ3x + σ6zσ5xσ4y + σ5zσ6yσ1x
(2)
where we have numbered the vertices of a plaquette as in Figure 2a.
■
THE MODEL The system consists of spin-1/2 particles arranged on vertices of a honeycomb lattice, interacting with nearest neighbors via a twobody interaction. Specifically, each spin is coupled to its neighbors by one of three links (x, y, z), depending on its orientation. For a link of orientation α, the spins at the ends of the link interact through Kα = σα• σα° where the subscripts refer to one or the other sublattice of the bipartite honeycomb system, as illustrated in Figure 1. The basic Hamiltonian for the system is given as H0 = −Jx
∑ x ‐links
σ•xσ x − Jy °
∑ y ‐links
σ•yσ y − Jz °
∑ z ‐links
σ•zσ z °
Figure 2. Sites of hexagonal plaquettes are numbered as depicted on the left. Note that the plaquettes along the defect line, not including the ends, all have six edges and so the sites of these plaquettes are numbered in the same way. The sites of the plaquettes at the top and bottom of the defect line are numbered as depicted in the center and on the right, respectively.
The defect is introduced to the honeycomb model by decoupling a string of spins from their nearest neighbors and then recoupling the neighbors with each other, as illustrated in Figure 1. The spins that are decoupled from the rest of the system are no longer considered to be part of the system and so their degrees of freedom are projected out of the Hilbert space and are said to have been deleted. The plaquettes along the defect line, not including the ends, all have six vertices coupled to each other the same way the vertices of plaquettes away from the defect are; therefore, their contribution to the time-reversal and paritybreaking potential will also be given by (2). The contribution to the potential from the plaquettes at the ends of the defect are given as
(1)
where the Ji are corresponding coupling constants.
Vp = σ2yσ1zσ8x + σ3zσ2xσ1y + σ4xσ3yσ2z + σ5yσ4zσ3x + σ6zσ5xσ4y + σ7yσ6xσ5z + σ8zσ7xσ6y + σ1xσ8yσ7z
(3)
where, if p is the top plaquette, the subscripts label the sites as in Figure 2b and if p is the bottom plaquette they label the sites as in Figure 2c. The full Hamiltonian for the system is then formally given as H = H0 + κ ∑ Vp p
Figure 1. Spins are coupled to their nearest neighbors by an interaction that depends on the orientation of the corresponding link, (x, y, z). The defect is introduced by decoupling a chain of spins (circled in red) and recoupling the remaining spins as shown. The length of the defect is the number of z-links deleted along this line.
(4)
where κ is a constant and the sum is over all plaquettes of the system. We can define a vortex operator Wp for each plaquette p of the lattice. If p is not one of the plaquettes at either end of the defect, then the plaquette is a hexagon and we define the vortex operator for p as
In this form, the system exhibits a simple phase diagram with three equivalent abelian topological phases, which can be mapped onto the toric code, and a gapless phase. Adding to this Hamiltonian a time-reversal and parity-breaking potential opens a spectral gap in the originally gapless phase and turns it into the non-abelian topological phase of the Ising type. The time-reversal and parity-breaking terms are given as a sum of three-body spin terms. The contribution from each plaquette will take the following form
Wp = σ1zσ2xσ3yσ4zσ5xσ6y
(5)
However, if p is one of the plaquettes at either end of the defect, then it has eight edges and we define the vortex operator as Wp = σ1zσ2xσ3xσ4xσ5yσ6zσ7xσ8y
(6)
The vortex operators Wp commute mutually and with the full Hamiltonian, including the time-reversal and parity-breaking potential terms. Consequently, the Hilbert space can be written B
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A as / = ⊕ /{wp}, where /{wp} is the common eigenspace of
H0 = −Jx
{wp}
∑
(bq†1 + bq1)τqx2(bq†2 + bq2)
x ‐links
the Wp operators corresponding to the particular configuration of eigenvalues {wp}. We say that a vortex occupies the plaquette p if wp = −1.
− Jy
∑
iτqz1(bq†1 − bq1)τqy2(bq†2 + bq2)
y ‐links
■
− Jz ∑ (I − 2Nq)
SOLUTION OF THE MODEL WITH DEFECTS We use the exact solution of the model described in ref 15 to transform our honeycomb lattice model with a defect into a square lattice of spins with hardcore bosons living on the vertices, also with a defect, as depicted in Figure 3. During this
(8)
q
b†q
where the operators and bq are creation and annihilation operators for the hardcore boson at the lattice site q, Nq is the bosonic number operator at site q and the operators ταq (α = x, y, z) are the Pauli spin operators for the effective spin at site q. In this basis, the vortex operators become ⎧(I − 2Na)(I − 2Nb)(τazτ yτczτdxτey) b ⎪ ⎪ if p is at the top of the defect, ⎪ ⎪(I − 2Na)(I − 2Nc)(τazτbxτcyτdzτey) Wp = ⎨ ⎪ if p is at the bottom of the defect, ⎪ z y z y ⎪(I − 2Na)(I − 2Nb)(τa τb τc τd ) ⎪ ⎩ otherwise,
Figure 3. Mapping the honeycomb lattice can be visualized as shrinking the z-links of the lattice to a point, creating a square lattice while the degrees of freedom of each z-link are mapped to the degrees of freedom of the square lattice according to eq 7. The resultant defect in the square lattice is equivalent to the toric code defect introduced in ref 10 via a unitary transformation of the vortex operators described by Watts et al.9
(9)
where we label the sites of the plaquettes on the square lattice as in Figure 4. Moreover, in this basis, the potential terms for the plaquettes at the ends of the defect line are Vp = τbyτazτex(I − 2Na)(bb† + bb)(be† + be) + iτbzτbxτaz(I − 2Nb)(bb† + bb)(ba† − ba)
transformation, the spin degrees of freedom of each z-link of the honeycomb are relabeled as the effective spin and hardcore boson degrees of freedom of the square lattice as follows:
+ iτcx(bc† + bc)(bb† − bb) + iτdzτcz(bd† − bd)(bb† + bb) + τdzτcy(bd† + bd)(bc† + bc) + iτezτdxτdz(be† − be)(bd† + bd)(I − 2Nd) + τezτdy(be† + be)(bd† + bd)
where 0 and 1 refers to the boson occupation number. This induces a change of basis for the operators on the Hilbert space in which the Hamiltonian becomes
+ τeyτez(ba† + ba)(be† + be)(I − 2Ne)
(10)
Figure 4. After mapping the z-links of the honeycomb to the sites of the square lattice, we letter the sites of the plaquettes as above. C
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A and the expression for Vp when p is not at either end of the defect can be found in Kells et al.15 We can now fermionize the bosons by defining a Jordan− Wigner type string operator for each site of the lattice.15 The composition of these string operators with the boson creation and annihilation operators will be fermionic creation and annihilation operators. Then changing basis again so that the Hamiltonian is expressed in terms of these new operators will effectively transform the hardcore bosons of the model into fermions. To define a string operator for a site q, we first pick a particular site of the lattice as a reference point, labeled O. Then we decide on a procedure to construct a set of operators which act on sites of the lattice forming a string connecting O to q, the string operator will be the composition of all these operators. We use the same reference point and procedure to define the string operator for each site on the lattice. To provide an example, we can consider our model on a torus by taking a Nx × Ny lattice with the defect situated away from the boundary and then identifying the opposite edges. Following ref 15, we choose the lower left-hand corner of the square lattice as our reference point O and give each site a pair of coordinates (qx, qy), qx being the distance in the x-direction, in units of lattice spacing, between O and q, and qy being the distance in the ydirection. Note, in this coordinate system, the point O is the origin (0, 0). Also note, although every site on the lattice has a pair of coordinates, not all pairs of coordinates (qx, qy), with qx ∈ [0, Nx] and qy ∈ [0, Ny], correspond to a site on the lattice because some sites have been deleted along the defect. To write down the string operator for a site q = (qx, qy), we use the following procedure: beginning at the origin, we move along the x-direction to the site (qx − 1, 0) and consider the set of operators {−(I − 2Ni)τxi } which contains a single operator for each site i that we cross along the way. Then, at the corner of the y string (qx, 0) we add the operator −τ(q to our set of operators. x,0) Now, moving along the y-direction, we add the operator τxi for each site we cross to reach the site q, and for q itself we add the operator τyq. The string operator for q is then defined as the composition of these operators. As an example, the definition of the string operator for the site (2, 6) is given as
Figure 5. String operator S(2,6) for the site (2, 6) is shown above. Note, the composition of the string operator with the boson creation/ annihilation operator for a given site will anticommute with that of other sites because either the string of operators connecting the sites to the origin will have anticommuting Pauli operators at exactly one site, or one of the creation/annihilation operators will anticommute with one of the factors (I − 2Ni) of a string operator.
iτzq1Sq1τyq2Sq2. Note, the string operators square to the identity and so the number operator Nq = b†q bq = c†q cq is unchanged when expressed in fermion operators and hence the vortex operators are also left unchanged. The basic Hamiltonian is now quadratic in fermionic operators and can be written as 1 H0 = 2
Δq , q ′ = Jx Xq , q ′(δq(,xq)′ − δq(′x,)q) + Jy Yq , q ′(δq(,yq)′ − δq(′y,)q) (15)
(11)
cq = bqSq
(x) Here, δq,q′ is the usual kronecker delta and δq,q′ is defined to be 1 if
q and q′ are the sites on the left- and right-hand side of an x-link, respectively, and zero otherwise. Similarly, δ(y) q,q′ is 1 if q and q′ are the sites on the bottom and top of a y-link, respectively, and zero otherwise. The potential also becomes quadratic in fermionic operators in this basis and can be written as
(12)
Expressing the basic Hamiltonian in terms of these fermionic creation and annihilation operators yields H0Jx
∑
Xq1, q2(cq†1
−
cq1)(cq†2
1 V = ∑ Vp = 2 p
+ cq2)
x ‐links
+ Jy
∑ y ‐links
(14)
ξq , q ′ = Jz δq , q ′ + Jx Xq , q ′(δq(,xq)′ + δq(′x,)q) + Jy Yq , q ′(δq(,yq)′ + δq(′y,)q)
and is illustrated graphically in Figure 5. We can now define fermionic creation and annihilation operators for a site q by composing the string operator for q with the boson operators associated with q as follows: cq† = bq†Sq
q,q′
⎛ ξq , q ′ Δq , q ′ ⎞⎛ cq ′ ⎞ ⎟⎜ ⎟ cq )⎜ ⎜ Δ† −ξ T ⎟⎜ c † ⎟ q , q ′ ⎠⎝ q ′ ⎠ ⎝ q,q′
where the sum is over the sites q of the lattice and
x x S(2,6) = −(I − 2N(0,0))τ(0,0) [−(I − 2N(1,0))τ(1,0) ] y y x x x x × [−τ(2,0) ]τ(2,1) τ(2,2) τ(2,4) τ(2,5) τ(2,6)
∑
( cq†
∑ q,q′
( cq†
⎛ ξ̅ Δ̅ q , q ′ ⎞⎛ cq ′ ⎞ ⎜ q,q′ ⎟⎜ ⎟ cq )⎜ † T ⎟⎜ c † ⎟ ⎝ Δ̅ q , q ′ −ξ q̅ , q ′⎠⎝ q ′⎠
(16)
where
Yq1, q2(cq†1 − cq1)(cq†2 + cq2) + Jz ∑ (2Nq − I )
ξq̅ , q ′ = i ∑ Xq , ρYρ , q ′( −δq(,xρ)δq(′y,)ρ + δq(′x,)ρδq(,yρ) + δρ(x, q)′δρ(y, q)
q
(13)
ρ
− δρ(x, q)δρ(y, q)′)
where, if q1 and q2 are sites on the left- and right-hand side of a xlink, respectively, Xq1,q2 = −(I − 2Nq1)Sq1τxq2Sq2 and, if q1 and q2 are sites at the bottom and top of a y-link, respectively, Yq1,q2 =
(17)
and D
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
Figure 6. (a) Formation of zero-energy modes as the system crosses the phase transition from the abelian to the non-abelian phase. The results for the case with the defect and the case with two vortices replacing the defect are essentially indistinguishable. (b) Fermion probability density at J = 0.7 with the Majorana modes localized at the defect end points.
Δ̅ q , q ′ = i ∑ Xq , ρYρ , q ′(δq(,xρ)δq(′y,)ρ − δq(′x,)ρδq(,yρ) + δρ(x, q)′δρ(y, q) − δρ(x, q)δρ(y, q)′)
Here, the quantum numbers lx and ly are the eigenvalues of the two loop symmetries Lx and Ly, which are given by homologically nontrivial loops on the torus.15 Now, if we consider a particular vortex configuration, i.e., a configuration of eigenvalues {wp}, we can replace the Xq,q′ and Yq,q′ operators in (15), (17), and (18) by their eigenvalues and diagonalize the matrix in (19) numerically to find the fermionic spectrum in the Hilbert space /{w p}.
ρ
− 2iXq , q ′(δq(,xq)′ − δq(′x,)q) + 2iYq , q ′(δq(,yq)′ − δq(′y,)q)
(18)
Hence, the full Hamiltonian is 1 H= 2
∑
( cq†
q,q′
⎛ ξ + κξ ̅ Δq , q ′ + κ Δ̅ q , q ′ ⎞⎛ cq ′ ⎞ q,q′ ⎜ q,q′ ⎟⎜ ⎟ cq )⎜ † T ⎟⎜ c † ⎟ † T ⎝ Δq , q ′ + κ Δ̅ q , q ′ −ξq , q ′ − κξ q̅ , q ′⎠⎝ q ′⎠
■
NUMERICAL CALCULATIONS The analytical representation of the model with defects described above opens up the possibility of further investigation of this system. This has to proceed by numerical means as the defect breaks the translational invariance of the model. We first investigate how the presence of the defect affects the spectrum of the model in its non-abelian phase. Specifically, we wish to confirm the expectation that the defect leads to an emergence of zero modes similar to those we observe in the case of vortex excitations.16 We also expect that the number of zero modes will be the same as the number of the defect end points like in the case of vortices.21 We then proceed to demonstrate the numerical capabilities of our methodology on the system consisting of four zero-energy modes attached to two end points of a defect and a pair of vortices. We study the energy spectrum of this system for various configurations of vortices. We then simulate the dynamics of the fusion of a vortex zero mode with a zero mode at the defect end point in real time. This process is equivalent to fusion of the non-abelian Ising anyons. Basic Setup. We consider the full Hamiltonian of the model with the defect on the torus. It includes the time-reversal and parity-breaking potential and thus the corresponding quantum phase diagram exhibits both the abelian +(2 ) toric code phase and non-abelian Ising phase. On the torus, homologically nontrivial loop symmetry operators Lx and Ly lead to four distinct symmetry sectors labeled by the eigenvalues lx = ±1 and ly = ±1, respectively. In our calculations, we focus only on the sector where both these eigenvalues are negative. The fermionic representation of the model permits favorable scaling of the numerical diagonalization with the system size. Our numerical computations are readily carried out for a square lattice with 80 × 80 plaquettes. The defect, if present, is positioned diagonally across the center with its end points
(19)
The operators Xq,q′ and Yq,q′, appearing in the basic Hamiltonian, can be expressed as products of vortex operators. As discussed in ref 15, the product of two string operators associated with neighboring sites can result in a sequence of operators acting on a closed loop of sites that will be equal to the product of vortex operators associated with the plaquettes enclosed by the loop. If we identify each plaquette p with the site on the bottom left corner q (so now Wp = W(qx,qy)) and identify each x-link with the site on its left-hand side and each y-link with the site on the bottom (so now Xq,q′ = X(qx,qy) and Yq,q′ = Y(qx,qy)), we find ⎧ qy − 1 ⎪∏ W (qx , qi) ⎪ ⎪ qi = 0 ⎪1 ⎪ X(q , q ) = ⎨ x y qy − 1 ⎪ ⎪−lx ∏ W(q , q ) x i ⎪ qi = 0 ⎪ ⎪−lx ⎩
if qy ≠ 0 and qx ≠ Nx − 1 if qy = 0 and qx ≠ Nx − 1 if qy ≠ 0 and qx = Nx − 1 if qy = 0 and qx = Nx − 1 (20)
⎧1 ⎪ ⎪ Y(q , q ) = ⎨ x y ⎪− l y ⎪ ⎩
if qy ≠ Ny − 1 Ny − 1 qx − 1
∏ ∏ W(q , q ) j
qi = 0 qj = 0
i
if qy = Ny − 1 (21) E
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A separated by 40 plaquettes in both the x and y directions. The defect end points are thus positioned at the plaquettes with the coordinates (20, 20) and (60, 60). For the purpose of the scaling study, described further below, we extended these diagonalizations of the system up to a size of 700 × 700 plaquettes. Formation of Majorana Zero Modes. Spectrum and Phase Transition. We first numerically diagonalize the Hamiltonian for the 80 × 80 lattice with the diagonal defect described above while driving the system across the transition from the abelian to the non-abelian phase. Specifically, we fixed the coupling Jz = 1 and varied equally the other two couplings J = Jx = Jy from the value 0.45−0.7. The time-reversal and paritybreaking potential was applied with the coupling constant κ = 0.2. The phase transition is exactly at J = 0.5 in the vortex-free case. In the presence of vortices, it is known22 that the exact phase boundary of the non-abelian phase lies between J = 1/2 and J = 1/ √2 with the latter value being relevant to the full vortex lattice. In the thermodynamic limit, we would expect to see the formation of zero modes exactly at the phase transition. The calculations of the energy spectrum in the vicinity of E = 0 reveal that zero modes start forming as the system enters the non-abelian phase. They correspond to two eigenstates represented by fermionic wave functions, each with two fermionic modes sharply localized at the end points of the defect. The isolated Majorana modes can be obtained by taking an appropriate superposition of both near zero-energy eigenstates. As in the case of vortex zero modes, the energy of the two eigenstates converges to zero with increasing value of the coupling J, which drives the phase transition. At the same time, the spectral gap in the non-abelian phase increases. This leads to shorter penetration depth of the Majorana modes and their better localization. This behavior is essentially identical to the situation where the Majorana zero modes are localized around vortices.16 To see this quantitatively, we performed the same computation for the system without the defect and with a pair of vortices replacing the defect end points. The spectral results in both cases are shown in Figure 6. Scaling at the Critical Point. We would expect energy states near E = 0 to become the zero-energy modes exactly when the system passes the transition to the non-abelian phase. This is indeed strictly speaking the case if the system is of infinite size. For any finite size, the phase transition is necessarily smeared out by the finite size effects. To see that our calculations correctly reproduce the expected behavior in the thermodynamic limit, we investigated the scaling with the system size of the lowest positive energy eigenvalue of the model at the phase transition. We fixed the model parameters to Jz = 1, J = 0.5, and κ = 0.2 and diagonalized the system for N × N lattices with the linear sizes ranging from N = 200 to 700 in steps of 20. When increasing the size of the lattice, we kept the length of the defect string constant with the end points separated by 40 plaquettes in both directions. At this separation, the contribution to the overall spectrum from the mutual interaction between the two defect zero modes can be considered negligible. We plot the results of these calculations in Figure 7. The lowest positive energy eigenvalue of the system converges to zero linearly with decreasing 1/N2. This behavior indicates conformal invariance of the quantum critical point. The exact location of this point is modified by the presence of vortices and is exactly at the minimal value of J = 0.5 only for the vortex-free situation.22 Similarly, in the presence of the defect, the phase transition is slightly shifted relative to the vortex-free value. This is confirmed
Figure 7. Scaling of the energy of the lowest positive state with the system size in the presence of the defect at the phase transition between the abelian and non-abelian topological phase.
by detailed examination of the scaling plot in doubly logarithmic scales for various values of J around 0.5, which reveals that the exact position of the phase transition is at J = 0.50005. This shift, though quite subtle, is observable within our methodology, which can access a very large lattice sizes. Interaction and Fusion of Zero Modes. We now study how vortex and defect zero-energy modes interact and fuse. We first investigate the spectrum and the eigenstates of the system with a defect and a pair of vortices. Second, we simulate the dynamics of the fusion process by propagating the zero-energy fermionic modes attached to both the vortices and defect end points under a time varying Hamiltonian that transports one of the vortices from its initial position to the location of one of the defect end points. The latter approach is equivalent to timedependent simulation of the fusion of two Ising anyons. Spectrum and Eigenstates. We calculate the spectrum of the system in the presence of two vortex zero-energy modes and two zero modes attached to the defect end points. The vortices on neighboring plaquettes are induced by flipping the value of the coupling on the lattice link that is shared by both plaquettes from J to −J. The defect is fixed as above. We diagonalize the system for various configurations of the vortices starting with the creation of vortices on the lattice at plaquettes (20, 59) and (20, 60). We then move the first vortex vertically downward by flipping the consecutive couplings toward the end of the defect at (20, 20), and each time we diagonalize the system. As the vortex moves away from the vortex fixed at the position (20, 60), the energy of the two corresponding eigenstates converges to zero as the interaction between vortices diminishes. The energy starts diverging from zero as the vortex starts approaching the defect at (20, 20). In both cases the energy splitting originates from the interaction between the zero-energy modes. The energy of both states remains at ±0.4 after the zero modes are fused. The vortex−vortex and vortex−defect interaction energies as a function of distance are shown in Figure 8. Dynamics of Zero Mode Fusion. The initial states for the simulation of the fusion dynamics consists of all four zero-energy eigenstates of the system with the defect located diagonally with the end points at the plaquettes (20, 20) and (60, 60), and a pair of vortices located at the plaquettes (20, 60) and (20, 50). The F
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
The vortex zero mode at the plaquette (20, 50) is transported, under an intrinsically time-dependent Hamiltonian, vertically from its initial position to its final destination at the plaquette (20, 20), where one of the defect end points is located. The basic system parameters of these simulations are Jz = 1, J = 0.7, and κ = 0.2. We again choose the periodic boundary conditions with lx = ly = −1. The plaquette-to-plaquette transport of the vortex zero mode is achieved by slowly varying the coupling strength in the Hamiltonian from J to −J on the lattice link that is shared by both plaquettes.17,18 The average speed of vortex motion is 0.003 lattice constants per unit time, and the time step is 0.05. This corresponds to the change of J on a given link in each propagation step by ΔJ = 3 × 10−4 and the total of 2 × 105 steps for the entire propagation from the initial state to the final state where both zero modes fuse. To refer to different stages of propagation, we introduce the parameter T, which is defined in such a way that it reaches T = 1 when the propagation is complete. The results of this simulation are documented in Figure 9, which shows the evolution of the four fermionic initial wave functions at different stages of propagation. The top two sets correspond to the dynamics of the initial states whose zero modes are localized around the vortices. The bottom two sets
Figure 8. Change of the energy eigenvalues of the six eigenstates closest to E = 0 for various vortex configurations. One of the vortices is kept fixed at its initial position at the plaquette (20, 60) while the other is moved toward the defect end point at the position (20, 20) .
four zero-energy eigenstates appear in two pairs; one pair contains zero modes localized at the vortices and the other at the defect end points.
Figure 9. Four zero-energy eigenstates appearing in two pairs: one pair contains zero modes localized at the vortices (ψ1 and ψ2) and the other at the defect end points (ψ3 and ψ4). Above we show the fermion wave functions at various stages of the propagation of the four initial zero-energy eigenstates. The vertical axis indicates the modulus of the wave function and the color shows its phase. The parameter T indicates the stage of the propagation and equals to one when the propagation is complete. G
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
■
refer to the initial states with the defect zero modes. We observe the translation of one of the initial vortex zero modes toward the defect end point. At T = 0.6, we observe doubling of the fermionic peaks that comes from the interaction of fermions attached to the mobile vortex and the stationary defect end point. Both peaks eventually fuse into a single peak and eventually disappear as localized zero modes while being emitted into the bulk fermionic wave function in circular waves. The other two zero modes, attached to the stationary vortex and to the other defect end point remain unchanged throughout the process. It is interesting to look at the final stages of the fusion in detail. For this purpose, we show the annihilation of a single Majorana fermion located at the same defect end point as above. The initial state is taken to be the appropriate superposition of the zeroenergy eigenstates that show the zero modes localized at the defect end points (ψ3 and ψ4 in Figure 9). The final stages of the life of the Majorana fermion are illustrated in Figure 10. While
Article
CONCLUSIONS
We introduced lattice defects into the Ising topological phase of the Kitaev honeycomb model using its underlying toric code structure, and we accordingly generalized our analytical solution of this model for the situation including these defects. The crucial step in this effort was the formulation of an appropriate Jordan− Wigner fermionization procedure that enabled us to solve the system with defects as a quadratic fermionic problem using techniques developed in the theory of superconductivity. We believe that introducing dislocation defects into the nonabelian topological phase opens up interesting opportunities. For example, they enable a more complex fractional statistics as suggested by Kitaev1 and shown by Petrova et al.14 in the abelian phase of the honeycomb model. They are also related to interesting algebraic structures that generalize the braid group like those in ref 23 for example. And, on a more abstract level, they are also connected to extended topological quantum field theories.24 One particularly promising application we envision in the context of defects in topological lattice models stems from the possibility of realizing topological phases on higher genus surfaces. Though a honeycomb lattice can only realize a genus g = 1 surface like a torus, appropriately introducing defects may change the Euler characteristic of the lattice and permit the generalization to surfaces of different genii. Our analytical representation of the honeycomb model with defects moreover opens up the possibility of further investigation of this system by numerical means. We have demonstrated the numerical capabilities in the context of both stationary diagonalizations and also in dynamical calculations. We first focused on showing that the defect end points carry fermionic zero modes which are equivalent to zero modes attached to vortices. These composite particles are known to exhibit nonabelian fractional statistics of the Ising anyons. Our scaling study demonstrated that the solution of the model permits diagonalization of this system up to the size 700 × 700 plaquettes using readily available computers. A further increase of the scale of these computations is possible using highperformance computing resources. The presented numerical effort culminated with the real-time simulation of the vortex motion on the lattice and fusion of the vortex zero-energy modes with the zero-energy modes associated with the defect. We believe that these are the first numerical simulations of the dynamics of these quasiparticles in real time. We specifically focused this effort on the process of fusion of the vortex and defect Majorana fermions, which is equivalent to simulating the fusion of non-abelian Ising anyons. The results of our dynamical calculations also demonstrates that we can effectively manipulate the Majorana modes attached to static defects, essentially being able to turn them on and off in real time simulations. We believe that all these dynamical processes open up interesting applications in topological quantum information processing.20 A physical realization of the Kitaev honeycomb lattice model was proposed using neutral atoms trapped in an optical honeycomb superlattice.5 In this system, the trapped atoms carry spin-like degrees of freedom whose anisotropic interaction is mediated by suitably tuned polarized laser fields. Another implementation of the Kitaev model is also possible using polar molecules trapped in an optical lattice via a proposed toolbox for the simulation of condensed matter systems.6 This toolbox can also be augmented to include three-body interactions, which are
Figure 10. Annihilation of the Majorana zero mode at the defect end point (20, 20) during the last stages of the propagation.
the Majorana fermion gradually broadens, the initial superposition of zero modes, which has constituted its initial state, splits while opening an excitation gap as shown in Figure 11. The expectation value of the energy in each step of the process remains essentially zero (within numerical accuracy). Eventually the entire process dramatically terminates with the Majorana zero mode “melting” into the bulk fermionic population by emitting circular waves.
Figure 11. Transition probability to the eigenstates of the timedependent Hamiltonian at various stages of propagation of the isolated defect Majorana zero mode. The initial state (black), which is an equal superposition of two zero-energy eigenstates, splits in the final phase of the propagation (T = 0.965) (blue) and eventually “melts” into the fermionic bulk population (red). H
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A required to break time-reversal symmetry.25 There are also several types of solid-state materials that may exhibit topological phases that are in the same universality class as the non-abelian phase of the Kitaev model. These are represented by superconductors with p-wave pairing, like the experimentally studied strontium ruthenate.7,8
■
(19) Kosloff, R. Propagation Methods for Quantum Molecular Dynamics. Annu. Rev. Phys. Chem. 1994, 45, 145. (20) Nayak, C.; Simon, S. H.; Stern, A.; Freedman, M. H.; Das Sarma, S. Non-Abelian anyons and topological quantum computation. Rev. Mod. Phys. 2008, 80, 1083. (21) Ivanov, D. A. Non-Abelian Statistics of Half-Quantum Vortices in p-Wave Superconductors. Phys. Rev. Lett. 2001, 86, 268. (22) Lahtinen, V.; Kells, G.; Carollo, A.; Stitt, T.; Vala, J.; Pachos, J. K. Spectrum of the non-abelian phase in Kitaev’s honeycomb lattice model. Ann. Phys. 2008, 323, 2286. (23) Burella, G.; Watts, P.; Pasquier, V.; Vala, J. Graphical Calculus for the Double Affine Q-Dependent Braid Group. Ann. Henri Poincare 2014, 15, 2177. (24) Kapustin, A. Topological Field Theory, Higher Categories, and Their Applications, In Proceedings of the International Congress of Mathematicians; Hindustan Book Agency: New Delhi, 2010; Vol III, pp 2021−2043. (25) Büchler, H. P.; Micheli, A.; Zoller, P. Three-body interactions with cold polar molecules. Nat. Phys. 2007, 3, 726.
AUTHOR INFORMATION
Corresponding Author
*J. Vala. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Graham Kells and Paul Watts for inspiring discussions and critical comments to the manuscript. Our research is supported by the Science Foundation Ireland under Principal Investigator Award No. 10/IN.1/I3013.
■
REFERENCES
(1) Kitaev, A. Anyons in an exactly solved model and beyond. Ann. Phys. 2006, 321, 2. (2) Kitaev, A. Fault-tolerant quantum computation by anyons. Ann. Phys. 2003, 303, 2. (3) Kitaev, A.; Shen, A. H.; Vyalyi, M. N. Classical and Quantum Computation; Graduate Studies in Mathematics; American Mathematical Society: Providence, Rhode Island, 2002. (4) Elliott, S. R.; Franz, M. Majorana fermions in nuclear, particle, and solid-state physics. Rev. Mod. Phys. 2015, 87, 137. (5) Duan, L.-M.; Demler, E.; Lukin, M. D. Controlling Spin Exchange Interactions of Ultracold Atoms in Optical Lattices. Phys. Rev. Lett. 2003, 91, 090402. (6) Micheli, A.; Brennen, G. K.; Zoller, P. A toolbox for lattice-spin models with polar molecules. Nat. Phys. 2006, 2, 341. (7) Maeno, Y.; Rice, T. M.; Sigrist, M. The Intriguing Superconductivity of Strontium Ruthenate. Phys. Today 2001, 54, 42. (8) Nelson, K. D.; Mao, Z. Q.; Maeno, Y.; Liu, Y. Odd-Parity Superconductivity in Sr2RuO4. Science 2004, 306, 1151. (9) Watts, P.; Kells, G.; Vala, J. From Topological Quantum Field Theories to Topological Materials. In Quantum Information and Computation for Chemistry; Kais, S., Ed.; Advances in Chemical Physics; Rice, S. A., Dinner, A. R., Series Eds.; John Wiley and Sons: New York, 2014; Vol. 154. (10) Kitaev, A.; Kong, L. Models for Gapped Boundaries and Domain Walls. Commun. Math. Phys. 2012, 313, 351. (11) Bravyi, S. B.; Kitaev, A. Quantum codes on a lattice with boundary, arXiv:quant-ph/9811052. (12) Bombin, H.; Martin-Delgado, M. A. Family of non-Abelian Kitaev models on a lattice: Topological condensation and confinement. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 115421. (13) Beigi, S.; Shore, P. S.; Whalen, D. The quantum double model with boundary: condensations and symmetries. Commun. Math. Phys. 2011, 306, 663. (14) Petrova, O.; Mellado, P.; Tchernyshyov, O. Unpaired Majorana modes on dislocations and string defects in Kitaev’s honeycomb model. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 134404. (15) Kells, G.; Slingerland, J. K.; Vala, J. Description of Kitaev’s honeycomb model with toric-code stabilizers. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 125415. (16) Kells, G.; Vala, J. Zero energy and chiral edge modes in a p-wave magnetic spin model. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 125122. (17) Bolukbasi, A.; Vala, J. Rigorous calculations of non-Abelian statistics in the Kitaev honeycomb model. New J. Phys. 2012, 14, 045007. (18) Lahtinen, V.; Pachos, J. K. Non-Abelian statistics as a Berry phase in exactly solvable models. New J. Phys. 2009, 11, 093027. I
DOI: 10.1021/acs.jpca.6b00149 J. Phys. Chem. A XXXX, XXX, XXX−XXX
