Kitaev model on the honeycomb and diamond lattices Majorana hopping and vortex excitations in 2D and 3D

(1)

University of Amsterdam

MSc Physics

Theoretical physics

Master Thesis

Kitaev model on the honeycomb and diamond lattices

Majorana hopping and vortex excitations in 2D and 3D

by

Ludo Nieuwenhuizen

5970326

December 2015

60 ECTS

Research carried out between February 2014 and December 2015

Supervisor:

Prof. Dr. Kareljan Schoutens

Dr. Philippe Corboz

Examiner:

(2)

(3)

Abstract

The Kitaev model on the 2D honeycomb lattice, as introduced by Alexei Kitaev in [1], is explained and through the exact diagonalisation of the Kitaev Hamiltonian on a periodic lattice the influence of fluxes penetrating the lattice plaquettes on the low energy excitations is determined. It is argued that the model has a non-Abelian topological phase that shows zero modes that can be identified as belonging to the Ising anyon model. This work was already done in [2] and the current thesis reproduces the same results, based on own computations.

The Kitaev model on a 3D diamond lattice, as proposed by Shinsei Ryu in [3], is also explained. This ‘Ryu-Kitaev model’ respects both time-reversal and particle-hole symmetries, putting it in the topological symmetry class DIII of free fermion models, introduced in [4, 5]. Both analytic and numerical results of energy gap-values are presented for the periodic system in the vortex free phase.

The new work in this thesis consists of deriving explicit Hamiltonian terms that fully generate the Kitaev model on the periodic diamond lattice for arbitrary large unit cell for arbitrary vortex sector. It is shown that upon opening one boundary, the energy gap closes for thick systems, and opens when the length in the open lattice direction is minimalised. It is also shown that applying parallel vortex loops to the periodic lattice lowers the energy gap, but the existence of degenerate zero modes due to the vortex loops cannot be concluded from our current data.

(4)

(5)

Introduction and summary

for a non-physics audience

Theoretical physics is concerned with the rationalisation, explanation and prediction of physical phenomena. A theoretical physicist uses mathematical models and the abstraction of physical phenomena in its work to obtain a better understanding of the laws of nature.

A physical theory can have the property that it gives an explanation to a phenomenon that is observed in experiments, thus embedding the phe-nomenon in a framework that is valid also outside the laboratory where the phenomenon is observed. Theoretical physics can also provide models and ideas that challenge the existing views and can lead to new experiments and the finding of new phenomena.

This thesis is about the research that I have done for my master project on topological models in theoretical condensed matter physics. I will explain what these terms mean and explain what I have done during my research.

Condensed matter and phases

Condensed matter physics is the field that is concerned with solid and liquid phases of matter and their properties. These liquid and solid phases are typically found when particles have little energy to move around. While the classical phases solid and liquid are well known from high school, many other phases can be identified in the domain of condensed matter physics that, for instance, result from temperature changes, quantum properties that find their origin at tiny length scales, or magnetic fields. In their properties, these phases can differ substantially and in a sense understanding these differences between phases are what drives condensed matter phycists. To give some idea as to what these phases in condensed matter might be, I give two examples.

The first example is the phase of ordinary bar magnets: the ferronetic phase. This phase is characterised by the collective behavior of mag-netic spins1 that are present in a material (iron for instance). In the

(6)

(a) Ferromagnetic spin configuration (b) Anti-ferromagnetic spin configuration

Figure 1: The spin configurations that are behind the ferro- and antiferro-magnetic phase.

magnetic phase all the spins point in the same direction, which makes the tiny magnetic fields of the individual spins add up to a total magnetic field that you can feel on large scales. The magnetic field of a bar magnet or horse shoe-magnet is due to this ferromagnetic phase.

The spins that we just described can also be in an anti -parallel config-uration, such that every spin has its magnetic field pointing opposite to its neighbours. This phase of anti-parallel spins is called the anti-ferromagnetic phase and has no net magnetic field, since all tiny fields are cancelled by their neighbours. The ferromagnetic and anti-ferromagnetic phases are depicted in Figure [1].

A second example is that of the superconducting phase. Ordinary con-ducting materials have a resistance to concon-ducting electricity. This resistance absorbs energy and limits the amount of current the material can transport. By cooling certain materials down to below a critical temperature (which is usually a little above absolute zero, 0 degrees Kelvin, or −273◦C), these materials can suddenly become superconducting. When a material is super-conducting, it no longer has any resistance, so electricity can flow without energy loss. This superconducting phase is induced by the drop in temper-ature, but also relies on material properties (a piece of copper string will never become superconducting, not even when you cool it down to absolute zero). The topic of superconductivity is a very active field in physics and researchers hope to find superconductors that can work in room tempera-ture.

The phases in these two examples are called quantum phases, as they find their description in quantum theories that describe things on a very small scale. Now that we have seen an example of several types of phases in condensed matter, I will explain the meaning of the next word that was in my description: topology.

(7)

Topology

Topology is a field of mathematics that is concerned with the shape of things. A well known joke about topology is that a topologist2 cannot distinguish between a coffee mug with one ear and a doughnut. The explanation behind this is that to the topologist all three dimensional objects that have one hole are identical3_{. This persiflage of the viewpoint of a topologist signals}

that all three dimensional objects with one hole are topologically the same and we say that all these objects belong to the same topological phase. Two different topological phases in this context then correspond to two collections of shapes that each have a different number of holes.

One can extend this concept of topological classes to physical systems, in which physical models are the same if they can be deformed into each other without changing ‘the number of holes’. Talking about the ‘number of holes’ of a theory is not really correct, but it shows that the mathematical concept of topology is used in physics to characterise different phases, which is what condensed matter is all about. It turns out that these ‘topological phases’ have quite remarkable properties, such as hosting very exotic particles that might be used for making a quantum computer.

Thesis work

To end this layman’s introduction, I will give a short outline of the work I have done.

For the first part of my work I have been reading papers and textbooks to gain insight into the field of topological phases. This was not only in the beginning, because I also kept track of new literature4 _{and needed to learn}

new things as my research progressed.

The second part of my work consisted of understanding a certain topo-logical model in two dimensions5 _{and trying to reproduce the results that}

several authors have published about the properties of this model. See Fig-ure [3.1] on page [12] for the lattice. For reproducing these results, I have used the mathematical programming software Matlab. First I had to teach myself how to use the software, after which I was able to write computer

2_{A topologist is a mathematician that works in the field of topology.}

3

More precisely, to a topologist, all objects that can be deformed into each other without cutting or gluing parts together are in the same topological class.

4

The physics community has a freely accesible website where nearly every new physics paper is published before it appears in a journal. It is called the Arxiv (pronounced as the english word ‘archive’). Although papers on the Arxiv are not peer reviewed, it is a valuable resource for each physicist to keep track of the new publications in his field. It can be found on arxiv.org.

5

In condensed matter physics, models are investigated that live in either 1, 2 or 3 dimensions. A one dimensional model is defined on a line, a two dimensional model on a plane and a three dimensional model is defined in the three dimensional space that we seem to live in.

(8)

scripts that simulate the physical model of interest. With these scripts I ob-tained full control over the parameters6 of the model and by choosing those parameters the same as in the literature, I could reproduce the published results. This reproduction was important, since it confirmed that I really understood the model and it also helped me to improve my proficiency in Matlab, which I needed for the next part of my research.

The third part of my work concerned another model in three dimensions, that is like the previous model in that it can be in a topological phase. See Figure [5.1] on page [38] and Figure [6.3] on page [59] for figures showing the lattice. Gaining insight into the properties of the model and the way it is formulated (the lqmathematical’ language that it is written in), allowed me to translate it to a Matlab script and eventually to let me compute the basic properties of the model that were already published.

Taking this model further, I was able to write a more general code in Matlab, with which I had more control over the model than was reported in the literature. I then calculated some new properties of the model and that was it.

The last part of the work was writing it all down in this thesis, which was a substantial undertaking.

I would like to make clear that the parts I just described did not at all follow linearly. The actual day-to-day experience of doing research in physics is that most of these processes are happening at the same time: while I was writing parts of my thesis, I often had to revisit the literature and also the writing improved my knowledge, which affected the Matlab code I was writing. Also there were quite a few moments at which not much seemed to be happening.

If this introduction has peaked your interest, there are many popular blogs and articles online that can shed light on the fantastic domain of condensed matter physics. And if you have a background in physics, I of course heartily invite you to read this thesis.

6

A parameter is characteristic of a system that takes certain values. If you listen to music, the volume is a parameter that indicates the loudness of the sound waves.

(9)

Dankwoord

Graag wil ik Kareljan Schoutens bedanken voor zijn begeleiding en de er-varing die ik via hem heb opgedaan met het onderzoek in de theoretische natuurkunde. Een onderdeel hierin was het mede faciliteren van het bezoek aan een workshop in Dublin. Echter, vooral het feit dat hij me zelf heeft laten zoeken naar een interessant onderwerp en ik er zodoende mijn eigen onderzoeksproject van kon maken, hebben hier sterk aan bijgedragen.

Tijdens mijn onderzoek heb ik af en toe contact gehad met Ville Lahti-nen, welke ik tevens wil bedanken voor de uitleg van zijn werk.

Jan van Maarseveen, zonder wiens moleculen-bouwdoos ik niet ver zou zijn gekomen met het begrijpen en opdelen van het diamantrooster, wil ik van harte bedanken, evenals Adri Bon voor de goede ict-dienstverlening waardoor ik onder andere met Matlab aan de slag kon.

Als laatste wil ik mijn liefste Rosa danken voor haar steun en de kritische spiegel die ze me soms heeft voorgehouden.

(10)

(11)

Introduction

1.1 Topology and condensed matter

In mathematics, the domain of topology is concerned with the properties of a space that stay the same under continuous transformations, where the operations of ’cutting’ and ’gluing’ spaces are discontinuous, so not allowed for the condition of continuity.

To elucidate this abstract definition, consider a coffee mug with one ear. If we deform the mug without cutting it or pasting parts together, we can change it into a doughnut (or punctured solid disk). And vice versa, a doughnut can be deformed to a coffee mug with one ear. The space that is considered in this example is the collection of 3D shapes with a condition that forbids cutting and gluing. The parts that are in the same continuously connected space or phase, are characterised by their number of holes.

The space that is considered by topologists can be linked to the space of quantum systems. For the space of all quantum systems, we can impose the condition that there must be a gap1 in the energy spectrum and look at which parts of the space of all quantum systems are continuously connected under this energy gap-condition. This gap condition yields an interesting distinction between quantum systems, whose phases are subsequently called topological.

Topological phases of matter became an active research field when Von Klitzing et al. reported in [6] on quantised Hall conductances on a two di-mensional plane . After the publication2, a great research effort in explaining the reported effect started and since then, the field has been continuously developing.

1_{A gap in the spectrum is the region above the minimal energy that has no energy}

modes. The value of this gap is the energy of the first excited state above the ground state.

2

The effect that Von Klitzing et al. reported on, was later named the integer quantum Hall effect or IQHE.

(14)

In the search for topological models, classifications of topological phases were made. Before topological phases of matter were first found, ordered phases often were due to a spontaneous breaking of symmetry and could generally be understood through the values of an order parameter. Only on continuously shifting the order parameter, the system’s phase could be changed. However, with topological models no order parameter can be de-fined and the condensed matter community had to come up with something else.

It turned out that the different topological phases can be classified by regarding certain discrete symmetries of the system. In a paper by Altland and Zirnbauer, [7], a total of ten symmetry classes were distinguished for random matrices, which were linked to Hamiltonians of gapped free fermion models in [4, 5]. The relevant symmetries used in this classification are time-reversal, particle-hole, and chiral symmetries.

Some topological phases were reported to be able host non-Abelian anyons which, upon actual physical realisation, might be of use in the error free performing of quantum computation [8].

In this thesis we investigate models in 2 and 3 spatial dimensions3. The two dimensional model has a non-Abelian phase that can be linked to the Ising anyon model and the three dimensionsal model, while being in a topo-logical phase, cannot support anyons.

1.2 Structure of this thesis

In Chapter [2], we present a brief outline of anyons and braiding. Also, in this chapter, the discrete symmetries are introduced that are used to fully classify topological phases of free fermion models.

In Chapter [3] the Kitaev model on the diamond lattice is presented, which was invented by Alexei Kitaev in [1]. The model has a non-Abelian topological phase with ν = ±1 in symmetry class D and is presented as a lattice model with spin-1/2 particles on each site that couple along three directions on the honeycomb lattice. We will show how this spin model can be re-expressed as a Majorana hopping model and a lattice gauge theory interpretation is presented.

In Chapter [4] the Kitaev model on the honeycomb lattice is presented with a unit cell that can host any amount of vortices on the lattice plaque-ttes. For several configurations of these vortices, the low-energy properties are calculated and results are obtained that exactly agree with those in the literature ([2]). An identification of these vortices as Ising anyons is made and from the data this claim is made feasible.

3_{This is a condensed matter thesis, so when me mention a model to be 3D we mean}

that it has three spatial dimensions, which in high energy physics would be written as 3+1 D

(15)

In Chapter [5], we look at a different topological model, namely the Kitaev model on the diamond lattice, as presented by Shinsei Ryu in [3]. This model has not 2 but 4 degrees of freedom per lattice site. Following the literature, it is transformed into a Majorana hopping problem with static Z2

gauge field on the diamond lattice bonds. The model shows a topological phase with winding number ν = 1 in symmetry class DIII. The gap opening of the model is presented.

In Chapter [6] we present out own work on the Ryu-Kitaev diamond model. We present the explicit Hamiltonian for first a 4-site unit cell and then, after showing that it gives identical results to Ryu’s, we increase its size to hold 4M N Z sites per unit cell, where M , N and Z can be any positive integer. Finally, we present the effect that introducing vortices on the lattice has and present what happens if we cut one of the periodic boundaries open, resulting in a system with two boundaries.

(16)

(17)

2

Preliminary pieces

In this first part we will give an overview of the topics related to topological phases that we will need in later chapters. We will start with a few words on anyons and braiding. Then we go on to explain how we can classify topological models by the use of discrete symmetries.

2.1 Anyons and braiding

In the field of topological quantum computation, the notion of anyons and the braiding of their spacetime paths is central. Although in this thesis we will not address topological quantum computation explicitly, the approach we have used is based at least in part on the ideas of quantum computation. We will see in a later chapter that the Kitaev honeycomb model can be identified as an Ising anyon model when choosing suitable parameter values. While the Kitaev model on the diamond lattice, which we discuss in the last two chapters, cannot support anyons, we were inspired to look for flux lines or loops that can be braided around each other, providing some parallel with the anyons on the honeycomb model.

Presently, we give a flavour of anyons and their braiding, for which we have used [9, chapter 4]. The notation we use agrees with this reference.

In arbitrary dimensions, if we have a two-particle state of indistinguish-able bosons or fermions 1 and 2, |ψ1(r)ψ2(r0)i, we can interchange them:

|ψ1(r)ψ2(r0)i = ± |ψ2(r)ψ1(r0)i , (2.1)

where + refers to the particles being bosons and − refers to them being fermions. In two dimensions, however, there is another possibility, namely for anyonic particles to obtain an arbitrary complex phase upon interchang-ing:

(18)

The term ’anyon’ was first used by Frank Wilczek in [10], where he indicated particles that were not bosonic of fermionic, but could support any phase, hence the name.

In physical situations, anyons can be realised as quasiparticles that are interesting precisely because of their non-trivial exchange statistics. When we want to address the statistics of anyons, we typically look at interchanging particles, which can also be represented as the braiding of their world lines in the 2+1 dimensional spacetime, hence the word ‘braiding’.

In three dimensional space, point excitations have no topological prop-erties, since the path of twice exchanging two point excitations is reducible to an identity operation and is therefore topologically trivial. In [11], the braiding statistics of 3D excitations is treated. In three dimensions, one can look at the statistics of point-like excitations and line/loop excitations. We can associate a flux to a loop and a charge to a point excitation and following the Aharonov-Bohm effect [12], braiding the charge through the loop yields a statistical phase

θ = q · φ. (2.3)

A more involved braiding operation is that of two loops that braid around each other. For this we take two vortex loops α, β with fluxes φα, φβ

carrying charges qα, qβ respectively. We thread one loop through the other

and back to its initial position, which is equivalent to first threading loop α through loop β and then threading loop β through α. To simplify, we turn-wise contract the loops. Following the Aharanov-Bohm effect that we just mentioned, we can write the statistical phase as

θαβ = qα· φβ+ qβ· φα. (2.4)

The authors of [11] note that statistical phases obtained from braiding can be used to identify the topological phase of a system. They argue that θ and θαβ are not sufficient if we want to fully distinguish different 3D

topological phases by braiding. To fully determine the topological phase, it was argued that the phases obtained from braiding two loops around each other that are threaded by a third loop provide a sufficient characterisation to make the distinction.

Thus, when using braiding to determine the topological phase of a 3D model the braiding properties of threading two loops around each other, while being threaded themselves by a third loop is needed. While in this thesis we will not perform these braiding operations, they might be done with the explicit model that is presented in the last chapter.

(19)

2.2 Symmetries

In 1931, Eugene Wigner proved a theorem that specifies the form of any symmetry operation. The theorem, called Wigner’s theorem, postulates that the general form of an arbitrary symmetry operation S acting on a state x of Hilbert space H can be written as

Sx = φ(x)U x (2.5)

where U is either a unitary or anti-unitary matrix and φ(x) is a phase with modulus one. Since this theorem is valid for all symmetries, we know that any operator generating symmetries is either unitary or anti-unitary1_,

which gives us clear instructions to find the explicit form of any symmetry operator.

For materials that show topological properties, it turns out that discrete symmetries are important [7]. The two symmetries that are important in this respect, are particle-hole (PH) and time-reversal (TR). PH and TR operators are denominated by P and T respectively. We call an operator A symmetric under PH (TR) if it commutes with P (T):

[A, P ] = 0 or A = P AP−1 (particle-hole symmetry) (2.6)

[A, T ] = 0 or A = T AT−1 (time-reversal symmetry) (2.7) A third symmetry is chiral symmetry, which is automatically present when both TR and PH symmetries are present. There is also a possibility to have a chiral symmetry when TRS and PHS are not on their own obeyed.

2.2.1 Time-reversal symmetry

A system that is time-reversal invariant, is invariant under the transforma-tion

T : t → −t (2.8)

Much of the physics we know explicitly breaks time-reversal. Think, for instance, of a charged particle in a magnetic field. If we were to send t to −t, the particle would not take the same path back, but would follow another

1_{Anti-unitary operators A fullfill A(φα + ψβ) = (Aφ)α}?_{+ (Aψ)β}?_{, while they are a}

bijection on the full Hilbert space and they preserve the norm of the states. The most common example of an anti-unitary operator is the complex conjugate operator K.

(20)

trajectory due to the Lorentz force having changed direction. When consid-ering topological phases, it however turns out that TR symmetry appears quite often.

To find out what a time-reversal operator looks like, we focus on the action of T on the commutator of ˆx and ˆp, the position and momentum operators. Position is obviously unaffected by time-reversal, but momentum is, since p contains a velocity, ~p = m~v = m d~x/dt, and sending t to −t will flip the sign of p. Hence, we can write,

T i~T−1= T [ˆx, ˆp]T−1 = −[ˆx, ˆp]T T−1 = −i~. (2.9) This shows that the TRS operator acts as T iT−1 = −i, so T must be proportional to the complex conjugation operator K, which is an anti-unitary operator.

Discrete symmetries for spin-1₂

In the introduction we gave some definitions of discrete symmetry operators for time-reversal (TR) and particle-hole (PH) operations. For spin-1₂, we are now in a position to give an explicit expression of TRS operator.

Since spin is a form of momentum, it must change sign upon TR:

T ST−1= −S, (2.10)

Upon time reversal spin changes direction: S → −S. We can interpret this as a rotation by π around an arbitrary axis. The convention is to pick the y-axis [13]. This fixes the proportionality that we needed to find the form of a TR operator and we can write it down as

T = e−iπSy_K. _(2.11)

We can expand the exponent and simplify, remembering that (σ2)2 = 1,

e−iπSy _{= e}−iπ₂σ2 _{= cos(}π

2σ2) − i sin( π 2σ2). (2.12) cos(π 2σ2) = 1 − (π₂σ2)2 2! + (π₂σ2)4 4! − ... = cos( π 2) = 0, sin(π 2σ2) = π 2σ2− (π₂σ2)3 3! + ... = σ2sin( π 2) = σ2,

which ultimately yields a representation for the two-dimensional time-reversal operator:

(21)

2.2.2 Particle-hole symmetry

A particle-hole symmetry exists when a system is indifferent to whether it is described in terms of particles or in terms of holes (absence of particles). In the language of second quantisation the effect of a system that has particle-hole symmetry is easiest to see. If one would replace the creation operators of a second-quantised Hamiltonian into annihilation operators, the Hamil-tonian must remain the same. A direct consequence of a system with PHS is that its energy is symmetric around E = 0: counting minus the energy of all holes must yield the same as counting the energy of all particles. Act-ing with a particle-hole conjugation operator on a Hamiltonian must yield minus the same Hamiltonian.

To see that a PHS operator is anti-unitary, we follow the same logic as we did for the TRS operator in Section [2.2.1]. Instead of starting with {ˆx, ˆp}, we start with {t, E} = i~, the anti-commutator of time with energy. We take P to be the particle-hole symmetry operator. Acting with P leaves t invariant, but, as we mentioned, it tranforms E to −E. This yields

P i~P−1 = −i~, (2.14)

and so the anti-unitarity of P is shown.

For Hamiltonians that are of Bogoliubov- de Gennes (BdG) type, we can find an eplicit form for P . To do this, we take a BdG Hamiltonian

H =Ψ† ΨH Ψ

†

Ψ† !

, (2.15)

where Ψ†and Ψ are the fermion creation and annihilation operators respec-tively and

H = H0 ∆ ∆? −H0

!

(2.16)

is the BdG Hamiltonian. Acting with σx on the state Ψ

†

Ψ† !

exchanges Ψ with Ψ†, which is exchanging particles with holes. And hence, for the BdG Hamiltonian, we can write

P = iσxK. (2.17)

2.2.3 Symmetry classes

With these symmetry operators, we can give the classification of topolog-ical phases for free fermions. Each class has only topologtopolog-ical phases in a restricted number of dimensions. The phases that each symmetry class can have, are indicated by the topological invariant that takes a certain value.

(22)

Figure 2.1: The tenfold way of classifying topological phases in n mod 8 dimensions as discussed in [4, 5]. The two entries encircled in red are the two phases that we will encounter in this thesis. Symmetry class D in 2D is the topological phase of the Kitaev honeycomb model and we will encounter the DIII class in the 3D Ryu-Kitaev model on the diamond lattice.

In most cases this is an integer number that can be either 1 or 0 (Z2)2, and

in some cases it takes values in all integers of Z.

In [7] an exhaustive classification was presented of ten different symmetry classes for random matrices. These ten random matrix classes were later identified to correspond to ten symmetry classes for topological phases of free fermion theories in [4, 5], which led to a much-used framework in the search for new topological phases.

For the 2D Kitaev model on the honeycomb lattice that we introduce in the next chapter, we will encounter symmetry class D that has only particle-hole symmetry. in 2D, class D can host any integer topological invariant. The 3D Ryu-Kitaev model on the diamond lattice has both TRS and PHS, and is in symmetry class DIII. Just like class D in 2D, class DIII in 3D can host any integer topological invariant.

2

(23)

3

Kitaev honeycomb model in

two dimensions

The Kitaev honeycomb model was introduced by Alexei Kitaev in his sem-inal paper [1]. It is a spin-1₂ system on a honeycomb lattice, with nearest-neighbour interactions that are of XX, YY and ZZ type, depending on the direction of the link joining the nearest neighbours.

The model makes it possible to put vortices on the lattice, which subse-quently can be used for braiding in their 2+1 dimensional spacetime. Braid-ing these vortices is interestBraid-ing, because we will later see that the model can enter a non-Abelian phase, which hosts non-Abelian anyons that can be used for topological quantum computation.

In this chapter we will explain the model. First we introduce the spin-1₂ Hamiltonian and discuss some of its features. Next, we identify the spin operators as products of Majorana fermions which simplify the model to a free fermion model with a fixed Z2 gauge field.

The content of this chapter is based largely on references [1] and [2]. [9] Was also used, which relies in part on the work done in the earlier mentioned references. The notation we use, agrees with [2] and we note that [9] has different notation.

3.1 The Kitaev model on the honeycomb lattice

The Kitaev model consist of spin-1₂ particles residing on vertices of the honeycomb lattice. The interaction between nearest neighbours is of XX, YY and ZZ type for link in directions x, y and z, respectively. The lattice can be presented as the union of two triangular sublattices ΛAand ΛB. We

will call ΛAthe odd sublattice and represent it pictorially with filled circles.

ΛB is the even sublattice and will be drawn as empty circles.

From each site there are three directions to a neighbouring site, x, y and z. In Figure [3.1], the picture is made clear. Each of the three nearest

(24)

z x y σz σz σxσ x σy σy 1 4 2 3 6 5 p ΛA ΛB 1 2 3 σy σx σz (a) (b) (c) (d)

Figure 3.1: The bipartite honeycomb lattice with spins on every vertex. (a) The link directions x, y and z assigned. (b) The nearest neighbour and next-nearest neighbour terms. (c) The directionality of the next-nearest-neighbour hopping. (d) A plaquette p is defined through the loop around a hexagon. The lattice is numbered such that each A(B)-sublattice site is drawn as a filled (empty) dot.

neighbour terms has its own coupling stength Jα, α = x, y, z. The model

also has three-spin terms with coupling strength K, given by the product of spin matrices of three consecutive sites.

The Kitaev honeycomb Hamiltonian with nearest and next-nearest neigh-bour terms is given by [2]:

H = −Jx X x links σx_iσ_jx−Jy X y links σ_iyσy_j−Jz X z links σ_izσ_jz−K X (i,j,k) σx_iσ_jyσz_k. (3.1) The couplings Jx, Jy , Jz and K are non-negative. We can tune these

couplings separately, making the model explicitly anisotropic. Three-spin terms, involving next-nearest neighbours i and k, connected through site j, can be seen as an effective magnetic field with field strength K [9]. We will see that this term explicitly breaks time-reversal symmetry (TRS)1.

To be precise about the three-spin term, we remark that the summation of the i, j, k and the permutation of the indices x, y, z are site-specific. In Figure [3.1] we have depicted the three-spin term concerning sites 1, 2 and 3. The Hamiltonian term for this three-spin case will be2 σ₁yσx₂σz₃.

1

From now on we will only write TRS to denominate time-reversal symmetry.

2_{In the general case, the term σ}α iσ

β jσ

γ

(25)

The shortest loop that can be made on the honeycomb lattice is along a hexagon, see Figure [3.1](d). Such a hexagonal loop circumferes a plaquette that we will indicate with integer index p. We introduce the plaquette operator or Wilson loop operator Wp:

Wp = σz1σx2σ y 3σ z 4σ5xσ y 6. (3.2)

Wp squares to the identity3, which ensures its eigenvalues to be wp = ±1.

It can be shown that the operators Wp commute with the Hamiltonian,

[Wp, H] = 0, (3.3)

which makes the plaquette operator a local symmetry. To name some other properties, the plaquette operators commute with each other and their prod-uct over the whole lattice is 1:

[Wp, Wp0] = 0,

Y

p

Wp = 1. (3.4)

3.2 Majorana fermionization

Now that we have provided the initial set-up of the Kitaev model, we will introduce a transformation from spins to Majorana fermions. This transfor-mation will turn the Kitaev model into a Majorana fermion hopping problem that is easier to solve. The treatment of the Majorana substitution that is presented here follows [1] and [2].

For the spin on a given site i, we can introduce two fermionic modes a1

and a2. From these fermion modes we construct four modes, bαi (α = x, y, z)

and ci by taking their real and complex values4:

bx_i = i(a†_1,i− a_1,i), by_i = a_2,i+ a†_2,i, ci= a1,i+ a † 1,i, bzi = i(a † 2,i− a2,i). (3.5)

These operators obey the Majorana condition

(bα_i)†= bα_i, c†_i = ci, (3.6)

as well as the anti-commutation relation

{bα_i, bβ_j} = 2δα,βδi,j, {c_i, bα_j} = 0. (3.7)

associated with links from i (k) to j. β is the index of the link connecting j with the link pointing away from the three-spin path.

3

This is true because spins on different sites commute and on the same site {σα, σβ} =

2δα,β 4

The notation is the same as in [1] and [2]. The reason of naming the modes bα and c

(26)

The Majorana condition and fermion anti-commutation relation mean that bα_i and ci are Majorana fermion operators for any lattice site i.

We can write the spin operators of the Hamiltonian in terms of these Majorana operators [1]:

˜

σα_i = ibα_ici, (3.8)

where the subscript i is the lattice index and α is x, y or z.

This representation in terms of Majorana fermions obeys the Pauli ma-trix algebra of the spin matrices. We would like to write the Kitaev Hamil-tonian in terms of Majoranas, but we must first address the fact that we have introduced too many degrees of freedom in the identification of Equa-tion (3.8).

3.2.1 Degrees of freedom

With the introduction of two fermions on each site, we obtain a local 4-dimensional Hilbert space per spin (both fermion modes empty, either one occupied, or both occupied), but since a spin-1₂ particle only has a two-dimensional space (from its two possible projections), we must find a re-striction to the number of fermion modes we allow to represent the model. To clarify this, we say that we are in the extended or full space ˜L when addressing the full, 2-fermion case. This extended space contains extra, un-physical, states and the original physical states that are described in the intial spin context. In order to restrict the number of degrees of freedom, we will project out states ˜|Ψi from the extended space ˜L onto the physical subspace L ⊂ ˜L that contains the physical states that are required for a truthful representation of spin operators as Majorana fermion products.

We remind the reader of the mathematical identity

σxσyσz = i. (3.9) Using the Majorana representation, this becomes

e

σxσ_ey_eσz= −ibxc byc bzc = ibxbybzc = iD, (3.10) where we identified D = bxbybzc. We can add a site-index i to D to make it act on a specific site: Di.

We find that in order to satisfy Equation (3.9), our physical states |Ψi must obey

Di|Ψi = |Ψi ∀ |Ψi ∈ L. (3.11)

So the Di operator acts as the identity on the physical subspace.

To see which fermion states we should choose to represent a single spin, we address the two complex fermion modes we started out with, rewritten in terms of Majorana modes:

(27)

a1,i= 1 2(ci+ ib x i), a2,i = 1 2(b y i + ib z i). (3.12)

We rewrite the physical requirement, Equation (3.11), as

D |Ψi = (1 − 2a†₁a1)(1 − 2a†2a2) |Ψi . (3.13)

Since we want D to act as the identity operator, we see that this is only possible for the states |0i and a†₁a†₂|0i, which consequentially lie in the physical space L. The other two states, a†₁|0i and a†₂|0i, lie only in the extended, unphysical space ˜L. Note that for these last states, acting with D yields a minus sign. Using D we can remove two of the four states, reducing the number of states with a factor of 2. The physical subspace L is exactly half of the extended space ˜L.

We can now construct a projector P to project out the non-physical states: PL= Y j 1 + Dj 2 , (3.14)

where the product runs over all site indices j. Acting with P on an unphys-ical state yields zero.

We have introduced a way to project out the unphysical states that were introduced by the initial fermionization and can conclude that the representation of Equation (3.8) is truthful for the case where we restrict ourselves to the physical subspace L. This means that the Hamiltonian can be written in terms of Majorana fermions, without changing the physics.

3.2.2 Final Hamiltonian form

Using Majorana operators, we can express a two-spin Hamiltonian term as follows:

σ_ixσ_jx= ibx_ici ibxjcj

= −i(ibx_ibx_j)cicj

= −iˆux_ijcicj,

(3.15)

where we have identified the link operator ˆuα_ij = ibα_ibα_j (α = x, y, z). This link operator will pop up quite often. For now, we remark that uα_jk = −uα_kj.

(28)

Figure 3.2: the positive direction of the link operators is indicated. The operator uij is +1 when j is on a black (odd) site and i on a white (even)

site. With this one can assign a positive-direction arrow as indicated. The operator uij can be intepreted as “ the directional link from site j to site i”

A three-spin term yields, remembering (bα)2 = c2 = 1, σ₁yσ₂xσ₃z= iby₁c1 ibx2c2 ibz3c3 = i3by₁bx₂bz₃c2c1c3 = i3by₁(by₂)2bx₂(bz₂)2bz₃c2c1c3 = −i ûy₁₂uˆz₂₃bx₂by₂bz₂c2c1c3 = −i ûy₁₂uˆz₂₃D2c1c3 = −i ûy₁₂uˆz₂₃c1c3, (3.16)

where from unity we introduced the squares and identified the operator D. Since we are only interested in states in the physical space, the D operator acts as unity and can be taken out.

We will now give the full Hamiltonian in terms of Majorana fermions:

H = i 2   X (i,j) Jαijuˆ αij ij cicj + K X (i,k,j) ˆ uαik ik uˆ αkj kj cicj  . (3.17)

Here, the first sum is over nearest neighbours i and j and the second is over next-nearest neighbours i and k connected through the site at j. αij is x, y

or z, depending on the directionality of the link between sites i and j. Because of the antisymmetry of the link operators ˆuij, we must fix an

orientation. We choose the orientation such that hopping from an odd to an even sublattice site is in the positive direction and hopping the other way yields a minus sign5. This is visualized in Figure [3.2]. Note that uij is the

link from j to i and not the other way round.

Since the direction αij of the link between sites i and j is implicit from

the sites themselves, we will from now on refrain from writing the directions α in the link operators uij. We can further simplify the Hamiltonian, by

noting that the ˆuij link operators square to 1. The link operators commute

with H, so we can replace them with the set of values {uij} with each

uij = ±1. Note that the factor 1₂ with respect to Equation (3.1) is due to

counting the hoppings (i, j) and (j, i) separately.

(29)

by by bz bz c bx bx

Figure 3.3: Identification of spins as Majorana particles. The bα are con-nected with each other and the c’s remain on the original lattice site. The Kitaev model is expressed as a hopping model between the Majoranas c on the honeycomb lattice sites.

Effecting the above remarks, the Hamiltonian takes its final form

H = i 2   X (i,j) Jαuijcicj+ K X (i,k,j) uikukjcicj  , (3.18)

where the first terms sums over nearest neighbors and the second term over next-nearest neighbors i and j connected through k.

It now remains to make a choice for the configuration of the set {uij} of

all link operators, effectively choosing a flux configuration {wp}, and then

to compute the spectrum. However, in the following we will argue that the link operators uij are outside the physical space L and present a Z2 gauge

field interpretation.

3.3 Emerging lattice gauge theory

One can verify that the Hamiltonian commutes with the D projection op-erators, [H, Di] = 0, which means that we can compute the states of the

model, while at the same time restricting them to the physical space. Aside from the projection operators, one can verify that the link opera-tors are also local symmetries. This is shown by proving that [H, ˆuij] = 0:

ˆ

(30)

operator through the Hamiltonian yields a plus sign, resulting in the zero-commutator.

From this one might conclude that it then is straightforward to fix the oriented ˆuij’s to ±1 on all links and obtain the Hamiltonian that is quadratic

in ˆuij’s, as done in the previous section. However, when we remind ourselves

of the requirement of Equation (3.11), we must be a little more careful. Upon computing, one will find that the link operators ˆuij do not all

commute with Di6. In fact,

{ˆuij, Di} = 0. (3.19)

The anti-commutativity with D indicates that sectors labeled by their con-figuration of uij operators do not per se lie in the physical subspace L.

Other than focussing on separate link operators, we note that the loop operator Wp, which in the Majorana fermionization translates to a

prod-uct of link operators going around a hexagon, does belong to the physical subspace L: Wp = Y i,j∈p ˆ uij, (i ∈ ΛB, j ∈ ΛA)7, [Wp, Di] = 0. (3.20)

We can now identify the link operators ˆuij as local Z2 gauge fields that

live on the bonds between sites [1]. The local gauge transformations that act on them are the Di operators.

Di flips the sign of the three link operators that are connected to site i.

This does not change the plaquette values, since it adds a factor of (−1)2 to the loop. Going around a hexagonal loop on the honeycomb lattice always involves two links connected to the same site, so the plaquette operators Wp are the gauge invariant objects and can be identified as Wilson loop

operators. The different physical sectors of the Kitaev model are equivalent to the configuration of plaquettes on the lattice.

Plaquettes with an eigenvalue wp = −1 can be identified as having a

π-flux vortex living on their surface [9]. The eigenvalues uij = −1

(direc-tionalised) form an unphysical string from one vortex to another or form an unphysical loop. We call these strings unphysical, since their link operators uij = −1 do not commute with Di. To make sure that the description of our

model stays physical, we use the projector PL to project out the unphysical

states to stay safely in the physical world:

˜

|Ψi ∈ ˜L → PL|Ψi = |Ψi ∈ L,˜ PL= N Y i=1 I + Di 2 , (3.21) where the projector was already introduced in the previous section.

6

Still [ˆuij, Dk] = 0 for i, j 6= k due to the local nature of both operators.

(31)

3.3.1 Consideration of the unphysicality of gauge sectors

In the next chapter we will treat the properties of the Kitaev honeycomb model for different vortex sectors w that are identified by their configuration of vortex plaquettes, w = {wp}. In our analysis of these vortex sectors, we

will pick a single gauge choice to produce the desired vortex sector. However, as we have argued above, these link operators do not lie in the physical subspace.

It turns out that the energy spectrum of the model is invariant under the choice of Z2 gauge fields, as long as this choice corresponds to the same

vortex sector. This means that if we want to address the spectrum of the model in a certain vortex configuration, it suffices to use only one Z2 gauge

choice.

However, when one is interested in the wave functions belonging to the energy spectrum, the unphysicality of the Z2 gauge fields does need

con-sideration. To obtain the physical states |Ψwi of the system with a certain

vortex configuration w = {wp}, we perform an equal-weight sum over all Z2

gauge field configuration states |Ψui that yield the vortex configuration w.

The equal weighted sum of gauge sector states then does lie in the physical subspace, since it commutes with the projector D.

3.4 Concluding remarks

In this chapter we have introduced Kitaev’s honeycomb model. We have seen that the spin connections in three directions could be reformulated to a Majorana hopping problem. Important in the treatment was the fact that upon Majoranising the model, we had to prevent the introduction of addi-tional degrees of freedom by explicitly projecting onto the physical states. In the last section, we provided a gauge theory perspective. We saw that the plaquette operators act gauge invariantly, since at the site they act on, they flip the sign of all incoming connections, while leaving the netto vortex configuration intact.

In the next chapter, we will use what we learned here and compute the low-energy spectrum of the model and connect it to the Ising anyon model.

(32)

(33)

4

Spectrum of the Kitaev

honeycomb model

In the previous chapter we presented the Kitaev model on the honeycomb lattice. Now that we have reached a quadratic Hamiltonian in terms of Ma-jorana fermions ci, we can solve the spectrum. In this chapter we compute

the low-energy spectrum of the Kitaev honeycomb model, which was done by Lahtinen et al in [2].

We will follow the approach of first diagonalising the Hamiltonian, using Fourier transforms, and then choosing a vortex sector to obtain the low-energy spectrum of the Hamiltonian. After a short outline we will present the Hamiltonian with site-specific hopping terms. We then move on to treating the fluxless lattice and focus on its low-energy spectrum and the tolopogical phases that occur. In the last section, we look at the influence of flux plaquettes on the spectrum and finish with describing their connection to the Ising anyon model.

Additional to [2], the original work of Kitaev, [1], was also used. While the set-up and treatment is taken from these sources, the actual plots presented are ours. Since we later extend the current model to a three-dimensional model on the diamond lattice, we used the 2D honeycomb model to become familiar with the procedures and computations. Especially the hand work for 2D in Matlab was a good preparation for the Matlab approach in 3D.

4.1 Explicit Hamiltonian

To simplify the task ahead, we first deform the honeycomb lattice to a square one. This is done by setting the links in the z-direction to zero length, which makes a single site on the square lattice correspond with two sites, one on ΛA and the other on ΛB of the honeycomb lattice. This is

(34)

sy sx sy sx

sx

sy

Figure 4.1: The honeycomb lattice can be deformed by shrinking the z-bonds to zero length, rotating it and stretching the bonds to squares.

sx = (1, 0) and sy = (0, 1). In order to keep describing the same system, we

introduce an index λ ∈ {1, 2} to distinguish between ΛA (λ = 1) and ΛB

(λ = 2).

For this new square lattice, we take a unit cell composed of M N lattice-sites, where M is the number of sites in the sx-direction and N the number

of sites in the sy-direction as indicated in Figure [4.1]. A site in this unit

cell is reached by the vector k = (m, n); m ∈ {1, 2, ..., M }, n ∈ {1, 2, ..., N }. To be explicit, we rewrite the lattice index i, as done in [2], to i = (r, k, λ), where λ is the sublattice index that takes values in {1, 2}, k the unit cell vector that defines the location (m, n) on the unit cell (M, N ) and r the vector that indicates on which unit cell the lattice site is located. With this identification, the Hamiltonian is written as

H = i 4 X r,v (M,N ) X k,l 2 X λ,µ=1 Aλkµl(v)cλk(r)cµl(r + v), (4.1)

where A contains the nearest neighbour and next-nearest neighbour hopping terms, together with the link operator eigenvalues uij.

The Hamiltonian now sums over unit cells at r and v, which we exploit by taking the Fourier transform with respect to the unit cell vectors r.

The Fourier transform of the Majorana operators is given by [2]:

cλk(r) = (M,N )

Z

B.Z.

d2p eip.rcλk(p). (4.2)

After the Fourier transform, temporarily using indices a = (λ, k) and b = (µ, l) and suppresing the link operators, we encounter terms such as

(35)

X r,v ca(r) cb(r + v) = 2 C X r,v Z B.Z. d2p d2p0 eip·r eip0·(r+v)cacb = 2 C Z B.Z. d2p d2p0X r ei(p+p0)·reip0·v cacb = 2 Z B.Z. d2p d2p0 δ(p − p) eip0·vcacb (4.3) = 2X v Z B.Z. d2p eip·vcacb.

We see that we are only left with a summation over unit cell vectors v. Our hopping matrix only contains nearest- and next-nearest terms, so for unit cells of dimension greater than 1, we can only reach the nearest unit cells located at (M, 0), (0, N ), (M, N ) and their negative counterparts. For the rest of this thesis we, however, set the number of unit cells to one, C = 1, so that these nearest neigbouring unit cells are actually the same unit cell as the original one. As a direct effect, we will only encounter momentum-bearing exponents for hopping entries that go across the periodic boundary of the single unit cell.

Earlier, we remarked that shrinking the z-bonds to zero to obtain a square lattice with lattice vectors nx = (1, 0) and ny = (0, 1) does not

affect the physics of the system. We now briefly point out why. As can be concluded from the above, the only point at which the geometry of the lattice comes into play, is in the product of the unit cell vectors v with the reciprocal lattice vectors p. If we consider the honeycomb lattice, we have the unit cell lattice vectors ˜v1 = M (1, 0)1 and ˜v2 = N₂(1,

√

3) for which our lattice is periodic. From these lattice vectors we obtain the reciprocal lattice vectors ˜p1 = 2π_M(1, 1/

√

3) and ˜p2 = 2π_N(0, 2/

√

3). Taking the inner product yields ˜v · ˜p = (2π, 2π). If we do this for the square lattice vectors vx= (M, 0), vy = (0, N ) and their reciprocal counterparts p1 = (2π_M, 0) and

p2 = (0,2π_N), we get exactly the same result and hence making the lattice

square is allowed.

The first Brillouin Zone (B.Z.) of a square lattice with unit cell M nx+

N ny = (M, N ) is a rectangle in reciprocal space (px, py), bounded by the

lines at px = −_Mπ , px = _Mπ, py = −_Nπ and _Nπ. Therefore the integral over

the B.Z. becomes [2] (M,N ) Z B.Z. d2p = r 2 C π/M Z −π/M dpx p2π/M π/N Z −π/N dpy p2π/N, (4.4) 1

(36)

where C denotes the number of unit cells on the lattice. We obtain H = 1 2 π/M Z −π/M dpx p2π/M π/N Z −π/N dpy p2π/N c1(p) c2(p) !† A11(p) A12(p) A21(p) A22(p) ! c1(p) c2(p) ! , (4.5) where c†_λ = (c_λ(1,1), c_λ(1,2), ..., c_{λ(M,N )})† and the matrices Aλµ are M N ×

M N matrices, containing the nearest neighbour (λ 6= µ) and next-nearest neighbour (λ = µ) hopping amplitudes.

In [2], the nonzero elements of the above Hamiltonian were given. We copy them out in full.

c†₁A12c2= 2i (M,N ) X k=(1,1) ( + uk,k Jz c † 1,kc2,k + uk,k−nx Jxe iδ(m−1)p·vx _c† 1,kc2,k−nx + uk,k−ny Jye iδ(n−1)p·vy _c† 1,kc2,k−ny), (4.6) ˆ y ˆ x c†₂A21c1= 2i (Lx,Ly) X k=(1,1) ( − uk,k Jz c†_2,kc1,k − u_k,k+n_x Jxe−iδ(m−M )p·vxc†_2,kc1,k+nx − u_k,k+n_y Jye−iδ(n−N )p·vy c † 2,kc1,k+ny). (4.7) −ˆx −ˆy The next-nearest neighbour hoppings on ΛA yield the terms

c†₁A11c1 = 2iK (M,N ) X k=(1,1) ( + uk_k,k+n_y e−iδ(n−N )p·vy _c† 1,kc1,k+ny − uk−nx k,k−nx+ny e iδ(m−1)p·vx_e−iδ(n−N )p·vy _c† 1,kc1,k−nx+ny − uk k,k+nx e −iδ(m−M )p·vx _c† 1,kc1,k+nx (4.8) + uk−ny k,k+nx−ny e −iδ(m−M )p·vx_eiδ(n−1)p·vy _c† 1,kc1,k+nx−ny + uk−nx k,k−nx e iδ(m−1)p·vx _c† 1,kc1,k−nx − uk−ny k,k−ny e iδ(n−1)p·vy _c† 1,kc1,k−ny),

(37)

and next-nearest neighbour hopping terms for ΛB are c†₂A22c2 = 2iK (M,N ) X k=(1,1) ( − uk+ny k,k+ny e −iδ(n−N )p·vy _c† 2,kc2,k+ny + uk+ny k,k−nx+ny e iδ(m−1)p·vx_e−iδ(n−N )p·vy _c† 2,kc2,k−nx+ny + uk+nx k,k+nx e −iδ(m−M )p·vx _c† 2,kc2,k+nx (4.9) − uk−ny k,k+nx−ny e −iδ(m−M )p·vx_eiδ(n−1)p·vy _c† 2,kc2,k+nx−ny − uk_k,k−n_x eiδ(m−1)p·vx _c† 2,kc2,k−nx + uk_k,k−n_y eiδ(n−1)p·vy _c† 2,kc2,k−ny).

Note that ul_k,m is a short hand notation for uk,lul,m. We have made use

of the vector k = (m, n) to run over the lattice sites and M and N are the dimensions of the unit cell.

For the pair of sites inside the same unit cell, we see that v = 0 and we get terms without momenta. When the connected lattice sites are on different unit cells, say one at r and the other at r0= r + avx+ bvy, the corresponding

terms obtain a momentum-bearing complex phase of eip·(avx+bvy)_{with a, b =}

±1. The delta-functions in Equations (4.6) to (4.9) select the instances where the boundary is crossed and for those cases give the expected phase. We see the hermiticity of H from A12= A

† 21.

For all the above terms, the signs of the uk,l have already taken into

account the directionality of the link. This means that we arrive at the fluxless sector if we choose uij = 1 ∀ i, j, without taking into account the

orientation of the links. From now on, when we talk about a link uij having

value −1, this means that the oriented link has a negative sign. The next-nearest neighbour terms have additionally been suited with signs that follow the sublattice directionality indicated in Figure [3.1](c).

The spectrum that we gain from this Hamiltonian is double [1], in the sense that for every E(q), there is an energy −E(q) with the same modulus. This automatically leads to the identification that the negative part of the spectrum forms the valence band and the positive energies represent the levels in the conduction band. In this light, we speak of the energy gap ∆0

which we define as the energy of the first excited state. When the system is gapless, it means that there are two E = 0 states, so that the bands touch and the gap is closed.

The model obeys particle-hole symmetry since the eigenvalues of the Hamiltonian that are mirrored around E = 0. The next-nearest neighbour hopping terms can be seen as introducing a local magnetic field on the lattice, which breaks time-reversal symmetry. With only particle-hole symmetry,

(38)

the Kitaev model is in symmetry class D, that has phases with integer topological invariants ν = Z.

In what follows we will look at different vortex configurations. We adopt the denomination of [9] and [2] and speak of n-vortex sectors, by which a configuration of n vortices on the otherwise vortex-free lattice is meant.

4.2 Vortex Sectors

In this section, we present several vortex sectors. We remind of the plaquette operator Wp =Qi,j∈∂puij. Since an eigenstate with wp = −1 corresponds

to a vortex on plaquette p, we can make vortices by putting specific links uij

to −uij, which allows us to put vortices in the center of specific plaquettes.

In [14], it is remarked that instead of locally flipping the sign of the local gauge field uij, one can also change the sign of the local coupling strength Jij

and the nnn coupling parameter K. This identification allows for a physical manipulation of the vortex sectors, since in a an experiment one might be able to control these local couplings, but not the gauge fields on the links.

We will start with discussing the vortex-free sector. Second, we will look at the full-vortex sector and after this, we will compute the low-energy spec-trum of several sparse vortex configurations. The results shown here are of our own computations. As was done in [2], we look at unit cells of dimen-sions M, N = (20, 20). These dimendimen-sions are large enough to clearly see the spectrum for the inifite lattice and choosing these dimensions also allows us to scan over a fine-meshed Brillouin zone, which results in sufficiently accurate results for moderate desktop computer runtime.

4.2.1 Vortex-free sector

The Kitaev model on the honeycomb lattice is in the ground state when there are no vortices present, i.e. wp = 0 for all p [1]. The original

argu-ment for this was presented by Lieb in [15]. When we consider the fluxless case, all uij’s are replaced with +1, reminding that we already took the

orientation of the links into account in the formulation of the Hamiltonian in Equations (4.6) to (4.9). This yields the needed value of +1 for every plaquette wp.

It was shown in [1] that the vortex-free sector can be in four different phases, Ax, Ay, Az and B. The three A phases are gapped and support

Abelian anyons and it was shown that they have Chern number ν = 0. Kitaev makes the link between the A phase and the toric code model, for which the species of anyons have Abelian braiding statistics2.

The system is in the Aαphase when one of the following three inequalities

(39)

Jy Jx+ Jz Jz Jx+ Jy Jx Jy+ Jz Ax Ay Az B Jx≈ Jy≈ Jz Jz= Jx+ Jy J x = J y + J z _Jy= Jx+ Jz

Figure 4.2: Cross section of the positive quadrant (Jα≥ 0), with Jx+ Jy+

Jz = 1, of the phase diagram of the Kitaev honeycomb model showing the

Aα and B phases. The B-phase, indicated with a gray surface, obeys all

three relations of Equation (4.10). At the boundaries of the B phase one of the inequalities is changed into an equality, which is consequently destroyed if we move further into the Aα phase.

is broken: |J_x| ≤ Jy +|J_z| , Jy ≤ |J_x| +|J_z| , |J_z| ≤ |J_x| + Jy . (4.10) When all three inequalities are obeyed, the system is in a non-Abelian phase, which we denote with B. In the B-phase with K = 0, there are two points in momentum space where the bands touch and have linear disper-sion3, which makes the B-phase gapless. These two Dirac points disappear when for instance |Jx|+

Jy

= |J_z|, which constitutes one of the three bound-aries of the B-phase. The phase diagram is presented in Figure [4.2].

In Figure [4.3](left), we show the behavior of the fermion energy gap ∆0

when the system crosses from the Az to the B phase. This happens when

the coupling strength J (Jx = Jy = J ) increases past the J = 0.5 point.

Setting K > 0 opens a gap in the B phase and since the K-part is magnetic, it breaks TRS. The opening of the gap means that a stable phase is formed. The data in Figure [4.3] were obtained from an (M, N ) = (20, 20) unit cell. We set Jz = 1 and let Jx = Jy = J vary from 0 to 1 for certain values

of K. These values for Jα and K are the same as used in [2] and we will

keep following this reference until the end of this chapter.

4.2.2 Full-vortex sector

3

These points are often called Dirac points [13], which are defined though their char-acteristic cone of linear dispersion called the Dirac cone.

(40)

J 0.2 0.4 0.6 0.8 1 " 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 K=0 K=1/15 K=2/15 K=1/5 J 0 0.2 0.4 0.6 0.8 1 " fv 0 0.5 1 1.5 2 K=0 K=1/15 K=2/15 K=1/5

Figure 4.3: Fermionic energy gap in the vortex-free (left) and full-vortex (right) sector for different values of K. (left) Vortex free sector. The system is in the Az phase for 0 ≤ J < 1/2 and in the B phase for 1/2 < J ≤ 1. The

gap in the B-phase opens for nonzero K. (right) Full-vortex sector. The phase boundary is shifted to higher values of J than in the vortex-free case and the maximum gap value is lower than in the vortex-free case.

1 2

10 20

10 100

After the vortex-free sector we look at the proper-ties of the full-vortex sector, that is formed when we put a vortex on each lattice site. The smallest periodic unit cell that we need for a full vortex lat-tice is composed of two latlat-tice sites, where one of the links, say u2,20 in the figure to the right is set

to −1. The figure shows the original lattice sites in black and the periodic doubles of these in gray, additionally decorated with a prime 0.

In Figure [4.3](right), we have displayed the fermion gap ∆f vfor the

full-vortex sector. Note that the gap still closes, but that it is smaller than it was in the vortex-free sector. What also is clear is that the phase boundary (the point where the gap closes) is shifted to higher values of J .

In [2] the same phase boundary was found, but the dependence of the energy gap on J is different from what our data show. Apart from this difference, all other data that we obtained for the Kitaev honeycomb model are identical to those presented in [2] and we cannot find a cause for this difference at this time.

4.2.3 Arbitrary flux sectors

Now that we have seen the vortex-free and full-vortex sectors, we will discuss the spectrum when we introduce a set of vortices. To compare situations, the plots we provide, address both the Abelian A-phase and the non-Abelian B-phase. For the A-phase, we choose Jx = Jy = 1₃ and Jz = 1. For the

(41)

(1, 1) (2, 1) (1, 2) (1, N ) (1, 2) (d + 1, 1)(d + 2, 1) (M, 1) (2, 1) (2, 2)

Figure 4.4: The M × N lattice with the choice of link operator-signs pro-ducing two vortices separated by a distance d.

As was mentioned earlier, we look at lattices of dimension 20 × 20 and scan across the Brillouin zone. On this lattice, we can pick links and set them to −1 to make vortices. We can only make the vortices in pairs, since the unphysical string of negative-sign link operators must terminate at two separate vortices, or go around in a loop, where in the last case, we again obtain the vortex-free sector that we treated earlier.

For the remainder of this chapter, it is useful to mention our nota-tion. For a 2n-vortex sector we analyse the first few excitation gaps i, i = 1, 2, 3, ..., which we write as ∆2nv_i . The exponent of 2nv indicates the vortex sector and i the excitation number. We will mostly talk about the vortex gap which, we will see, originates from the vortex-vortex interactions. Also the fermion gap ∆2nv_f will be mentioned, which is the excitation gap for free fermions in the model.

Two-vortex sector

We make vortices separated by a distance d, by setting d consecutive y-links to −1. We choose to set the negative-valued links in such a way that we end up with a vortex on the lower left plaquette and the other d places to the right, as is drawn in the figure below. Note that with d = 1, we mean that the vortices are next to each other. See Figure [4.4] for a visual representation.

In Figure [4.5](left), we present the fermion gap for the two-vortex sector, where the vortices are a distance d apart. We show two cases: J = 1/3 and J = 1. When the system is in the A-phase (J = 1/3), the fermion gap seems insensitive to the vortex separation distance. However, for the two-vortex sector in the B-phase (J = 1), we see that the vortex gap is closed when K = 0 and that it also soon closes for higher values of K when the vortices

(42)

d 1 2 3 4 5 6 7 8 " 2v 0 0.5 1 1.5 2 "2v 1 "2v 2

Figure 4.5: (left) The fermion gap for distances d between the vortices for several values of K. (right) The first and second excited state in the B phase, with J = 1 and K = 1/5 for separation distances d.

are brought further apart. We see that when we extrapolate d to infinity, we get a zero-mode in the B-phase for the 2-vortex sector.

In Figure [4.5](right) we have plotted the first two excited states in the B-phase for K = 1/5, which is the same as the figure in the right figure minus the lower K-value plots. We see that the vortex gap ∆2v

1 goes to zero

and that the second excited state is stable at a value of ∆2v,2 = 2, which

is the same value as the one we found when we looked at the vortex-free sector.

Four-vortex sector

We produce the four-vortex sector by creating two vortex pairs a distance of 9 sites apart in the y-direction of the lattice on the 20 × 20 lattice. As was done with the two-vortex sector, we vary the separation distance of these vortices in the x-direction and investigate the low energy modes that we obtain. For the four vortices, we obtain from the direct diagonalisation a doubly degenerate energy level with energies 4v₁ and 4v₂ , with the same values as the one we found for the two-vortex sector, 2v₁ . This degeneracy in energy levels has a clear meaning, namely that it is the interaction energy (for small distances) of the two vortex pairs. As in the two-vortex case, this energy disappears when d becomes large.

We can describe the low-energy modes, using b†_i for the low-energy (vor-tex) creation operators, as [2]:

|0i = |g.si

|1i = b†₁(p0) |g.s.i , ∆4v1 = minp0|

4v 1 (p)|

|2i = b†₂(p0) |g.s.i , ∆4v2 = minp0|

4v

2 (p)| (4.11)

|3i = b†₁(p0)b†2(p0) |g.s.i , ∆4v3 = minp0|

4v

(43)

d 1 2 3 4 5 6 7 8 " 4v 0 0.5 1 1.5 2 "4v₁ "4v₂ "4v₃ "4v₄ d 1 2 3 4 5 6 7 8 " 6v 0 0.5 1 1.5 2 "6v 1 "6v₂ "6v₃ "6v 4 "6v 5 "6v₆ "6v 7 "6v 8

Figure 4.6: The excited states of two different vortex sectors showing the energy gap against the separation distance d. For both figures the additional zero mode coming from the ground state is not depicted. (left)The first four excited states for the four-vortex sector. (right) The first eight excited states of the six-vortex sector.

|4i = b†₃(p0) |g.s.i , ∆4v3 = minp0|

4v

3 (p)| = ∆4vf ,

where b†_i occupies the ith _{energy level above the ground state. These}

excited states are plotted in Figure [4.6](a). For large d we now have four zero-modes. The first comes from the ground state, two are the degenerate single-occupied modes and the last is the double-occupied mode.

Six-vortex sector

The last sector we treat in this chapter is the one where we introduce three vortex pairs which are set equally far apart (at y = 1, y = 7 and y = 13) and pull these pairs simultaneously apart in the x-direction as was also done for the other vortex sectors. We obtain a thrice degenerate mode, one mode for each vortex pair. We will give the first excited states, given in [2] and present the obtained plots.

|0i = |g.s.i

|1i = b†₁(p0) |g.s.i , ∆6v1 = minp0|