UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A study of topological insulators in three dimensions Rosenberg, Gilad Amir 2012

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2012_fall_rosenberg_gilad.pdf [ 4.84MB ]
Metadata
JSON: 24-1.0073186.json
JSON-LD: 24-1.0073186-ld.json
RDF/XML (Pretty): 24-1.0073186-rdf.xml
RDF/JSON: 24-1.0073186-rdf.json
Turtle: 24-1.0073186-turtle.txt
N-Triples: 24-1.0073186-rdf-ntriples.txt
Original Record: 24-1.0073186-source.json
Full Text
24-1.0073186-fulltext.txt
Citation
24-1.0073186.ris

Full Text

A Study of Topological Insulators in Three Dimensions  by Gilad Amir Rosenberg B.Sc., Ben-Gurion University, 2004  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Physics)  The University Of British Columbia (Vancouver) September 2012 c Gilad Amir Rosenberg, 2012  Abstract In this work we study four interesting effects in the field of topological insulators: the Witten effect, the Wormhole effect, topological Anderson insulators and the absence of bulk magnetic order in magnetically doped topological insulators. According to the Witten effect, a unit magnetic monopole placed inside a medium with θ = 0 is predicted to bind a fractional electric charge. We conduct a first test of the Witten effect based on the recently established fact that the electromagnetic response of a topological insulator is given by the axion term θ (e2 /2πh)B · E, and that θ = π for strong topological insulators.  We establish the Wormhole effect, in which a strong topological insulator, with  an infinitely thin solenoid of magnetic half flux quantum carries protected gapless and spin filtered one-dimensional fermionic modes, which represent a distinct bulk manifestation of the topologically non-trivial insulator. We demonstrate that not only are strong topological insulators robust to disorder but, remarkably, under certain conditions disorder can become fundamentally responsible for their existence. We show that disorder, when sufficiently strong, can transform an ordinary metal with strong spin-orbit coupling into a strong topological ‘Anderson’ insulator, a new topological phase of matter in three dimensions. Finally, we lay out the hypothesis that a temperature window exists in which the surface of magnetically doped topological insulators is magnetically ordered but the bulk is not. We present a simple and intuitive argument why this is so, and a mean-field calculation for two simple tight binding topological insulator models which shows that indeed a sizeable regime such as described above could exist. This indicates a possible physical explanation for the results seen in recent experiments. ii  Preface This thesis is based almost exclusively on extensive notes written by myself throughout my program. In addition, some short sections are based on notes written by my research supervisor Marcel Franz, and on publications authored by me, my supervisor and our collaborators Huaiming Guo (who was a postdoc at UBC at the time) and Gil Refael (Prof. at Caltech). Chapter 1 motivates the rest of the thesis by providing an overview of the research field of topological insulators. It is inspired and based on several review articles published on the subject. Chapter 2 provides an introduction to tight-binding Hamiltonians in second quantization and Chapter 3 presents the main model used in this work and its topological classification. They are based on my notes, but some of the work is covered in basic Condensed Matter textbooks and papers. The model itself was my supervisor’s idea, and was based at first on symmetry arguments and a paper by Fradkin, Dagotto and Boyanovsky [1], but later we discovered it was used as a “regularized Bi2 Se3 ” model in many papers around the same time or before. Section 2.7 is based on Ref. [2]. Chapter 4 is based on my notes while working on this project. The idea was my supervisor’s, I carried out the calculations and he wrote the paper which I then helped proof and edit, titled “Witten effect in a crystalline topological insulator”, published in Physical Review B, 82, 035105 (2010). Chapter 5 is based on my notes while working on this project. The idea was my supervisor’s, I carried out the calculations and he wrote the paper which I then helped proof and edit, titled “Wormhole effect in a strong topological insulator”, published in Physical Review B 82, 041104(R) (2010). Section 5.3 was done by iii  my supervisor and is included for completeness. Chapter 6 is based on my notes while working on this project. We published a paper titled “Topological Anderson insulator in three dimensions” in Physical Review Letters 105, 216601 (2010). The sections in the paper on the Born approximation and on the numerical simulation of the Witten effect in the quasi 1D system were based on my notes while working on the project, which are presented in this chapter. Section 6.3 follows Ref. [3] closely, filling in the gaps. Chapter 7 is based on my notes while working on this project. This project started with an idea of my supervisor’s of a possible explanation for some recent experiments. I carried out the calculations and wrote the paper (based on the notes), which was then edited and proofed by my supervisor: “Surface magnetic ordering in topological insulators with bulk magnetic dopants”, published in Physical Review B 85, 195119 (2012). Section 7.5 is based on notes written by my supervisor. The appendix is based on my notes. Sections A.1 and A.2 are known (but rarely explained in detail), Section A.3 is original, and Section A.4 follows [4] (Sections 5.3-5.4), filling in the gaps.  iv  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  viii  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ix  Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiii  1  2  Introduction to Topological Insulators  . . . . . . . . . . . . . . . .  1  1.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.2  A brief history of topological insulators . . . . . . . . . . . . . .  1  General Tight Binding Treatment . . . . . . . . . . . . . . . . . . .  13  2.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  13  2.2  Real space representation . . . . . . . . . . . . . . . . . . . . . .  13  2.3  Momentum space representation . . . . . . . . . . . . . . . . . .  14  2.4  Diagonalizing the Hamiltonian . . . . . . . . . . . . . . . . . . .  16  2.5  Adding degrees of freedom . . . . . . . . . . . . . . . . . . . . .  17  2.6  Common terms . . . . . . . . . . . . . . . . . . . . . . . . . . .  17  2.7  Topological classification . . . . . . . . . . . . . . . . . . . . . .  19  v  3  4  5  6  Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  3.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  3.2  Our cubic model . . . . . . . . . . . . . . . . . . . . . . . . . .  22  3.3  Explicit calculation of θ  . . . . . . . . . . . . . . . . . . . . . .  36  The Witten Effect in a Topological Insulator . . . . . . . . . . . . .  39  4.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  39  4.2  Axions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  40  4.3  Classical view of the Witten effect I . . . . . . . . . . . . . . . .  41  4.4  Classical view of the Witten effect II . . . . . . . . . . . . . . . .  45  4.5  Adding a single monopole . . . . . . . . . . . . . . . . . . . . .  46  4.6  Charge bound to a monopole . . . . . . . . . . . . . . . . . . . .  48  4.7  Proposal for experimental realization . . . . . . . . . . . . . . . .  52  4.8  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  Wormhole Effect in a Strong Topological Insulator  . . . . . . . . .  59  5.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  5.2  Laughlin’s argument . . . . . . . . . . . . . . . . . . . . . . . .  59  5.3  Analytical calculation . . . . . . . . . . . . . . . . . . . . . . . .  61  5.4  Numerical study . . . . . . . . . . . . . . . . . . . . . . . . . . .  63  5.5  Experimental detection . . . . . . . . . . . . . . . . . . . . . . .  65  5.6  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  Topological Anderson Insulators . . . . . . . . . . . . . . . . . . . .  68  6.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  68  6.2  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  68  6.3  TAI in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  69  6.4  TAI in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  74  6.4.1  General treatment . . . . . . . . . . . . . . . . . . . . . .  74  6.4.2  Iterative solution . . . . . . . . . . . . . . . . . . . . . .  79  6.4.3  Non iterative solution . . . . . . . . . . . . . . . . . . . .  81  6.4.4  Analytical zero order Born approximation . . . . . . . . .  84  6.4.5  Witten effect in quasi 1D system . . . . . . . . . . . . . .  87  vi  7  Topological Insulator with Magnetic Impurities  . . . . . . . . . . .  92  7.1  Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  92  7.2  Experimental results . . . . . . . . . . . . . . . . . . . . . . . .  92  7.3  Surface magnetic ordering in topological insulators . . . . . . . .  93  7.3.1  Surface doping . . . . . . . . . . . . . . . . . . . . . . .  93  7.3.2  Bulk doping . . . . . . . . . . . . . . . . . . . . . . . . .  94  7.4  Determining parameters for Bi2 Se3  . . . . . . . . . . . . . . . .  95  7.5  2D surface mean field calculation . . . . . . . . . . . . . . . . . . 100  7.6  3D bulk mean field calculation . . . . . . . . . . . . . . . . . . . 103  7.7  Adding the surfaces . . . . . . . . . . . . . . . . . . . . . . . . . 119  7.8  Perovskite lattice - bulk . . . . . . . . . . . . . . . . . . . . . . . 123  7.9  Perovskite lattice - z dependent case . . . . . . . . . . . . . . . . 126  7.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 8  Conclusions  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133  8.1  Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133  8.2  Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 A Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 A.1 Derivatives in different sets of coordinates . . . . . . . . . . . . . 147 A.1.1 Complex plane . . . . . . . . . . . . . . . . . . . . . . . 147 A.1.2 Polar coordinates . . . . . . . . . . . . . . . . . . . . . . 148 A.2 Exploiting rotational symmetry . . . . . . . . . . . . . . . . . . . 151 A.2.1 Spinless case . . . . . . . . . . . . . . . . . . . . . . . . 151 A.2.2 Adding the spin . . . . . . . . . . . . . . . . . . . . . . . 157 A.3 Adding monopoles in practice: calculating the Peierls factors . . . 163 A.3.1 Adding a lattice of monopoles . . . . . . . . . . . . . . . 163 A.3.2 Adding a single monopole . . . . . . . . . . . . . . . . . 167 A.3.3 Adding a planar monopole . . . . . . . . . . . . . . . . . 170 A.4 Derivation of the self consistent Born approximation . . . . . . . 172  vii  List of Tables Table 3.1  Time reversal invariant momenta and the parity eigenvalues for our model . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Table 7.1  27  Lattice constants for Bi2 Se3 . . . . . . . . . . . . . . . . . . . 100  viii  List of Figures Figure 1.1  Chiral edge state on the edge of a quantum Hall state . . . . .  3  Figure 1.2  Trivial and topological states . . . . . . . . . . . . . . . . . .  4  Figure 1.3  Quantum Hall state vs. Quantum Spin Hall State . . . . . . .  5  Figure 1.4  Destructive interference leads to a suppression of backscattering  6  Figure 1.5  Energy spectrum of massless and massive Dirac fermions . . .  9  Figure 1.6  ARPES measurements for Bi2 Se3 . . . . . . . . . . . . . . .  10  Figure 1.7  Local density of states: energy and momentum . . . . . . . .  11  Figure 3.1  Phase diagram for our 3D topological insulator (TI) . . . . . .  28  Figure 4.1  Charge density of a monopole embedded in a strong topological insulator (STI) . . . . . . . . . . . . . . . . . . . . . . .  49  Figure 4.2  Charge density of a monopole embedded in a STI vs. the radius 50  Figure 4.3  Log-log plot of the excess charge of a monopole embedded in a STI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Figure 4.4  3D charge density plot for the layer just below the planer monopole for weak disorder . . . . . . . . . . . . . . . . . . . . . . . .  Figure 4.5  51 54  2D charge density plot for the layer just below the planer monopole for weak disorder . . . . . . . . . . . . . . . . . . . . . . . .  55  Figure 4.6  Excess charge for different values of disorder strength . . . .  56  Figure 5.1  Laughlin’s argument for a TI . . . . . . . . . . . . . . . . . .  61  Figure 5.2  Additional charge for a flux tube in a TI . . . . . . . . . . . .  64  Figure 5.3  Additional charge for a flux tube in a TI . . . . . . . . . . . .  65  ix  Figure 6.1  Born approximation phase boundary lines with conductivity showing the 2D topological Anderson insulator (TAI) phase .  74  Figure 6.2  Born approximation phase boundary lines . . . . . . . . . . .  88  Figure 6.3  Born approximation phase boundary lines with conductivity showing the 3D TAI phase . . . . . . . . . . . . . . . . . . .  Figure 6.4  89  Quasi 1D system with a monopole and anti-monopole with periodic boundary conditions (PBC) . . . . . . . . . . . . . . .  90  Figure 7.1  Bulk energies at high symmetry points . . . . . . . . . . . . . 108  Figure 7.2  The magnetization of the impurities as a function of temperature 113  Figure 7.3  Magnetization as a function of the ”side of cube” in k-space . 114  Figure 7.4  Impurity spin expectation value as a function of the z coordinate, for different temperatures . . . . . . . . . . . . . . . . . 122  Figure 7.5  Impurity spin expectation value as a function of temperature for the discretized Bi2 Se3 lattice . . . . . . . . . . . . . . . . 123  Figure 7.6  Critical temperature as a function of the chemical potential µ for the discretized Bi2 Se3 lattice . . . . . . . . . . . . . . . . 124  Figure 7.7  A unit cell of the Perovskite lattice . . . . . . . . . . . . . . . 125  Figure 7.8  Bulk impurity magnetization on each basis site for the Perovskite lattice . . . . . . . . . . . . . . . . . . . . . . . . . . 127  Figure 7.9  Ferromagnetic order parameter as a function of temperature for the Perovskite lattice . . . . . . . . . . . . . . . . . . . . 130  Figure 7.10 Antiferromagnetic order parameter as a function of temperature for the Perovskite lattice . . . . . . . . . . . . . . . . . . 131  x  Glossary FCC  face centered cubic  TI  topological insulator  STI  strong topological insulator  WTI  weak topological insulator  TAI  topological Anderson insulator  STAI  strong topological Anderson insulator  TR  time reversal  TRS  time reversal symmetry  TRI  time reversal invariance  TRIM  time reversal invariant momenta  ARPES angle resolved spectroscopy STM  scanning tunnelling microscopy  QH  quantum Hall  QSH  quantum spin Hall  FQHL  fractional quantum Hall liquid  TKNN  Thouless-Kohmoto-Nightingale-Nijs xi  BZ  Brillouin zone  FT  Fourier transform  PBC  periodic boundary conditions  FM  ferromagnet  AF  anti-ferromagnetic  xii  Acknowledgments First and foremost, I would like to thank my research supervisor Prof. Marcel Franz, for coming up with excellent ideas for projects that would be doable but still provide a challenge, having an always open door, and his unwavering willingness to dive into the details when I needed assistance. I could not envision a better supervisor. Thanks also for funding me partially for all these years. I would like to thank my committee members, Prof. Ian Affleck, Prof. Josh Folk and Prof. Gordon Semenoff for the constructive feedback and guidance. I would also like to thank the examining committee at my final defence for reading my thesis and providing feedback on it: Prof. Ian Affleck, Prof. Rob Kiefl, and Prof. Roman Krems, and the chair Prof. Joel Feldman. To the University of British Columbia (UBC) for giving me the opportunity to discover Vancouver and British Columbia, and for funding me through the Graduate Entrance Scholarship, the International Partial Tuition Scholarship, the PhD Tuition Fee Award, the University Graduate Fellowship, the Pacific Century Graduate Fellowship and the Four Year Fellowship. To Prof. Doron Cohen, my supervisor at Ben-Gurion University (BGU) while I was working there for a year as a Research Assistant after my BA. Thanks so much for giving me the opportunity for my first experience in research, and instilling in me good practices, such as writing notes throughout my PhD instead of cramming at the end. To Prof. Ehud Meron from BGU, thanks for suggesting that I apply to UBC. To my parents, far away but close to my heart, for nurturing me throughout my childhood and believing in me, supporting me and allowing me to develop in my own way. To Maya, my best friend, and my beloved partner, for her love and encouragement and of course for baking my way to the PhD. xiii  Chapter 1  Introduction to Topological Insulators 1.1  Overview  In recent years, the field of topological insulators (TIs) has attracted much attention and research in the condensed matter community. The advance has been rapid, on both the theoretical and experimental fronts, and the excitement in the community is evident. In this chapter we present a short non-technical introduction to the field, outline and motivate the problems considered in the thesis, and place them in the context of the field. Several comprehensive review papers have been published on the subject, and we refer the interested reader to them for a more detailed introduction to the field [5–9].  1.2  A brief history of topological insulators  A topological insulator is a bulk insulator in two dimensions (2D) or three dimensions (3D) which is in a topological phase. The ordered phases (or states) we are most familiar with from our lives are trivial phases, which are characterized by symmetry breaking. For example, a crystal breaks the continuous translation symmetry and a ferromagnet breaks the continuous rotation symmetry of the vacuum. In contrast, a topological phase includes all states which can be deformed contin1  uously from one to another without closing the bulk gap. These states then share some fundamental properties (topological invariants) which can only be changed through a quantum phase transition, a process in which the bulk gap must close. Therefore, a topological state must have gapless (metallic) states on its surfaces, since that area is the border between a topological phase and a trivial phase (e.g. vacuum). The simplest example of a topological state is the integer quantum Hall (QH) state [10, 11]. This state can be created by confining electrons to two dimensions and subjecting them to a strong perpendicular magnetic field (which gives Landau levels) and an in-plane electric field. The conductivity in the direction perpendicular to the electric field is then quantized σxy = n  e2 h  (1.1)  The integer n is an example of an integer topological invariant (also called a first Chern invariant). It was later shown [12] that the topological invariant n (known as the TKNN invariant, after the authors of that paper: Thouless-Kohmoto-NightingaleNijs) can be calculated by a surface integral of the Berry flux over the Brillouin zone (BZ) n=  1 2π  ∑ j  BZ  ∇ × A j d2k  (1.2)  where A j = i u j |∇k |u j is the so called Berry connection, |u j (k) are the Bloch  wave functions, and the sum runs over the contributions from all occupied bands. This result was later generalized to include interactions and disorder [13]. The  QH state features gapless chiral edge states, either spin up or down, depending on the direction of the magnetic field. This can be viewed as a result of the classical cyclotron motion of the electrons around the edge (see Fig. 1.1a). The extremely high precision with which the quantum of conductivity has been measured (1 in 109 ) is characteristic of quantization resulting from a topological invariant [14]. The bulk topological state is manifest on the edge by a metallic edge state of electrons that propagate only in one direction (which depends on the direction of applied magnetic field), as in Fig. 1.1b. If we were to switch off 2  Conduction Band  E Insulator  n=0  (a)  (b) EF  Quantum Hall n=1 State Valence Band −π /a  0  k −π /a  Figure 1.1: Chiral edge mode on the edge of a QH state (a) A classical view of the edge state: electrons in a magnetic field make cyclotron orbits, and the chiral edge state can be attributed to this motion. (b) The energy spectrum for a semi infinite system based on Haldane’s model. This can be viewed as half a topological insulator: a single edge state connects the bulk bands. Figure from Ref. [5]. the magnetic field, the material would go into an ordinary (trivial) metallic state. Fig. 1.2 illustrates how a bulk property can be manifest on the surface, in this case the genus (number of holes). In this analogy, the orange is trivial, since it can be continuously deformed into nothing (vacuum), whereas the bagel is topological, since it cannot. It is impossible to deform the orange into a bagel without closing the hole, analogous to trying to deform a topological insulator into a trivial insulator without closing the energy gap. Topological insulators in 2D, also referred to as the quantum spin Hall (QSH) phase [15], were first discovered in a model consisting of two superimposed copies of Haldane’s model [16], with opposite spins moving in opposite directions on the edge [17–19]. This can be viewed as two copies of the QH state, with opposite spin and effective magnetic field (see Fig. 1.3). For a single copy of the QH state, time reversal invariance (TRI) is broken by the magnetic field. However, for the combined version, time reversal invariance is preserved. The reason is that the time reversal operator swaps both the direction of the magnetic field and flips the spin. In this model, there are two edge states propagating in opposite direction and with opposite spin. In order for time reversal invariance to be conserved, a spin 3  (a)  Insulating State  (b)  (c)  E  EG  −π /a  (d)  Quantum Hall State  B  0 k −π /a  (e)  (f)  E hωc  −π /a  0 k −π /a  Figure 1.2: (a)-(c) illustrate the trivial state and (d)-(f) sketch the QH state. (a) An atomic insulator. (b) An insulating band structure. (d) Electrons in a magnetic field move in circles. (e) The spectrum of electrons in a magnetic field: Landau levels. (c) and (f) Two surfaces with different genus g, the number of holes in the surface. g = 0 for the sphere (c) and g = 1 for the doughnut (f). The genus is topological invariant that distinguishes the two states, similar to the Chern number. Figure from Ref. [5]. up electron traveling clockwise (for example) must match up with a spin down electron traveling counterclockwise, which is indeed the case. It turns out that this type of effective magnetic field, opposite for opposite spins, can originate in spinorbit interactions. Hence, the known topological insulators have strong spin orbit coupling, and therefore tend to be heavier elements. The QSH state and the 3D TI have an odd/even topological invariant, known as a Z2 invariant. A physical intuition for this can be understood via the notion of a Kramers’ degeneracy (also known as a Kramers pair) [20]. Kramers showed that applying the time reversal operator T twice on a state of total spin S gives T 2 = (−1)2S  (1.3)  so that electrons pick up a minus sign, due to their half integer spin. As a con4  a  b  spinless 1D chain  spinful 1D chain  4=2+2  2=1+1  QH  QSH  Figure 1.3: The QH and QSH states can be viewed as arising from separation in space. (a) A spinless one-dimensional chain, with both right and left movers. These two movers are separated in a QH state, here the right mover is on the top edge and the left mover is on the bottom edge. (b) If we add spin to the one-dimensional chain, it now has two right movers and two left movers. These are separated in a quantum spin Hall state. Both edges now have two movers going in opposite directions, but the spin is opposite. Figure from Ref. [9]. sequence, for a half integer spin system, if the Hamiltonian is time reversal (TR) invariant, then all its eigen-states are at least two-fold degenerate. To illustrate this point, assume for a moment that there exists a non-degenerate state |Ψ . Since  the Hamiltonian is TR invariant, applying the time reversal operator should give us a|Ψ , where a is the wavefunction’s eigen-number (in general complex). Applying the time reversal operator twice should then give us |a|2 |Ψ (since the TR operator includes complex conjugation), but |a|2 = −1 which leads to a contradiction.  Given TRI, backscattering is not allowed. Consider low energy elastic scat-  tering off a non-magnetic impurity on the surface/edge of a system with a single Kramers pair with opposite spin. These two states are degenerate and connected by the time reversal operator. If a right mover with spin up scatters off a non-magnetic impurity, it can either go around it clockwise or counterclockwise (see Fig. 1.4). But in order to move left, it must have spin down, so its spin must rotate by an angle of π or −π. The difference between these two paths is 2π, which for spin 5  1/2 gives a relative minus sign between them. Therefore, when adding the contributions from these two paths they cancel (destructive interference). Therefore, given TRI and one Kramers pair, no backscattering is allowed. However, if we have two Kramers pairs, a right mover of the first pair could scatter into a left mover of the second pair without spin reversal, no longer giving destructive interference and leading to dissipation. Systems with an odd number of Kramers pairs are topologically distinct (in a different topological phase) from systems with an even number of Kramers pairs, since to go from one to the other requires closing the energy gap.  a  b  Figure 1.4: (a) Destructive interference on a lens with an antireflective coating. Rays reflected by the top surface interfere with the rays reflected by the bottom surface. They reflect destructively only if the thickness of the coating matches the wavelength. (b) A quantum spin Hall edge state scatters off a non-magnetic impurity. It can go clockwise or anticlockwise, giving a spin rotation of ±π, netting a difference of 2π. When spin 1/2 particles rotate by 2π they gain a −1 factor, which leads to destructive interference of the two paths. This is similar to the optical case in (a), but the important difference is the robustness - no fine tuning is necessary in this case. Figure from Ref. [9]. The gapless surface states are said to be topologically protected by the presence of time reversal invariance. If we add non-magnetic disorder to the surface, we expect the surface states to be preserved even for relatively high disorder. However, if we break time reversal invariance, for example by coating a topological 6  insulator with a ferromagnetic layer (effectively applying a magnetic field to the surface), the surface states would be gapped. This protection is of much interest for anticipated future applications, since decoherence has proven to be a serious challenge in quantum electronics. Based on the honeycomb lattice model, Graphene was initially thought to be a good candidate for the first realization of a 2D TI, however it soon became apparent that the native spin orbit coupling is too weak to support the topological order [21–25]. The search continued, focusing on heavier elements which typically have stronger spin orbit coupling. Two dimensional TIs were predicted [26] and observed [27] in HgTe/CdTe quantum wells (for a comprehensive review, see [28]). It was natural to then ask if there is a 3D counterpart to the QSH state. If we are interested in creating a 3D TI, we might naively try to stack many two dimensional topological insulators one on top of the other. However, this creates a weak topological insulator (WTI), which has an even number of Dirac cones on its surface and is not robust to (non-magnetic) disorder. In order to get a strong topological insulator, we require a spin-orbit coupling that mixes all the spin components, and hence it cannot be constructed by layering. In the two-dimensional QH state the topological invariant n is an integer topological invariant (sometimes called a Z invariant). For the 2D and 3D TIs the invariant ν is even/odd, also known as a Z2 invariant. An insulator with ν = 0 is trivial, and an insulator with ν = 1 is non-trivial, also called a strong topological insulator (STI). The electromagnetic response of a STI is characterized by the axion term [29, 30] e2 B·E 2πh  ∆Laxion = θ  (1.4)  This can be viewed as a magnetoelectric polarizability - if one applies an electric (magnetic) field, a magnetic (electric) moment dipole is generated [29, 30]. The axion field θ can generally vary in space and time and it has been shown that in periodic space and time θ is invariant under θ → θ + 2π [31]. However, if TRI is present, only θ = 0, π are allowed, since T takes E → E and B → −B. This led to  7  the realization that θ is in fact the Z2 invariant, θ = ν0 π  (1.5)  A constant and uniform θ does not affect the equations of motion, since it couples to a total derivative. As we discuss in Chapter 4, if θ is uniform and non-zero, Witten predicted[32, 33] that magnetic monopoles would bind fractional charge −eθ /2π. In this context, topological insulators present a platform for testing  fundamental physical paradigms. As in two dimensions (see Eq. 1.2), the threedimensional topological invariant can be calculated by an integral over the BZ of a Berry connection[29, 30] θ= where Ai  αβ  1 4π  BZ  ε i jk Tr Ai ∂ j Ak +  2i Ai A j Ak d 3 k 3  (1.6)  = −i uα |∂i |uβ are the elements of the Berry connection matrices  (here A is a vector of matrices), |uα (k) are the Bloch wavefunctions, ∂i = ∂ /∂ ki  and the trace runs over the occupied states. In practice, this expression is cumbersome to use. In many cases the model studied has an inversion symmetry and then a simpler expression can be used (in both 2D and 3D) [2] (−1)ν0 = ∏ ξ2m (Γi )  (1.7)  m,i  where Γi are the 4 (in 2D) or 8 (in 3D) time reversal invariant momenta (TRIM), ξ2m are the parity operator eigenvalues of the 2m-th occupied energy band (recall that at the TRIM they must be degenerate) and the product runs over the TRIM Γi and the Kramers pairs indexed by m. We use this method in Chapter 3 to topologically classify our model. The first model shown to be a TI in 3D was a toy model on the diamond lattice, followed soon after by a prediction and observation that the semiconducting alloy Bi1−x Sbx is a 3D TI [2]. The band gap of this material is rather small and the surface state is complicated, so the search continued for a more suitable candidate for experiments, which led to the discovery of the so called “second generation” of 3D TIs. These feature a single Dirac cone on the surface, observed in 8  angle resolved spectroscopy (ARPES) measurements, and larger band gaps. One of the most studied and most promising 3D TIs is the semiconducting thermoelectric Bi2 Se3 , with a relatively large band gap of ∼ 0.3eV (equivalent to about 3600K)  [34, 35]. The spectrum of Bi2 Se3 and other TIs has been studied using ARPES  [35–37] (see Fig. 1.6) and scanning tunnelling microscopy (STM) [38–40], showing that the surface states form an almost ideal Dirac cone, familiar from studies of graphene [41] (see Fig. 1.7). A simple tight-binding model on the cubic lattice has been used successfully to study Bi2 Se3 , and this is the model most of the numerical calculations in this work are based on [34, 42, 43]. We note that although TIs have been identified by studying the spectrum of the surface, the bulk spectrum should appear the same (qualitatively) for a topological insulator and ordinary insulator. The surface states of a topological insulator exhibit an odd number of gapless  1  Dirac cones, in other words, their surface states have a linear spectrum E = ±v|k|  (see Fig. 1.5a). This also implies a spin-momentum locking (in three dimensions)  - there is a spin direction associated with each momentum vector. If we cut the Dirac cone into circles, the spin is tangential to the circle at each point. For each surface electron with momentum k and spin up (down), we find an electron with momentum −k and spin down (up) on the opposite side of the circle. This shows  (a) reversal invariance is(b) that indeed time fulfilled.  (a)  2.5  (b)  �SF M �(µ B )  1.5 1.0  �SAF �(µ B )  Figure 1.5: (a) Energy spectrumsurf1 of a massless Dirac fermion. The bottom 2 surf2 fully occupied at half filling. The top T cbulk (red) cone is the valence band, z bulk T c (green) cone is the conduction band, which 1 is assumed to be vacant. (b) bulk Massive Dirac fermion. This figure shows how the opening of a gap 0 tends to lower the free energy. surf  2.0  Tc  ï1  0.5  0  2  T csurf  ï2  0.0 (a)  surf1 surf2 z avg bulk  20  40  T  60  80  avg2 avg4  9  0 (b)  2  20  40  T  60  80  0 40 60  NATURE PHYSICS DOI: 10.1038/NPHYS1274  LETTERS  Low  High  a  b 0.1  c  EF  0.1 Γ  M  K  M  0  Γ  K  0 kx (Ŭ1)  0.15  0  EB (eV)  ¬0.1  EB (eV)  ¬0.1  Intensity (a.u.)  SS  SS  ¬0.2  ¬0.2  ¬0.3  ¬0.3  ¬0.4  ¬0.4 ¬0.15  d  0 ky (Ŭ1)  kz  ¬0.15  0.15  f  (111)  ¬0.12  0 ky (Ŭ1)  0.12  Theoretical calculations  1.0 k K measurements of electronic dispersion on the [111] surFigure 1.6: ARPES Γ face of Bi2 Se3 . The incident photon energy was 22 eV, and data was Bulk M collectd near the Γ-point along 0.5 the Γ-MNo(a) SOC and Γ-K (b) lines in momenSOC tum space. (c) The two surface lines meet at a Dirac point at Γ. Figure k Surface No SOC Γ SOC and caption from Ref. [35]. Z F x  Energy (eV)  y  L  e  0  Breaking TRI, for example by adding magnetic dopants, is expected to open a  M ¬0.5 provided that the magnetic moments are gap in the spectrum of the surface states, 2  ky  M1  Γ  M3  ordered perpendicular to the surface. The resulting spectrum then resembles that of a ‘massive’ Dirac fermion (Fig. 1.5b). ¬1.0 There is considerable interest in having k x  K  Γ  M  a system with an odd number of massive Dirac fermions, since it is predicted to  Figure 1 | Strong spin–orbit interaction gives rise to a single SS Dirac cone. Theory (see the Methods section) versus experiments. a,b, High-resolution ARPES measurements of surface electronic band dispersion on Bi2 Se3 (111). Electron dispersion data measured with an incident photon energy of 22 eV near the -point along the –M (a) and –K (b) momentum-space cuts. c, The momentum distribution curves corresponding to a suggest that two 2 at . The V-shaped pure SS band pair observed in a–c is nearly isotropic in the momentum plane, forming surface bands converge into a single Dirac point a Dirac cone in the energy–kx –ky space (where kx and ky are in the –K and –M directions, respectively). The U-shaped broad continuum feature inside the V-shaped SS corresponds roughly to the bottom of the conduction band (see the text). d, A schematic diagram of the full bulk three-dimensional BZ of Bi2 Se3 and the two-dimensional BZ of the projected (111) surface. e, The surface Fermi surface (FS) of the two-dimensional SSs along the K– –M momentum-space cut is a single ring centred at if the chemical potential is inside the bulk bandgap. The band responsible for this ring is singly degenerate in theory. The TRIMs on the (111) surface BZ are located at and the three M points. The TRIMs are marked by the red dots. In the presence of strong spin–orbit coupling (SOC), the surface band crosses the Fermi level only once between two TRIMs, namely and M; this ensures the existence of a π Berry phase on the surface. f, The corresponding local density approximation (LDA) band structure (see the Methods section). Bulk band projections are represented by the shaded areas. The band-structure topology calculated in the presence of SOC is presented in blue and that without SOC is in green. No pure surface band is observed to lie within the insulating gap in the absence of SOC (black lines) in the theoretical calculation. One pure gapless surface band is observed between and M when SOC is included (red dotted lines).  exhibit many interesting topological phenomena, including the half quantum Hall effect on the surface (e /2h Hall conductance) [30], the image magnetic monopole  (an electric charge adjacent to a TI results in the field of a magnetic monopole  embedded in the TI) [44, 45], and a Kerr/Faraday angle quantization in units of the fine structure constant [46, 47]. A tunable gap would also allow the control  of the surface transport, and could in addition lead to unique practical applications associated with purely electric control of the surface magnetization [48, 49]. A  signature of the been observed recently using ARPES −1 experiment. The ‘V’ bands crossmassive EF at 0.09 ÅDirac along fermion –M and athasThe calculated band structure with and without SOC are overlaid 0.10 Å−1 along –K, and have nearly equal band velocities, approx- together for comparison. The bulk band projection continuum on magnetically Bi2 SeA3continuum-like , [50, 51] although interesting effects 5 imately 5in × 10 m s−1 , along thedoped two directions. the (111) the surface is represented by the associated shaded areas, blue with manifold of states—a filled U-shaped feature—is observed inside SOC and green without SOC. In the bulk, time-reversal symmetry with it have yet to be seen in a laboratory. A surprising feature of these experiments the V-shaped band pair. All of these experimentally observed demands E(k,↑) = E(−k,↓) whereas space inversion symmetry features can be identified, to first order, by a direct one-to-one demands E(k, ↑) = E(−k, ↑). Therefore, all the bulk bands are comparison with the LDA band calculations. Figure 1f shows the doubly degenerate. However, because space inversion symmetry theoretically calculated (see the Methods section) (111)-surface is broken at the terminated surface in the experiment, SSs are electronic structure of bulk Bi2 Se3 along the K– –M k-space cut.10 generally spin-split on the surface by spin–orbit interactions except NATURE PHYSICS | VOL 5 | JUNE 2009 | www.nature.com/naturephysics  © 2009 Macmillan Publishers Limited. All rights reserved.  399  NATURE PHYSICS DOI: 10.1038/NPHYS1270 Sb2Se3  8  Energy (eV)  0.05 0  ¬0.15  4  0.1  2  ¬4  Γ  M  0.3  4  0.1  2  0  0  ¬0.1  ¬2  ¬0.2  ¬4  ¬0.3 K  Γ  M  ¬6  ¬4 Γ  K  d  0.2  6  0.2  ¬2  ¬0.2  ¬8  8  Bi2Se3  0  0 ¬0.1  ¬6  0.4  Energy (eV)  6  0.2  ¬2  K  8  Sb2Te 3  4 0  ¬0.05  0.4 0.3  2  ¬0.10  c  b  6  Energy (eV)  0.15 0.10  M  ¬6 8  Bi2Te 3  6  0.1  Energy (eV)  a  ARTICLES  4 0  2  ¬0.1  0  ¬0.2  ¬2 ¬4  ¬0.3 K  Γ  M  Figure 4 | Surface states. a–d, Energy and momentum dependence of the LDOS for Sb2 Se3 (a), Sb2 Te3 (b), Bi2 Se3 (c) and Bi2 Te3 (d) on the [111] surface. Here, the warmer colours represent higher LDOS. The red regions indicate bulk energy bands and the blue regions indicate bulk energy gaps. The surface states can be clearly seen around the point as red lines dispersing in the bulk gap for Sb2 Te3 , Bi2 Se3 and Bi2 Te3 . No surface state exists for Sb2 Se3 .  Figure 1.7: Local density of states (warmer colours are higher values) asother it depends onfrom momentum andtoenergy fortheSblow-energy Sb2 Te3 (b),properties of be carried out on the three materials, which we see that characterize 2 Se3 (a),long-wavelength Sb2 Te3 and Bi2 Te3 are qualitatively the same as Bi2 Se3 , whereas the the system. Starting from the four low-lying states |P1+ Bi2 Se3enough (c), toBiinduce (d)anon the [111]andsurface. Red areas correspond to z , ↑ (↓) 2 Te3such SOC of Sb2 Te3 is not strong inversion. |P2− point, such a Hamiltonian can be z , ↑ (↓) at the 25 constructed by the theory of invariants for the finite wave bulk energy bands and blue regions correspond to bulk energy gaps. Topological surface states vector k. On the basis of the symmetries of the system, the The red lines near the Γ point show the surface states clearly for Sb Te , be written 2 3can The existence of topological surface states is one of the most generic form of the 4 × 4 effective Hamiltonian 2 important properties of the topological insulators. To see the down up to the order of O(k ), and the tunable parameters in Bi Se3 , Bi , but they do not exist for Sb2 Se3can . Figure from Ref. [34]. 2 Te3explicitly, topological features of2the four systems we calculate the the Hamiltonian be obtained by fitting the band structure surface states of these four systems on the basis of an ab initio of our ab initio calculation. The important symmetries of the calculation. First we construct the maximally localized Wannier system are time-reversal symmetry T , inversion symmetry I and function (MLWF) from the ab initio calculation using the method three-fold rotation symmetry C3 along the z axis. In the basis of is that the gap in the surface spectrum without bulk magnetic ordering, 21,22 − + − developed by Marzari and co-workers . We divide the semi-appears (|P1+ z , ↑ , |P2z , ↑ , |P1z , ↓ , |P2z , ↓ ), the representation of the infinite system into a surface slab with finite thickness and the symmetry operations is given by T = K · iσ y ⊗ I2×2 , I = I2×2 ⊗ τ3 even the areparameters uniformly everywhere in the sample. z remaining partthough as the bulk. Thedopants MLWF hopping for thedistributed and C3 = exp(i(π/3)σ ⊗I2×2 ), where K is3D the complex conjugation bulk part can be constructed from the bulk ab initio calculation, and operator, σ x,y,z and τ x,y,z denote the Pauli matrices in the spin and Wefordiscuss analyze a possible to this puzzle in Chapter the ones the surfaceand slab can be constructed from thesolution ab initio orbital space, respectively. By requiring7.these three symmetries and calculation of the slab, in which the surface correction to the lattice keeping only the terms up to quadratic order in k, we obtain the Another peculiar effect thatself-consistently occurs in a strong insulator is the ’wormconstants and band structure have been considered followingtopological generic form of the effective Hamiltonian: and the chemical potential is determined by the charge neutrality hole’ effect. Itand is surface well MLWF known thatparameters, an infinitely thin solenoidcarrying magnetic flux  condition. With these bulk hopping M(k) A1 kz 0 A2 k− we use an iterative method23,24 to obtain the surface Green’s Φ (a ‘Dirac string’) inserted into an band insulator has significant A1 kno −M(k) A2 k−effect0  z function of the semi-infinite system. The imaginary partordinary of the H (k) = 0 (k)I4×4 +  0 A2 k+ M(k) −A1 kz  surface Green’s function is the local density of states (LDOS), from on the spectrum of electrons. In a STI, such a solenoid carries protected gapless A2 k+ 0 −A1 kz −M(k) which we can obtain the dispersion of the surface states. The surface LDOS on the [111] surface for all four systems is shown in Fig. 4. For one-dimensional fermionic modes when Φ = hc/2e. These modes are spin-filtered 2 + o(k ) (1) Sb2 Te3 , Bi2 Se3 and Bi2 Te3 , one can clearly see the topological surface states that form a single Dirac cone at the point. In comparison, and represent a distinct bulk manifestation of the topologically non-trivial insulaSb2 Se3 has no surface state and is a topologically trivial insulator. Thus, the surface-state calculation agrees well with the bulk parity with k± = kx ± iky , 0 (k) = C + D1 kz2 + D2 k⊥2 and M(k) = M − B1 tor.andWe establish this the ‘wormhole’ effect by kboth general qualitative considerations 2 2 analysis, confirms conclusively topologically non-trivial z − B2 k⊥ . By fitting the energy spectrum of the effective nature of the three materials. For Bi2 Se3 , the Fermi velocity of the Hamiltonian with that of the ab initio calculation, the parameters and surface by numerical calculations within a minimal latticemodel model in determined. Chapter For 5 Bi2 Se3 , our fitting can be topological states is vF 5.0 × 105 m s−1 , which is similar to in the effective 2 leads to M = 0.28 eV, A1 = 2.2 eV Å, A2 = 4.1 eV Å, B1 = 10 eV Å , that of the other two materials. 2 2 2 Topological insulators survive under non-magnetic long as it is B2 = 56.6 eV Å , C =disorder, −0.0068 eV, Das 1 = 1.3 eV Å , D2 = 19.6 eV Å . Low-energy effective model Except for the identity term 0 (k), the Hamiltonian (1) is not too strong compared with the band gap. Once disorder is with of the order of along nothing but the the 3D Dirac model uniaxial anisotropy As the topological nature is determined by the physics near the point, it is possible to write down a simple effective Hamiltonian the z-direction and k-dependent mass terms. From the fact NATURE PHYSICS | VOL 5 | JUNE 2009 | www.nature.com/naturephysics  11  © 2009 Macmillan Publishers Limited. All rights reserved.  441  the band gap, we might expect the topological insulator to break down. However, recently, in a remarkable development, it has been noted first by numerical simulations [52] and shortly thereafter by analytical studies [3], that a phase similar to the two dimensional topological insulator [27, 53] can be brought about by introducing non-magnetic disorder into a 2D metal with strong spin orbit coupling. This new 2D topological phase, referred to as topological Anderson insulator (TAI), has a disordered insulating bulk with topologically protected gapless edge states that give rise to precisely quantized conductance e2 /h per edge. In a TAI, remarkably, conductance quantization owes its very existence to disorder. In Chapter 6 we show that this scenario is possible in three-dimensions as well.  12  Chapter 2  General Tight Binding Treatment 2.1  Overview  Our numerical work is based on the tight-binding model, an approximation which is often used in the study of weakly correlated insulators. We neglect interactions, which could be included to some degree as perturbations of this model. The tight binding model can provide a reasonable approximation if the interaction strength is much smaller than the bulk gap. When dealing with tight-binding Hamiltonians in second quantization, there are some standard steps which are typically undertaken to find the energies and wave functions. We outline these steps here, and when dealing with specific examples we will skip some of the details presented below. We start with the spinless particle case, and proceed to add spin at the end of this section.  2.2  Real space representation  We assume that the system is periodic over a lattice given by vectors Ri . The periodicity of the lattice motivates a change to the momentum space representation, which simplifies the Hamiltonian and allows us to diagonalize it efficiently. We write our Hamiltonian in real space as H = ∑ c†i Mi j c j i, j  13  (2.1)  where c†i creates a particle (typically an electron) on site i located at Ri , and M is an infinite discrete matrix which contains the physical information about the sites (on-site energy, hopping, interactions, etc). Note that Mi j can only depend on the relative position Ri − R j (and not the absolute position), due to the periodicity of the lattice.  The Hamiltonian must be Hermitian (H = H † ), so that H † = ∑ c†j Mi∗j ci = ∑ c†i M ∗ji c j  (2.2)  i, j  i, j  We see that Mi j = M ∗ji so that M is Hermitian too. Note that if Mi j is real then the Hermitian condition gives Mi j = M ji , i.e. it is symmetric under inversion, which means that there is no difference between hopping from i to j compared with hopping from j to i. This is the case, for example, for a model with on-site energy ε and hopping amplitude ti j with no magnetic field (so ti j is real). In this case Mi j = εδi j − ti j , which is indeed symmetric under inversion.  Often the basis contains more than one site, with each type of site sitting on a separate Bravais sublattice. For example, the honeycomb lattice can be described as two triangular sublattices, and the diamond lattice can be described as two face centered cubic (FCC) sublattices. In this case, we often define new creation operators for each sublattice, and rewrite the real space representation using the sublattice operators. The rest of the treatment is the same, so we do not treat this case separately here.  2.3  Momentum space representation  Next we convert to a Fourier representation c†k = c†i =  1 √ ∑ eik·Ri c†i N i 1 √ ∑ e−ik·Ri c†k N k  14  (2.3) (2.4)  where c†k creates a particle with momentum k and N is the number of lattice points.We substitute the Fourier representation into the Hamiltonian H = ∑ c†i Mi j c j = i, j  1  ∑ c†k ck N ∑ Mi j eik ·R −ik·R ∑ c†k ck ∑ Mi j eik·δ  k,k  =  i  (2.5)  i, j  k,k  =  j  δ  1 ei(k −k)·R j N∑ j  ∑ c†k ck ∑ Mi j eik·δ = ∑ c†k H(k)ck k  k  δ  where δ ≡ R j − Ri , and in the last line we have defined the Fourier transform of M as  H(k) ≡ ∑ Mi j eik·δ  (2.6)  δ  This is the momentum space analog of (2.1). Note that H(k) is invariant under translation (in k-space) by a reciprocal lattice vector K, i.e. H(k + K) = H(k). This can be shown easily H(k + K) = ∑ Mi j eik·δ eiK·(R j −Ri ) = ∑ Mi j eik·δ = H(k)  (2.7)  δ  δ  since for any lattice vector R and reciprocal lattice vector K we have eiK·R = 1. Using the definition (2.6) and the Hermitianity of M, we prove that H(k) is Hermitian too H † (k) = ∑ Mi∗j e−ik·δ = ∑ M ∗ji eik·δ = ∑ Mi j eik·δ = H(k) δ  δ  (2.8)  δ  In (2.5) we used the fact that δkk =  1 ei(k−k )·Ri N∑ i  15  (2.9)  This is intuitively true for k = k since then the exponent on the right gives just 1, and the sum counts the number of points on the lattice, giving N sites, and therefore N/N = 1. It can also be proved by substituting (2.4) into (2.3) 1 c†k = √ ∑ eik·Ri c†i = ∑ c†k N i k  1 ei(k−k )·Ri N∑ i  (2.10)  On the other hand, we can write c†k = ∑ c†k δkk  (2.11)  k  By comparing the last two equations we get (2.9). Similarly, by substituting (2.3) into (2.4) we get δi j =  1 eik·(R j −Ri ) N∑  (2.12)  k  2.4  Diagonalizing the Hamiltonian  In the next step, we wish to diagonalize the Hamiltonian in momentum space. The periodicity of the lattice guarantees that the Hamiltonian will be block diagonal, with each block corresponding to a value of k. This can be seen by writing H = ∑ Ψ†k H(k)Ψk  (2.13)  k  where Ψk is a vector containing the basis of creation operators. If the primitive cell of the direct lattice has N inequivalent sites, Ψk is a vector with N elements, and H(k) is an N × N matrix. Due to the block diagonal structure of the Hamiltonian,  we diagonalize the whole Hamiltonian by performing just one diagonalization - we diagonalize the matrix H(k) to find N energies Eki and N wavefunctions ϕki .  16  2.5  Adding degrees of freedom  It is possible to add to this treatment additional degrees of freedom such as spin or different orbitals on each site. For example, for spin s H=  c†iσ Mi jσ σ c jσ  ∑  (2.14)  i, j,σ ,σ  where M is now a four dimensional tensor, so that Mi j is a (2s + 1) × (2s + 1)  matrix, and σ , σ run from −s to s.  The momentum representation is written as in (2.3)-(2.4) with the additional  spin index c†k → c†kσ and c†i → c†iσ so that the momentum space representation of the Hamiltonian is  H=  c†kσ H(k)ckσ  (2.15)  H(k) ≡ ∑ Mi jσ σ eik·δ  (2.16)  ∑ k,σ ,σ  where  δ  Once again, we write the Hamiltonian as H = ∑k Ψ†k H(k)Ψk , but now Ψk is a 2s+1 vector, and H(k) is a (2s + 1)N × (2s + 1)N matrix.  2.6  Common terms  The most basic tight binding Hamiltonian we will consider is H0 = ε ∑ c†i ci − t i  ∑ c†i c j  (2.17)  i, j  where ε is an on-site energy, t is the real nearest neighbor hopping, and i, j stands for all pairs of sites i and j that are nearest neighbors. In general there could be different hopping amplitudes ti j for different pairs of sites. To add magnetic field we make ti j complex, where the Peierls phase factors θi j  17  encode the magnetic field [54] θi j =  2π Φ0  Rj Ri  A · dl  (2.18)  and H0 = ε ∑ c†i ci − t i  ∑ eiθ  ij  c†i c j  (2.19)  i, j  The hopping term is gauge invariant, as can be seen by making a gauge transformation ci →  eiχi ci  (2.20)  θi j → θi j + χi − χ j So that eiθi j c†i c j → eiθi j +iχi −iχ j c†i e−iχi c j eiχ j = eiθi j c†i c j  (2.21)  This Hamiltonian does not involve spin, so [H0 , Sz ] = 0. If we add spin and assume that the hopping term is spin independent, we can write the Hamiltonian as H0 = ε ∑ c†iσ ciσ − t i,σ  eiθi j c†iσ c jσ  ∑  (2.22)  i, j ,σ  We can decompose this Hamiltonian into separate Hamiltonians for each value of sz (the eigenvalues of Sz ), and H(k) will be block diagonal, with a block corresponding to each value of sz . For example, for spin  1 2  we can write H0 = H0,↑ + H0,↓ , and  H(k) can be written in block diagonal form, with two blocks. In general, we can write H0 = ∑σ H0,σ with H0,σ ≡ ε ∑ c†iσ ciσ − t i  ∑ eiθ  ij  c†iσ c jσ  (2.23)  i, j  Spin-orbit terms are common in tight-binding models of topological insulators. A term we will use (based on [17], for a derivation see [55]) which has been  18  modified here to allow for magnetic field is HSO = iα  ∑  e  i, j  where  i, j  iθik iθk j  e  c†i↑  T  c†i↓  c j↑ σ · (dˆik × dˆk j ) c j↓  (2.24)  stands for summation over all next nearest neighbors, dˆik is a unit  vector pointing from site i to its nearest neighbor k, σ = (σ1 , σ2 , σ3 ) is the vector of Pauli matrices and θi j are the Peierls phase factors.  2.7  Topological classification  This section is based on [2]. Motivated by the lattice we will use below, assume we have a bipartite lattice where the sublattices are exchanged by inversion, and a spin degree of freedom. The Pauli matrices in the sublattice subspace will be denoted by τµ and Pauli matrices in the spin subspace (assuming spin 1/2) will be denoted by σµ . Therefore, the parity operator is the operator that swaps sublattices but leaves the spin unchanged P ≡ τx ⊗ I  (2.25)  For example, we define the operator c†k,σ that creates an electron with momentum † creates an electron on sublattice k and spin σ on sublattice A, and similarly dk,σ  B. Now define a vector Ψk ≡ (ck↑ , ck↓ , dk↑ , dk↓ ), and check that P indeed swaps  sublattices    1 0    ck↑      dk↑          ck↓  dk↓   0 1   =   PΨk =   d   c  1 0   k↑   k↑   dk↓ ck↓ 0 1  (2.26)  Note that P2 = I, and therefore P = P−1 . Now define a time-reversal operator Θ ≡ i(I ⊗ σy )K  19  (2.27)  where K is the complex conjugation operator, and the inverse time-reversal operator Θ−1 ≡ iK(I ⊗ σy )  (2.28)  ΘΘ−1 = [i(I ⊗ σy )K] [iK(I ⊗ σy )] = −i2 (I ⊗ σy2 ) = I ⊗ I  (2.29)  So that  We verify that spin changes sign under the time-reversal transformation Θ(I ⊗ σµ )Θ−1 = [i(I ⊗ σy )K] (I ⊗ σµ ) [iK(I ⊗ σy )]  (2.30)  = (−1)µ+1 (I ⊗ σy σµ σy ) = −(I ⊗ σµ )  also note that Θ2 = [i(I ⊗ σy )K] [i(I ⊗ σy )K] = −I ⊗ I  (2.31)  Given a lattice Hamiltonian H = ∑i, j c†i Mi j c j , we define a Bloch Hamiltonian H(k), as in (2.6) H(k) ≡ ∑ Mi j eik·(Ri −R j )  (2.32)  i, j  We have already shown that H(k + K) = H(k), where K is a reciprocal lattice vector. Assuming that H is symmetric under inversion around the origin (Ri → −Ri for all i), we have [H, P] = 0, or equivalently H = PHP−1 . For time reversal, we find ΘH(k)Θ−1 = ∑ ΘMi j eik·(Ri −R j ) Θ−1 = ∑ Mi j e−ik·(Ri −R j ) = H(−k) i, j  (2.33)  i, j  At special points in the BZ, we have Γ = −Γ + G, so that H(Γ) = H(−Γ). These  points are called time reversal invariant momenta (TRIM), and can be written as Γi = (nx bx + ny by + nz bz )/2 20  (2.34)  with i = (nx ny nz ) and n j = 0, 1. We see that the number of TRIM in dimension d is 2d , in particular: 4 for 2D, and 8 for 3D. There are several ways of calculating the topological invariants. In many cases the model studied has an inversion symmetry (as in our case) and then a simple expression can be used (in both 2D and 3D) [2] for the strong topological invariant ν0 (−1)ν0 = ∏ δi  (2.35)  δi = ∏ ξ2m (Γi )  (2.36)  i  where  m  The index i = n1 n2 n3 (ni = 0, 1) runs over all the TRIM and ξ2m are the parity operator eigenvalues of the 2m-th occupied energy band (recall that at the TRIM they must be degenerate). At times, we might wish to also calculate the three weak topological invariants νi (with i = 1, 2, 3), which are given by (−1)νi =  ∏  ni =1,n j=i =0,1  δn1 n2 n3  We use this method in Chapter 3 to topologically classify our model.  21  (2.37)  Chapter 3  Model 3.1  Overview  We describe the basic cubic model with two orbitals per site used in most of the numerical work in this thesis and discuss its topological classification.  3.2  Our cubic model  Basic model The sites are on a cubic lattice of spacing a, with primitive vectors aµ = aµˆ (µ = x, y, z). We look at a nearest neighbor spin-orbit term which was suggested by Fradkin, Dagotto and Boyanovsky [1], but here we modify it so that the hopping amplitude is complex (the reason will become clear below) HSO = λ  ∑  iΨ†i σµ Ψi+µ + h.c.  (3.1)  µ=x,y,z i  with Ψi ≡ (ci↑ , ci↓ ). This can also be written as HSO = λ  ∑  (ic†i,σ σµ,σ σ ci+µ,σ + h.c.)  (3.2)  µ=x,y,z i,σ ,σ  where λ is the strength of the spin-orbit term, σµ are the Pauli matrices, c†iσ cre22  ates an electron on site i with spin σ , and site i + µ is taken to mean the nearest ˆ Note that the hopping amplitude is iλ for bonds pointing neighbor in direction µ. in directions x, ˆ y, ˆ zˆ, and −iλ for −x, ˆ −y, ˆ −ˆz. Therefore, the sum of all phase factors  (±i) on each face is zero. We rewrite the Hamiltonian in momentum space c†k↑  HSO = iλ ∑  T  c†k↓  k  = −2λ ∑ k  = −2λ ∑ k  c†k↑  σx ( eikx a − e−ikx a ) + σy ( eiky a − e−iky a ) + σz ( eikz a − e−ikz a ) T  [σx sin(kx a) + σy sin(ky a) + σz sin(kz a)]  c†k↓ c†k↑  T  c†k↓  ∑ σµ sin(kµ a) µ  ck↑  ck↑ ck↓ (3.3)  ck↓  ck↑ ck↓  Now it becomes obvious why we added an i to the hopping amplitude - without it, we would have gotten cosines instead of sines, which does not give a Dirac fermion at k = (0, 0, 0). We define Ψk ≡ (ck↑ , ck↓ ), so that HSO = ∑k Ψ†k HSO (k)Ψk , and HSO (k) = −2λ ∑ σµ sin(kµ a) = −2λ µ  sin(kz a) sin(kx a) + i sin(ky a)  sin(kx a) − i sin(ky a) − sin(kz a)  This is a sum of anti-commuting matrices, so that the energies can be written immediately εk = ±2λ  sin2 (kx a) + sin2 (ky a) + sin2 (kz a)  (3.5)  The two energy bands are degenerate at Dirac points located at K=  π (nx , ny , nz ) a  (3.6)  where nµ ∈ Z. The reciprocal lattice primitive vectors are bµ = (2π/a)µˆ (µ =  x, y, z), so that there are 8 inequivalent Dirac points, defined by nµ = (0, 1) for µ = x, y, z. Note that these correspond exactly to the time reversal invariant momenta (TRIM). In order to get a topological insulator, the spectrum must be gapped over the  23  (3.4)  whole Brillouin zone. However, a general 2 × 2 Hermitian matrix can be written  using the three Pauli matrices and the identity, where the identity shifts both bands and therefore cannot give a gap. We have already used up the three Pauli matrices, so we cannot open a gap without modifying the Hamiltonian. To get a gap, we can enlarge the Hamiltonian by adding a degree of freedom, and then add a gap inducing term. One way of enlarging the Hamiltonian is to break translational symmetry, for example by adding a staggered potential which would make the Hamiltonian 4 × 4 (this is the rock-salt structure), which gives a trivial insulator.  We choose instead to add an orbital to each site, and show below that it gives a topological insulator (for a range of parameters). Adding an orbital Each site now has two orbitals, so that we can have a total of four electrons † on the two orbitals, and on each site. If we define creation operators c†iσ and diσ  Ψk ≡ (ck↑ , ck↓ , dk↑ , dk↓ ), then the spin-orbit term is HSO (k) = −2λ ∑ sin(kµ a)(τz ⊗ σµ )  (3.7)  µ  In four dimensions there are five anti-commuting matrices. Three are already present here, the other two are τx ⊗ I and τy ⊗ I. The first term is a hopping term be-  tween orbitals, which does not depend on spin or direction of hopping. The most simple terms of this type are an on-site orbital hopping, and a nearest neighbor orbital hopping term. In real space these two terms can be written as Ht = ε ∑(c†iσ diσ + h.c.) − t i,σ  ∑  (c†iσ d jσ + h.c.)  (3.8)  i, j ,σ  In momentum space Ht = ∑{c†kσ dkσ (ε − 2t[cos(kx a) + cos(ky a) + cos(kz a)]) + h.c.}  (3.9)  Ht (k) = mk (τx ⊗ I)  (3.10)  k,σ  so that  24  where mk ≡ ε − 2t ∑µ cos(kµ a). Therefore, the full Hamiltonian H ≡ Ht + HSO  gives  H(k) = mk (τx ⊗ I) − 2λ ∑ sin(kµ a)(τz ⊗ σµ )  (3.11)  µ  As before, this is a sum of anti-commuting matrices, so that the energies can be written immediately Ek = ± εk2 + m2k  (3.12)  At half filling, the lower band will be filled and the upper band will be empty, so that the Fermi energy will be positioned precisely between the two bands EF = 0. Low energy excitations can be described by expanding the Hamiltonian around the Dirac points. We take k = K + q, and expand the Hamiltonian to first order in q. For example, choose K0 = (0, 0, 0), so that sin(kµ a)  qµ a and cos(kµ a)  1 so  that H(K0 + q)  [ε − 6t] (τx ⊗ I) − 2λ a ∑ qµ (τz ⊗ σµ )  (3.13)  µ  so that the energies are EK0 +q = ± (2λ a)2 (q2x + q2y + q2z ) + (ε − 6t)2  (3.14)  We see that near K0 the spectrum is gapped with mass mK0 = ε − 6t. Near one of the Dirac points K, the spectrum will have a mass  mK = ε − 2t ∑ cos(Kµ a) = ε − 2t ∑(−1)nµ = µ  µ     ε − 6t     ε − 2t   ε + 2t     ε + 6t  (3.15)  Topological classification We would like to answer the question, is this a topological insulator, and under  25  what conditions will it be an ordinary insulator? Here we use a simple method for topologically classifying inversion symmetric TIs [2]. By inspecting (3.11), we identify d0 = 0  (3.16)  d1 = ε − 2t[cos(kx a) + cos(ky a) + cos(kz a)]  d2 = 0  d3 = −2λ sin(kx a) d4 = −2λ sin(ky a) d5 = −2λ sin(kz a)  In order to classify this insulator, we must find δi = −sgn (d1 (k = Γi )) for each  of the 8 TRIM, defined by Γi = 21 (nx bx + ny by + nz bz ) with i = (nx ny nz ) and nµ ∈  {0, 1}. For a cubic lattice of spacing a, the reciprocal lattice vectors are bµ = ˆ so that we find (2π/a)µ, Γi=(nx ny nz ) =  π (nx , ny , nz ) a  (3.17)  Note that for zero gap (ε = 6t) there is a Dirac point at each of the TRIM (see (3.6)). We see that d1 (k = Γi ) = ε − 2t ∑ cos(kµ a) = ε − 2t ∑ cos(nµ π) µ  (3.18)  µ  = ε − 2t ∑(−1)nµ µ  This expression can be evaluated easily for each Γi . We summarize the results in a table below, assuming t > 0 (for t < 0, take t → −t and ε → −ε). There are five  columns for δi , ordered by increasing ε: ε < −6t, −6t < ε < −2t, −2t < ε < 2t, 2t < ε < 6t and 6t < ε. The four Z2 invariants are  26  i 1 2 3 4 5 6 7 8  nx 0 0 0 1 0 1 1 1  ny 0 0 1 0 1 0 1 1  nz 0 1 0 0 1 1 0 1  ∑µ (−1)nµ 3 1 1 1 -1 -1 -1 -3  d1 (Γi ) ε − 6t ε − 2t ε − 2t ε − 2t ε + 2t ε + 2t ε + 2t ε + 6t  δi + + + + + + + +  δi< + + + + + + + −  δi + + + + − − − −  δi> + − − − − − − −  δi − − − − − − − −  Table 3.1: Time reversal invariant momenta and the parity eigenvalues for our model (−1)ν1  =  ∏  δi = δ4 δ6 δ7 δ8  ∏  δi = δ3 δ5 δ7 δ8  ∏  δi = δ2 δ5 δ6 δ8  nx =1 ny ,nz =0,1  (−1)ν2  =  ny =1 nx ,nz =0,1  (−1)ν3  =  nz =1 nx ,ny =0,1  (−1)ν0  =  (3.19)  8  ∏ δi i=1  So that    (0; 000) ε < −6t or ε > 6t     (1; 111) −6t < ε < −2t (ν0 ; ν1 ν2 ν3 ) =   (0; 111) −2t < ε < 2t     (1; 000) 2t < ε < 6t  (3.20)  Note that the special values ε = 2t, 6t are excluded, since for those values the mass of at least one Dirac point is zero - the spectrum is not gapped over the whole Brillouin zone (see (3.15). To summarize, the on-site orbital hopping amplitude ε classifies this model: for  27  Ε Ε  Ε  6t  Ε 6t  Ε 2t  2t 1;111  0;000  0;111  1;000 0;111  2t  Figure 3.1: Phase diagram for this topological insulator. Shaded regions are in the STI phase, unshaded regions are in the WTI phase. Each part of the diagram is labelled by the four Z2 invariants (ν0 ; ν1 , ν2 , ν3 ). The labels for the ε > 0 half plane are indicated, the other half plane is simply a mirror image. Here ε and t are in units of λ , and we assume λ = 0 (i.e. there is non-zero spin orbit coupling). 2t < |ε| < 6t this is a topological insulator, for 6t < |ε|, |ε| < 2t it is an ordinary  insulator, and for ε = 2t, 6t it is a semi-metal. Edge states - continuum limit  To study the surface states of the topological insulator, we study the simplest 3D example - an infinite system consisting of a half space topological insulator (z < 0) and a half space ordinary insulator (z > 0). Therefore, the plane z = 0 is the domain wall, and we would like to study the states (“edge states”) on this plane. In this section we consider a low energy approximation to the Hamiltonian which leads to the continuum limit. We are basically reducing the problem to an effective 4 × 4 Hamiltonian which captures the essential physics of the problem.  In our present model, this construction can be achieved by varying ε in space,  so that ε = ε(z). We notice that if we have 2t < ε < 6t our model will be a topological insulator and for 6t < ε it will be an ordinary insulator. If we vary ε(z) continuously, it must cross the value ε = 6t on a plane, and we will have a semi28  metal on this plane. This is the domain wall. We define m(z) ≡ ε(z) − 6t and choose a solitonic profile with asymptotic values  lim m(z) = +2t  (3.21)  z→∞  lim m(z) = −2t  z→−∞  and m(0) = 0. By looking at (3.11), we see that in the orbital subspace, only two Pauli matrices are present: τx and τz . Therefore, it is possible to rotate the Hamiltonian (in this subspace), so that (τx , τz ) → (τx , τy ), which will make the Hamiltonian block off-  diagonal. This is achieved by a rotation by −π/4 around τx , which will therefore leave τx unchanged. If we define U ≡ e−i(π/4)τx ⊗ I, then we find  H = U † HU = mk (τx ⊗ I) − 2λ ∑ sin(kµ a)(τy ⊗ σµ )  (3.22)  µ  We find the basis of operators in the rotated system is      1 0 i 0 ck↑ ck↑ + idk↑       ck↓  ck↓ + idk↓  0 1 0 i 1  1 †   = √   Ψk → U Ψk = √     2 2  i 0 1 0 dk↑  ick↑ + dk↑  0 i 0 1 dk↓ ick↓ + dk↓  (3.23)  Looking at Table 3.1 we see that the only δi that changes when we take ε past the phase change at ε = 6t, is δ1 . Therefore, we expect the interesting physics to happen near the associated TRIM, K0 = (0, 0, 0). Consider the low energy expansion of H(k) around k = K0 as done in (3.13), but now using the off-diagonal block form of the Hamiltonian. In addition, we insert the z dependence of ε, and change the font of the Hamiltonian to signify that we are in the continuum limit H (q) ≡ H (K0 + q)  vF m(z)(τ ˜ x ⊗ I) − ∑ qµ (τy ⊗ σµ )  (3.24)  µ  where we have defined vF ≡ 2λ a and m(z) ˜ ≡ m(z)/vF . Since the system is periodic  in x and y but not in z, it is natural to leave x and y in the momentum representa29  tion, but to Fourier transform z back to real space, and to define a perpendicular momentum vector q ≡ (qx , qy ) H (q, z)  vF [m(z)(τ ˜ x ⊗ I) + i∂z (τy ⊗ σz ) − qx (τy ⊗ σx ) − qy (τy ⊗ σy )]  (3.25)  We now separate the Hamiltonian into two parts H (q, z) = hz + Hq  (3.26)  where hz ≡ m(z)(τx ⊗ I) + ivF ∂z (τy ⊗ σz )  (3.27)  Hq ≡ −vF [qx (τy ⊗ σx ) + qy (τy ⊗ σy )] Note that this basis is especially convenient, since it is block off-diagonal, and in addition hz is diagonal in each block, and Hq is off-diagonal in each block. We now wish to find the wavefunctions Ψ of the full Hamiltonian H (q, z). It is obvious that the two parts of the Hamiltonian anti-commute {hz , Hq } = 0. If the  two parts had commuted, we could perform a standard separation of variables, and Ψ would be a product of the wavefunctions of each part of the Hamiltonian. We start by looking for zero energy states of hz , i.e. solutions to the equation hz ϕ(z) = 0  (3.28)  Since hz is block off-diagonal, we look at each block separately. We label the top block h+ z and write h+ z ϕ+ (z) = vF  m(z) ˜ + ∂z  0  a(z)  0  m(z) ˜ − ∂z  b(z)  =0  (3.29)  so that (m(z) ˜ + ∂z )a(z) = 0 (m(z) ˜ − ∂z )b(z) = 0  30  (3.30)  By rearranging and integrating both sides, we find a(z) =  e−  z ˜ 0 m(z  )dz  b(z) =  +  z ˜ 0 m(z  )dz  e  (3.31)  However, when the argument of the exponent is positive, the resulting wave function will be non-normalizable for the asymptotic values of m(z) given by (3.21), therefore the second equation has only the trivial solution b(z) = 0, and the single wavefunction for this block is ϕ+ (z) =  1 0  e−  z ˜ 0 m(z  )dz  e−  z ˜ 0 m(z  )dz  (3.32)  Similarly for the lower block of hz , ϕ− (z) =  0 1  (3.33)  So that the two wavefunctions of the full matrix hz are   1   0 −  ϕ+ (z) = N  0 e   0   0   0 −  ϕ− (z) = N  0 e   1  z ˜ 0 m(z  )dz  z ˜ 0 m(z  )dz  (3.34)  where N is the normalization constant which cannot be written explicitly without knowledge of the profile of the domain wall m(z). Note that the normalization constant is the same for the two wavefunctions. In the 1D system we have Hq = 0, so that H = hz . In this case ϕ± are the exact  zero modes of the complete Hamiltonian. In the 3D (and 2D) case, they are only zero modes of hz and a further step must be carried out to find the wavefunctions 31  of H . We now wish to project the full Hamiltonian onto the basis ϕ± . We define a  projection matrix P which has these vectors as its columns, and drop the exponential (and bring it back at the end). Since we are projecting onto zero states of hz , we have hz ϕ± = 0, so that in effect we are projecting only Hq .  1  0 P≡ 0  0  0     0  0  1  (3.35)  So that the projected Hamiltonian is  h± = P Hq P = −vF †  = vF   1 0 0 0  −iqx + qy   0 0 0 1  iqx + qy iqx − qy  0  qy + iqx  qy − iqx  0  −iqx − qy  = vF (qy σx − qx σy )  1 0 0 0 0 0 0 1    cq↑ + idq↑   cq↓ + idq↓  1    = √1 √  2 icq↑ + dq↑  2  icq↓ + dq↓  cq↑ + idq↑ icq↓ + dq↓  (3.37)  So that the full projected Hamiltonian, which we interpret as the effective Hamiltonian for the edge states, can be written as HB =  ∑ Ψ˜ †q h± Ψ˜ q  (3.38)  q  =  † c†q↑ − idq↑ vF † 2∑ −ic†q↓ + dq↓ q  T  0  qy + iqx  cq↑ + idq↑  qy − iqx  0  icq↓ + dq↓  32  1 0       0 0     0 0    0 1 (3.36)  We project the operator basis vector too  ˜ q = P † Ψq = Ψ    which indeed contains only one Dirac cone, as expected for the edge states of a topological insulator. The edge Hamiltonian h± is a sum of Pauli matrices so that  the energies can be written down immediately  q2x + q2y  ε± = ±vF  (3.39)  We define the phase eiα ≡ (qy + iqx )/ q2x + q2y , so that h± = ε  0  eiα  e−iα  0  (3.40)  We find the wavefunctions of h± φ+ =  1 √ 2  φ− =  1 √ 2  1  (3.41)  e−iα − eiα 1  We construct the unitary matrix U that has φ± as its columns 1 U=√ 2  1 e−iα  − eiα  (3.42)  1  which is guaranteed to diagonalize h± U † h±U =  ε 2  1 − e−iα  eiα  0  1  e−iα  eiα  1  0  e−iα  − eiα 1  =ε  1  0  0 −1  (3.43)  and the operator basis vectors ˜q = γq = U † Ψ  1 2  1  eiα  cq↑ + idq↑  − e−iα  1  icq↓ + dq↓   =  (3.44)  ei 2 e−i 2 (cq↑ + idq↑ ) + ei 2 (icq↓ + dq↓ ) α  α  α    1  2 e−i α2 − e−i α2 (cq↑ + idq↑ ) + ei α2 (icq↓ + dq↓ )  33  If we define γk ≡  γk+  then  γk− HB =  ∑  εi γqi† γqi  (3.45)  q,i=± † The γq± operators create an excitation along the domain wall. By calculating  the anticommutation relations of these operators, we find that they are ordinary Fermions † {γq+ , γq+ } = † {γq− , γq− } =  1 † † {cq↑ , c†q↑ } + {cq↓ , c†q↓ } + {dq↑ , dq↑ } + {dq↓ , dq↓ } =1 4 1 † † {cq↑ , c†q↑ } + {cq↓ , c†q↓ } + {dq↑ , dq↑ } + {dq↓ , dq↓ } =1 4 (3.46)  In conclusion, the two bound states are given by a linear combination of ϕ± with  the coefficients given by the elements of φ±   Ψ+ (q, z) =  Ψ− (q, z) =  1       0  − z m(z )dz   0 ˜ (3.47)  0 e   e−iα   − eiα    0  − z m(z 1  e 0 ˜ )dz √ − eiα ϕ+ (z) + ϕ− (z) = N    2  0  1 1 √ ϕ+ (z) + e−iα ϕ− (z) = N 2  (where the normalization constant has been redefined N →  √ 2N )  To find the energies, we note that the Hamiltonian is a sum of anticommuting matrices H = hz +Hq . Therefore, it is convenient to square it, so H 2 = h2z +Hq2 . If we now operate with H 2 on Ψ± and recall that ϕ± are zero modes of hz , i.e. hz ϕ± = 0, we see that  E± = ε± = ±vF 34  q2x + q2y  (3.48)  Note that these bound states are written in the block off-diagonal basis given by (3.23)   ck↑ + idk↑    c + id 1  k↓ k↓  Ψk = √   2 ick↑ + dk↑   ick↓ + dk↓  (3.49)  For the special case of a step function domain wall (zero width wall) m(z) = 2t[Θ(z) − Θ(−z)]  (3.50)  So that Ψ± (q, z) ∝ e−  z ˜ 0 m(z  )dz  − v2t |z|  =e  (3.51)  F  Giving an exponential decay, as expected. Recall that in the lattice model vF ≡ 2λ a.  As an example of a wall with a width, we choose m(z) = 2t tanh  z ξ  (3.52)  where ξ is a length scale for the width of the domain wall. This gives Ψ± (q, z) ∝ e  z tanh zξ F 0  − v2t  dz  =e  − v2t ξ ln cosh(z ) F  z/ξ 0  z = cosh ξ  − v2t ξ F  (3.53)  Another similar function is m(z) = 2t arctan(z/ξ ). We can also check the limit of an ordinary insulator or topological insulator in all space, i.e. m(z) = ±2t for all z, which gives  Ψ± (q, z) ∝ e−  z ˜ 0 m(z  )dz  =e  ∓ v2t z F  (3.54)  which is non-normalizable, since for z → ∓∞ the wave function will blow up. This  is consistent, since there is no reason for a bound state to exist in these cases. Another limiting case is m(z) = 0 for all z, i.e. the spectrum is ungapped over the  35  whole system. Phrased differently, we expand the boundary to include all space there is nothing but boundary. The exponent above gives the identity in this limit the “bound states”’ are no longer localized around z = 0 but rather extend over all z. Similar domain walls in Graphene were studied in [56] in a different context: the domain walls were calculated for Graphene gapped by a staggered on-site energy, which is not a topological insulator, the system is 2D, no spin is involved, and the 4 × 4 Hamiltonian is formed by unifying the two Dirac points (each is represented by a 2 × 2 matrix) in one matrix, whereas in our case it is formed  by considering only the Dirac point (represented by a 4 × 4 matrix) whose mass  changes in the phase change.  3.3  Explicit calculation of θ  We start from the full Hamiltonian in momentum space (3.11) H(k) = mk (τx ⊗ I) − 2λ ∑ sin(kµ a)(τz ⊗ σµ )  (3.55)  µ  and expand around each of the Dirac points K, so that k = K + q. We find H (q) = mK (τx ⊗ I) − 2λ ∑(−1)nµ qµ (τz ⊗ σµ )  (3.56)  µ  with (3.15)  mK = ε − 2t ∑ cos(Kµ a) = ε − 2t ∑(−1)nµ = µ  µ     ε − 6t     ε − 2t   ε + 2t     ε + 6t  (3.57)  where the indices nµ = 0, 1 index the Dirac points K=  π (nx , ny , nz ) a  36  (3.58)  We define the momentum space non-abelian gauge field, in this case a set of three matrices A j of dimensions 4 × 4 Aj  µν  = −i uµ |∂ j |uν  (3.59)  where ∂ j ≡ ∂ /∂ k j and |uν are the momentum space wavefunctions. Note that this is the 3D analogue of the 2D Berry connection scalar, defined as A ≡ −i u|∇k |u . Then we calculate θ (see [30]) θ=  1 4π  BZ  2 d 3 k ε i jk Tr Ai ∂ j Ak + iAi A j Ak 3  (3.60)  where the trace is performed only over the occupied bands. This formula bears some resemblance to the Thouless-Kohmoto-Nightingale-Nijs (TKNN) formula for 2D systems. To evaluate this integral, we start by finding the wavefunctions of the linearized Hamiltonian. We then calculate the three matrices A j and substitute into this equation. For each Dirac point K, we find θK =  Λ 0  d3q  2mK v3F = tan−1 (m2K + q2 v2F )2  vF Λ mK vF Λ − 2 mK mK + v2F Λ2  (3.61)  The integrand is sharply peaked around q = 0. We take the limit Λ → ∞ and find  that the contribution of each Dirac point is θK =  π sign(mK ) 2  (3.62)  and θ = ∑ θK = K  π sign(mK ) 2∑ K  (3.63)  where the result is understood to be mod 2π, so that 0 ≤ θ < 2π. Since the lin-  earized Hamiltonian is of the standard form of a Dirac Hamiltonian, we expect this result to hold for any similar model. In this particular model, mK is given by (3.57),  37  so that θ = ∑ θK = K  π sign ε − 2t ∑(−1)nµ 2∑ µ K  (3.64)  We now vary ε and find that  0 θ= π  |ε| < 2t or |ε| > 6t 2t < |ε| < 6t  which verifies our earlier findings (3.20) if we identify ν0 = θ /π.  38  (3.65)  Chapter 4  The Witten Effect in a Topological Insulator 4.1  Overview  It was noted a long time ago that a term of the form θ (e2 /2πh)B · E may be added to the standard Maxwell Lagrangian without modifying the familiar laws of elec-  tricity and magnetism. θ is known to particle physicists as the ‘axion’ field and whether or not it has a nonzero expectation value in vacuum remains a fundamental open question of the Standard Model. A key manifestation of the axion term is the Witten effect: a unit magnetic monopole placed inside a medium with θ = 0 is predicted to bind a (generally fractional) electric charge −e(θ /2π + n) with n  integer. Here we conduct a first test of the Witten effect based on the recently established fact that the axion term with θ = π emerges naturally in the description of the electromagnetic response of a new class of crystalline solids called topological insulators – materials distinguished by strong spin-orbit coupling and non-trivial band structures. Using a simple physical model for a topological insulator we demonstrate the existence of a fractional charge bound to a monopole by an explicit numerical calculation. We also propose a scheme for generating an ‘artificial’ magnetic monopole in a topological insulator film, that may be used to facilitate the first experimental test of Witten’s prediction.  39  4.2  Axions  The idea of the axion was introduced in 1977 by Peccei and Quinn [57, 58] as a means to solve what is known as the ‘strong CP problem’ in the physics of strong interactions. The strong CP problem, the details of which are quite subtle, has to do with the vacuum structure of Quantum Chromodynamics. In simple physical terms it can be stated as a question: Why is the electric dipole moment of the neutron (currently unobserved) so small? The Standard Model predicts a value for the neutron dipole moment |dn | ∼ 10−16 θ e cm, with θ of order unity, that should  be readily measurable. Peccei-Quinn’s solution promotes θ to a dynamical field  describing a new elementary particle, the axion, whose vacuum expectation value has relaxed to a very small value, explaining the smallness of |dn |. The actual value  of θ , and the validity of the Peccei-Quinn solution and its variants [59, 60] remain open questions of considerable importance to fundamental physics. The axion is also believed to be a viable candidate for the elusive dark matter that comprises the majority of matter in our universe [61] and is subject to active experimental searches [62, 63]. In a remarkable development axion electrodynamics has recently emerged as a key tool in the description of STIs. It has been realized [29, 30], that the electromagnetic response of a STI is characterized by the axion term ∆Laxion = θ  e2 B·E 2πh  (4.1)  with θ = π, the only non-zero value permitted by the time-reversal symmetry. When the time-reversal symmetry is broken, e.g. in a crystal showing weak magnetism, θ can acquire an arbitrary value. Fluctuations in the magnetic order parameter then act as a dynamical axion field and can be thought of as emergent axion particles [43]. Thus, aside from possible practical applications, crystalline solids with topologically non-trivial band structures have the potential to provide tabletop laboratories for the testing and exploration of fundamental physical paradigms. A fundamental property of the axion medium is the Witten effect [32]: in the quantum theory, a magnetic monopole of unit strength (i.e. projecting magnetic flux Φ0 = hc/e) immersed in an axion medium must carry electric charge  40  −e(θ /2π + n) with n integer. This effect, although theoretically well established,  has never been experimentally tested because until now both a suitable axion medium and the means to produce a magnetic monopole have been lacking. In this section we demonstrate how the connection between the axion response [29, 30] and strong topological insulators[2, 34, 35, 37, 64–66] may serve to overcome both obstacles. We remark that a 1D realization of the Witten effect in antiferromagnetic spin chains was proposed a long time ago [67]. Here we furnish the first concrete physical example of the Witten effect in 3D by modeling a STI with a magnetic monopole inserted in its bulk. We show that the monopole binds a fractional charge ±e/2 consistent with Witten’s prediction [32]. We then discuss possible ways to overcome the second obstacle by creating an emergent magnetic monopole in a  topological insulator. This can be achieved by exploiting the degrees of freedom associated with a vortex in the exciton condensate that may emerge in a thin film topological insulator under external bias [68]. We conclude that the prospects for experimental verification of the Witten effect in a tabletop experiment using a STI appear promising.  4.3  Classical view of the Witten effect I  Although the Witten effect is quantum-mechanical in nature its essence can be understood by studying the classical Maxwell equations modified in the presence of the axion term ∆Laxion = θ  e2 B·E 2πh  (4.2)  We recall Maxwell’s equations with magnetic monopoles (and currents) and θ = 0 (see page 273 in Jackson [69], for example) ∇ · E = ρe  ∇ · B = ρm  ∂B − Jm ∂t ∂E ∇×B = + Je ∂t  ∇×E = −  41  (4.3)  and compare with Maxwell’s equations in the presence of a θ field but with no magnetic charges and currents (ρm = 0, Jm = 0), as found by Wilczek [70] ∇ · E = ρe −  e ∇θ · B 4π 2  (4.4)  ∇·B = 0  ∂B ∂t e ∂E + Je + 2 ∇×B = ∂t 4π ∇×E = −  ∂θ B + ∇θ × E ∂t  (where we are working in units where the flux quantum is dimensionless hc/e = 1, so that the fine structure constant α = e, so that in Wilczek’s notation κ = α/4π 2 = e/4π 2 ). We infer that Maxwell’s equations in the presence of a non-zero θ and magnetic charges and currents should be ∇ · E = ρe −  e ∇θ · B 4π 2  (4.5)  ∇ · B = ρm  ∂B − Jm ∂t ∂E e ∇×B = + Je + 2 ∂t 4π ∇×E = −  ∂θ B + ∇θ × E ∂t  We observe that for uniform, constant θ these equations revert to the familiar Maxwell’s equations, consistent with the notion that ∆Laxion can be written as a total derivative in this case. An important related property [31] is the periodicity under θ → θ + 2πn of the axion action Saxion , implying that θ can be chosen from the interval [0, 2π).  The axion term (4.2) is invariant under the duality transformation that swaps magnetic and electric fields (or in the more general case, mixes the fields via a rotation), but the modified Maxwell equations for axion electrodynamics (4.4) and (4.5) are not invariant under this transformation. However, it appears to be well known that adding an axion term to the Lagrangian breaks the duality transformation and several attempts have been made to resolve this issue in the past [71–75]. 42  However, at present, it does not appear that there is a consensus on how to solve this problem. Now consider a monopole of strength g, ∇ · B = 4πgδ (3) (r), placed at the ori-  gin, in a medium initially characterized by θ = 0. We wish to understand what happens when we turn on θ as a function of time (but keep it uniform in space). To this end we set ∇θ = 0, Je = Jm = 0 (no currents in vacuum) and take the divergence of the last equation to obtain ∇·  ∂E e ∂θ + 2 ∇·B = 0 ∂t 4π ∂t  (4.6)  We see that an electric field is generated in this process. Integrating Eq. (4.6) over space and time, we find that this field can be thought of as originating from a point electric charge Q located at the origin with magnitude Q = −ge  ∆θ π  (4.7)  where ∆θ is the net change in θ and we assumed that there was no initial electrical charge bound to the monopole. In a topological insulator ∆θ = π, thus one expects a unit magnetic monopole (g = 1/2) to bind fractional charge Q = −e  1 +n 2  (4.8)  The integer n accounts for the possibility of binding extra electrons, which can always occur – only the fractional part of Q is non-trivial. The question arises, how does the inclusion of the dielectric constant for a specific TI affect this calculation. We write down Maxwell’s equations for this  43  case (see Eq. 104 of [29] , for example) ∇ · D = ρe −  e ∇θ · B 4π 2  (4.9)  ∇ · B = ρm  ∂B − Jm ∂t ∂D e ∇×H = + Je + 2 ∂t 4π  ∇×E = −  ∂θ B + ∇θ × E ∂t  As before, we take the divergence of the last equation, finding ∇·  ∂D e ∂θ + 2 ∇·B = 0 ∂t 4π ∂t  (4.10)  and upon integration in space and time, we find the same result as before, Eq. (4.8). We conclude that the Witten effect holds for dielectric materials. One may ask, how does the localization length of the charge localized at the monopole depend on the dielectric constant. To derive a simple result, we assume that the charge density around the monopole is uniform. Then one can derive the energy required to localize charge Q within radius R, which is U = 3Q2 /20πεR. If we then assume that the work done by the system to localize the charge is the same, regardless of the dielectric constant, we find R=  ε0 R0 ε  (4.11)  where R0 and R are the localization length in vacuum and in the dielectric material (respectively). Since ε > ε0 , we find that under these assumptions, the localization length is shorter in a dielectric. This can, perhaps, be understood as a result of the higher capacitance of a dielectric under potential - the dielectric can store more charge within a given volume. Since the charge is the same in both cases, it is possible to store it in a smaller volume in the dielectric. This result could be made more realistic, for example by using an exponential density profile, but we leave this for further studies.  44  4.4  Classical view of the Witten effect II  There is another way to use the classical Maxwell equations to derive a similar result. Consider a magnetic point charge (a monopole) of charge g, located at the origin ρm = 4πgδ (3) (r)  (4.12)  with no electric or magnetic currents (Je = Jm = 0), and no electric charges (ρe = 0) embedded within a topologically insulating sphere of radius R. Due to a divergence of integrals at 0, we are forced to consider the case of a shell with inner radius ε and outer radius R, and at the end take ε → 0, giving the same result. Then θ (r) = π[Θ(R − r) − Θ(ε − r)]  (4.13)  where Θ is the Heaviside function. We note that ∇θ (r) =  dθ dΘ(R − r) dΘ(ε − r) rˆ = π − dr r r = −π[δ (r − R) − δ (r − ε)]ˆr  (4.14)  where we have used the fact that the derivative of the Heaviside function is the Dirac delta function dΘ(r)/dr = δ (r), and the Dirac delta function is symmetric δ (r) = δ (−r). We use Gauss’ law B · ds =  ρm d 3 r  (4.15)  to find the magnetic field B(r) =  g rˆ r2  (4.16)  So that the first equation of (4.5) gives the effective charge density ∇·E = −  e δ (r − R) δ (r − ε) ∇θ · B = ge − ≡ ρe f (r) 2 4π 4πR2 4πε 2  45  (4.17)  By integrating over a narrow region around each sphere we find the charge on the outer and inner surfaces R+δ  Qo =  R−δ ε+δ  Qi =  ε−δ  Therefore, for g =  1 2  ρe f (r)4πr2 dr = ge  (4.18)  ρe f (r)4πr2 dr = −ge  (4.19)  and in the limit of ε → 0 we find fractional charge −e/2  bound to the monopole (as in the previous section), and +e/2 bound to the outside surface of the topological insulator.  4.5  Adding a single monopole  We first consider a single magnetic monopole, in other words a magnetic point charge. Therefore, we expect ∇ · Bg = 4πgδ (3) (r)  (4.20)  This equation is fulfilled by a Coulombic magnetic field Bg = g  r g = 2 rˆ 3 r r  (4.21)  However, if the magnetic field can be written using a vector potential B = ∇ × A, we immediately reach a contradiction  ∇ · B = ∇ · (∇ × A) = 0  (4.22)  since the divergence of a curl is identically zero, but ∇ · Bg = 0. Therefore, we write the magnetic field as a sum of the magnetic field of the magnetic monopole and an additional singular part B ≡ Bg + Bs =  g rˆ − 4πgθ (z)δ (x)δ (y)ˆz r2  (4.23)  The singular part Bs can be thought of as an infinitely long infinitesimal solenoid leading flux lines in the −ˆz direction. It is often referred to as a “Dirac string”. 46  Equation (4.22) is now fulfilled since ∇ · Bs = −4πgδ (x)δ (y)  dθ (z) = −4πgδ (3) (r) dz  (4.24)  where we have used the fact that the derivative of the heaviside function is the delta function. We see that by adding the Dirac string we have gotten rid of the contradiction, and therefore a vector potential can be written for this magnetic field. We expect the total flux through a sphere centered on the magnetic monopole to be zero, and indeed Φ=  rˆ · ds − 4π r2  B · ds = g 1 R2  = g  ds − 4π  θ (z)δ (x)δ (y)ˆz · ds  (4.25)  θ (z)δ (x)δ (y)ds = 0  Note that the reason we want to write the magnetic field of the monopole using a vector potential, is that if we can do this, we can conveniently add the magnetic field to the Hamiltonian via Peierls phase factors. The magnetic field of the monopole can be rewritten as a gradient Bg =  g 1 rˆ = −g∇ 2 r r  (4.26)  since in general, for a function of the radial coordinate only, we have ∇ f (r) =  df rˆ dr  (4.27)  Therefore, ∇ × Bg = −g∇ × ∇  1 r  =0  (4.28)  since the curl of a gradient is identically zero. If we now take the curl of the full magnetic field, we find that only the singular part remains ∇ × B = ∇ × (Bg + Bs ) = ∇ × Bs  47  (4.29)  On the other hand ∇ × B = ∇ × (∇ × A)  (4.30)  = ∇(∇ · A) − ∇ A = −∇ A 2  2  where we have chosen the Coulomb gauge ∇ · A = 0. Finally, we find that the vector potential is related only to the singular part of the magnetic field ∇ × Bs = −∇2 A  (4.31)  For a periodic system, it is convenient to Fourier transform this expression, and after rearranging we find (for k = 0) Ak = i  k × Bsk k2  (4.32)  Adding a monopole to the Hamiltonian in practice involves calculating the Peierls factors for a given configuration. We cover this calculation for a lattice of monopoles, a single monopole and a planar monopole in Appendices A.3.1, A.3.2, and A.3.3 (respectively).  4.6  Charge bound to a monopole  We now consider a magnetic monopole in the interior of a STI. We model this situation by the Hamiltonian H presented in Section 3.2 with a monopole positioned at the center of a cubic unit cell. The magnetic field of the monopole couples to both the electron charge and the electron spin through the orbital and Zeeman couplings, respectively. The form of the orbital coupling is dictated by gauge invariance and is thus universal; in our lattice model it is implemented by the Peierls substitution, which attaches factors eiθi j to all hopping terms connecting sites i and j. Here θi j = (2π/Φ0 )  j i  A · d l and A is the magnetic vector potential. The Zeeman cou-  pling is of the form −gµB B · S/¯h where µB = e¯h/2me c is the Bohr magneton and S  denotes the electron spin. For free electrons g is close to 2 but in solids the effective g can be substantially larger. The Zeeman coupling thus leads to an additional  48  term in the Hamiltonian, HZ = −gµB  1 B j · Ψ†j σ Ψ j , 2∑ j  (4.33)  where B j is the magnetic field at site j of the lattice. This term is non-universal and its importance will depend on the ratio of gµB |B| to the other relevant energy  scales in the Hamiltonian set by λ , ε and t.  We solve the Hamiltonian H = HSO + Hcd + HZ in a cube containing L3 sites by exact numerical diagonalization. The monopole is placed inside the central unit cell (at the origin), so that the magnetic field of the monopole is B = (Φ0 /4πr2 )ˆr. We choose a gauge in which the system retains the four-fold rotational symmetry around the z axis [76], A = −Φ0 (1 + cos θ )∇ϕ, with (θ , ϕ) the spherical angles. Exploiting this symmetry we are able to simulate system sizes up to L = 20, which  requires diagonalizing a complex valued Hermitian matrix of size 41 (4 × 203 ) =  8, 000. In order to calculate the charge density at half filling we require knowledge  ln (δρ + c)  of all the occupied eigenstates.  10 8 6 4 2 0 ï2ï4 ï6ï8 10 ï10 ï10ï8 ï6 ï4 ï2 0 2 4 6 8 y  x  Figure 4.1: Charge density δ ρ of the three closest layers below the monopole, for g = 0. We use our model TI on the cube-shaped lattice with 203 sites with a unit monopole at its center, with parameters t = λ , ε = 4t, leading to a bulk gap ∆ = 4t.  49  g=0 g=2 g=10  0.5  δQ(r)  0.4  0.3  0.2  0.1  0.0  0  2  4  6  8  r  10  12  14  16  Figure 4.2: The excess charge δ Q(r) (in units of e) for different Zeeman coupling g. The knee feature seen at r = 10 corresponds to the radius at which the sphere used to calculate δ Q(r) first touches the system boundary. The parameters are as in Fig. 4.1. We diagonalize the Hamiltonian, once with the magnetic monopole and once without, obtaining charge densities ρ1 and ρ0 , respectively. The monopole-induced charge density δ ρ = ρ0 − ρ1 is plotted in Fig. 4.1. To determine the total charge  bound to the monopole we calculate the excess accumulated charge in a sphere of radius r centered on the monopole, δ Q(r) = ∑|ri |<r δ ρ(ri ). We find (Fig. 4.2) that  it saturates at −e/2 to within 4 significant digits, comparable to the accuracy of our numerics. For g = 0 we find two localized zero modes, one at the monopole and  one on the surface. Fractional charge bound to the monopole can be understood in this case by appealing to the standard arguments [77, 78] developed originally to describe charge fractionalization in polyacetylene [79]. Briefly, when a topological defect (such as a domain wall in polyacetylene) produces a localized zero mode inside the gap in a particle-hole symmetric system, one can show that the spectral weight of the state contains equal contributions from the valence and the conduction bands. Thus, the valence band shows a net deficit of half a state in the vicinity 50  ï1.6  g=2 g=6 g=10  log |δQ(r) − δQ0 (r)|  ï1.8 ï2.0 ï2.2 ï2.4 ï2.6 ï2.8 ï3.0 ï3.2  ï3.4 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90  log r  Figure 4.3: Log-log plot of δ Qg − δ Q0 showing the power-law approach ∼ r−α of the accumulated charge to its asymptotic value of 1/2. The least-square fit yields exponents α = 2.85, 3.04, 2.79 for g = 2, 6, 10, respectively. We attribute the deviations of the numerically determined exponent α from the expected value of 3 to the finite size effect. The parameters are as in Fig. 4.1. of the defect. This translates into the defect carrying fractional charge ±e/2, the  sign depending on whether the zero mode is empty or occupied.  Like in polyacetylene we find the saturation of charge to be exponential ∼  exp (−r/ξ ), where ξ ∼ 1/∆ and ∆ is the bulk gap.  When g > 0, the Zeeman coupling causes changes in the charge distribution  near the monopole but the total accumulated charge remains quantized at −e/2. In  this case there are no exact zero modes in the spectrum and δ Q(r) approaches e/2 as a power law with the exponent close to −3 (Fig. 4.3).  The power-law dependence can be understood as follows. The Zeeman term  acts as an additional time-reversal breaking field which modifies the value of axion θ away from π close to the monopole. This causes a non-vanishing ∇θ and thus, according to Eq. (4.5), an additional contribution to the effective charge density.  51  The simplest assumption, δ θ ∼ B2 , gives δ ρ ∼ ∇θ · B ∼ r−7 and δ Q ∼ r−4 , a power  law but with the exponent not quite in agreement with our numerical simulation. On further reflection one realizes that in our model δ θ cannot be proportional to B2 but rather must be proportional to its gradients. This is because in the presence of a uniform Zeeman term the system remains inversion-symmetric. Inversion symmetry dictates a quantized value of θ = 0, π even when T is explicitly broken [29, 30]. Thus, non-vanishing δ θ requires a spatially varying Zeeman field. The simplest assumption satisfying these requirements is ∇θ ∼ ∇2 B. In the vicinity of the monopole one finds δ ρ ∼ ∇θ · B ∼ r−6 and δ Q ∼ r−3 in agreement with our  numerical results.  We note that the above considerations are based on the effective axion action (4.1) and apply on length scales large compared to ξ . The power law tail in the fractional charge distribution for g > 0 appears on top of a short-lengthscale structure with a roughly exponential profile that is controlled by the properties of the microscopic Hamiltonian and is thus non-universal. At the intermediate lengthscales the interplay of the two contributions can give rise to interesting structures such as the peak in δ Q(r) at r  2.5 seen in Fig. 4.2 for g = 10.  By the same method described above we have investigated the spin density induced by the monopole. We find that there is no net spin S attached to the monopole. Thus, in addition to charge fractionalization, a magnetic monopole inserted in a STI constitutes an example of spin-charge separation in three spatial dimensions. This is perhaps not surprising in view of the fact that spin-orbit coupling present in the Hamiltonian breaks the SU(2) spin symmetry and, as a result, electron spin is not a good quantum number in the model describing our system.  4.7  Proposal for experimental realization  Although there is no known theoretical principle that prohibits the existence of fundamental magnetic monopoles in nature [76], none have been observed to date despite extensive searches [80]. This null observation has led to a consensus that fundamental monopoles either do not exist for some heretofore unknown reason or they are very rare in our part of the universe. In either case the observed absence of fundamental monopoles poses a challenge to the idea of experimental verification  52  of the Witten effect using a STI. At best, one could conceive of a new scheme for possible detection of magnetic monopoles using a STI in the role of a sensor if a convenient way to detect the fractional charge could be found. A much more promising avenue for the verification of the Witten effect is suggested by exploiting emergent instead of fundamental monopoles. A classic example of such an emergent behavior in a crystalline solid is the 2007 theoretical prediction [81] and the subsequent experimental observation [82–84] of monopoles in frustrated magnetic systems called ‘spin ice’, realized in certain magnetic pyrochlore compounds such as Dy2 Ti2 O7 or Ho2 Ti2 O7 . Magnetic monopoles in these systems arise as elementary excitations above the collective ground state of spins and the monopole-like magnetic field configuration originates from the magnetic moments of the constituent spins. In principle, the emergent monopoles in the spin ice could be used to test the Witten effect if a compound that is simultaneously a STI and a spin ice could be identified. Unfortunately no such material is known at present although we note that STI behavior has been theoretically predicted to occur in crystals with the same pyrochlore structure [85, 86] that underlies the spin ice behavior. It is thus possible that a suitable material will be discovered in the future. Here we focus on a different type of emergent magnetic monopole that can arise in a thin film STI placed in a uniform external electric field. The basic idea and the feasibility of its experimental realization have been discussed in Ref. [68]. Following that work we envision the simplest STI with just one gapless Dirac state per surface and the chemical potential initially tuned to the neutral point. When a strong enough electric field is applied perpendicular to the plane of the film the chemical potential undergoes a shift that is opposite in the two surfaces. This creates a small electron Fermi surface in one surface and a small hole Fermi surface in the other. The essence of the proposal [68] lies in the observation that an arbitrarily weak Coulomb interaction between the surface states produces an exciton condensate, which may be viewed as a coherent fluid of electron-hole pairs drawn from the opposite surfaces. Such an exciton condensate is characterized by a complex scalar order parameter Φ, which can fluctuate in space and time. In particular Φ = Φ0 eiχ can contain vortices – point-like topological defects with the phase χ winding by ±2π around a vortex. It has been pointed out in Ref. [68] that to electrons in a STI 53  such a vortex is indistinguishable from a ‘planar monopole’, i.e. a monopole with magnetic field radiating in the plane of the surface. 0.04 0.03  2 − ρ1  0.02 0.01 0.00  ï0.01 7  5  3  1  ï1  y  ï3  ï5  ï7 ï7  ï5  ï3  ï1  1  3  5  7  x  Figure 4.4: A cubic sample of a TI including disorder with a planar unit monopole at its center, size L = 14 and parameters as in Fig. 4.1. The plot shows the charge density 2 − ρ1 for the layer just below the planar monopole for weak disorder µ = 0.05∆. A planar monopole can be viewed as an adiabatic deformation of an ordinary monopole achieved by flattening the field lines in a cylindrically symmetric fashion. One expects that the total charge bound to the monopole via the Witten effect should be insensitive to such an adiabatic deformation and therefore a vortex in the exciton condensate should bind fractional charge ±e/2. This indeed has been ar-  gued to happen in Ref. [68] based on the Dirac equation describing the low-energy physics of the surface states in the presence of the exciton condensate. Here, taking a more general point of view, we establish the existence of the fractional charge in such a condensate by studying a planar monopole embedded inside a STI. Our calculation below does not rely on the low-energy approximation for the surface states and is insensitive to the detailed microscopic structure of the condensate. Rather, it exploits only the most fundamental property of the STI given by its nontrivial  54  7  7  0.04  5  0.06  5  0.03  0.04  3  3 0.02  0.02  y  1  y  1 0.01  ï1 ï3  0.00 ï1 ï0.02  ï3  0.00  ï5  ï0.04  ï5 ï0.01  ï7 ï7  ï5  ï3  ï1  x  1  3  5  ï7 ï7  7  ï5  ï3  ï1  (a)  x  1  3  5  7  ï0.06  (b)  7  7  0.04  5  0.06  5  0.03  0.04  3  3 0.02  0.02  y  1  y  1 0.01  ï1 ï3  0.00 ï1 ï0.02  ï3  0.00  ï5  ï0.04  ï5 ï0.01  ï7 ï7  ï5  ï3  ï1  x  1  3  5  ï7 ï7  7  (c)  ï5  ï3  ï1  x  1  3  5  7  ï0.06  (d)  Figure 4.5: A cubic sample of a TI including disorder with a planar unit monopole at its center, size L = 14 and parameters as in Fig. 4.1. The plots show the charge density 2 − ρ1 for the layer just below the planar monopole for: (a) weak disorder µ = 0.05∆, and (b) stronger disorder µ = 0.20∆. Panels (c) and (d) show the difference in charge density δ ρ = ρ0 − ρ1 for µ/∆ = 0.05, 0.20, respectively, for the same layer. axion response. In general, the fractional charge is expected to be robust against weak disorder that does not break TRI. Such a disorder will be present in a real sample and we  55  0.5 µ=0.00 µ=0.02 µ=0.05 µ=0.20  δQ(r)  0.4 0.500  0.3  0.498  0.2  0.496 0.1 0.494 0.0  0  4  2  5 4  6  r  6  8  10  12  Figure 4.6: The excess charge δ Q(r) (in units of e) for different values of disorder strength µ. The inset shows a close up of the saturation. At this scale a small deviation from the expected asymptotic value 1/2 that increases with the disorder strength becomes visible. We attribute this deviation to the finite-size effect in our numerical calculation. This identification is supported by the fact that the deviations grow more pronounced for smaller system sizes and close to the surface. Also, it is consistent with the notion that the bound charge is localized on the length scale ξ ∼ 1/∆ which increases as the disorder reduces the spectral gap. model it here by adding a term † − † HD = ∑ µ + j Ψ j Ψ j + ∑ µ j Ψ j τz Ψ j , j  (4.34)  j  to the Hamiltonian. The first term represents a parity-preserving on-site disorder (independent of the orbital), while the second term is a parity breaking disorder. As before we solve the Hamiltonian H = HSO + Hcd + HD , in a cube containing L3  sites, by exact numerical diagonalization. The planar monopole projects an  effective magnetic field Beff = (Φ0 /2πr)δ (z)ˆr (in cylindrical coordinates) and the 56  vector potential can be chosen as A = (Φ0 /2π)ϕδ (z)ˆz. The effective field does not couple to electron spin [68], so there is no Zeeman term in this case. The disorder coefficients µ ± j are chosen from a Gaussian distribution with standard deviation µ. Note that the disorder breaks the four-fold rotational symmetry of the system, so we cannot exploit this symmetry in this case to efficiently diagonalize the Hamiltonian. Consequently we are limited to system sizes up to L = 14. For weak disorder ∆ the charge bound to the planar monopole remains −e/2 (see Fig. 4.5), and  µ  for strong disorder µ  ∆ the charge bound is zero. Remarkably, even for fairly  significant disorder (such that it generates charge density fluctuations comparable to the charge density induced by the monopole) the difference in charge density δ ρ shown in Fig. 4.5d is only weakly affected. In the framework of the current proposal the key ingredient required to produce a monopole-like configuration is the exciton condensate. As explained in Ref. [68] it is difficult to reliably estimate the critical temperature TEC for the formation of the exciton condensate, but under optimal conditions it should be higher than it is in bilayer graphene, where the occurrence of this effect is hotly debated. Once the exciton condensate is formed, vortices can be nucleated by applying an inplane magnetic field. Since the exciton condensate is itself insulating, the main conduction channel in this situation will be through vortices, each carrying −e/2  charge. Fractional charge of the carriers then can be detected using established techniques [87, 88].  4.8  Conclusions  Predicted more than 30 years ago in the context of high-energy physics, but never before observed in a real or numerical experiment, the Witten effect is realized in a strong topological insulator. A unit magnetic monopole inserted in a model STI binds electric charge −e/2 in accordance with the prediction [32] and fur-  nishes a rare example of charge fractionalization and spin-charge separation in 3 spatial dimensions. In the special case when the underlying system possesses  particle-hole symmetry and when the Zeeman term (4.33) can be neglected, the appearance of fractional charge follows from the same ‘zero-mode’ arguments that underlie charge fractionalization in a one-dimensional system of fermions coupled  57  to a scalar field with a soliton profile [77, 78] as realized in dimerized polyacetylene [79]. In the more general case when the Zeeman term or weak disorder are present, there exist no exact zero modes in the spectrum of electrons yet the fractional charge remains precisely quantized. This reflects the more subtle topological order that underlies the axion response of a STI, which is robust against any weak perturbation that respects time reversal symmetry [29, 30, 64, 65]. An interesting open question is the fate of the Witten effect in the presence of magnetic disorder. Experimentally, this can be implemented by adding a small concentration of magnetic ions (such as Fe or Mn) into a topological insulator. Although at the microscopic level θ is no longer quantized in the presence of time reversal symmetry (TRS)-breaking perturbations, we expect its effective value, relevant to the physics at long lengthscales, to remain pinned at π as long as the magnetic moments stay disordered. This is because the net magnetic moment in a macroscopic region containing many impurities will effectively vanish. Qualitatively, this suggests that the Witten effect may survive inclusion of a moderate concentration of magnetic dopants with randomly oriented moments. At low temperatures moments may order ferromagnetically [89]. In this case both TRS and the inversion symmetry are broken (the latter due to the random position of magnetic dopants) and the effective θ can acquire an arbitrary value. In this situation we expect a monopole to still bind fractional charge according to Eq. (4.8) but we leave a detailed study of this case to future investigation. Can the Witten effect be observed experimentally in the near future? We believe that the answer is affirmative. One essential ingredient, the axion medium, is now widely available in any of the recently discovered STIs [35, 37, 66]. If an emergent monopole can be realized, exploiting the proposed exciton condensate [68], the spin ice-type physics [81–84], or by some other means, then the experimental challenge is reduced to designing a suitable method for the detection of fractional charge bound to the monopole. The fractional charge of elementary excitations in fractional quantum Hall fluids has been previously detected [87, 88] and it should be possible to adapt these methods to topological insulators. In this way, studies of crystalline quantum matter with non-trivial topological properties could help settle one of the enduring challenges of fundamental physics and provide new insights into the behaviour of electrons placed in unusual situations. 58  Chapter 5  Wormhole Effect in a Strong Topological Insulator 5.1  Overview  An infinitely thin solenoid carrying magnetic flux Φ (a ‘Dirac string’) inserted into an ordinary band insulator has no significant effect on the spectrum of electrons. In a strong topological insulator, remarkably, such a solenoid carries protected gapless one-dimensional fermionic modes when Φ = hc/2e. These modes are spin-filtered and represent a distinct bulk manifestation of the topologically non-trivial insulator. We establish this ‘wormhole’ effect by both general qualitative considerations and by numerical calculations within a minimal lattice model. We also discuss the possibility of experimental observation of a closely related effect in artificially engineered nanostructures.  5.2  Laughlin’s argument  Surface electrons in a STI [2, 34, 35, 37, 64–66, 90] form a gapless helical liquid, protected by TRS through the topological invariants that characterize the bulk band structure. When TRS is broken, which may be accomplished by coating the surface with a ferromagnetic film, the helical liquid transforms into an exotic insulating  59  state characterized by a precisely quantized Hall conductivity σxy = n +  1 2  e2 , h  (5.1)  with n integer. This result follows from the microscopic theory of the surface state [2] and also from the effective electromagnetic action describing the bulk of a topological insulator, which contains the axion term [29, 30]. Although it might not be possible to measure this ‘fractional’ quantum Hall effect in a transport experiment [91], Eq. (5.1) is predicted to have observable physical consequences, such as the low-frequency Faraday rotation [29] and the image magnetic monopole effect [44]. It is instructive to apply Laughlin’s flux insertion argument [92] to the STI surface described by Eq. (5.1). This argument was devised to establish the fractional charge of quasiparticles in fractional quantum Hall liquid (FQHL) [93] and involves the adiabatic insertion of an infinitely thin solenoid carrying magnetic flux Φ(t) into the system, as illustrated in Fig. 5.1a. As the flux is ramped up from 0 to Φ0 = hc/e, a circumferential electric field is generated in accord with Faraday’s law ∇ × E = −(1/c)(∂ B/∂t). This induces a Hall current on the STI surface j = σxy (E × zˆ) which brings electric charge δ Q = σxy  Φ0 1 = n+ e c 2  (5.2)  to the solenoid. Since the flux tube carrying a full flux quantum Φ0 can be removed from the electronic Hamiltonian by a gauge transformation, one concludes, as in FQHL, that an excitation with fractional charge (5.2) must exist. This finding stands in contradiction to the well established microscopic theory of these surface states given by an odd number of massive Dirac Hamiltonians [2]. Elementary excitations of a massive Dirac Hamiltonian are particle-hole pairs which are charge neutral. Yet, this same Dirac Hamiltonian exhibits Hall conductivity (5.1), which, through Laughlin’s argument outlined above, implies fractionally charged quasiparticles. The resolution to this paradox comes from the realization that the quantum Hall state realized on the surface of a STI is inextricably linked to the bulk of the STI. 60  a)  b)  \  E  FM  \  STI  Figure 5.1: a) Topological insulator coated with a FM film. The flux tube employed in Laughlin’s argument and the induced electric field are indicated. b) Flux tube threading a cylindrical hole in a STI. Arrows illustrate the helical spin state for upward moving electrons (for downmovers the arrows are reversed). Laughlin’s argument fails because the flux tube inserted into the bulk of the STI is not inert. We demonstrate below that when Φ = (s + 1/2)Φ0 , with s integer, the flux tube carries topologically protected gapless fermionic modes and forms a conducting quantum wire – a ‘wormhole’ – along which the accumulated surface charge can escape to another surface of the sample. In the end, no net fractional charge is accumulated at the surface and Laughlin’s argument instead predicts, indirectly, a new effect associated with a Dirac string in the bulk of a STI that we propose to call a ‘wormhole effect’.  5.3  Analytical calculation  The analytical calculation was performed by Prof. Marcel Franz and is included here for completeness. We begin by considering a bulk STI with a cylindrical hole of radius R threaded 61  by magnetic flux Φ = ηΦ0 with 0 ≤ η < 1 as illustrated in Fig. 5.1b. By solving the  Dirac equation for the surface electrons we show that a gapless state exists when η=  1 2  and persists in the limit R → 0. According to Ref. [94] electron states on a  curved surface of a STI, characterized by a normal unit vector n, ˆ are described by a Dirac Hamiltonian of the form 1 H = v h¯ ∇ · nˆ + nˆ · (p × σ ) + (p × σ ) · nˆ 2  (5.3)  where v is the Dirac velocity, p = −i¯h∇ is the momentum operator and σ =  (σ1 , σ2 , σ3 ) is the vector of Pauli spin matrices. The magnetic flux is included  by replacing p with π = p − (e/c)A, where A = ηΦ0 (ˆz × r)/2πr2 is the vector potential. For a cylindrical inner surface nˆ = −(cos ϕ, sin ϕ, 0), the Hamiltonian (5.3)  becomes, in cylindrical coordinates and taking v = h¯ = 1, Hk = −  1 + 2R  1 R (i∂ϕ  + η)  −ike−iϕ  − R1 (i∂ϕ + η)  ikeiϕ  .  (5.4)  We assumed a plane-wave solution eikz along the cylinder axis and replaced −i∂z → k.  The eigenstates of Hk are of the form Ψk (ϕ) =  fk iϕ e g  eiϕl  (5.5)  k  ˜ k = ( fk , gk )T is an eigenstate of H˜kl = σ2 k − σ3 (l + with l integer. The spinor Ψ 1 2  − η)/R with an energy eigenvalue  Ekl = ±v¯h k2 +  (l + 12 − η)2 . R2  (5.6)  For a generic strength of the magnetic flux the spectrum of electrons along the cylindrical surface shows a gap ∆=  2v¯h 1 −η . R 2 62  (5.7)  When η = 12 , i.e. at half flux quantum, the l = 0 mode becomes gapless, Ek0  η=1/2  =  ±v¯h|k|, independent of the hole radius R. This is the wormhole effect introduced  above. Physically, the necessity of the flux for the gapless state to occur stems from the Berry’s phase π acquired by electron spins in the helical state depicted in Fig. 5.1b. The gapless state occurs at half flux quantum when the Aharonov-Bohm phase exactly cancels the spin Berry’s phase.  5.4  Numerical study  We now study the wormhole effect using a concrete lattice model of a topological insulator which we solve by exact numerical diagonalization. We use the model presented in Section 3.2. At half filling, depending on the values of the parameters λ , t, ε, the system can be a trivial insulator, as well as a STI and WTI [29, 33]. Below, unless stated otherwise, we work with parameters 2t < ε < 6t, corresponding to a STI phase characterized by the Z2 invariant (1; 000). All energies are expressed in units of λ which we take equal to 1. Using our model with a single flux tube in a geometry with open boundary conditions it is possible to visualize the flow of charge at the intermediate steps of Laughlin’s flux insertion argument. To this end we consider a cube of size L3 and supplement the Hamiltonian with a surface magnetization term HS = −ΩS  ∑ j∈surf  rˆ j · Ψ†j σ Ψ j .  (5.8)  Here rˆ j represents the unit vector pointing outward from the origin located at the cube’s center and ΩS is the surface magnetization strength. HS breaks TRI at the sample surface and a gap of size ∼ 2|ΩS | opens up in the spectrum of the surface  states. Figure 5.2 shows the evolution of charge δ Q accumulated near the intersection of the flux tube with the magnetized surface as a function of η = Φ/Φ0 .  For small η we observe δ Q = 12 eη, consistent with the fractional Hall conductivity σxy = e2 /2h expected on the basis of Eq. (5.1). At η =  1 2  a charge e/2 travels  along the flux tube and combines with the negative charge that has built up on the opposite surface. For η >  1 2  the charge δ Q grows again at the rate controlled by  σxy until it reaches δ Q = 0 at η = 1. As already mentioned above, a Dirac string 63  carrying a full flux quantum Φ0 can be removed by a gauge transformation and the above evolution is thus consistent with the expectation that this weakly interacting system returns to the original configuration at the end of a full cycle.  Figure 5.2: Charge δ Q (in units of e) induced in the lower half of the sample by a flux tube carrying flux Φ. Insets show the charge density at indicated values of flux. A cube-shaped sample with 203 lattice sites and parameters t = 1, ε = 4 is used, ΩS = 1. δ Q is defined relative to the η = 0 situation. Figure 5.3 displays a modified arrangement with the flux tube terminated by a magnetic monopole located at the center of the sample. This furnishes a realization of the Witten effect [32] in a STI [33]. As a function of increasing η, the charge first accumulates at the intersection of the flux tube and the surface. At η =  1 2  a  charge e/2 travels along the wormhole to the monopole, the corresponding charge density clearly visible in the inset to Fig. 5.3. At η = 1 the flux tube becomes invisible but the e/2 charge remains bound to the monopole as expected on the basis of general arguments [32, 33].  64  Figure 5.3: Charge δ Q in the sphere of radius r0 centered at the monopole as a function of monopole strength. We have chosen r0 = 9.097 so that the sphere is large enough to enclose the expected charge e/2 for the unit monopole projecting flux Φ0 . A cube-shaped sample with 203 lattice sites and parameters t = 1, ε = 4 is used, ΩS = 0. δ Q is defined relative to the η = 0 situation.  5.5  Experimental detection  We now address the possibility of experimental detection of the wormhole effect predicted in this section. In a real physical system it is not possible to confine magnetic flux to an area of size comparable to the crystal lattice spacing as would be necessary to probe the wormhole effect in its pure form. However, it should be possible to observe a closely related effect in a nanoscale hole fabricated in a STI crystal with a uniform magnetic field applied parallel to its axis, Fig. 5.1b. Sweeping the magnetic field strength will result in a periodic variation of the conductance along the hole with minima at nΦ0 and maxima at (n + 1/2)Φ0 as the excitation spectrum oscillates between insulating and metallic. Such variations should be observable experimentally if certain conditions are met. First, the hole radius R must be sufficiently large so that several oscillations can be observed in the available 65  range of the laboratory field B. This gives R  (NΦ0 /πB)1/2 for N oscillations.  Second, R must be sufficiently small so that the maximum spectral gap Eq. (5.7) is large compared to kB T , or otherwise the oscillations in the conductance will be √ washed out by thermal broadening. This gives R h¯ v/ 2kB T . Taking typical values B = 10T, N = 10, v = 5 × 105 m/s and T = 1K yields 36 nm R 450 nm.  Thus, the experimental challenge would lie in fabricating a sub-micron size hole (or an array of holes) in a STI crystal or a thick film and measuring the conductance along the holes. In reality a magnetic field can not, at present, be localized within such a small distance. Therefore, we assume that the magnetic field is uniform, so that it exists both in the holes and in the bulk of the TI. This magnetic field breaks TRI in the bulk, so we may ask, what is the maximum magnetic field allowed? As a rough estimate, we require that the spacing between the Landau levels caused by the magnetic field be much smaller than the bulk gap. If the spacing approaches the size of the bulk gap, then it becomes difficult to discern the bulk gap from the spacing between the levels in the conduction and valence bands, and we expect the effect  to break down. This condition can be written as h¯ ωc ∆, where ωc = eB/m∗ is the cyclotron frequency. Therefore, we require B m∗ ∆/e¯h 6T , by substituting the bulk gap ∆  0.3 eV and the effective mass m∗  0.14me , for Bi2 Se3 [95]. The  condition for the radius (derived above) shows that as we decrease the magnetic field, the minimum hole size required increases, but the maximum hole size is not affected. These two conditions can be used in conjunction, taking into account the radius of the holes that can be fabricated, to find the maximum magnetic field. For example, for B = 0.5T we find that the above minimum hole size increases from 36nm to 160nm, and B  6T , so that both conditions could conceivably be fulfilled  at once.  5.6  Conclusions  The wormhole effect studied in this section represents a distinct bulk manifestation of the unusual electron properties in a strong topological insulator. Its existence resolves a conceptual dichotomy that arises when Laughlin’s argument is applied to the magnetized STI surface and exemplifies a unique bulk-surface correspondence  66  inherent to STIs. A closely related counterpart of the wormhole effect should be observable in artificially engineered nanostructures fabricated from available STIs.  67  Chapter 6  Topological Anderson Insulators 6.1  Overview  Disorder, ubiquitously present in solids, is normally detrimental to the stability of ordered states of matter. In this section we demonstrate that not only is STI robust to disorder but, remarkably, under certain conditions disorder can become fundamentally responsible for its existence. We show that disorder, when sufficiently strong, can transform an ordinary metal with strong spin-orbit coupling into a strong topological ‘Anderson’ insulator, a new topological phase of quantum matter in three dimensions.  6.2  Introduction  Disorder is well known to play a fundamental role in low-dimensional electronic systems, leading to electron localization and consequent insulating behavior in the time-reversal invariant systems [96]. Disorder also underlies much of the phenomenology of the integer quantum Hall effect [97]. Recently, in a remarkable development, it has been noted first by numerical simulations [52] and shortly thereafter by analytical studies [3], that a phase similar to the two dimensional topological insulator (also known as the quantum spin-Hall insulator [27, 53]) can be brought about by introducing non-magnetic disorder into a 2D metal with strong spin orbit coupling. This new 2D topological phase, referred to as TAI, has a dis-  68  ordered insulating bulk with topologically protected gapless edge states that give rise to precisely quantized conductance e2 /h per edge. In TAI, remarkably, conductance quantization owes its very existence to disorder. A question naturally arises whether such behavior can occur in three spatial dimensions. More precisely, one may inquire whether an inherently 3D topological phase analogous to the STI [64, 65, 90] could be reached by disordering a clean system that is initially in a topologically trivial phase. This is a nontrivial question because just as the 3D STI cannot be reduced to the set of 2D topological insulators, the existence of a 3D ‘strong’ TAI presumably cannot be deduced from the physics of 2D TAI. Below, we show the answer to the above question to be affirmative. We construct an explicit example of a disorder-induced topological phase in three dimensions with physical properties analogous to those of the STI. We propose to call this new phase a strong topological Anderson insulator (STAI). We argue that some of the topologically trivial compounds with strong spin-orbit coupling discussed in the recent literature, such as e.g. Sb2 Se3 , could become STAI upon introducing disorder. In other compounds that already are STIs in their clean form, disorder can reinforce this behavior by rendering the bulk truly insulating. We note that the authors of Ref. [3] anticipated the existence of a 3D disorder-induced topological phase.  6.3  TAI in 2D  This section follows [3] closely, filling in the gaps. We start from the two dimensional low energy effective Hamiltonian of a HgTe quantum well H = α(kx σx − ky σy ) + (m + β k2 )σz + [γk2 +U(r)]σ0  (6.1)  The full Hamiltonian is a 4 × 4 block diagonal matrix, with a block H and a block  H ∗ , but since there is no spin coupling, we can treat each block separately. This is a Dirac Hamiltonian, with the important addition of the two quadratic terms. The parameters α, β , γ are material dependent and U is the site specific disorder. We define the ‘clean’ Hamiltonian, H0 which is defined as H with no disorder (U = 0), which has gap m. 69  The question which we wish to answer is: if we start with an insulator - noninverted gap (m > 0) and Fermi energy in the gap |εF | < m (here εF is measured from the location of the Dirac point for m = 0), and switch on disorder, can the disorder  make this insulator into a topological insulator? To achieve this, the renormalized mass would have to be inverted m¯ < 0, with the renormalized Fermi energy in the gap |ε¯F | < |m|. ¯ Phrased differently, we will show that disorder moves the phase  transition from m = 0 to some positive value of m.  We start from the self consistent Born approximation with point-like impurities. The self energy is given by (see derivation in Appendix A.4) Σ(k) = Nu2  2  a 2π  d 2 k [εF − H0 (k) − Σ(k) + i0+ ]−1  (6.2)  where u2 is the variance of a random variable uniformly distributed in the range [−U0 /2,U0 /2], and is given by 1 u2 = NU0  U0 2  −  U0 2  2 u3 u2 du = NU0 3  U0 2  =  0  U02 12N  (6.3)  and the infinitesimal imaginary part is required in order to make the integral converge. From now on we suppress it, remembering that wherever εF appears inside the integral, it includes an infinitesimal imaginary part. For convenience we write the self energy as Σ = ∑µ Σµ σµ and similarly, the Hamiltonian is H0 = ∑µ αµ σµ with α0 ≡ γk2 , α1 = αkx , α2 = −αky , α3 = m + β k2 . We substitute these into the  argument of (6.2)  εF − H0 (k) − Σ(k) = εF − ∑(αµ + Σµ )σµ  (6.4)  µ  ≡ ε˜F − ∑ α˜ i σi i  where we have defined ε˜F ≡ εF − (α0 + Σ0 ), and α˜i ≡ αi + Σi .  However, we need the inverse of this matrix, which is inconvenient to work  with. To get rid of it, we now take the inverse, but then multiply by the same matrix but with a plus, and its inverse. This is similar to the trick used when we have the inverse of a complex number, and multiply and divide by the complex 70  conjugate, to get a real number in the denominator. −1  [εF − H0 (k) − Σ(k)]  =  −1  ε˜F − ∑ α˜ i σi  ε˜F + ∑ α˜ i σi i  i  =  ε˜F + ∑i α˜ i σi ε˜F2 − ∑i α˜ i2  −1  ε˜F + ∑ α˜ i σi i  (6.5)  We substitute this expression into (6.2), to find Σ=  U02 a 12 2π  2  ε˜F + ∑i α˜ i σi ε˜F2 − ∑i α˜ i2  (6.6)  εF − Σ0 − γk2  (6.7)  d2k BZ  which we can separate into four equations Σ0 = A  d2k BZ  Σ1 = A  d2k BZ  Σ2 = A  d2k BZ  Σ3 = A  d2k  D(k) αkx + Σ1  D(k) −αky + Σ2  D(k) m + Σ3 + β k2 D(k)  BZ  where we have defined the prefactor A ≡ ε˜F2 − ∑i α˜ i2 .  U02 12  a 2 2π  and the denominator D(k) ≡  Note that the form of the self energy shows that the Fermi energy and mass  are renormalized by Σ0 and Σ3 (respectively), so we define the renormalized values ε¯F ≡ εF − Re Σ0 and m¯ ≡ m + Re Σ3 , and we have ε¯F = εF − A Re  d2k BZ  m¯ = m + A Re  d2k BZ  71  ε˜F − γk2 D(k)  m¯ + β k2 D(k)  (6.8) (6.9)  We write out the denominator D(k) ≡ ε˜F2 − ∑ α˜ i2  =  i  (ε¯F − γk2 )2 − (αkx + Σ1 )2 − (αky − Σ2 )2 − (m¯ + β k2 )2  → (εF − γk2 )2 − α 2 k2 − (m + β k2 )2 =  (6.10)  (εF2 − m2 ) − [α 2 + 2(εF γ + mβ )]k2 + (γ 2 − β 2 )k4  Where the limit is the zero order Born approximation, for which we take Σµ → 0 and m¯ → m and ε¯F → εF on the right side of the above expressions for the self  energy and renormalized mass and Fermi energy. In this limit Σ1 = Σ2 = 0, since we are left with an odd integrand over even limits. Furthermore, in this limit we can evaluate the integrals in (6.8) and (6.9). We note that the integrand depends only on k2 so we can change to polar coordinates, d 2 k = 2πkdk, and impose a cutoff Λ, to find  2πAγ = εF − 2 Re  γ −β2  ε¯F  dk 0  εF2 −m2 γ 2 −β 2    εF 3 γ k−k  Λ  −α  2 +2(ε  F γ+mβ ) γ 2 +β 2  k2 + k4    (6.11)  The integral we need to solve is of the form Λ  dk 0  Ak − k3 B +Ck2 + k4  1 Λ4 − ln 4 B  (6.12)  where we have kept only the part that diverges logarithmically for Λ → ∞. The logarithmic divergence is expected, since for k → ∞ the integrand is approximately −1/k. We substitute, and take Λ = π/a, to find  U02 a 2 γ 1 γ2 − β 2 π ln 12 2π β 2 − γ2 4 εF2 − m2 a 1 U02 a2 γ β 2 − γ2 π 4 = εF − ln 2 48π β 2 − γ 2 εF2 − m2 a  ε¯F = εF − 2π  4  (6.13)  and similarly m¯ = m −  1 U02 a2 β β 2 − γ2 π ln 2 48π β 2 − γ 2 εF2 − m2 a 72  4  (6.14)  Note that the correction has two relative minus signs compared with the correction to the Fermi energy. The first comes from −k3 → k3 in the integrand, which gives  an overall minus in the result of the integral. The second comes from the fact that the pre-factor of the integral is negative for the Fermi energy but positive for the  mass (see (6.8) and (6.9)). We note that the correction to the mass is negative for β > γ which is a necessary but insufficient condition to get a TAI. We can use these equations to find the phase boundaries. We start by looking for the critical mass inversion line, defined by m¯ = 0. We find m =  U02 a2 β β 2 − γ2 π ln 48π β 2 − γ 2 εF2 − m2 a  4  (6.15)  which can be solved to give π |εF − m | = a 2  2  4  2  2  |β − γ | e  − 48πm 2 a  β 2 −γ 2 1 β U02  (6.16)  We see that if we impose m¯ = 0 and keep m = const, then as we increase U0 , we must also increase εF . The curve described by this equation separates the phase space into two regions: m¯ > 0 and m¯ < 0, and if a TAI phase exists, it will be in the second region. If we start from the curve m¯ = 0 and increase U0 or decrease εF , we decrease the effective mass, so this would bring us into the region m¯ < 0. In addition to the mass inversion, the renormalized Fermi energy must be inside the renormalized mass. The critical condition for this is ε¯F = ±|m|. ¯ For m¯ ≥ 0, we  have ε¯F = ±m, ¯ so  εF = ±m −  U02 a2 1 β 2 − γ2 π ln 2 48π γ ± β εF − m2 a  4  (6.17)  For m¯ ≤ 0, we have ε¯F = ∓m, ¯ so the signs are reversed from the above, but the  curves are the same. This corresponds to starting from U0 = 0 and increasing the  disorder. When we cross the mass inversion line, the mass inverts, so that the curves swap. The TAI phase is to the right of the first line, and in between the two other lines. The three lines: the mass inversion (critical) line and the two others intersect at a point. This can be seen by looking for an intersection for m¯ ≥ 0, so εF = m¯ = −m¯ 73  week ending  P H Y S I C A L soRthat E Vm¯I = E 0Wwhich L EisTthe T Econdition RS for the critical line. 6 NOVEMBER 2009  tic term p2 ¼ À@2 r2 in H, ate / eÀx= , adds a negative ogical mass. The renormalized refore have the opposite sign as gical mass renormalization by change in the phase diagram, without the terms quadratic in " and m then remains the of m nnot appear. " as well zed topological mass m, cal potential , " from the selfveraged effective medium. To mputer simulations [2,3], we attice (lattice constant a) and sorder potential U, uniformly ÀU0 =2, U0 =2). We denote by an of the clean quantum well in FIG. 1 (color online). Computer simulation of a HgTe quan5,7]. [10],conductance showing theG average conductance hGi strength as a U0 Figuretum 6.1: well Average as a function of disorder by function of disorder strength U (logarithmic scale) and Fermi  0 and the Fermi energy EF in a 2D strip of dimensions 400a × 100a. , in a disordered strip of widthfrom W ¼the 100a length Born energy E F ¼ hðEF À (2) Curves A and B are the phase boundaries selfand consistent L ¼ 400a. The TAI phase is indicated. Curves A and B are the approximation. Figure from Ref. [3]. phase boundaries resulting from the effective medium theory. age, is a 2  2 matrix which we Curve A separates regions with positive and negative renormalatrices: Æ ¼ Æ0 0 þ Æx x þ " while curve B marks the crossing of the ized topological mass m, malized topological mass and renormalized chemical potential  " with the band edge (jj " ¼ given by 6.4 TAI 3Dcurves have been calculated without any adjustable " inBoth Àm). parameter. The phase boundary of the TAI at strong disorder is  " ¼ EF À lim ReÆ0 : (3) k!0 outside oftreatment the regime of validity of the effective medium theory. 6.4.1 General  À1  HÞÀ1 i;  " ¼ topological insulator is at m After nters the (negative) band gap considering the 2D case, it is natural to ask whether the TAI phase can be much Equationa generalization (5) implies that theHamiltonian lower branch of previous found in 3D also.[10]. We consider of the in the " < 0) is then curve B (defined by  " ¼Hamiltonian m fixed EFpresented % n approximation, Æ is givensubsection by to three dimensions. The we will useatwas in m > 0 independent of U0 . This explains the puzzling abSection 3.2 and is the same one we used previously to study the Witten effect sence of the TAI phase in the valence band (EF < 0), and the Wormhole with an additional constant term to shift the chemical observed effects, in the computer simulations [2,3]. ½EF þ i0þ À H0 ðkÞ À ƊÀ1 : To quantitatively test the phase diagram resulting from (4) the effective medium theory, we performed computer simulations similar to those reported in Refs. [2,3]. The first Brillouin zone.) The selfconductance G is calculated from the lattice Hamiltonian omentum and diagonal (so there 74 [5,7] in a strip geometry, using the method of recursive the parameters , , 
). By Green functions. The strip consists of a rectangular disornction of EF and U0 we obtain dered region (width W, length L), connected to semiFig. 1. infinite, heavily doped, clean leads [14]. Theory and simuapproximate solution in closed lation are compared in Figs. 1 and 2. Figure 1 shows the phase diagram. The weak-disorder    2 À 
2 @4      boundary of the TAI phase observed in the simulations is  ln ; (5a)     2  2  described quite well by the self-consistent Born approxi-  potential, and a hopping term Hγ = 6γ ∑ c†iσ ciσ − γ i,σ  ∑  c†iσ c jσ  (6.18)  i, j ,σ  which is diagonal in spin and orbit, so the additional term in momentum space multiplies the identity. When looking at the low energy theory, it is obvious that this term is necessary in order to get the quadratic term that multiplies the identity which is present in the 2D case (see Section 6.3). The clean Hamiltonian (no disorder) can be written as H0 = ∑k Ψ†k H0 (k)Ψk , with H0 (k) = −2γ ∑ cos(ki a) − 2λ ∑ sin(ki a)τ3 ⊗ σi + mk τ1 ⊗ σ0 i  (6.19)  i  = d4 + ∑ d µ γ µ µ  where we have chosen the gamma matrices so that γ0 = τ1 ⊗σ0 and γi = τ3 ⊗σi , and  define d0 ≡ mk = ε − 2t ∑i cos(ki a), di ≡ −2λ sin(ki a), d4 = 2γ[3 − ∑i cos(ki a)].  The Hamiltonian is a sum of anti-commuting matrices (such that {γµ , γν } =  δµν ) and the identity. Therefore, the spectrum can be found easily. We subtract the part that multiplies the identity and square H0 (k) − d4  2  2  =  ∑ dµ γµ  =  µ  ∑ dµ dν {γµ , γν } = ∑ dµ2  µ>ν  (6.20)  µ  So that the energies are E0 (k) = d4 ±  ∑ dµ2 = d4 ±  m2k + εk2  (6.21)  µ  where εk = ±2λ  sin2 (kx a) + sin2 (ky a) + sin2 (kz a)  (6.22)  are the energies of the system with no ordinary hopping ε = t = 0, so only the spin orbit term remains. 75  To find the bandwidth ∆ we note that E0 (k) is maximum and minimum for k = (π, π, π), where it has the values 12γ ± |ε + 6t|, so the bandwidth is ∆ = 2(ε + 6t) = 2(m + 12t)  (6.23)  (where in the last line we have assumed ε > 0 and t > 0, and used m = ε − 6t)  As in the two dimensional case, we start from the self consistent Born approx-  imation with point-like impurities. The self energy is given by (see derivation in Appendix A.4) Σ(k) = Nu2  a 2π  3  d 3 k [εF − H0 (k) − Σ(k) + i0+ ]−1  (6.24)  where u2 is the variance of a random variable uniformly distributed in the range [−U0 /2,U0 /2], and is given by 1 u = U0 2  U0 2  −  U0 2  2 u3 u du = U0 3 2  U0 2  0  =  U02 12  (6.25)  and the infinitesimal imaginary part is required in order to make the integral converge. From now on we suppress it, remembering that wherever εF appears inside the integral, it includes an imaginary infinitesimal part. For convenience we write the self energy as Σ = ∑µ Σµ γµ +Σ4 . We substitute the Hamiltonian and self energy into the argument of (6.24) εF − H0 (k) − Σ(k) = (εF − d4 − Σ4 ) − ∑(dµ + Σµ )γµ  (6.26)  µ  ≡ ε˜F − ∑ d˜µ γµ µ  where we have defined ε˜F ≡ εF − d4 − Σ4 , and d˜µ ≡ dµ + Σµ . We note that d4 = 2γ[3 − ∑i cos(ki a)], so that at the Γ point, d4 does not shift the energies. Back to  the calculation, we need the inverse of this matrix, which is inconvenient to work with. To get rid of it, we take the inverse, but then multiply by the same matrix but with a plus, and its inverse. This is similar to the trick used when we have the inverse of a complex number, and multiply and divide by the complex conjugate.  76  We rely on the fact that the gamma matrices anti-commute with one another. −1  [εF − H0 (k) − Σ(k)]  =  −1  ε˜F − ∑ d˜µ γµ  ε˜F + ∑ d˜µ γµ  ε˜F + ∑ d˜µ γµ µ  µ  µ  =  −1  ε˜F + ∑µ d˜µ γµ ε˜F2 − ∑µ d˜µ2  (6.27)  We substitute this expression into (6.24), to find Σ=  U02 a 12 2π  3  d3k BZ  ε˜F + ∑µ d˜µ γµ ε˜F2 − ∑µ d˜µ2  (6.28)  εF − Σ4 − d4  (6.29)  which we can separate into five equations Σ4 = A  d3k BZ  Σ0 = A  d3k BZ  Σi = A  d3k BZ  where we have defined the prefactor A ≡ ε˜F2 − ∑µ d˜µ2 .  D(k) mk + Σ0 D(k) di + Σi D(k)  U02 12  a 3 2π  and the denominator D(k) ≡  Note that the form of the self energy shows that the Fermi energy and mass  are renormalized by Σ4 and Σ0 (respectively), so we define the renormalized values ε¯F ≡ εF − Re Σ4 and m¯ ≡ m + Re Σ0 , and m is the minimal value of |mk |. We find ε¯F = εF − Re Σ4 = εF − A Re  d3k BZ  ε¯F − γck D(k)  (6.30)  and m¯ = m + Re Σ0 = m + A Re  d3k BZ  77  m¯ + tck D(k)  (6.31)  where we have defined ck ≡ 2 3 − ∑ cos(ki a)  (6.32)  i  and we have assumed that the smallest gap is at the Γ point, so that m = ε − 6t and then we rewrite  mk = ε − 2t ∑ cos(ki a) = m + tck  (6.33)  i  We consider first the three equations for Σi . In the zero order Born approximation we have Σi = A  d3k BZ  di  (6.34)  D(k)|Σi =0  Write out the denominator D(k) ≡ ε˜F2 − ∑ d˜µ2 µ  = (εF − d4 − Σ4 )2 − ∑(dµ + Σµ )2  (6.35)  µ  = (εF − d4 − Σ4 )2 − (mk + Σ0 )2 − ∑(−2λ sin(ki a) + Σi )2 i  We see that for Σi = 0 the denominator is even D(k)|Σi =0 = D(−k)|Σi =0 , since mk and d4 depend on k only through cos, and sin appears squared. However, di ∝ sin(ki a) and so it is odd in k. Therefore, we have an odd integrand over even limits, and the integral gives zero, so we get Σi = 0 from (6.34). Therefore, in this case the zero order solution Σi = 0 is the final self consistent solution. The denominator is 2  D(k) =  2  − m¯ + 2t 3 − ∑ cos(ki a)  ε¯F − 2γ 3 − ∑ cos(ki a) i  i  = (ε¯F − γck ) − (m¯ + tck ) + λ sk 2  2  2  78  − (2λ )2 ∑ sin2 (ki a) i  (6.36)  where we have defined sk ≡ −4 ∑ sin2 (ki a)  (6.37)  i  6.4.2  Iterative solution  The equations for the self energy can be solved numerically by self consistent iteration. As explained above, Σi = 0 is the self consistent solution for the last three equations of (6.29). So we are left with two equations, for Σ0 and Σ4 . In addition, we wish to probe only the phase boundaries - probing the whole phase space would be a lengthy procedure. Therefore, we use the definition of the renormalized mass m¯ ≡ m + Re Σ0 and Fermi energy ε¯F ≡ εF − Re Σ4 and the phase boundary condition ε¯F = ±|m| ¯ to write the three self consistency equations  εF = Re Σ4 ± |m| ¯ = Re Σ4 ± |m + Re Σ0 | 2 U0 εF + i0+ − γck − Σ4 Σ4 = ∑ D(k)| 12N 3 k∈ BZ +  (6.38)  εF →εF +i0  Σ0 =  U02 m + tck + Σ0 ∑ 3 12N k∈BZ D(k)| + εF →εF +i0  For the second and third equations, we notice that we are summing over expressions that are even in ki - they are invariant under the three reflections ki → −ki .  Therefore, instead of calculating the sum for the full Brillouin zone, we can calculate it just for the positive eighth of the BZ and multiply by 8 (while being careful not to over count). The algorithm for the calculation is as follows: 1. For each U0 in some range (U0 ≥ 0) 2. While we haven’t reached self consistency (a) I0 = I4 = 0 (sums) (b) If this is the first run, set Σ0 , Σ4 (complex seed values), otherwise use Σ0 and Σ4 from previous run as seeds 79  (c) For each ki = 2π nNi , ni = 1, . . . , N and i = x, y, z and N = 50 (for example) i. Calculate D(k)−1 ii. Add new terms to sums I0 = I0 + (m + tck + Σ0 )D(k)−1 , I4 = I4 + (εF + i0+ − γck − Σ4 )D(k)−1  (d) Calculate the new values Σ˜ 0 = (U02 /12N 3 )I0 , Σ˜ 4 = (U02 /12N 3 )I4 , ε˜F = Re Σ˜ 4 ± |m + Re Σ˜ 0 | ˜ 0 −Re Σ0 ˜ 0 −Im Σ0 < tol, Re ΣRe < tol and Im ΣIm < tol (and the Σ0 Σ0 same for Σ4 ), we have reached self consistency - save these values  (e) If  ε˜F −εF εF  (f) Else if number of iterations > 20 (for example) then the exact value will remain undetermined Using an imaginary seed and adding the infinitesimal imaginary part to εF are critical, or else the sums will not converge. In a certain part of the phase diagram, the iterative procedure fails to reach self consistency. When using the correct values as seed values, it does reach self consistency. We conclude that the iterative method is consistent, however it is not strong enough for these particular parameters. Instead we use the Matlab built in function fsolve, which solves a system of nonlinear equations. We take the equations (6.38) and rearrange them so that they are homogeneous, as required for this method Re Σ4 ± |m + Re Σ0 | − εF = 0 U02 εF + i0+ − γck − Σ4 − Σ4 = 0 ∑ D(k)| 12N 3 k∈ BZ + εF →εF +i0  U02 12N 3  ∑ k∈BZ  m + tck + Σ0 D(k)|ε  + F →εF +i0  80  − Σ4 = 0  (6.39)  6.4.3  Non iterative solution  We can also probe the phase transition lines directly. We start by rewriting the denominator D(k) using the identity (a + b)2 − (c + d)2 = [(a + c) + (b + d)][(a − c) + (b − d)]  (6.40)  So D(k) = (ε¯F − γck )2 − (m¯ + tck )2 + λ 2 sk  (6.41)  = [ε¯F + m¯ + (t − γ)ck ] [ε¯F − m¯ − (t + γ)ck ] + λ 2 sk  This expression is more convenient to work with when looking for the phase boundaries, since it contains the combinations ε¯F ± m, ¯ of which at least one will be zero on the phase boundaries.  Recall (6.30) and (6.31) ε¯F = εF − A  d3k BZ  m¯ = m + A  d3k BZ  ε¯F − γck  D(k) m¯ + tck  (6.42)  D(k)  and we must remember to take only the real part of the integrals. Add and subtract m¯ + ε¯F = m + εF + A  d3k BZ  m¯ − ε¯F = m − εF + A  d3k  (m¯ − ε¯F ) + (t + γ)ck  D(k) (m¯ + ε¯F ) + (t − γ)ck  BZ  (6.43)  D(k)  Now look for critical lines |ε¯F | = −m¯ > 0 (since m¯ < 0). We break this down into two cases and consider first the case m¯ + ε¯F = 0. Subsitute into the last two equations to find 0 = m + εF + A  d3k BZ  2m¯ = m − εF + A  d3k BZ  81  2m¯ + (t + γ)ck D+ (k) (t − γ)ck D+ (k)  (6.44)  and D+ (k) = −(t − γ)ck [2m¯ + (t + γ)ck ] + λ 2 sk  (6.45)  Now define the integrals I1+ ≡  BZ  I2+ ≡  BZ  d3k d3k  1 D+ (k) ck  (6.46)  D+ (k)  Rewrite the two equations in terms of these integrals 0 = m + εF + A 2mI ¯ 1+ + (t + γ)I2+  (6.47)  2m¯ = m − εF + A(t − γ)I2+ Solve for A and εF (in terms of A) to find A =  m¯ − m mI ¯ 1+ + tI2+  (6.48)  εF = −m − A 2mI ¯ 1+ + (t + γ)I2+ Similarly, for the second case, ε¯F − m¯ = 0 we have 2m¯ = m + εF + A  d3k BZ  0 = m − εF + A  d3k BZ  (t + γ)ck D− (k) 2m¯ + (t − γ)ck  (6.49)  D− (k)  We continue as for the previous case, by defining the integrals I1− and I2− , and solving for A. The final result can be written as A± =  m¯ − m mI ¯ 1± + tI2±  εF± = ∓m ∓ A 2mI ¯ 1± + (t ± γ)I2±  82  (6.50)  and the denominator (we write both cases together now) D± (k) = − [(t ∓ γ)ck ] [2m¯ + (t ± γ)ck ] + λ 2 sk  (6.51)  for the phase boundary lines m¯ ± ε¯F = 0.  Note that for m¯ = m there is no renormalization of the mass, and indeed we get  A±  = 0 so U0 = 0 - no disorder, and εF± = ∓m, since the phase boundary occurs  when the Fermi energy is at the edge of the gap, which is just m in this case.  Let us now discuss the meaning of zeros of the denominator. Recall that the spectrum is given by (6.21) [E0 (k) − d4 ]2 = ∑ dµ2  (6.52)  µ  and the Fermi surface is defined by E0 (k) = εF . Now write the denominator in the zero order Born approximation, and equate it to zero D(k) = (εF − d4 )2 − ∑ dµ2 = 0  (6.53)  µ  We see that this equation is equivalent to (6.52) with E0 (k) = εF , so that the zeros of the denominator define the Fermi surface. If we start with a metal close to a phase transition to an insulator, the Fermi surface shrinks as we bring the Fermi energy closer to the edge of the gap. Once the Fermi energy crosses the gap, we have an insulator and the Fermi surface does not exist. Hence, at the phase transition we expect to get a Fermi surface of zero volume, in other words a point. Indeed, by inspecting (6.51) we see that it is zero only at the Γ point (0,0,0). The algorithm for finding the phase boundary lines is 1. For each m¯ (a) Calculate the integrals I1± and I2± numerically (b) Calculate U0± =  12A± (2π/a)3 and εF±  The advantages of this method are the elimination of the iterative process, probing only the zone that interests us (the boundaries), so we scan a 1D (m) ¯ phase space 83  instead of a 2D phase space (εF and U0 ). To check our equations, we can evaluate (6.50) analytically in the low energy limit. We recover the zero order equations for the phase boundaries, which are calculated in the next section.  6.4.4  Analytical zero order Born approximation  In order to study the phase boundaries analytically, we consider the zero order Born approximation, by plugging in Σ = 0 on the right side of (6.29). As argued above, Σi = 0 is a self consistent solution. We are left with equations for Σ0 and Σ4 (in the limit Σ = 0 on the right side) Σ4 = A  d3k BZ  Σ0 = A  d3k BZ  εF − γck  D(k)|Σ=0 m + tck  (6.54)  D(k)|Σ=0  and the denominator D(k)|Σ=0 = (εF − γck )2 − (m + tck )2 + λ 2 sk  (6.55)  These integrals are probably too complicated to solve exactly analytically, so we try a low energy approximation. If we take ε > 6t > 0 then the mass at all Dirac points (which are located at the TRIM) is positive, since ε ± 2t > 0 and ε ± 6t > 0.  In this case, the gap at the Γ point (0, 0, 0) is the smallest. To see this, set ε = 6t +δ . Then the gap at the Γ point is δ and the gap at the other points is 4t + δ , 8t + δ , 12t + δ . So, as long as t > 0 and δ  t we may approximate by linearizing around  the Γ point. This makes sense, since if the gap at the Γ point is by far the smallest, and we turn on disorder, it is reasonable to expect this gap to be the first to invert. If only the mass at the Γ point inverts, we will have seven Dirac points with positive mass and one with a negative mass, and therefore a STI. For very high disorder, it is possible that the other masses will also invert, but we consider only the simplest  84  case here. We approximate to second order in k 1 2 3 − 3 + (ak)2 = (ak)2 2  ck = 2 3 − ∑ cos(ki a) i  sk = −4 ∑ sin2 (ki a)  (6.56)  −4(ak)2  i  Substitute in the self energy equations Σ4 = A  d3k Λ  Σ0 = A  d3k  εF − γ(ak)2  (6.57)  D(k)|Σ=0 m + t(ak)2 D(k)|Σ=0  Λ  and in the denominator D(k)|Σ=0 =  εF − γ(ak)2  2  − m + t(ak)2  2  − v2 (ak)2  We see that the integral only depends on k2 , therefore we can integrate over a sphere d 3 k = 4πk2 dk, so that Σ4 = 4πA Σ0 = 4πA  Λ  dk 0 Λ  dk  [εF − γ(ak)2 ]k2  (6.58)  D(k)|Σ=0 [m + t(ak)2 ]k2  0  D(k)|Σ=0  We rearrange the integrals Σ4 =  4πA γ a2 γ 2 − t 2  Σ0 =  4πA t a2 γ 2 − t 2  Λ  dk 0 Λ  dk 0  εF2 −m2 (γ 2 −t 2 )a4  εF 2 k − k4 γa2  2  − 2(ε(γF γ+mt)−v k2 + k4 2 −t 2 )a2  (6.59)  m 2 k + k4 ta2 2 εF2 −m2 − 2(ε(γF γ+mt)−v k2 + k4 2 −t 2 )a2 (γ 2 −t 2 )a4  The integrals we need to solve are of the form Λ  dk 0  Ak2 ± k4 B +Ck2 + k4 85  ±Λ  (6.60)  Where we have kept only the part that diverges linearly for Λ → ∞. The additional  terms are constant in the limit Λ → ∞. We set Λ = π/a, substitute the result of the integral and the prefactor A to find Σ4 = −  4π U02 a a2 12 2π  3  U02 γ π γ = 2 2 2 γ −t a 24π t − γ 2  (6.61)  and similarly Σ0 = −  U02 t 2 24π t − γ 2  (6.62)  Therefore, the renormalized Fermi energy and mass, in the zero order Born approximation are given by U02 γ 24π t 2 − γ 2 U2 t m¯ ≡ m + Re Σ0 = m − 0 2 24π t − γ 2  ε¯F ≡ εF − Re Σ4 = εF −  (6.63)  In this rough approximation, the correction to the Fermi energy and the mass does not depend on εF , and also not on m. Recall that we first took the zero order Born approximation by taking Σ = 0 in the integrals, thereby replacing m¯ → m, and ε¯F →  εF , and then took the low energy approximation and imposed a cutoff. However,  had we considered the full Born expressions in the low energy approximation, the results would have been exactly the same. Therefore, this result is not actually a zero order Born approximation, but rather a low energy approximation. To get an inverted renormalized mass, the correction to the mass must be negative (assuming the initial mass is positive m > 0), so we must have γ < t (assuming t, γ > 0). We can use these equations to find the phase boundaries. We start by looking for the critical mass inversion line, defined by m¯ = 0. We find m−  U02 t =0 24π t 2 − γ 2  86  (6.64)  or U0 =  24πm  t2 − γ2 ≡ Uc t  (6.65)  Assuming t > γ > 0 and noting that the disorder strength U0 ≥ 0 by definition.  This is the critical value of disorder Uc necessary for mass inversion. Unlike the 2D case, the critical value in this approximation does not depend on the Fermi energy εF . In addition to the mass inversion, the renormalized Fermi energy must be inside the renormalized mass. The critical condition for this is ε¯F = ±|m|. ¯ For m¯ ≥ 0, we have ε¯F = ±m, ¯ so  εF = ±m −  U02 1 24π γ ± t  (6.66)  and for m¯ ≤ 0, we have ε¯F = ∓m, ¯ so the signs are reversed from the above, but  the curves are the same. This corresponds to starting from U0 = 0 and increasing the disorder. When we cross the mass inversion line, the mass inverts, so that the curves swap. The TAI phase is to the right of the first line, and in between the two other lines. These two lines are parabolas in the (εF ,U0 ) phase space. The three lines: the mass inversion (critical) line and the two others intersect at a point, by definition. This can be seen by looking for an intersection for m¯ ≥ 0, so  εF = m¯ = −m¯ so that m¯ = 0 which is the condition for the critical line.  6.4.5  Witten effect in quasi 1D system  To further confirm the 3D nature of the observed topological phase we probed for the Witten effect [32] in our model STAI. According to Refs.[29, 30] the effective electromagnetic Lagrangian of a 3D strong topological insulator contains an unusual ‘axion’ term ∼ θ E · B with θ = π. According to Witten [32] a  magnetic monopole inserted into a medium with non-zero θ binds electric charge −e(θ /2π + n) with n integer.  Specifically, as for the ordinary TI, it is expected that on the insertion of a unit  magnetic monopole into a TAI, fractional charge −e/2 will bind to the monopole.  At first we used the set up used previously, namely a cube of side 14. However, it 87  60 50 40  EF  30 20 10 0 ï10 0  20  40  60  80  100 U0  120  140  160  180  200  Figure 6.2: Plots of the Born approximation phase boundary lines. The two solid black lines are equations (6.66), which tell us when the Fermi energy crosses the gap. The vertical dashed black line is the critical disorder strength for mass inversion. These three lines are for the zero order Born approximation for which we have a closed form expression. The coloured curves were obtained via the self consistent Born approximation. Parameters: t = 24 λ = 20, ε = 6t + δ , δ = t/24, γ = 16. The critical disorder strength Uc 31.7 and the bandwidth ∆/m 139. For U0 = 0 the difference between the curves is 2m. In that region the Fermi energy is in the gap and the gap is positive (since the disorder is smaller than the critical disorder), so the material is an insulator. Above and below this, it is a metal. In the triangle on the right, the material is a TAI, since the mass is inverted (disorder larger than the critical disorder) and the Fermi energy is in the gap since we are between the two solid black lines. . was found that the disorder is strong enough in this case to enlarge the localization length so that it is comparable to the system size and the effect cannot be observed (with the computational resources at our disposal). Therefore, we had to come up with a geometry that would allow us to probe a much longer localization length. The solution is a quasi 1D system containing a monopole and anti-monopole with periodic boundary conditions. Periodic boundary conditions are preferable since then the low energy edge states are eliminated and therefore they do not mix  88  Bulk+Surfaces Bulk  c)  ST  AI  STAI  0  100  150  200  250  300  U0 (meV) d) Figure 6.3: False color plot of conductance as a function of disorder strength U0 and the Fermi energy EF . Each data point corresponds to a single realization of the disorder potential. The color scale is chosen to emphasize the effect of fluctuations in G. The dashed line marks the band inversion boundary defined as Re(m) ¯ = 0 based on the Born approximation and the solid lines represent the Born approximation phase boundaries separating a band insulator and a metal defined by |Re(E¯F )| = |Re(m)|. ¯ We use parameters ε = 145meV, t = 24meV, λ = 20meV, and γ = 16meV, corresponding to m = 1meV. The calculation of G is outside the scope of this dissertation and is included in Ref. [98].  35 30 25 20  z  15 10 5 0  . x  with the monopole and anti-monopole low lying states, therefore doubling the localization length that we can probe. Periodic boundary conditions can only be used for a configuration which has zero flux through all faces. This is achieved by altering the shape of the magnetic field lines, so that the field is uniform in the z direction in between the monopole and anti-monopole and zero elsewhere. It is expected that the Witten effect will survive for any field configuration which can be adiabatically deformed to a magnetic point charge. The reason is that the charge localized at the monopole/anti-monopole is embedded in an insulator. In the absence of metallic modes (such as the worm hole effect, see Chapter 5), this charge has  89  a)  δρ  b) 6.4: 0.20 δρ Figure A prism-shaped sample with monopole (+), anti-monopole (−), δ ρ¯ allows us to and field lines indicated. The asymmetric field distribution use 0.10 PBC in all directions and thus avoid difficulties associated with the surface states. 0.00  nowhere to go. In addition, in order to maximize the localization length probed, ï0.10 we lengthen the z axis compared with the x and y axes. We tune the parameters ε  and t so that we have an ordinary insulator with a very small gap, and then switch ï0.20  1 on a strong disorder to get a25TAI.  50 z  75  100  We measured the induced fractional charge in a configuration containing a  c)  4997  0.50 monopole and an anti-monopole depicted in Fig. 6.4. Our results 4999presented in Fig. 5000 6.5 and Fig. 6.6 clearly indicate fractional charge ±e/2 bound to the monopole, 0.25  5003  δQ  confirming the expected value of θ = π. This result lends additional support to our identification0.00 of STAI as a genuinely 3D topological phase characterized by the bulk axion term. ï0.25 ï0.50 1  25  50 z  90  75  100  δρ δ ρ¯  0.20  δρ  0.10 0.00 ï0.10 ï0.20 1  25  50 z  75  100  Figure 6.5: Electric charge density δ ρ(z) = ∑x,y [ρ(r)−ρ0 (r)] induced by the monopole/anti-monopole pair in a 5 × 5 × 100 sample at half filling. Here ρ and ρ0 represent the charge density with and without the pair, ¯ respectively. A smoothed charge density function δ ρ(z), obtained by convolving δ ρ(z) with a Gaussian of width σ = 1.5, is also plotted to emphasize the charge bound to the monopole. 4997 4999 5000 5003  0.50  δQ  0.25 0.00 ï0.25 ï0.50 1  25  50 z  75  100  Figure 6.6: Integrated charge δ Q(z) = ∑z ≤z δ ρ(z ) in units of e for EF in the range 19 − 25meV, corresponding to the STAI phase. The number of filled electron states is indicated in the legend (5000 represents the halffilling). The steps of magnitude ±e/2 at the monopole/anti-monopole locations show the expected localized fractional charge due to the Witten effect. Steps with magnitude ±e located elsewhere are due to the shift of some bound states in applied magnetic field and can be viewed as a finite-size effect which should diminish for a system with a larger cross section. In both panels we use U0 = 150meV and other parameters as in Fig. 6.3.  91  Chapter 7  Topological Insulator with Magnetic Impurities 7.1  Overview  In this section, we propose that a temperature window exists in which the surface of a magnetically doped TI is magnetically ordered but the bulk is not. We present a simple and intuitive argument why this is so, and back it up by a meanfield calculation for two simple tight binding TI models: a cubic-lattice regularized Bi2 Se3 and a model on the Perovskite lattice. Our results [99] show that indeed a sizeable regime such as described above could exist in real TIs, and this indicates a possible physical explanation for the results seen in experiments [50, 51].  7.2  Experimental results  A signature of the massive Dirac fermion has been observed recently using ARPES in magnetically doped Bi2 Se3 [50, 51], although the interesting effects associated with it have yet to be seen in a laboratory. A surprising feature of these experiments is that the gap in the surface spectrum appears without bulk magnetic ordering, even though the dopants are uniformly distributed everywhere in the 3D sample. These findings raise several important questions concerning the precise conditions under which TRI-breaking perturbations open up a gap. Can a gap open in the surface 92  state of a TI in a TRI-broken phase which however lacks global magnetic ordering? Although we know of no systematic study of this problem, simple arguments suggest that unordered magnetic moments do not open a gap. Consider creating such a disordered state from a uniform 2D FM in the surface of a TI by introducing domains with opposite magnetization (taken to point in the direction perpendicular to the surface). It is well known that the resulting domain walls carry topologically protected gapless fermionic modes [100]. As the number of the domains grows so does the density of the low-energy fermion modes, ultimately presumably recovering the 2D gapless state characteristic of the system with unbroken TRI. The above argument thus suggests that uniform magnetic ordering over large domains is necessary to gap out the surface modes in a TI.  7.3 7.3.1  Surface magnetic ordering in topological insulators Surface doping  The most natural way to attempt to open up a gap in the surface state of a TI is to coat the surface with a ferromagnetic material, with magnetization perpendicular to the surface. Theoretically, this causes a gap to open up, proportional to the magnetization of the FM coating [29, 100]. To illustrate this point, consider the effective low energy Hamiltonian for electrons on the surface of a 3D TI that lies parallel to the x-y plane H0 = v(kx σy − ky σx )  (7.1)  where v is the Fermi velocity, and σi are Pauli matrices in the spin subspace. If we coat this surface with a ferromagnetic coating with magnetization M = M zˆ then we get an additional term in the Hamiltonian H = H0 + JM  93  σz , 2  (7.2)  where J is the exchange coupling strength. Since H is a sum of anti-commuting matrices we can write the spectrum down immediately Ek = ± v2 k2 + (JM/2)2 ,  (7.3)  where k2 = kx2 + ky2 . We see that a gap of size JM has opened up. However, even though from a theoretical standpoint this proposition seems promising, experimentally it has proven very difficult to fabricate a sample with the requisite properties. Two key challenges need to be overcome: first, to observe most of the interesting surface phenomena, one requires the surface to remain insulating; however most ferromagnets in nature are metallic. Second, for a ferromagnet in a thin-film geometry, the magnetization vector usually lies in the plane, whereas a perpendicular magnetization is required to open up a gap in the TI surface state. To the best of our knowledge, this has yet to be achieved in an experiment, although some theoretical work has been done on this topic [101].  7.3.2  Bulk doping  If surface doping with magnetic impurities fails, it is natural to try bulk doping. In ARPES experiments[50, 51] it was found that doping the bulk with non-magnetic impurities (such as Ca, Sn and Tl) did not result in a gap in the Dirac cone, as expected since they do not break TRI. Conversely, doping with magnetic impurities, for example Bi2−x Fex Se3 , resulted in a spectral gap that increased with the concentration of magnetic dopants x, with a gap of 60meV for x = 0.25 (the bulk gap for Bi2 Se3 is ∼ 0.3eV). For the magnetic dopants Fe and Mn it was found that, at least  for small x, the bulk was paramagnetic, while for the undoped samples the bulk was found to be diamagnetic. The magnetization measurements were not sensitive to the surface. This raises the question of magnetic ordering in the bulk versus the surface. In general, ordered phenomena in lower dimensions are more fragile (Tc3D > Tc2D ), for example, the XY model and Heisenberg models - in 1D they do not order at any temperature, in 2D they order only for T = 0 and in 3D they order for T < Tc . This is also the case for superconducting order and general stability of lattices. 94  However, in the case of a TI, we will argue that it is possible that Tcbulk < Tcsurf (for another example of this, see Ref. [102, 103]). Therefore, there is a regime Tcbulk < T < Tcsurf in which the bulk is unordered (paramagnetic) and the surface is ordered (for example, ferromagnetic). To illustrate why this could be the case, first recall that magnetic ordering with the magnetization perpendicular to the TI surface implies opening of a gap in the spectrum of the surface states. This can be seen directly from Eqs. (7.2,7.3). Now consider the ungapped surface spectrum, assuming half filling, so that the surface valence band is fully occupied and the conduction band is empty (Fig. 1.5a). Gapping the surface states causes the occupied states to move down in energy (Fig. 1.5b), so the total kinetic energy decreases. Therefore, the formation of a surface gap is favourable. If the chemical potential is shifted either up or down then the net gain in kinetic energy is diminished and we expect Tcsurf to decrease. Contrast this with the situation in the insulating bulk, which is gapped to begin with. In an ordinary insulator with negligible spin-orbit coupling it is not possible to generate magnetization in the initially spin-degenerate bands without first closing the gap. Equivalently, one may recall that the spin susceptibility of an ordinary insulator with a negligible spin-orbit coupling vanishes. In the bulk of a topological insulator the situation is more complicated as a result of the strong spin-orbit coupling that is necessary for the occurrence of the topological phase. In this case, magnetic susceptibility can be significant [89], and can lead to bulk magnetic states at non-zero temperatures. Nevertheless, in this study we find that generically the surface critical temperature for the formation of magnetic order exceeds the bulk critical temperature.  7.4  Determining parameters for Bi2 Se3  Several papers have considered a Hamiltonian which is a so called “discretized” model of Bi2 Se3 . We understand the discretization here as projecting the Bi2 Se3 lattice on a tetragonal (elongated cubic) lattice. The tetragonal model’s low energy approximation agrees with the low energy approximation of the Bi2 Se3 model, and in this respect it is an approximation of it. We use it since it is significantly simpler, and we expect the physical properties of the Bi2 Se3 model to be well captured by  95  this approximation although specific details may vary. Here we briefly survey the different Hamiltonians in the literature, some of which appear to be incorrect. We note also that the tetragonal model is in fact the same as the Hamiltonian presented in Chapter 3.2 up to some additional terms, and the different lattice constant and hopping parameters for the z axis and the x − y plane.  Liang Fu alerted us to two errors in [34], one of the first theoretical papers on  Bi2 Se3 . The first is that the kz term of that Hamiltonian can be written (in our notation) as A1 kz (σz ⊗ τx ). This term is odd under mirror symmetry about the x-  axis, since kz → kz (unaffected), τx → τx (orbital unaffected) and σz → −σz (for  spin). Therefore, this term breaks the crystal symmetry and cannot be correct. In addition, the spin term has the form (kx , ky ) · (σx , σy ), which should be replaced by  (kx , ky ) · (σy , −σx ), as has been observed (spin tangent to Dirac cone).  The same group published a later paper [43] that corrected the first mistake but  not the second. Unfortunately, the Hamiltonian from the first paper has propagated with its mistakes. The correct Hamiltonian appears in [42]. Our attempt To determine realistic parameters for Bi2 Se3 we compare the low energy approximation of our model, with the low energy approximation from the literature [34]   M(k)  0  A1 kz  A 2 k−       A1 kz −M(k) A2 k− 0    + O(k2 ) H(k) = ε0 (k) +   0 A k M( k) −A k 2 + 1 z  A2 k+ 0 −A2 kz −M(k)  (7.4)  2 and M(k) = M − B k2 − B k2 . The with k± = kx ± iky , ε0 (k) = C + D1 kz2 + D2 k⊥ 1 z 2 ⊥  time reversal and inversion operators are given by T = K iσy ⊗ I and P = I ⊗ τz .  We notice two minor differences with our model: the orbital and spin subspaces appear in reverse order, and the inversion operator is different (ours contains τx ). Note that a slightly different representation appears in a more recent paper  96  ([89], see supplement)   M(k)  0  Axy k−  A z kz      Axy k+ −M(k) −Az kz  0  + O(k2 ) H(k) = ε0 (k) +   0  −A k M( k) A k z z xy +   Az kz 0 Axy k− −M(k)  (7.5)  which can be obtained from the previous matrix by swapping rows and columns 2 + D k2 and M(k) = M + B k2 + B k2 . However, (2 ↔ 4), and ε0 (k) = C + Dxy k⊥ z z z z 0 xy ⊥  they also have M0 , Bz ≡ B1 and Bxy ≡ B2 that appear with the opposite sign. They  claim that changing the sign of M(k) does not affect any observables, which appears to be true. A more significant difference is that the parameters Bz and Bxy (and similarly for D and A) are switched. Since the parameters of the first model appear to give an ordinary insulator, we choose to work with the parameters from the second paper, and assume for now that the first paper contains some typos in the parameters. Therefore, we can go from our model to their model (the first representation above) by swapping subspaces and then applying the unitary transformation U = e−i 4 σ2 σx π  (7.6)  So that our model gives   d0   −d1 − id2 H(k) = d4 − µ +   0  −d3  −d1 + id2  0  −d0  d3  d3  d0  0  −d1 + id2  −d3       −d1 − id2   −d0 0  (7.7)  We now look at the low energy approximation around the Γ point (0, 0, 0) sin(ki ai )  (7.8)  ki ai  cos(ki ai )  1−  97  (ki ai 2  )2  so that d0 = ε − 2 ∑ ti cos(ki ai ) i  ε − 2 ∑ ti + ∑ ti (ki ai )2 γ0 − 2 ∑ γi + ∑ γi (ki ai )2  d4 = γ0 − 2 ∑ γi cos(ki ai )  i  i  di = −2λi sin(ki ai )  (7.9)  i  i  i  −2λi ki ai = −vi ki  By comparing we find that the x and y axes are equivalent, so we set γx = γy ≡ γ⊥  and similarly for t and λ . We also set ax = ay = a and az = c. We find the following four independent equations ε0 (k) = d4  (7.10)  M(k) = d0 Az kz = −d3  Axy k− = −d1 + id2 which yield simple equations for the parameters of our model. The first equation gives γ0 = C + 2 γz = γ⊥ =  Dxy Dz +4 2 2 c a  (7.11)  Bxy Bz +4 2 c2 a  (7.12)  Dz c2 Dxy a2  The second equation gives ε = M0 + 2 tz = t⊥ =  Bz c2 Bxy a2  98  and the third and fourth equations give Az 2c Axy 2a  λz = λ⊥ =  (7.13)  The parameters were obtained in [34] by fitting to an ab initio calculation of the ˚ Az = spectrum of Bi2 Se3 and are given by [89] as M0 = −0.28 eV, Axy = 2.2 eVA,  ˚ Bxy = 10 eVA˚ 2 , Bz = 56.6 eVA˚ 2 , C = −0.0068 eV, Dxy = 1.3 eVA˚ 2 , Dz = 4.1 eVA, 19.6 eVA˚ 2 . In order to determine the lattice constants a and c we need two equations to equate the unit cell of our model (cubic) with the unit cell of Bi2 Se3 . We choose to equate the volume and the ratio between the base area and the height. The lattice constants of our model are labeled a and c, and the lattice constants of Bi2 Se3 (obtained experimentally) are labelled a and c . The two equations are √ 3 2 a c 2√ 3 2 2 3a c  a c = 2  a2 c  =  (7.14)  and the solutions 3 a 2  a =  (7.15)  c √ 3  c = For the lattice constants of Bi2 Se3 we use a  = 4.1387A˚  c  = 28.629A˚  99  (7.16)  x 0 0.0023 0.0044 0 0.03  ˚ a (A) 4.1387 4.1381 4.1365 4.1396 4.106  ˚ c (A) 28.629 28.627 28.625 28.636 28.562  Table 7.1: Lattice constants for Bi2 Se3  So that we have a = 5.0689A˚  (7.17)  c = 16.5290A˚ The lattice constants we found in the literature are tabulated in Table 7.1. These results were obtained from X-ray powder diffraction in two separate experiments [104, 105]. The researchers observe the decrease in lattice constants and argue that this is due to the smaller atomic radius of Mn compared with Bi. In addition, they argue that the decrease argues strongly that the Mn atoms replace the Bi atoms, as opposed to being located interstitially.  7.5  2D surface mean field calculation  This section follows notes written by my supervisor, Prof. Marcel Franz. Consider the 2D surface of a 3D TI. We define the axes so that the 2D surface is positioned perpendicular to the z axis. The effective Hamiltonian for electrons on the surface is He (k) = vk · σ  (7.18)  where v is the Fermi velocity, which is assumed the same for the x and y directions (as is reasonable for Bi2 Se3 ), k = (kx , ky ), and σ = (σy , −σx ) are Pauli matrices in  the spin subspace.  100  On substitution into (7.46) we find HeMF (k) = He (k) + JM  σz JM − µ = v(kx σy − ky σx ) + σz − µ 2 2  (7.19)  so that the energies are Ee = α  v2 k2 + M˜ 2  (7.20)  where k2 = kx2 + ky2 , M˜ = JM/2 and α = ±1.  We minimize the free energy with respect to m and M. minimizing by m does  not involve the electron Hamiltonian (HeMF does not depend on m), and we find M = −2SxBS (β JmS)  (7.21)  To minimize with respect to M we take the derivative of the energy ∂ Ee ∂M  ∂ v2 k2 + M˜ 2 ∂M MJ 2 2 2 ˜ 2 v k +M = α 4 = α  where we have defined uk ≡  (7.22) − 21  =  MJ 2 α 4 uk  v2 k2 + M˜ 2 , so that Ee = αuk − µ  Substitute into (7.64) to find  ∂ Ee (k, α) 1 f (Ee ) NJ ∑ ∂M kα JM α = f (αuk − µ) ∑ 4N kα uk JM 1 = − [1 − f (uk + µ)) − f (uk − µ)] ∑ 4N k uk  m =  (7.23)  In the limit µ = 0, which is the case for half filling at T = 0, we find tanh( 2 k ) JM 1 − 2 f (uk ) JM m=− =− ∑ ∑ 4N k uk 4N k uk βu  (7.24)  We now investigate our result near the transition with the objective of finding 101  the critical temperature Tc . In this limit m → 0 as we take T → Tc so that β JmS  1.  We can approximate the Brillouin function for small arguments, as done in the next section (see (7.88)), and we find −  m  3 M 2S(S + 1) β Jx  (7.25)  Which we can now substitute into the expression for m (in the limit µ = 0 for all T ), to find tanh( 2 k ) 3 M JM =− 2S(S + 1) β Jx 4N ∑ uk k βu  −  (7.26)  After dividing by M and an additional factor on both sides, we find an equation for Tc 1=  tanh( 21 βc v|k|) S(S + 1) 1 βc J 2 x ∑ 6 N k v|k|  where we have set M = 0 in the energies uk ≡  v2 k2 + M˜ 2  (7.27)  v|k|, since this term  would just give higher order terms in M (we expanded both sides to first order in M). We integrate tanh( 21 βc v|k|) 1 N∑ v|k| k  = 2π = = =  a 2π a 2π a 2π  a 2π  2  Λ  k dk 0  tanh( 12 βc vk) vk  (7.28)  1  2 β vΛ 4π dz tanh z βc v2 0 2 4π 1 ln cosh βc vΛ βc v2 2 2 2π β vΛ βc v2  2  In the last step we assume that the cutoff is large βc vΛ  1 so that ln cosh x  x.  Substituting back into the gap equation, we find Tcsur f  S(S + 1) π 3kB 102  Λa 2π  2  J2x Λv  (7.29)  For S = 5/2, J = 0.5eV, x = 0.05, v = 2λ⊥ a and a cutoff Λa/(2π)  Tcsurf  1/5 we find  100K. This equation is identical to a result in [101] (middle of 3rd page)  and in [106] (equation 19). Note that mean field treatments usually overestimate Tc so in reality we might expect a significantly lower Tc . Of course, this result also depends strongly (linearly) on the cutoff. A simple result can also be found for the case β µ  1 which is relevant for  most values of µ inside the bulk gap, Tcsurf  S(S + 1) π 3kB  Λa 2π  2  J2x (Λv − |µ|). (Λv)2  (7.30)  As expected, Tcsurf is seen to decrease away from half filling. To complete the argument we would now like to give a similar simple estimate for Tcbulk . Unfortunately, the bulk critical temperature is not easy to estimate, since unlike the topologically protected surface state, whose physics is simple and universal, the bulk of a TI can be complicated and Tcbulk will generally depend on the details of the band structure and other factors. For Bi2 Se3 with 5% concentration of Cr dopants, Ref. [89] estimates Tcbulk  70K using first-principles numerical cal-  culations. Comparing with our rough estimate for Tcsurf given above we see a clear indication that a Tcbulk < Tcsurf regime can be easily obtained.  7.6  3D bulk mean field calculation  We consider a 3D topological insulator with magnetic impurities at random sites. The idea is to find a simple model for the material Bi2−x Mnx Se3 (and similar materials), which has been experimented on recently. We are interested in finding the dependence of the magnetization on temperature. We assume that the density of impurities is low enough that we can neglect impurity-impurity interactions and that there is no clustering. This is supported in experiments on doped Bi2 Te3 [107] (see FIG 11 in that paper). Therefore, we add only an on-site electron-impurity spin-spin interaction term H = ∑ Ψ†k He (k) − µ Ψk + J ∑ SI · sI I  k  103  (7.31)  where He is the Hamiltonian for the electrons, J is a coupling constant, the sum runs over all impurity sites I, SI is the spin of the impurity on site I, sI ≡ Ψ†I σ2 ΨI  is the spin of the electrons on site I and Ψ†I is a four component vector containing the creation operators (to be defined explicitly below). We rewrite the interaction term, as is usual for a mean field treatment  ∑ SI · s  I  =  ∑ ηi Si · s  (7.32)  i  i  I  =  ∑(ηi Si − M) · (si − m) + ∑ M · si + ∑ ηi Si · m − ∑ M · m i  i  i  i  where  1 i ∈ I ηi ≡ 0 i ∈ /I  (7.33)  and the sum over i runs over all sites, M and m are the average magnetization of the impurities and electrons (respectively), which do not depend on the spatial coordinates, in this treatment. The impurity magnetization is averaged over all impurities and over all sites, often referred to as the virtual crystal approximation (VCA) [108]. This is like “spreading out” the total spin of the impurities, so that there is some impurity spin on all sites, not only impurity sites, so M≡  ∑I SI N  (7.34)  The quantity of interest here is the average impurity spin, given by S ≡  ∑I SI NI  (7.35)  Similarly, the electron magnetization m is averaged over all electrons and all space, in other words the total spin of the electrons is spread out evenly on each site m=  ∑i si N  104  (7.36)  The quantity of interest is the average electron spin, given by s =  m ∑i si = 2N 2  (7.37)  where we have assumed half filling, so that 2N electrons fill half the available states. In a mean field calculation, we assume that the fluctuations around the mean are small, and hence neglect the term which is second order in the fluctuations. We are left with J ∑ SI · sI I  J M · ∑ si + Jm · ∑ SI − NJ M · m i  (7.38)  I  where N is the number of sites. The last equation shows that the interaction term has now been decoupled. We started with a term that coupled the spin of the electrons to the spin of the impurities. Now we have three terms: the first involves only the spin on the impurities, the second involves only the spin of the electrons, and the third depends on neither - it depends only on the fields. Therefore we can rewrite our Hamiltonian, in its approximate mean-field form, as H MF = ∑ Ψ†k He (k) + J M · k  σ − µ Ψk + Jm · ∑ SI − NJ M · m 2 I  (7.39)  where we have used σ  σ  ∑ si = ∑ Ψ†i 2 Ψi = ∑ Ψ†k 2 Ψk i  i  (7.40)  k  Conside the case in which the magnetization is in the z direction, so M = M zˆ and m = mˆz. We have H MF  =  ∑ Ψ†k k  He (k) + JM  σz − µ Ψk + Jm ∑ SIz − NJMm 2 I  = HeMF + HIMF − NJMm  105  (7.41)  The electron mean field Hamiltonian is just the Hamiltonian for electrons in the “clean” system (no impurities) with an additional term. The “clean” electron Hamiltonian is given by He (k) = d4 + ∑ dµ γµ  (7.42)  µ  where we have chosen the gamma matrices so that γ0 = τ1 ⊗σ0 , γ1 = −τ3 ⊗σ2 , γ2 =  τ3 ⊗ σ1 , γ3 = τ2 ⊗ σ0 (so that γ4 = τ3 ⊗ σ3 and define d0 ≡ mk = ε − 2 ∑i ti cos(ki ai ),  di ≡ −2λi sin(ki ai ), d4 = γ0 − 2 ∑i γi cos(ki ai ). This is the same model presented in Section 3.2 but with an additional hopping term γ as in Chapter 6, allowing for  anisotropic hopping, as expected for example in Bi2 Te3 and Bi2 Se3 , and with two corrections based on Liang Fu’s comments (see Section 7.4). We can now write the energies for the mean field Hamiltonian E(k, λI ) = Ee + EI − NJMm  (7.43)  The first term is the energies of the electron Hamiltonian, which are computed below. The second term involves a sum over the expectation value of the spin operator in the z direction at each impurity site EI = Jm ∑ λI  (7.44)  I  where λI is the component of the spin in the z direction on impurity I. The last term is a scalar, so it gives a shift in energy. We now diagonalize the mean field Hamiltonian for the electrons, which is given by σz − µ Ψk = ∑ Ψ†k HeMF (k)Ψk 2  (7.45)  σz σz − µ = d4 + ∑ dµ γµ + JM − µ 2 2 µ  (7.46)  HeMF = ∑ Ψ†k He (k) + JM k  k  and we can write explicitly HeMF (k) = He (k) + JM  106  and the energies are Ee = d4 + γ  ∑ dµ2 + M˜ 2 + 2δ M˜ µ  d02 + d32 − µ  (7.47)  2  = d4 + γ  d12 + d22 +  d02 + d32 + δ M˜  −µ  where M˜ ≡ JM/2, and γ, δ = ±1. We see that if we take away the impurities by setting J = 0 we get the “clean” doubly degenerate electron spectrum Ee = d4 ±  ∑ dµ2 − µ  (7.48)  µ  The two degenerate energies at the Γ point, for J = 0 are Ee (k = 0) = d4 + γ|m0 | − µ  (7.49)  So that the bulk gap is 2|m0 |, where m0 = ε − 4t⊥ − 2tz . For J = 0 we find that the ˜ = |JM| degeneracy is broken by a splitting of 2|M| ˜ −µ Ee (k = 0) = d4 + γ(|m0 | + δ M)  (7.50)  ˜ and the bulk gap is smaller 2|(|m0 | − |M|)|.  We want to find equations for M and m. We do this by minimizing the free  energy F with respect to both M and m. Based on the energy, we can write the free energy as F = Fe + FI − NJMm  (7.51)  In general, F =−  1 ln Z β  (7.52)  where Z is the many particle partition function, β = 1/(kB T ) and kB is the Boltz-  107  5 4 3 2  E  1 0 ï1 ï2 ï3 ï4 Γ  X  Γ  M  R  M  X  R  Figure 7.1: Bulk energies at high symmetry points for S = 5/2, x = 0.05, T =0 mann factor. For the electrons the partition function is given by Ze = ∏ Zkδ γ = ∏  ∑  kδ γ ni =0,1  kδ γ  e−β ni Ee (δ ,γ,k) = ∏ 1 + e−β Ee (δ ,γ,k)  (7.53)  kδ γ  So that the electron free energy is Fe = −  1 1 ln Ze = − β β  ∑ ln  1 + e−β Ee (δ ,γ,k)  (7.54)  kδ γ  The partition function for the NI impurities is ZI =  NI  s  ∑  −β Jmλ  e  (7.55)  λ =−s  so that the impurity free energy is FI = −  s NI 1 ln ZI = − ln ∑ e−β Jmλ β β λ =−s  (7.56)  We now minimize the free energy by requiring the derivative with respect to M and 108  m to be zero. Start with m ∂F ∂ FI = − NJM = 0 ∂m ∂m  (7.57)  So M=  1 ∂ FI NJ ∂ m  =  S NI 1 ∂ − ln ∑ e−β Jmλ NJ ∂ m β λ =−S  = −2x  (7.58)  S ∂ ln ∑ e−αλ ∂α λ =−S  where α ≡ β Jm so that ∂ /∂ m = β J∂ /∂ α and x = NI /(2N) is the ratio of number  of impurities to number of Bismuth atoms. This expression occurs frequently in calculations involving magnetic moments (see for example [109], page 655). We evaluate the sum separately S  ∑  e−αλ  =  λ =−S  = =  e−αS 1 + eα + e2α + · · · + e2αS  (7.59)  e(S+ 2 )α − e−α(S+ 2 ) e(2S+1)α = e α α 1 − eα e 2 − e− 2 sinh[(S + 12 )α] sinh α2 1  −αS 1 −  1  Substitute into the magnetization, to find M = −2x  sinh[(S + 21 )α] ∂ ln ∂α sinh α2  = −x (2S + 1) coth (2S + 1)  (7.60) α α − coth = −2SxBS (αS) 2 2  where BS (y) is the Brillouin function BS (y) =  2S + 1 2S + 1 1 1 coth y − coth y 2S 2S 2S 2S  109  (7.61)  So that our final equation for M is M = −2SxBS (β JmS)  (7.62)  In the zero temperature limit kB T → 0 we expect the impurity magnetic moments  to be aligned and at their extremum value ±S, so that the total magnetization is  ±NI S/N = ±2Sx. Indeed, in this limit we get BS → 1 (assuming m > 0) so that  −2Sx. In the limit of infinite temperature kB T → ∞ we expect the impurity  M  moments to be random and hence their average magnetization should be zero. In this limit BS → 0 and we indeed find M → 0. Similarly, for M  ∂F ∂ Fe = − NJm = 0 ∂M ∂M  (7.63)  So m=  1 ∂ Fe NJ ∂ M  =  1 ∂ 1 − NJ ∂ M β  ∑ ln  1 + e−β Ee (k,δ ,γ)  (7.64)  kδ γ  =  e−β Ee (k,δ ,γ) ∂ Ee (k, δ , γ) 1 ∑ 1 + e−β Ee (k,δ ,γ) ∂ M NJ kδ γ  =  1 ∂ Ee (k, δ , γ) f (Ee ) ∑ NJ kδ γ ∂M  where f (E) = (1 + eβ E )−1 is the Fermi-Dirac distribution. We take the derivative of the energy ∂ Ee ∂M  = γ =  ∂ ∂M  δ γJ 2  2  d12 + d22 +  d02 + d32 + δ M˜  d02 + d32 + δ M˜  d12 + d22 +  110  (7.65) d02 + d32 + δ M˜  2 −1/2  =  δ γJ εkδ 2 ukδ  where we have defined εkδ  =  d02 + d32 + δ M˜  ukδ  =  2 d12 + d22 + εkδ  (7.66)  So that Ee = γukδ − µ˜  (7.67)  and µ˜ ≡ µ − d4 . We substitute into the magnetization m to find m = =  1 ε δ γ kδ f (Ee ) ∑ 2N kδ γ ukδ 1 ε ˜ − f (−ukδ − µ)] ˜ δ kδ [ f (ukδ − µ) ∑ 2N kδ ukδ  = − =  (7.68)  1 ε ˜ − f (ukδ + µ)] ˜ δ kδ [1 − f (ukδ − µ) ∑ 2N kδ ukδ  1 εk− ε ˜ − f (uk− − µ)] ˜ − k+ [1 − f (uk+ + µ) ˜ − f (uk+ − µ)] ˜ [1 − f (uk− + µ) ∑ 2N k uk− uk+  where we have used f (−E) =  1 1 = 1− = 1 − f (E) −β E 1+ e 1 + eβ E  (7.69)  In the limit µ = 0 and γ0 = γi = 0 we find µ˜ = 0 so that m = − = − =  1 ε δ kδ [1 − 2 f (ukδ )] 2N ∑ u kδ kδ  (7.70)  1 ε β ukδ δ kδ tanh 2N ∑ u 2 kδ kδ  1 εk− β uk− tanh ∑ 2N k uk− 2  111  −  εk+ β uk+ tanh uk+ 2  where we have used 1 − 2 f (x) = 1 −  2 eβ x − 1 βx = = tanh β ε β x 2 1+ e e +1  (7.71)  The chemical potential µ can be determined by summing over the average occupation for each energy (given by the Fermi-Dirac distribution) and equating to the total number of states. Since in this case the system is infinite, it is more appropriate to discuss the density n=  1 1 ˜ f (Ee ) = ∑ f (γukδ − µ) ∑ N kδ γ N kδ γ  (7.72)  For µ = 0 (half filling) at T = 0 we have n = 2 since there are four states per site, but only half of those are occupied. The equations for m and M are coupled non-linear equations which can be solved self consistently by an iterative procedure. The algorithm is 1. For each temperature T 2. Initialize M and m to those obtained for previous T (or otherwise to a rough guess of the values) 3. While we haven’t reached self consistency (a) Calculate M from (7.62) (b) Calculate m from (7.68) (c) Calculate µ by solving (7.72) for µ Alternatively and equivalently, the value of m can be written as a sum over expectation values of the spin operator m = ∑ Ψi | k,i  σz σz |Ψi = ∑(U † U)ii f (Ei ) 2 2  (7.73)  k,i  where U is the matrix that diagonalizes the Hamiltonian. In the z dependent case (next section), it is more convenient to use this formalism. The reason is that in that 112  2.5  S (µ B )  2  1.5  1  0.5  0 0  10  20  30  40  50  T (K)  60  70  80  90  Figure 7.2: The magnetization of the impurities as a function of temperature for S = 5/2, x = 0.05 case we cannot diagonalize the Hamiltonian by hand, and therefore do not have a closed form for the energies. Taking the required derivative of the energies would then involve a numerical derivative which would be messy. To test whether the physics of our model is dominated by the TRIM, we compute M(T ) by summing only over a subset of the BZ shaped like a cube and centred on the Γ point (see Fig. 7.3). We find that for reasonable cutoffs (small compared with the BZ), the agreement is not good. Seeing as we have set our parameters by the Bi2 Se3 low energy model (which should be correct near the Γ point), we can expect our model to agree with this model and other models for this regime. However, there is no reason to believe that our model would agree with that model for the rest of the Brillouin zone (although it is possible). In particular, if we believe that the physics is dominated by the TRIM, then our model will grossly overestimate Tc . It is interesting to look at two limits: low temperatures, and close to the transition. We start with the first limit. At T = 0 we find M = −2Sx (assuming M < 0) so that M˜ = −JSx. We assume d4 = 0 for simplicity. The energies are Ee = γukδ .  113  5  4 8 12 16 20  4.5 4  M/x(µ B )  3.5 3 2.5 2 1.5 1 0.5 0 0  10  20  30  40  50  T (K)  60  70  80  90  100  Figure 7.3: M(T) as a function of the ”side of cube” in k-space, for a cube centred around the Γ point. Here L=20, so the ”side of cube” 20 result is the full result. The other parameters are J=0.5, S=5/2, x=0.05. At T = 0, for half filling, only the γ = −1 states will be filled, so m(T = 0) =  1 ε ∑ δ γ ukδkδ f (Ee ) 2N kδ γ  = −  (7.74)  1 ε δ kδ ∑ 2N kδ ukδ  Now consider a low energy approximation around the Γ point εkδ  =  d02 + d32 + δ M˜  m20 + v2z kz2 − δ JSx  ukδ  =  2 d12 + d22 + εkδ  2 + ε2 v2⊥ k⊥ kδ  (7.75)  Notice that εkδ depends on kz but not on kx and ky . We impose a cutoff and convert the sums to integrals. It is most convenient to first do the integral over the kx -ky  114  plane and then the integral over the kz axis, so we split the sums m(T = 0) = −  1 εkδ 1 δ ∑ ∑ Nz k ,δ 2N⊥ k⊥ ukδ  (7.76)  z  To check the factors and sign we can take the “surface limit” Nz = 1, kz = 0, δ = 1, m0 = 0 (no gap on the surface) whereupon we find that εk = −JSx and uk = v2 k2 + εk2 , so that  msur f (T = 0)  1  JSx 2N ∑ k  v2 k2 + (JSx)2  (7.77)  Consider first just the term in square parenthesis. We substitute the low energy approximation, define Ω = εkδ /v⊥ , convert the sum into an integral with a cutoff,  evaluate the integral, and assume that the cutoff is the largest energy scale in the problem 1 εkδ ∑ 2N⊥ k⊥ ukδ  1 2N⊥ ∑ k⊥  εkδ 2 + ε2 v2⊥ k⊥ kδ 2  Ω a 2 2π πΩ  Λ⊥  2πk⊥ dk⊥  2 + Ω2 k⊥  0  a 2π  (7.78)  2  Λ⊥  where the integral is the same as before, where it was evaluated explicitly. We plug this back into the expression for m to find m(T = 0) = −  π a v⊥ 2π  115  2  Λ⊥ ∑ δ δ  1 εk Nz ∑ kz  (7.79)  Once again, consider just the expression in the square brackets 1 εkδ Nz ∑ kz  c 2π  Λz 0  m20 + (vz kz )2 − δ JSx  dkz   c  vz 2π  =  m0 vz  Λz  dkz  0  (7.80)   2  + kz2 − δ JSxΛz   The integral gives  0  2  m0 vz  Λz  dkz  + kz2  1 2 Λ 2 z  (7.81)  so that we find 1 εkδ Nz ∑ kz  c 1 vz Λ2z − δ JSxΛz 2π 2  (7.82)  So m(T = 0)  −  π a v⊥ 2π  2π  JSx Λ⊥ v⊥  2  Λ⊥ ∑ δ δ  Λ⊥ a 2π  c 1 vz Λ2z − δ JSxΛz 2π 2 2  (7.83)  Λz c 2π  This result has a quadratic dependence on the cutoff, as expected. Perhaps surprisingly, it doesn’t depend on vz . We now investigate our result near the transition with the objective of finding the critical temperature Tc . In this limit m → 0 as we take T → Tc so that β JmS We do this by expanding both sides of the equation for m to first order in M˜ m =  ε 1 δ γ kδ f (Ee ) ∑ 2N kδ γ ukδ  1.  (7.84)  We start with the left hand side. Consider the limit of the Brillouin function for  116  small arguments. Since 1 x + + O(x3 ) x 3  coth x  (7.85)  we can approximate the Brillouin function to first order BS (y) =  2S + 1 2S + 1 1 1 coth y − coth y 2S 2S 2S 2S  S+1 y 3S  (7.86)  So that M = −2SxBS (β JmS)  −  2S(S + 1) β Jmx 3  (7.87)  and therefore m  −  M˜ 3 S(S + 1) β J 2 x  (7.88)  For the right hand side, focus on the argument of the sum. To find the first order approximation, we take the first derivative with respect to M˜ and take M˜ = 0 (higher order terms in M˜ disappear in this procedure) ∂ε ∂ f εkδ ∂ ukδ εkδ ∂ εkδ f f = kδ + − f ˜ ˜ u u ∂ M kδ ∂ M kδ ∂ M˜ ukδ ∂ M˜ u2kδ  (7.89)  We take the required derivatives, assuming M˜ > 0 ∂ εkδ = ∂ M˜  d02 + d32 + δ  (7.90)  and ∂ ukδ ∂u ε ε = kδ kδ = kδ ∂ εkδ ∂ M˜ ukδ ∂ M˜  d02 + d32 + δ  (7.91)  and finally ∂ f ∂ Ee ε ∂f = = −γβ f 2 eβ Ee kδ ∂ Ee ∂ M˜ ukδ ∂ M˜  117  d02 + d32 + δ  (7.92)  We substitute the last three equations into the previous equation, to find ∂ εkδ f = ∂ M˜ ukδ  d02 + d32 + δ  f ukδ  1 − γβ f eβ Ee  2 εkδ ε2 − kδ ukδ u2kδ  (7.93)  Substitute this equation and the derivative with respect to M˜ of (7.88) into (7.84) to find a closed equation for Tc 1= −βc J 2 x  (7.94) S(S + 1) 1 ∑ δγ 6 N kδ γ  d02 + d32 + δ  f ukδ  1 − γβc f e  2 β Ee εkδ  ukδ  −  2 εkδ u2kδ  ˜ M=0  We can take the “surface limit” to check this result. This involves taking all hopping parameters along the z direction to zero, taking µ˜ = 0, and taking the low energy approximation. This gives d12 + d22 = v2 k2 , d0 = d3 = d4 = 0. So δ = 1, ˜ uk = v2 k2 + M˜ 2 and therefore the energies are Ee = γ v2 k2 + M˜ 2 . On εk = M, substitution into the previous equation, we find 1= −βc J 2 x  (7.95) S(S + 1) 1 γ 6 N∑ kγ  f v2 k2 + M˜ 2  1 − γβc f eβ Ee  M˜ 2 v2 k2 + M˜ 2  −  M˜ 2 v2 k2 + M˜ 2  where the sum over k is now over a circle defined by |k| < Λ (instead of over the whole Brillouin Zone). On taking the M˜ = 0 limit we find 1 = −βc J 2 x = βc J 2 x  f (γv|k|) S(S + 1) 1 γ ∑ 6 N kγ v|k|  (7.96)  tanh 12 βc v|k| S(S + 1) 1 ∑ 6 N k v|k|  where we have used 1  ∑ γ f (γv|k|) = f (v|k|) − f (−v|k|) = 2 f (v|k|) − 1 = − tanh( 2 β v|k|)  (7.97)  γ  The last equation is identical to the equation derived for the surface case (7.27). 118  ˜ M=0  7.7  Adding the surfaces  Consider now a system with open boundary conditions in the z direction and periodic boundary conditions in the x and y directions (as before). The advantage of this setup is that in addition to the bulk, we can also investigate the magnetization near the edges of the sample. In the previous section, we looked at the bulk, and had periodic boundary conditions in all directions. Therefore, it was useful to Fourier transform in all three directions. We now take the real space Hamiltonian, and Fourier transform it in x and y, keeping the dependence on z in real space. The electron Hamiltonian is He = HSO + Hcd + Hγ where HSO = i ∑ λµ Ψ†j τz σµ Ψ j+µ + h.c.  (7.98)  j,µ  Hcd = ε ∑ Ψ†j τx Ψ j − j  Hγ  = γ0 ∑ j  Ψ†j Ψ j −  ∑  tµ Ψ†i τx Ψ j + h.c.  <i, j>  ∑  γµ Ψ†i Ψ j + h.c.  <i, j>  where Ψ j = (c j↑ , c j↓ , d j↑ , d j↓ ), and τµ and σµ are Pauli matrices in orbital and spin subspaces (respectively). As preparation for the next stage, we define ϕ = (Ψ1 , Ψ2 , . . . ΨLz ), a vector of length 4Lz and write the Hamiltonian in the expanded 4Lz space, rewriting it by separating into an in-plane part and an out of plane part. Then we Fourier transform the in-plane part in x and y. The result has a momentum dependent in-plane part, which is similar to the result for the bulk, except it does not contain terms in the z direction, and has a momentum limited to the plane. In addition, we have the out of plane real space terms which simply come from the above real space Hamiltonian, but are limited to the z direction. He (k⊥ ) = − 2λ⊥ τz (σx sin(kx a) + σy sin(ky a) + iλz T˜ τz σz ˜ x + [ε − 2t⊥ (cos(kx a) + cos(ky a))]τx − tz Rτ  + [γ0 − 2γ⊥ (cos(kx a) + cos(ky a)] − γz R˜  119  (7.99)  where He = ∑ ϕ † (k⊥ )He (k⊥ )ϕ(k⊥ )  (7.100)  k⊥  and R˜ and T˜ are matrices of dimension 4Lz given by R˜ i j = δi, j−1 + δi, j+1 and T˜i j = δi, j−1 − δi, j+1 , and Ψ j = ∑ Ψ(k⊥ ) eik⊥ ·R j  (7.101)  k⊥  (the dependence on z is implicit since it does not change in the transform) Now that we have the form of the electron Hamiltonian, we perform the mean field approximation as done in the previous section. We rewrite the interaction term, define the magnetization as a z dependent quantity, and neglect the second order term in the fluctuations, to find J ∑ ∑ SI · sI j I∈z  ∑ j  J M j · ∑ si + Jm j · ∑ SI − N⊥ J M j · m j i∈ j  (7.102)  I∈ j  where M j and m j are the magnetization of the impurities and electrons (respectively) for layer j (meaning z = j, in units of the lattice constant). The sum over j runs from 1 to Lz and the sum over i ∈ j runs over sites i in the x-y plane defined by  z = i. Since we have periodic boundary conditions only in the x and y directions, we Fourier transform only in those directions, and keep the dependence on z. So that the Hamiltonian, in its approximate mean-field form is H MF =  ∑ ϕ † (k⊥ ) k⊥  (7.103) He (k⊥ ) + J M ⊗  σ − µ ϕ(k⊥ ) + J ∑ m j · SI − N⊥ J ∑ M j · m j 2 j,I∈ j j  where Ψk⊥ is a vector of length 4Lz , and He (k⊥ ) is a matrix of dimension 4Lz . Here M is a vector of length Lz , and by M ⊗ σ we mean that M is expanded to a  diagonal matrix of dimension Lz and then we perform an external product with σ . We consider the case in which the magnetization is in the z direction, so M = M zˆ  120  and m = mˆz (M and m are both vectors of length Lz ). We have H MF  = =  ∑ ϕ † (k⊥ ) k⊥ HeMF  He (k⊥ ) + JM ⊗  σz − µ ϕ(k⊥ ) + J ∑ m j · SI − N⊥ J ∑ M j · m j 2 j j,I∈ j  + HIMF − N⊥ J ∑ M j · m j  (7.104)  j  We now diagonalize the full electron Hamiltonian He (k⊥ )+JM ⊗ σ2z − µ, and find a matrix U which diagonalizes it. A given vector ϕ of length 4Lz can be transformed into the basis of eigenvectors by ϕ˜ = U † ϕ which can be inverted to give ϕ = U ϕ˜  We want to calculate the electron magnetization of each layer, which is given by the sum over the expectation values of the spin operator mi = ∑ Ψi | k⊥  σz |Ψi 2  (7.105)  We define ϕi as a vector of length 4Lz which contains mostly zeros and four non-zero elements from Ψi in the corresponding location. For example, ϕ3 = (0, 0, Ψ3 , 0, . . . ) (the zeros are zero matrices of dimension 4) so that mi = ∑ ϕi | k⊥  σz σz |ϕi = ∑ ϕ˜ i |U † U|ϕ˜ i = 2 2 k⊥  4i  ∑  U†  j=4i−3  σz U 2  jj  f (E j ) (7.106)  The chemical potential can be determined as before, by solving n=  1 N  ∑ k⊥ ,α  f (Eα (k⊥ ))  (7.107)  for µ, with n = 2 for half filling and α runs from 1 to 4Lz . The results for the magnetization in this model are presented in Fig 7.4 and Fig 7.5. The bulk is ferromagnetically ordered up to Tcbulk face remains FM ordered up to  Tcsurf  73K, and the sur-  102K for the two surfaces (see Fig 7.4).  Therefore, the window in which the surface is ordered and the bulk is unordered (paramagnetic) is  29K.  We plot the magnetization in the z direction in Fig 7.5. Here the effect can be seen clearly. If we ramp up the temperature from T = 0 at first all spins, regardless 121  0 40 60 80 90 100 120  2.5 2.0  S (µ B )  1.5 1.0 0.5 0.0 ï0.5 ï1.0 2  4  6  8  10  z  12  14  16  18  20  Figure 7.4: Impurity spin expectation value as a function of the z coordinate, for different temperatures. The parameters are Lx = Ly = 40, Lz = 20, S = 5/2, J = 0.5eV, x = 0.05, and the rest of the parameters are chosen to fit the Bi2 Se3 dispersion close to the Γ point based on first principles calculations [34, 42, 43]: γ0 = 0.3391 eV, γ⊥ = 0.0506 eV, γz = 0.0717 eV, ε = 1.6912 eV,t⊥ = 0.3892 eV,tz = 0.2072 eV, λ⊥ = 0.2170 eV and λz = 0.1240 eV. of if they are in the bulk or surface, are fully polarized. As we increase the temperature, the magnetization of the bulk drops faster than the magnetization of the surface. Eventually we cross Tcbulk at which stage the magnetization of the bulk is zero, but the magnetization of the surface is finite. If we increase the temperature further, the magnetization of the surface drops, but the magnetization of the bulk remains at zero. Once we cross Tcsurf , the magnetization of both the bulk and the surface vanish – thermal fluctuations have broken the ordered phases. In samples with open surfaces we observed spatial fluctuations in the bulk magnetization as the temperature approached Tcbulk from above, e.g. the T = 80K curve in Fig 7.4. We originally attributed these to the large magnetic susceptibility of the bulk, which diverges at T → Tcbulk , and the proximity of the ordered surfaces, and  found Tcbulk from the bulk calculation with periodic boundary conditions, the results of which are plotted in Fig. 7.5. It was recently shown (after publication of our paper [99]) that these fluctuations are a spin density wave (SDW) phase that occurs due to a peak in the magnetic susceptibility [110], which could explain the  122  2.5  surf avg bulk  S (µ B )  2.0 1.5 1.0  0.0 0  T csurf  T cbulk  0.5  20  40  60  T  80  100  120  Figure 7.5: Impurity spin expectation value as a function of temperature for the discretized Bi2 Se3 lattice. We plot the order parameter on the surfaces (“surf”), the averaged bulk result for the 14 middle sites (“avg”), and the results of the separate bulk calculation (“bulk”). The parameters and curves are as in Fig. 7.4 source of the fluctuations that we observed. We plot the critical temperature for the surface and the bulk as a function of the chemical potential in Fig. 7.6. The surface critical temperature is maximal when the chemical potential intersects the surface Dirac point, which happens for µ  0.09eV, and falls on both sides, as expected. The result agrees well with  the predicted linear dependence of the surface critical temperature on µ predicted in Eq. (7.30). The bulk critical temperature is approximately constant within the bulk gap, as expected – since there are no bulk states within the gap, changing the chemical potential should not change the critical temperature.  7.8  Perovskite lattice - bulk  Due to the problems we encountered in the simulation of Bi2 Se3 (namely, the fluctuations in magnetization), we decided to try another TI model, the so called Perovskite lattice [111]. This lattice can be described as a simple cubic lattice with additional sites on the center of all edges. The nearest neighbor tight binding Hamil-  123  surf surfa bulk  100 95  Tc  90 85 80 75 70 ï0.10 ï0.05  0.00  0.05 µ 0.10  0.15  0.20  0.25  Figure 7.6: Critical temperature as a function of the chemical potential µ for the discretized Bi2 Se3 lattice. We plot the surface critical temperature (squares, blue), and the approximate surface temperature from Eq. (7.30) (solid line, black). In addition, we plot the bulk critical temperature as obtained from an average over the bulk sites (circles, red). The error bars were obtained by fitting the magnetization below the critical temperature to the known behaviour for a second order phase transition, M ∝ (Tc − T )1/2 . The parameters are as in Fig. 7.4, except for Lz = 10. The cutoff was chosen for best fitting by eye Λa/2π = 0.215. The vertical dashed lines mark the bulk gap. tonian is    cos(kx a) cos(ky a) cos(kz a)   cos(kx a)  0 0 0  H(k) = −2t  cos(k a) 0 0 0  y   cos(kz a) 0 0 0 0  124  (7.108)  Figure 7.7: A unit cell of the Perovskite lattice, which can be described as edge centered cubic. The four basis sites are labelled. To add spin-orbit interaction we must expand the Hamiltonian by adding a spin subspace, giving an 8 × 8 matrix    0    0  HSO (k) = 4λ  σx ⊗ 0   0  0  0 σy ⊗  0  0  0  0 σz ⊗  0  0  0  0  0  0  0  0  0 −i 0  0  0  0  0  0  −i 0 0  0  0  i  −i 0 0  0   0  0  sin(ky a) sin(kz a) + i  0  0  i  sin(kx a) sin(kz a) + 0  0   0    0  sin(kx a) sin(ky a)   0  0  (7.109)  The remainder of the process of determining the magnetization self consistently is the same, except that here we have a four point basis which is separated spatially, so we determine the magnetization separately for each of the basis points on each 125  site. For the electron magnetization, we have mi = ∑ Ψ j | k, j  σz σz 1 (U † U) j j f (E j ) ⊗ Oi |Ψ j = ∑ 2 N/4 2  (7.110)  k, j  where Oi is a 4 × 4 matrix that is zero everywhere except for Oii = 1. Here i  runs from 1 to 4 and corresponds to the basis site for which we are calculating the magnetization. For the impurity magnetization, we have Mi = −2SxBS (β Jmi S)  (7.111)  As before, these equations and the constraint on the chemical potential are used to successively calculate new values from the old ones, until reaching self consistency. In this case we see that as we ramp up T from zero, the system is initially in a FM state, characterized by the magnetization on all basis sites being equal. This persists up to a critical temperature T ∼ 5K where a phase transition occurs to an  AF state. This phase is characterized by a total magnetization of 0 with the sites 1 and 4 being polarized in the one direction and 2 and 3 in the opposite direction. This phase persists until Tcbulk  60K, after which the system is paramagnetic  (thermal fluctuations cause the breakdown of the ordered state).  7.9  Perovskite lattice - z dependent case  To study a system containing both edges and bulk, we now Fourier transform (FT) the Hamiltonian only in the x and y directions, keeping the z dependence, and using open boundary conditions. In the previous section we had a Hamiltonian matrix of dimension 8 = 4 × 2, due to the 4 point basis and the spin degree of freedom. Now we must include an additional subspace for the dependence on z,  so the Hamiltonian is of dimension 8Lz , where Lz is the number of sites in the z direction. The nearest neighbour hopping in the x and y directions remains as before, but with an additional identity matrix in the z subspace, and the hopping in the z direction is given by the hopping constant multiplied by a suitable matrix for the one dimensional hopping in the z axis. We have  126  2  1 2 3 4  S (µ B )  1 0 ï1 ï2 0  10  20  30  T (K)  40  50  2.5  AF FM  2 S (µ B )  60  1.5 1 0.5 0 0  10  20  30  T (K)  40  50  60  Figure 7.8: Bulk impurity magnetization on each basis site as a function of temperature (top) and FM and AF order parameters as a function of temperature (bottom) for S = 5/2, x = 0.05, L = 10,J = 0.25, λ = 0.2t     cos(kx a) cos(ky a) 0    cos(kx a) 0 0 0  − t R˜ ⊗ I2 H(k) = ILz ⊗ −2t  cos(k a) 0 0 0 y   0 0 0 0 0  (7.112)  where the matrix R˜ is a matrix of dimension 4Lz , acting in the z subspace and the 4-  127  point basis subspaces (with an identity in the spin subspace). Note that the nearest neighbour hopping in the z direction is only between the basis sites 1 and 4, so the matrix R˜ is all zero except for 1 for those hoppings. This can be written explicitly as  Ru + Rd Rd 0 0  Ru + Rd Rd 0  Ru   ˜ R= 0 Ru Ru + Rd . . .  ... ... ...  ... 0  0  0  Ru  0    0          0 Rd Ru + Rd  (7.113)  where Riuj = δi1 δ j4 , Ridj = δi4 δ j1 and 0 is understood to be a zero matrix of dimension 4. The spin orbit hopping is given by  HSO (k) = 2λ σx ⊗ T˜x sin(ky a) + σy ⊗ T˜y sin(kx a) sin(kz a)    0 0 0 0    0 0 i 0    + σz ⊗  0 −i 0 0 sin(kx a) sin(ky a)    0 0 0 0  (7.114)  Notice that the hopping between the 2 and 3 sites (last term) happens in the x −  y plane, so it is the same as in the bulk. The two other terms involve diagonal hopping: the hopping along x or y gives the corresponding sin function and the hopping in the z direction gives the Tx and Ty matrices. The matrix T˜x codes the hopping for sites 2 and 4 (site 2 is on the x axis)   Tu + Td    −Tu  ˜ Tx =   0   ... 0  −Td  Tu + Td  0  0  0  −Td  0  0  −Tu  Tu + Td  ...  ...  ...  ...  0  0  −Tu       0    −Td  Tu + Td  (7.115)  where Tui j = δi2 δ j4 , Tdi j = δi4 δ j2 and 0 is understood to be a zero matrix of dimen128  sion 4. Similarly, the matrix T˜y codes the hopping for sites 3 and 4 (site 3 is on the y axis)   Tu + Td    −Tu  T˜y =   0   ... 0  −Td  Tu + Td  0  0  0  −Td  0  0  −Tu  Tu + Td  ...  ...  ...  ...  0  0  −Tu       0    −Td  Tu + Td  (7.116)  where Tui j = δi3 δ j4 , Tdi j = δi4 δ j3 and 0 is understood to be a zero matrix of dimension 4. Note that in this case the magnetization as a function of z is not expected to have a mirror symmetry around the middle of the system. The reason is that the 1 site and 4 site are not equivalent, and such a symmetry would map the one on to the other. The results for magnetization in this model are presented in Fig 7.9 and Fig 7.10. In this case we observe that for a large range of temperatures, the magnetization on basis sites 1 and 4 is similar and opposite in sign from the magnetization on basis sites 2 and 3 (which are equal due to a rotational symmetry around the z axis). This motivates the definition of a FM order parameter SiFM = (Si,1 + Si,2 + Si,3 + Si,4 )/4 and an AF order parameter SiAF = (Si,1 − Si,2 − Si,3 + Si,4 )/4 (where Si, is the ex-  pectation value of the impurity spin on basis site of lattice site i). We see that the system is ferromagnetically ordered up to T The bulk remains AF ordered until dered until  Tcsurf  Tcbulk  4.5K, where it becomes AF ordered. 60.5K, and the surface remains AF or-  72K and 78K for the two surfaces. Therefore, the maximum  window in which the surface is ordered and the bulk is unordered (paramagnetic) is  17.5K. Note that the difference in surfaces results from the fact that the unit  cell is not symmetric under reflection along z - the “top” of the system ends with basis site 4 and the bottom ends with basis site 1 (see Fig. 7.7). In addition, we see that the surfaces undergo an additional partial phase change signified by discrete jumps of the order parameter seen at T  4.5, 5.5, 13.5, 35.5K, which we claim as  additional evidence that the bulk and surface differ.  129  2.5  surf1 surf2 z bulk  SF M (µ B )  2.0  T cbulk  1.5 1.0  T csurf  0.5 0.0 0  20  40  T  60  80  Figure 7.9: Ferromagnetic order parameter as a function of temperature for the Perovskite lattice, with parameters Lx = Ly = 20, Lz = 10, S = 5/2, J = 0.25 eV, t = 1 eV, λ = 0.2 eV, x = 0.05. We plot the order parameter on each of the surfaces (“surf1” and “surf2”), the result for the site at z = Lz /2 (“z”), and the results of the separate bulk calculation (“bulk”).  7.10  Conclusions  We have demonstrated that magnetically doped topological insulators can have a sizeable window where the bulk is paramagnetic and the surface is magnetically ordered. Our conclusions are based on general arguments that involve only universal properties of the topologically protected surface states and on numerical calculations performed on simple lattice models of 3D topological insulators. Physically, these results are in accord with the intuition that the metallic state on the surface should be more susceptible to magnetic ordering than the insulating bulk, although this expectation is only partially borne out in systems with strong spin-orbit coupling. The results reported in this study rely on two key approximations: (i) the mean-field decoupling of the exchange interaction between electrons and magnetic dopants indicated in Eq. (7.38), and (ii) the ‘virtual crystal’ approximation which replaces the localized magnetic moments of dopant atoms by their average over all lattice sites. We have attempted to bypass the latter approximation by solving 130  SAF (µ B )  2  T cbulk  1  surf1 surf2 z avg bulk  0 ï1  T csurf  ï2 0  20  40  T  60  80  Figure 7.10: Antiferromagnetic order parameter as a function of temperature for the Perovskite lattice, with parameters as in (a). Here we also plot the average over the four central sites (“avg”). the problem in real space without averaging over sites. This procedure gives reasonable results for Tcbulk compared to those obtained in the virtual crystal (thereby verifying that this approximation is reasonable) but is numerically more costly. The system sizes that we could simulate did not allow us to unambiguously determine the surface critical temperature. Nevertheless based on these results we feel that, in a large crystal, (ii) will not lead to a significant error in the determination of the critical temperatures. Since fluctuations around the mean field result are typically stronger in 2D than in 3D, they will likely reduce the size of the temperature window between Tcbulk and Tcsurf . On the other hand, our mean field calculation did not include electronelectron interactions, which would tend to stabilize the long range magnetic order and thereby strengthen the ordered phases. It has been shown recently that electronelectron interactions can lead to spontaneous breaking of TRI on the surface of a TI [112], even in the absence of magnetic dopants. Thus the interactions might in fact strengthen the effect found in our study, although a more detailed investigation would be needed in order to obtain a quantitative result. Overall, in our opinion it is very likely that the combined effect of the magnetic dopants and the electron-  131  electron interactions can account for the experimentally observed surface excitation gap without bulk magnetic order in Mn and Fe doped Bi2 Se3 . We note that the critical temperature for the bulk magnetic ordering of the magnetically doped topological insulator Bi2−x Mnx Te3 was recently measured [107] to be 9 − 12K for x = 0.04 − 0.09. This critical temperature is smaller than our results  and those of Ref. [89], which estimates Tcbulk  70K for Cr-doped Bi2 Se3 using  first-principles numerical calculations. We chose our coupling constant J based on the first-principles results, so we expect that for a reduced coupling constant the temperature window would shift to lower temperatures. As noted above, our arguments are universal, hence the exact details of the material and coupling might alter the result quantitatively but should not change it qualitatively. Taking the above finding for Bi2−x Mnx Te3 as a guideline, one may surmise that Tcbulk for Fe and Mn doped Bi2 Se3 lies in a similar range of temperatures. It is then entirely possible that the ARPES experiments [50, 51], performed at ∼ 20K, detected surface  magnetic ordering without bulk magnetism as advocated in this paper. Careful surface-sensitive magnetic measurements, as proposed recently [113], might be able to probe this intriguing phenomenon directly. The possible presence of a Kondo effect is an interesting question. Many years ago [114] it was shown that in a Dirac metal the Kondo effect takes place at T = 0K only if the Kondo coupling is above a critical value, in contrast with the normal metal case, where the transition happens at zero coupling. Several much more recent theoretical studies considered an Anderson impurity on the surface of a topological insulator [115, 116], and a recent experimental study on Bi2 Se3 nanoribbons found a Kondo temperature of ∼ 30K [117]. These results  indicate that the Kondo temperature might be on the order of the surface critical temperature for Bi2 Se3 , but further study would be necessary to see whether it is indeed so, and whether this would affect our results qualitatively.  132  Chapter 8  Conclusions 8.1  Summary  In this work we study exotic effects that are predicted to occur in topological insulators (the Witten Effect and the Wormhole Effect), establish a new phase of matter - the novel Topological Anderson Insulator, and offer a possible explanation for recent experiments on magnetically doped topological insulators. We summarize the key results from each chapter below. In Chapter 4 we observe that the long predicted Witten Effect in particle physics should occur in a topological insulators as well - a unit magnetic monopole in a medium with θ = 0 is predicted to bind a fractional electric charge −eθ /2π. Using  a simple physical model for a topological insulator we demonstrate the existence of a fractional charge bound to a monopole by an explicit numerical calculation. We also propose a scheme for generating an artificial magnetic monopole in a topological insulator film, that may be used to facilitate the first experimental test of Wittens prediction. In Chapter 5, we demonstrate that a Dirac string with a half flux quantum flux inserted into a STI carries protected gapless one-dimensional fermionic modes. These modes are spin-filtered and represent a distinct bulk manifestation of the topologically non-trivial insulator. We establish this wormhole effect by both general qualitative considerations and by numerical calculations within a minimal lattice model. We also discuss the possibility of experimental observation of a closely 133  related effect in artificially engineered nanostructures. In Chapter 6, we show that remarkably, under certain conditions disorder is responsible for the existence of a STI referred to as a strong topological ’Anderson’ insulator. We study the phase diagram using the Born approximation. In Chapter 7 we lay out the hypothesis that a temperature window exists in which the surface of a magnetically doped TI is magnetically ordered but the bulk is not. We present a simple and intuitive argument why this is so, and back it up by a mean-field calculation for two simple tight binding TI models: a cubic-lattice regularized Bi2 Se3 and a model on the Perovskite lattice. Our results show that indeed a sizeable regime such as described above could exist in real TIs, and this indicates a possible physical explanation for the results seen in experiments.  8.2  Future work  The field of TIs is a young field, and as such there are many basic questions which have yet to be fully addressed. The effect of interactions and the discovery of additional topological phases are some of the challenges that await. In Chapter 4 we proposed an experiment that could possibly allow for the first observation of the Witten effect. To the best of our knowledge, this experiment has yet to be done. The challenges are creating an artificial magnetic monopole and measuring the fractional charge. In Chapter 5 we suggest an experiment to observe the wormhole effect, by drilling sub-micron holes in a STI and sweeping the magnetic field, resulting in a periodic variation of the conductance along the hole(s). The main experimental difficulty appears to be the fabrication of sub-micron holes and measuring the conductance or the charge passing along the holes. In Chapter 6 we discuss the possibility of constructing a topological insulator from a trivial insulator (with a small gap) by adding disorder. There are some good candidates, such as Sb2 Se3 , bulk HgTe and the half Heusler compounds, but their low energy physics is not well described by our four-band model, so further theoretical work is necessary to determine whether disorder could drive the transition into the topological phase. Chapter 7 we used the virtual crystal approximation, a mean field approxima-  134  tion and neglected interactions to derive a result which we think is to be trusted qualitatively. It would be useful to derive a more accurate result, which could then, perhaps, be trusted quantitatively, and tested in an experiment.  135  Bibliography [1] Eduardo Fradkin, Elbio Dagotto, and Daniel Boyanovsky. Physical realization of the parity anomaly in condensed matter physics. Phys. Rev. Lett., 57:2967–2970, Dec 1986. → pages iii, 22 [2] Liang Fu and C. L. Kane. Topological insulators with inversion symmetry. Phys. Rev. B, 76:045302, Jul 2007. → pages iii, 8, 19, 21, 26, 41, 59, 60 [3] C. W. Groth, M. Wimmer, A. R. Akhmerov, J. Tworzydło, and C. W. J. Beenakker. Theory of the topological Anderson insulator. Phys. Rev. Lett., 103:196805, Nov 2009. → pages iv, 12, 68, 69, 74 [4] S. Doniach and E.H. Sondheimer. Green’s functions for solid state physicists. Imperial College Press, 1998. → pages iv, 172 [5] M. Z. Hasan and C. L. Kane. Colloquium: topological insulators. Rev. Mod. Phys., 82:3045–3067, Nov 2010. → pages 1, 3, 4 [6] J. E. Moore. The birth of topological insulators. Nature, 464:194, 2010. → pages [7] Xiao-Liang Qi and Shou-Cheng Zhang. Topological insulators and superconductors. Rev. Mod. Phys., 83:1057–1110, Oct 2011. → pages [8] M. Zahid Hasan and Joel E. Moore. Three-dimensional topological insulators. Annual Review of Condensed Matter Physics, 2(1):55–78, 2011. → pages [9] Xiao-Liang Qi and Shou-Cheng Zhang. The quantum spin Hall effect and topological insulators. Physics Today, 63(1):33–38, 2010. → pages 1, 5, 6 [10] R.E. Prange and S.M. Girvin. The quantum Hall effect:. Springer-Verlag, 1987. → pages 2 136  [11] K. v. Klitzing, G. Dorda, and M. Pepper. New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance. Phys. Rev. Lett., 45:494–497, Aug 1980. → pages 2 [12] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs. Quantized Hall conductance in a two-dimensional periodic potential. Phys. Rev. Lett., 49:405–408, Aug 1982. → pages 2 [13] Qian Niu, D. J. Thouless, and Yong-Shi Wu. Quantized Hall conductance as a topological invariant. Phys. Rev. B, 31:3372–3377, Mar 1985. → pages 2 [14] Klaus von Klitzing. Developments in the quantum Hall effect. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 363(1834):2203–2219, 2005. → pages 2 [15] Shuichi Murakami, Naoto Nagaosa, and Shou-Cheng Zhang. Spin-Hall insulator. Phys. Rev. Lett., 93:156804, Oct 2004. → pages 3 [16] F. D. M. Haldane. Model for a quantum Hall effect without landau levels: condensed-matter realization of the ”parity anomaly”. Phys. Rev. Lett., 61:2015–2018, Oct 1988. → pages 3 [17] C. L. Kane and E. J. Mele. Quantum spin Hall effect in graphene. Phys. Rev. Lett., 95:226801, Nov 2005. → pages 3, 18 [18] C. L. Kane and E. J. Mele. Z2 topological order and the quantum spin Hall effect. Phys. Rev. Lett., 95:146802, Sep 2005. → pages [19] B. Andrei Bernevig and Shou-Cheng Zhang. Quantum spin Hall effect. Phys. Rev. Lett., 96:106802, Mar 2006. → pages 3 [20] H.A. Kramers. Thorie gnrale de la rotation paramagntique dans les cristaux. Proceedings Koninklijke Akademie van Wetenschappen, 33:959–972, 1930. → pages 4 [21] Daniel Huertas-Hernando, F. Guinea, and Arne Brataas. Spin-orbit coupling in curved graphene, fullerenes, nanotubes, and nanotube caps. Phys. Rev. B, 74:155426, Oct 2006. → pages 7 [22] Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Kleinman, and A. H. MacDonald. Intrinsic and Rashba spin-orbit interactions in graphene sheets. Phys. Rev. B, 74:165310, Oct 2006. → pages 137  [23] Yugui Yao, Fei Ye, Xiao-Liang Qi, Shou-Cheng Zhang, and Zhong Fang. Spin-orbit gap of graphene: First-principles calculations. Phys. Rev. B, 75:041401, Jan 2007. → pages [24] J. C. Boettger and S. B. Trickey. First-principles calculation of the spin-orbit splitting in graphene. Phys. Rev. B, 75:121402, Mar 2007. → pages [25] M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl, and J. Fabian. Band-structure topologies of graphene: spin-orbit coupling effects from first principles. Phys. Rev. B, 80:235431, Dec 2009. → pages 7 [26] B. Andrei Bernevig, Taylor L. Hughes, and Shou-Cheng Zhang. Quantum spin Hall effect and topological phase transition in HgTe quantum wells. Science, 314(5806):1757–1761, 2006. → pages 7 [27] Markus Konig, Steffen Wiedmann, Christoph Brune, Andreas Roth, Hartmut Buhmann, Laurens W. Molenkamp, Xiao-Liang Qi, and Shou-Cheng Zhang. Quantum spin Hall insulator state in HgTe quantum wells. Science, 318(5851):766–770, 2007. → pages 7, 12, 68 [28] Markus K¨onig, Hartmut Buhmann, Laurens W. Molenkamp, Taylor Hughes, Chao-Xing Liu, Xiao-Liang Qi, and Shou-Cheng Zhang. The quantum spin Hall effect: theory and experiment. Journal of the Physical Society of Japan, 77(3):031007, 2008. → pages 7 [29] Xiao-Liang Qi, Taylor L. Hughes, and Shou-Cheng Zhang. Topological field theory of time-reversal invariant insulators. Phys. Rev. B, 78:195424, Nov 2008. → pages 7, 8, 40, 41, 44, 52, 58, 60, 63, 87, 93 [30] Andrew M. Essin, Joel E. Moore, and David Vanderbilt. Magnetoelectric polarizability and axion electrodynamics in crystalline insulators. Phys. Rev. Lett., 102:146805, Apr 2009. → pages 7, 8, 10, 37, 40, 41, 52, 58, 60, 87 [31] M. M. Vazifeh and M. Franz. Quantization and 2π periodicity of the axion action in topological insulators. Phys. Rev. B, 82:233103, Dec 2010. → pages 7, 42 [32] E. Witten. Dyons of charge eθ /2π. Phys. Rev. B, 86(34):283 – 287, 1979. → pages 8, 40, 41, 57, 64, 87 [33] G. Rosenberg and M. Franz. Witten effect in a crystalline topological insulator. Phys. Rev. B, 82:035105, Jul 2010. → pages 8, 63, 64 138  [34] Haijun Zhang, Chao-Xing Liu, Xiao-Liang Qi, Xi Dai, Zhong Fang, and Shou-Cheng Zhang. Topological insulators in Bi2 Se3 , Bi2 Te3 and Sb2 Te3 with a single Dirac cone on the surface. Nature Physics, 5(6):438–442, 06 2009. → pages 9, 11, 41, 59, 96, 99, 122 [35] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan. Discovery (theoretical prediction and experimental observation) of a large-gap topological-insulator class with spin-polarized single-Dirac-cone on the surface. Nature Physics, 5:398, 2009. → pages 9, 10, 41, 58, 59 [36] D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier, J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan. A tunable topological insulator in the spin helical Dirac transport regime. Nature, 460(7259):1101–1105, 08 2009. → pages [37] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen. Experimental realization of a three-dimensional topological insulator, Bi2 Te3 . Science, 325(5937):178–181, 2009. → pages 9, 41, 58, 59 [38] T. Hanaguri, K. Igarashi, M. Kawamura, H. Takagi, and T. Sasagawa. Momentum-resolved landau-level spectroscopy of Dirac surface state in Bi2 Se3 . Phys. Rev. B, 82:081305, Aug 2010. → pages 9 [39] Pedram Roushan, Jungpil Seo, Colin V. Parker, Y. S. Hor, D. Hsieh, Dong Qian, Anthony Richardella, M. Z. Hasan, R. J. Cava, and Ali Yazdani. Topological surface states protected from backscattering by chiral spin texture. Nature, 460(7259):1106–1109, 08 2009. → pages [40] Tong Zhang, Peng Cheng, Xi Chen, Jin-Feng Jia, Xucun Ma, Ke He, Lili Wang, Haijun Zhang, Xi Dai, Zhong Fang, Xincheng Xie, and Qi-Kun Xue. Experimental demonstration of topological surface states protected by time-reversal symmetry. Phys. Rev. Lett., 103:266803, Dec 2009. → pages 9 [41] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim. The electronic properties of graphene. Rev. Mod. Phys., 81:109–162, Jan 2009. → pages 9  139  [42] Liang Fu and Erez Berg. Odd-parity topological superconductors: theory and application to Cux Bi2 Se3 . Phys. Rev. Lett., 105:097001, Aug 2010. → pages 9, 96, 122 [43] Rundong Li, Jing Wang, Xiao-Liang Qi, and Shou-Cheng Zhang. Dynamical axion field in topological magnetic insulators. Nature Physics, 6(4):284–288, 04 2010. → pages 9, 40, 96, 122 [44] Xiao-Liang Qi, Rundong Li, Jiadong Zang, and Shou-Cheng Zhang. Inducing a magnetic monopole with topological surface states. Science, 323(5918):1184–1187, 2009. → pages 10, 60 [45] Jiadong Zang and Naoto Nagaosa. Monopole current and unconventional Hall response on a topological insulator. Phys. Rev. B, 81:245125, Jun 2010. → pages 10 [46] Wang-Kong Tse and A. H. MacDonald. Giant magneto-optical Kerr effect and universal Faraday effect in thin-film topological insulators. Phys. Rev. Lett., 105:057401, Jul 2010. → pages 10 [47] Joseph Maciejko, Xiao-Liang Qi, H. Dennis Drew, and Shou-Cheng Zhang. Topological quantization in units of the fine structure constant. Phys. Rev. Lett., 105:166803, Oct 2010. → pages 10 [48] Takehito Yokoyama, Yukio Tanaka, and Naoto Nagaosa. Anomalous magnetoresistance of a two-dimensional ferromagnet/ferromagnet junction on the surface of a topological insulator. Phys. Rev. B, 81:121401, Mar 2010. → pages 10 [49] Ion Garate and M. Franz. Inverse spin-galvanic effect in the interface between a topological insulator and a ferromagnet. Phys. Rev. Lett., 104:146802, Apr 2010. → pages 10 [50] Y. L. Chen, J.-H. Chu, J. G. Analytis, Z. K. Liu, K. Igarashi, H.-H. Kuo, X. L. Qi, S. K. Mo, R. G. Moore, D. H. Lu, M. Hashimoto, T. Sasagawa, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z. X. Shen. Massive Dirac Fermion on the surface of a magnetically doped topological insulator. Science, 329(5992):659–662, 2010. → pages 10, 92, 94, 132 [51] L. Andrew Wray, Su-Yang Xu, Yuqi Xia, David Hsieh, Alexei V. Fedorov, Yew San Hor, Robert J. Cava, Arun Bansil, Hsin Lin, and M. Zahid Hasan. A topological insulator surface under strong coulomb, magnetic and disorder perturbations. Nature Physics, 7(1):32–37, 01 2011. → pages 10, 92, 94, 132 140  [52] Jian Li, Rui-Lin Chu, J. K. Jain, and Shun-Qing Shen. Topological Anderson insulator. Phys. Rev. Lett., 102:136806, Apr 2009. → pages 12, 68 [53] B. Andrei Bernevig, Taylor L. Hughes, and Shou-Cheng Zhang. Quantum spin Hall effect and topological phase transition in HgTe quantum wells. Science, 314(5806):1757–1761, 2006. → pages 12, 68 [54] R. Peierls. Zur theorie des diamagnetismus von leitungselektronen. Zeitschrift fur Physik, 80:763–791, November 1933. → pages 18 [55] Hongki Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, Leonard Kleinman, and A. H. MacDonald. Intrinsic and Rashba spin-orbit interactions in graphene sheets. Phys. Rev. B, 74:165310, Oct 2006. → pages 18 [56] G. W. Semenoff, V. Semenoff, and Fei Zhou. Domain walls in gapped graphene. Phys. Rev. Lett., 101:087204, Aug 2008. → pages 36 [57] R. D. Peccei and Helen R. Quinn. CP conservation in the presence of pseudoparticles. Phys. Rev. Lett., 38:1440–1443, Jun 1977. → pages 40 [58] R. D. Peccei and Helen R. Quinn. Constraints imposed by CP conservation in the presence of pseudoparticles. Phys. Rev. D, 16:1791–1797, Sep 1977. → pages 40 [59] Steven Weinberg. A new light boson? Phys. Rev. Lett., 40:223–226, Jan 1978. → pages 40 [60] F. Wilczek. Problem of strong P and T invariance in the presence of instantons. Phys. Rev. Lett., 40:279–282, Jan 1978. → pages 40 [61] Leanne D Duffy and Karl van Bibber. Axions as dark matter particles. New Journal of Physics, 11(10):105008, 2009. → pages 40 [62] L. D. Duffy, P. Sikivie, D. B. Tanner, S. J. Asztalos, C. Hagmann, D. Kinion, L. J Rosenberg, K. van Bibber, D. B. Yu, and R. F. Bradley. High resolution search for dark-matter axions. Phys. Rev. D, 74:012006, Jul 2006. → pages 40 [63] M. Kuster, B. Beltr´an, G.G. Raffelt, and European Organization for Nuclear Research. Axions: theory, cosmology, and experimental searches. Lecture Notes in Physics. Springer, 2007. → pages 40 [64] Liang Fu, C. L. Kane, and E. J. Mele. Topological insulators in three dimensions. Phys. Rev. Lett., 98:106803, Mar 2007. → pages 41, 58, 59, 69 141  [65] J. E. Moore and L. Balents. Topological invariants of time-reversal-invariant band structures. Phys. Rev. B, 75:121306, Mar 2007. → pages 58, 69 [66] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan. A topological Dirac insulator in a quantum spin Hall phase. Nature, 452(7190):970–974, 04 2008. → pages 41, 58, 59 [67] Ian Affleck. Dyon analogs in antiferromagnetic chains. Phys. Rev. Lett., 57:1048–1051, Aug 1986. → pages 41 [68] B. Seradjeh, J. E. Moore, and M. Franz. Exciton condensation and charge fractionalization in a topological insulator film. Phys. Rev. Lett., 103:066402, Aug 2009. → pages 41, 53, 54, 57, 58 [69] J.D. Jackson. Classical electrodynamics. Wiley, 1999. → pages 41 [70] Frank Wilczek. Two applications of axion electrodynamics. Phys. Rev. Lett., 58:1799–1802, May 1987. → pages 42 [71] N. ANDERSON and A. M. ARTHURS. A variational principle for maxwell’s equations. International Journal of Electronics, 45(3):333–334, 1978. → pages 42 [72] Joe Rosen. Redundancy and superfluity for electromagnetic fields and potentials. American Journal of Physics, 48(12):1071–1073, 1980. → pages [73] A Sudbery. A vector lagrangian for the electromagnetic field. Journal of Physics A: Mathematical and General, 19(2):L33, 1986. → pages [74] S. C. Tiwari. On local duality invariance in electromagnetism. ArXiv e-prints, October 2011. → pages [75] L. Visinelli. Dual Axion Electrodynamics. ArXiv e-prints, November 2011. → pages 42 [76] Y.M. Shnir. Magnetic monopoles. Texts and Monographs in Physics. Springer, 2005. → pages 49, 52 [77] R. Jackiw and C. Rebbi. Solitons with Fermion number 21 . Phys. Rev. D, 13:3398–3409, Jun 1976. → pages 50, 58 [78] Jeffrey Goldstone and Frank Wilczek. Fractional quantum numbers on solitons. Phys. Rev. Lett., 47:986–989, Oct 1981. → pages 50, 58 142  [79] W. P. Su, J. R. Schrieffer, and A. J. Heeger. Solitons in polyacetylene. Phys. Rev. Lett., 42:1698–1701, Jun 1979. → pages 50, 58 [80] Kimball A Milton. Theoretical and experimental status of magnetic monopoles. Reports on Progress in Physics, 69(6):1637, 2006. → pages 52 [81] C. Castelnovo, R. Moessner, and S. L. Sondhi. Magnetic monopoles in spin ice. Nature, 451(7174):42–45, 01 2008. → pages 53, 58 [82] T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell. Magnetic coulomb phase in the spin ice Ho2 Ti2 O7 . Science, 326(5951):415–417, 2009. → pages 53 [83] D. J. P. Morris, D. A. Tennant, S. A. Grigera, B. Klemke, C. Castelnovo, R. Moessner, C. Czternasty, M. Meissner, K. C. Rule, J.-U. Hoffmann, K. Kiefer, S. Gerischer, D. Slobinsky, and R. S. Perry. Dirac strings and magnetic monopoles in the spin ice Dy2 Ti2 O7 . Science, 326(5951):411–414, 2009. → pages [84] S. T. Bramwell, S. R. Giblin, S. Calder, R. Aldus, D. Prabhakaran, and T. Fennell. Measurement of the charge and current of magnetic monopoles in spin ice. Nature, 461(7266):956–959, 10 2009. → pages 53, 58 [85] H.-M. Guo and M. Franz. Three-dimensional topological insulators on the pyrochlore lattice. Phys. Rev. Lett., 103:206805, Nov 2009. → pages 53 [86] Dmytro Pesin and Leon Balents. Mott physics and band topology in materials with strong spin-orbit interaction. Nature Physics, 6(5):376–381, 05 2010. → pages 53 [87] V. J. Goldman and B. Su. Resonant tunneling in the quantum Hall regime: measurement of fractional charge. Science, 267(5200):1010–1012, 1995. → pages 57, 58 [88] R. de Picciotto, M. Reznikov, M. Heiblum, V. Umansky, G. Bunin, and D. Mahalu. Direct observation of a fractional charge. Nature, 389(6647):162–164, 09 1997. → pages 57, 58 [89] Rui Yu, Wei Zhang, Hai-Jun Zhang, Shou-Cheng Zhang, Xi Dai, and Zhong Fang. Quantized anomalous Hall effect in magnetic topological insulators. Science, 329(5987):61–64, 2010. → pages 58, 95, 97, 99, 103, 132 143  [90] Rahul Roy. Topological phases and the quantum spin Hall effect in three dimensions. Phys. Rev. B, 79:195322, May 2009. → pages 59, 69 [91] Dung-Hai Lee. Surface states of topological insulators: the Dirac Fermion in curved two-dimensional spaces. Phys. Rev. Lett., 103:196804, Nov 2009. → pages 60 [92] R. B. Laughlin. Quantized Hall conductivity in two dimensions. Phys. Rev. B, 23:5632–5633, May 1981. → pages 60 [93] R. B. Laughlin. Anomalous quantum Hall effect: an incompressible quantum fluid with fractionally charged excitations. Phys. Rev. Lett., 50:1395–1398, May 1983. → pages 60 [94] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin. Interaction-induced criticality in Z2 topological insulators. Phys. Rev. Lett., 105:036803, Jul 2010. → pages 62 [95] James G. Analytis, Jiun-Haw Chu, Yulin Chen, Felipe Corredor, Ross D. McDonald, Z. X. Shen, and Ian R. Fisher. Bulk Fermi surface coexistence with Dirac surface state in Bi2 Se3 : A comparison of photoemission and Shubnikov de Haas measurements. Phys. Rev. B, 81:205407, May 2010. → pages 66 [96] Patrick A. Lee and T. V. Ramakrishnan. Disordered electronic systems. Rev. Mod. Phys., 57:287–337, Apr 1985. → pages 68 [97] Ferdinand Evers and Alexander D. Mirlin. Anderson transitions. Rev. Mod. Phys., 80:1355–1417, Oct 2008. → pages 68 [98] H.-M. Guo, G. Rosenberg, G. Refael, and M. Franz. Topological Anderson insulator in three dimensions. Phys. Rev. Lett., 105:216601, Nov 2010. → pages 89 [99] G. Rosenberg and M. Franz. Surface magnetic ordering in topological insulators with bulk magnetic dopants. Phys. Rev. B, 85:195119, May 2012. → pages 92, 122 [100] Liang Fu and C. L. Kane. Superconducting proximity effect and Majorana Fermions at the surface of a topological insulator. Phys. Rev. Lett., 100:096407, Mar 2008. → pages 93 [101] Qin Liu, Chao-Xing Liu, Cenke Xu, Xiao-Liang Qi, and Shou-Cheng Zhang. Magnetic impurities on the surface of a topological insulator. Phys. Rev. Lett., 102:156603, Apr 2009. → pages 94, 103 144  [102] T. C. Lubensky and Morton H. Rubin. Critical phenomena in semi-infinite systems. ii. mean-field theory. Phys. Rev. B, 12:3885–3901, Nov 1975. → pages 95 [103] Kaoru Ohno and Yutaka Okabe. Critical behavior of surface-layer magnetization at bulk Tc : Extraordinary transition. Phys. Rev. B, 39:9764–9767, May 1989. → pages 95 [104] P. Janicek, C. Drasar, P. Lostak, J. Vejpravova, and V. Sechovsky. Transport, magnetic, optical and thermodynamic properties of Bi2−x Mnx Se3 single crystals. Physica B: Condensed Matter, 403(1920):3553 – 3558, 2008. → pages 100 [105] Jeongyong Choi, Hee-Woong Lee, Bong-Seo Kim, Hyoyeol Park, Sungyoul Choi, S.C. Hong, and Sunglae Cho. Magnetic and transport properties of mn-doped Bi2 Se3 and Sb2 Se3 . Journal of Magnetism and Magnetic Materials, 304(1):e164 – e166, 2006. → pages 100 [106] T. Jungwirth, Jairo Sinova, J. Maˇsek, J. Kuˇcera, and A. H. MacDonald. Theory of ferromagnetic (III,Mn)V semiconductors. Rev. Mod. Phys., 78:809–864, Aug 2006. → pages 103 [107] Y. S. Hor, P. Roushan, H. Beidenkopf, J. Seo, D. Qu, J. G. Checkelsky, L. A. Wray, D. Hsieh, Y. Xia, S.-Y. Xu, D. Qian, M. Z. Hasan, N. P. Ong, A. Yazdani, and R. J. Cava. Development of ferromagnetism in the doped topological insulator Bi2−x Mnx Te3 . Phys. Rev. B, 81:195203, May 2010. → pages 103, 132 [108] L. Nordheim. Ann. Phys. (Leipzig), 9:607, 1931. → pages 104 [109] Neil W. Ashcroft and David N. Mermin. Solid state physics. Thomson Learning, Toronto, 1st edition, January 1976. → pages 109 [110] Martha Lasia and Luis Brey. Temperature-induced spin density wave in a magnetically doped topological insulator Bi2 Se3 . Phys. Rev. B, 86:045317, Jul 2012. → pages 122 [111] C. Weeks and M. Franz. Topological insulators on the lieb and perovskite lattices. Phys. Rev. B, 82:085310, Aug 2010. → pages 123 [112] Yuval Baum and Ady Stern. Magnetic instability on the surface of topological insulators. Phys. Rev. B, 85:121105, Mar 2012. → pages 131  145  [113] M. M. Vazifeh and M. Franz. Spin response of electrons on the surface of a topological insulator. Phys. Rev. B, 86:045451, Jul 2012. → pages 132 [114] David Withoff and Eduardo Fradkin. Phase transitions in gapless Fermi systems with magnetic impurities. Phys. Rev. Lett., 64:1835–1838, Apr 1990. → pages 132 ˇ [115] Rok Zitko. Quantum impurity on the surface of a topological insulator. Phys. Rev. B, 81:241414, Jun 2010. → pages 132 [116] Xiao-Yong Feng, Wei-Qiang Chen, Jin-Hua Gao, Qiang-Hua Wang, and Fu-Chun Zhang. Anderson impurity in a helical metal. Phys. Rev. B, 81:235411, Jun 2010. → pages 132 [117] Judy J. Cha, James R. Williams, Desheng Kong, Stefan Meister, Hailin Peng, Andrew J. Bestwick, Patrick Gallagher, David Goldhaber-Gordon, and Yi Cui. Magnetic doping and Kondo effect in Bi2 Se3 nanoribbons. Nano Letters, 10(3):1076–1081, 2010. PMID: 20131918. → pages 132  146  Appendix A  Appendix A.1  Derivatives in different sets of coordinates  A.1.1  Complex plane  We want to transform from (x, y) to (z, z∗ ) where z = x + iy, so that z∗ = x − iy, and z + z∗ 2 z − z∗ 2i  x = y =  (A.1)  Therefore ∂x ∂ z∗ ∂y ∂ z∗  =  1 2  = −  147  (A.2) 1 2i  If we have a function f (x, y), we can use the chain rule to write df dz  ∂ f ∂x ∂ f ∂y + ∂x ∂z ∂y ∂z 1 ∂f 1∂f + 2 ∂x i ∂y ∂f 1 ∂f −i 2 ∂x ∂y  = = =  (A.3)  so 1 d = dz 2  ∂ ∂ −i ∂x ∂y  (A.4)  d 1 = ∗ dz 2  ∂ ∂ +i ∂x ∂y  (A.5)  Similarly,  This can also be seen by taking the complex conjugate of the previous equation. Using these two derivatives, we write d2 dzdz∗  = =  1 4 1 4  ∂ ∂ −i ∂x ∂y 2 ∂2 ∂ + ∂ x 2 ∂ y2  ∂ ∂ +i ∂x ∂y 1 = ∇2 4  (A.6)  So that the Laplacian ∇2 = 4  A.1.2  d2 dzdz∗  (A.7)  Polar coordinates  We want to transform from (z, z∗ ) to (r, θ ), where z = r eiθ z∗ = r e−iθ  148  (A.8)  and √ zz∗ 1 z = ln ∗ 2i z  r = θ  (A.9)  Using the chain rule d dz  ∂r ∂ ∂θ ∂ + ∂z ∂r ∂z ∂θ 1 −iθ ∂ i ∂ e − 2 ∂r r ∂θ  = =  (A.10)  and similarly, d dz∗  =  1 iθ e 2  i ∂ ∂ + ∂r r ∂θ  (A.11)  Using (A.4) and (A.5) we find the useful result ∂ ∂ −i ∂x ∂y ∂ ∂ +i ∂x ∂y  =  e−iθ  =  eiθ  i ∂ ∂ − ∂r r ∂θ i ∂ ∂ + ∂r r ∂θ  (A.12)  This can also be obtained by going directly from Cartesian coordinates to polar coordinates x = r cos θ  (A.13)  y = r sin θ so x2 + y2 y = tan−1 x  r = θ  (A.14)  We will need the derivative of tan−1 which can be found by using a standard trick for finding derivatives of inverse trigonometric functions. Define b = tan−1 a, so 149  a = tan b, and we can write d tan−1 a db = = da da  da db  −1  d tan b db  =  −1  =  1 1 = 2 1 + tan b 1 + a2  (A.15)  Now calculate the derivatives we will need for the chain rule ∂r ∂x ∂r ∂y  x = cos θ r y = sin θ r  = =  (A.16)  and ∂θ ∂x ∂θ ∂y  = =  sin θ y y =− 2 =− x r r x cos θ y = 2= x r r  ∂ tan−1 ∂x ∂ tan−1 ∂y  (A.17)  Using the chain rule, we find ∂r ∂ ∂θ ∂ + ∂x ∂r ∂x ∂θ ∂ sin θ ∂ = cos θ − ∂r r ∂θ  ∂ ∂x  =  ∂ ∂y  = sin θ  (A.18)  and similarly ∂ cos θ ∂ + ∂r r ∂θ  By combining the last two equations and rewriting using eiθ = cos θ + i sin θ we find ∂ ∂ −i ∂x ∂y ∂ ∂ +i ∂x ∂y  =  e−iθ  =  eiθ  which is identical to (A.12).  150  ∂ i ∂ − ∂r r ∂θ ∂ i ∂ + ∂r r ∂θ  (A.19)  A.2  Exploiting rotational symmetry  A.2.1  Spinless case  In general, when a Hamiltonian has an n-fold rotational symmetry around some axis, it can be block diagonalized into n blocks. This arises from the fact that the Hamiltonian commutes with the rotation operator U (i.e. [H,U] = 0) that generates the above rotation. By block diagonalizing the Hamiltonian, we break the problem down into several smaller diagonalization problems. Since direct diagonalization is O(N 3 ), for large matrices this has the potential of drastically reducing the time and memory needed to diagonalize a given Hamiltonian. In particular, 3D problems often involve large matrices. For example, in our model, a cubic system of side N0 requires a matrix of side 4(N0 )3 (the factor of 4 comes from spin and orbit degrees of freedom). Even for N0 = 20 we must diagonalize matrices of side 32, 000, which is beyond the computation power, memory and time for most small computers. Exploiting the four-fold rotational symmetry of the problem, reduces this to diagonalization of four matrices of side 8, 000, which can be done easily. In addition, the memory necessary to store such matrices is large. A complex matrix of side 4(N0 )3 with double precision, requires 28 (N0 )6 bytes of memory (the additional factor of 24 comes from the 8 bytes necessary and the real and imaginary parts). Therefore, for N0 = 20 we need 214 = 8192 MB’s of memory, i.e. 8.2 GB, once again a significant amount of memory. To store four matrices of side (N0 )3 we need a quarter of that memory, i.e. just over 2 GB, which is within reason. In our problem, the Hamiltonian is 4-fold rotational symmetric around an axis parallel to the Dirac string. We treat this problem in detail, but other rotational symmetries could be treated similarly by followingt the same steps. We label the quadrants from 1 to 4 going clockwise around the system, and choose the numbering within each quadrant so that a rotation of π/2 around the axis of symmetry keeps the numbering intact. If we look at a particular wavefunction Ψ we can  151  separate the it by quadrant, and label the part in quadrant i by Ψi so that   Ψ1       Ψ2   Ψ= Ψ   3 Ψ4 and the Hamiltonian will have the form  h11 h12  h21 h22 H = h  31 h32 h41 h42  (A.20)  h13 h14     h23 h24   h33 h34   h43 h44  (A.21)  where hi j are matrices of side N03 , with elements connecting quadrants i and j. Using the fourfold rotational symmetry, we find h11 = h22 = h33 = h44 ≡ h  (A.22)  h12 = h23 = h34 = h41 ≡ g  h13 = h24 = h31 = h42 ≡  f  We require h to be Hermitian, giving hi j = h†ji so that h and f are Hermitian, but g is in general not Hermitian   h  g  f   † g h g H =  f g† h  g f g†  g†     f  g  h  (A.23)  We separate g into a Hermitian part a (so a = a† ) and an anti-Hermitian part b (so  152  b† = −b), therefore g = a + b and g† = a − b, and   h  a+b  f  a−b      a − b  h a + b f  = I4 ⊗ h + Ma ⊗ a + Ma ⊗ b + M f ⊗ f H =  f a−b h a + b   a+b f a−b h (A.24) where  Ma  Mb  Mf   0 1 0  1 0 1 =  0 1 0  1 0 1  0 1  −1 0 =   0 −1  1 0  1   0 =  1 0  0 1   1  0  = (I2 + σ1 ) ⊗ σ1 1  0  0 −1  1 0  = i(I2 − σ1 ) ⊗ σ2 0 1  −1 0  0  1  = σ1 ⊗ I2    (A.25)  Note that (I2 + σ1 )(I2 − σ1 ) = 0, so that [Ma , Mb ] = 0. In addition, M f commutes with the other two [M f , Ma ] = [M f , Mb ] = 0. Therefore these three matrices can  be diagonalized simultaneously. We do this in two steps. First we note that these matrices contain only the identity and σ1 in the left subspace, so that we can rotate  153  around σ2 by π/2, by using u = e−iπ/4σ2 ⊗ I2 to get M˜ a = (I2 + σ3 ) ⊗ σ1 =  2σ1 0 0  0  M˜ b = i(I2 − σ3 ) ⊗ σ2 =  0  0  M˜ f  = σ3 ⊗ I2 =  (A.26)  0 2iσ2  I2  0  0  −I2  Note that these matrices are all block diagonal, and the third is diagonal. The nonzero blocks of the first two matrices can be diagonalized separately: 2σ1 is diagonalized by 1 U1 = √ 2  −1 1  (A.27)  1 U2 = √ 2  i  −i  (A.28)  1  1  and 2iσ2 is diagonalized by  1  1  We combine these two matrices into a block diagonal matrix that diagonalizes M˜ a and M˜ b . Note that M˜ f is composed of the identity and minus the identity in those two blocks, so it is not affected by this rotation. Explicitly, U1†  0  0  U2†  I  0  0 −I  0  U1 0  U2  =  U1†U1  0  0  −U2†U2  I  =  0  0 −I  (A.29)  Putting the two rotations together we find  U =u  U1  0  0  U2  −i π4 σ2  = e  ⊗ I2  −1 1 0   1 1  √   2 0 0  154  0     0  0 i −i  0 1 1 1 0  (A.30)  and U block diagonalizes H so that    H˜ = U HU =     h − 2a + f  †   h + 2a + f h − 2ib − f  h + 2ib − f     (A.31)    Where h can be obtained by considering only hopping within a certain quadrant, g contains only the hopping elements between sites in quadrant one and two (for example) and f contains only hopping elements between sites in quadrants one and three (for example). Constructing h is trivial, but for the other matrices, it is helpful to consider the rotation operator R˜ which has the property ˜ 1 = Ψ2 RΨ  (A.32)  R˜ 2 Ψ1 = Ψ3 R˜ 3 Ψ1 = Ψ4 and R˜ 4 = I4 . Therefore, we can write ˜ 1) Ψ†1 gΨ2 = Ψ†1 g(RΨ Ψ†1 f Ψ3  =  (A.33)  Ψ†1 f (R˜ 2 Ψ1 )  which shows that we must make an additional transformation in this case. This transformation can also be thought of differently. If our full Hamiltonian is 4N ×  4N, then h, g and f are all N × N, and the vector Ψ is of length 4N. To form g we are connecting elements in the first quarter of Ψ with elements in the second quarter of Ψ. Using our normal ordering with no transformations (i.e.,Ψ†1 gΨ2 ) , this would necessarily result in a matrix of size 2N × 2N. The transformation  brings the indices from the second quadrant back into the first quadrant of indices. We now diagonalize each of the four blocks separately. We find a list of 4N03 energies. The energies do not change under rotation, so these are also the energies of the original Hamiltonian H. We also find four sets of N03 wavefunctions, and ˜i we label the four matrices which have these wavefunctions as their columns as Ψ 155  (i = 1, . . . , 4), and construct the matrix of wavefunctions of H˜  ˜1 Ψ  ˜2  Ψ ˜ = Ψ  ˜3 Ψ     ˜4 Ψ        (A.34)  To find the wavefunctions of the original Hamiltonian H, we rotate the wavefunc˜ (thereby rotating each of the wavefunctions separately), so that the tion matrix Ψ matrix Ψ contains the 4N03 wavefunctions of H. ˜ Ψ = UΨ  (A.35)  This can be seen by taking HΨ = EΨ, inserting the identity UU † in between H and Ψ, and multiplying both sides by U † on the left, so that (U † HU)(U † Ψ) = E(U † Ψ), ˜ = U † Ψ, so that Ψ = U Ψ. ˜ and we identify Ψ An equivalent but faster method of finding the matrix U is to diagonalize the explicit form of the n-fold rotational operator. In this case it is easy to see that  0  0 R= 0  1   1 0 0  0 1 0  0 0 1  0 0 0  (A.36)  which rotates the system by π/2, thereby taking quadrant one to quadrant two (for example), and in general     Ψ1 Ψ2     Ψ2  Ψ3     R Ψ  = Ψ   3  4 Ψ4 Ψ1  (A.37)  Also note that [H, R] = 0, or equivalently H = R† HR. If two matrices commute, then if one of them is diagonal in a particular basis, the other will be also be diagonal in this basis. Therefore, if we find a basis that diagonalizes R, we can 156  diagonalize H by rotating to this same basis. Therefore, we diagonalize R and construct the matrix whose columns are the normalized wavefunctions of R. This matrix is identical to the matrix U that was found above, up to the order of the columns. Again, we find that U † HU is diagonal.  A.2.2  Adding the spin  If the particles in the system are not spinless (or spin polarized), the rotation operator must also rotate spin. So we write the Hilbert space as an external product of a spin subspace and a spatial subspace. To rotate spin 1/2 around the z axis by an angle of π/2, we use the rotation operator π 1 1 Rs = e−i 4 σ3 = √ (I − iσz ) = √ 2 2  1−i 0  0 1+i  (A.38)  We can multiply a rotation by a phase with no effect on the result of the rotation, since when rotating we multiply by Rs on one side and R†s on the other side, so the √ phase cancels. For convenience we multiply the above operator by (1 + i)/ 2 to get 1 0  Rs =  0  i  (A.39)  We now impose [Rs ⊗ R, H] = 0, where R is the operator that rotates the spatial  subspace, as given in (A.36). We seek the form of H which will fulfill this equation.  The spin subspace of the Hamiltonian can be written as a linear combination of Pauli matrices and the identity 4  H = ∑ σi ⊗ Hi  (A.40)  i=0  where σ0 ≡ I2 , the identity matrix of rank two. Using this decomposition, we impose the above commutation condition, to find  ∑[Rs ⊗ R, σi ⊗ Hi ] = 0 i  157  (A.41)  Note that Rs can be written as I2 − iσz (up to a constant phase), so for i = 0, 3 the spin subspace commutes with the spin rotation operator [Rs , σi ] = 0, so that it is  sufficient to solve for H that commutes with the spatial rotation matrix [R, Hi ] = 0. This is the same equation we solved for the spinless, therefore we can use the same structure for H0 and H3 which is given by (A.23). So that for i = 0, 3 we have   hi  ai + bi  ai − bi  0      ai − bi hi ai + bi 0   = I4 ⊗ hi + Ma ⊗ ai + Mb ⊗ bi (A.42)  Hi =  ai − bi hi ai + bi    0 ai + bi  0  ai − bi  hi  Note that we have set fi = 0, as appropriate for our model, since there are no hopping elements that connect the quadrants one and three, for example. From now on we assume this is the case, purely for convenience. For a model that does include such elements, there will be additional matrices fi that will have to be constructed. For i = 1, 2 the spin subspace anticommutes with the spin rotation operator {Rs , σi } = 0, so one solution would be to look for H that anticommutes with the  spatial rotation matrix {Hi , R} = 0. In order to do this, it is convenient to write R as an external product of Pauli matrices,  1 R = [I2 ⊗ (σ1 + iσ2 ) + σ1 ⊗ (σ1 − iσ2 )] 2  (A.43)  For example, it is obvious that I ⊗σ3 and σ1 ⊗σ3 anticommute with R. The solution  for Hi would then be a linear combination of matrices that anticommute with R.  However, in our case Rs takes σ1 to σ2 so that we expect H1 and H2 to be coupled. This can also be written as RH1 R† = H2  (A.44)  We can postulate a form for H1 , and then find the form for H2 using the above  158  simple relation. After some trial and error, we find  H1  H2  g2  0   †  g2 h2 =   0 g†1  −g1 0  h2 g1  † g1 −h1 =   0 −g†  2 g2 0  g1  h1  −h1 −g†2 0 −g2 −h2 −g†1  −g†1     0   −g2   −h2  g†2  0   −g1   h1  (A.45)  (A.46)  Note that we can go from H1 to H2 by taking h1 → h2 , h2 → −h1 , g1 → −g2 , and  g2 → g1 .  In practice, we find hi by considering only one quadrant, and constructing the  matrix h which contains all the matrix elements for hoppings within that quadrant. Then we observe that h = ∑ σi ⊗ hi  (A.47)  i  In order to proceed, we must make sure to order our sites first by spin, so that the top of the spinor has only spin up and the bottom of the spinor has only spin down (or vice versa) Ψ=  Ψ↑  (A.48)  Ψ↓ so that we can label the quadrants of h by spin h=  h↑↑ h↑↓ h↓↑ h↓↓  159  (A.49)  and then we can find hi by decomposing h h = σ0 ⊗  h↑↑ + h↓↓ 2  + σ1 ⊗  h↑↓ + h↓↑ 2  + σ2 ⊗  h↓↑ − h↑↓ 2i  + σ3 ⊗  h↑↑ − h↓↓ 2 (A.50)  So that h0 = h1 = h2 = h3 =  h↑↑ + h↓↓ 2 h↑↓ + h↓↑ 2 h↓↑ − h↑↓ 2i h↑↑ − h↓↓ 2  (A.51)  Similarly, we construct the matrix g which includes all matrix elements of hopping from quadrant 1 to quadrant 2, and then find gi by the same recipe as above. When constructing the matrix g there is a complication, which is resolved as in the previous section, in (A.33). However, we must remember to now include also the spin ˜ rotation in the operator R. To find the matrix that block diagonalizes the Hamiltonian, we diagonalize the full rotation matrix RR ≡ Rs ⊗ R  160  (A.52)  and form the matrix that has its eigenvectors as its columns  1    0 1 0 −1 0 −1 0 1    0 −1 0 −i 0 i 0 1      1 0 1 0 1 0 1 1 0  U=  2 i 0 1 0 −1 0 −i 0     −1 0 1 0 1 0−1 0    −i 0 1 0 −1 0 i 0   1 0 1 0 1 0 1 0   0  −1 0  i  0  −i  0  (A.53)  Usually when forming a rotation matrix by putting the eigenvectors as the columns, the order of the eigenvectors does not matter. However, here it does. The reason is that this matrix block diagonalizes the Hamiltonian, as opposed to the more commonly encountered case (in Physics at least) where U fully diagonalizes the Hamiltonian. Swapping eigenvectors in the full diagonalization problem simply swaps the order of the eigenvalues in the diagonal matrix. However, in our problem (and any problem where U does not diagonalize H), swapping the order of the eigenvectors will in general change the form of the block diagonal matrix. In our case the eigenvectors are guaranteed to come in degenerate pairs, and swapping the order of two pairs will just change the order of the blocks in the block diagonal matrix. Therefore, we must make sure to order the eigenvectors by their eigenvalues. For example, in the above case, the order (−1, −1, i, i, −i, −i, 1, 1) was used. As mentioned, the matrix U block diagonalizes the Hamiltonian, so that    B1   B2 U † HU =   B3        B4  161  (A.54)  where B1 = B2 = B3 = B4 =  h0 − h3 + 2i(b0 − b3 ) ih1 − h2 + a˜1 + a˜∗2 + b˜ ∗1 − b˜ 2 −ih1 − h2 + a˜∗1 + a˜2 − b˜ 1 + b˜ ∗2 h0 + h3 − 2(a0 + a3 )  h0 − h3 + 2(a0 − a3 ) ih1 − h2 − a˜1 − a˜∗2 + b˜ ∗1 − b˜ 2 −ih1 − h2 − a˜∗1 − a˜2 − b˜ 1 + b˜ ∗2 h0 + h3 + 2i(b0 + b3 )  h0 − h3 − 2(a0 − a3 ) ih1 − h2 + a˜1 + a˜∗2 − b˜ ∗1 + b˜ 2 −ih1 − h2 + a˜∗1 + a˜2 + b˜ 1 − b˜ ∗2 h0 + h3 − 2i(b0 + b3 ) h0 − h3 − 2i(b0 − b3 ) ih1 − h2 − a˜1 − a˜∗2 − b˜ ∗1 + b˜ 2 −ih1 − h2 − a˜∗1 − a˜2 + b˜ 1 − b˜ ∗2 h0 + h3 + 2(a0 + a3 )  (A.55)  where b˜ i = (1 + i)bi , a˜i = (1 + i)ai , we have decomposed gi = ai + bi into Hermitian and anti Hermitian parts (respectively), as before, and this structure of Bi is dependent on the order of eigenvectors chosen in U. As done for the spinless case, we now diagonalize each of the four blocks Bi separately. We find a list of 4N03 energies. The energies do not change under rotation, so these are also the energies of the original Hamiltonian H. We also find four sets of N03 wavefunctions, and we label the four matrices which have these ˜ i (i = 1, . . . , 4), and construct the matrix of wavefunctions as their columns as Ψ wavefunctions of H˜    ˜ = Ψ      ˜1 Ψ ˜2 Ψ ˜3 Ψ ˜4 Ψ        (A.56)  To find the wavefunctions of the original Hamiltonian H, we rotate the wavefunc˜ (thereby rotating each of the wavefunctions separately), so that the tion matrix Ψ matrix Ψ contains the 4N03 wavefunctions of H. ˜ Ψ = UΨ  (A.57)  In summary, here is a possible algorithm for diagonalizing such a system: 162  1. Construct h and g and find hi and gi using (A.51) 2. Construct the four blocks Bi and diagonalize each separately (A.55) 3. Form a list of energies and the matrix of wavefunctions and sort them simultaneously by energy ˜ 4. Rotate the matrix of wavefunctions Ψ = U Ψ  A.3  Adding monopoles in practice: calculating the Peierls factors  A.3.1  Adding a lattice of monopoles  If we wish to construct a system with a monopole and study it numerically, it is most convenient to work with periodic boundary conditions. We construct a cubic unit cell that has a monopole at its center and eight anti-monopoles at its corners. Since each corner is shared by eight cubes, the total magnetic charge per unit cell is zero g(1 − 8 × 18 ) = 0. Each unit cell also includes a Dirac string going from an  anti-monopole at a corner to the monopole at the center. This setup can be viewed as a bar magnet - the field lines go out of the monopole, into the anti-monopole, and through the Dirac string back to the monopole. We define our coordinates so that the monopole is at (0, 0, 0) and the Dirac string goes to the anti-monopole at a 2 (1, 1, 1).  The singular part of the magnetic field of one monopole-antimonopole  pair is Bs,0 (r) = −4πgδ (x − y)δ (y − z) θ (z) − θ z −  a 2  lˆ  (A.58)  √ where lˆ ≡ 1/ 3(1, 1, 1) is the vector pointing from the monopole in the direction  of the Dirac string. The two delta functions enforce x = y = z and the heaviside functions set the ends of the string. For a periodic array of cells with direct lattice vectors R we have Bs (r) = ∑ Bs,0 (r − R) R  163  (A.59)  where R = a(nx , ny , nz ) and nα ∈ Z. The magnetic field is obviously translationally invariant Bs (r + R) = Bs (r) for all direct lattice vectors R. We now take the Fourier transform to find Bsk Bsk = F Bs (r)  Bs (r) e−ik·r d 3 r =  =  ∑ Bs,0 (r − R) e−ik·r d 3 r (A.60) R  =  ∑ e−ik·R  Bs,0 (r) e−ik·r d 3 r  R  =  ∑ e−ik·R F  Bs,0 (r) = ∑ δk,G F Bs,0 (r)  R  G  We first find the Fourier transform of Bs,0 . Two of the three integrals are trivial due to the two delta functions, and the third integral (over the heaviside functions) is just an integral of an exponent over the section (0, a/2) Bs,0k = F Bs,0 (r)  Bs,0 (r) e−ik·r d 3 r  =  = −4πglˆ = −4πglˆ = −4πglˆ  (A.61)  δ (x − y)δ (y − z) θ (z) − θ z − ∞ −∞ a/2 0  θ (z) − θ z −  a 2  a 2  e−ik·r d 3 r  ˜  e−ikz dz  ˜  ˜a  e−ikz dz = −2πga e−ik 4 sinc  ˜ ka lˆ 4  where k˜ ≡ kx + ky + kz and sinc(x) ≡ sin(x)/x, so that Bsk =  ∑ δk,G F  Bs,0 (r)  (A.62)  G  = −2πga ∑ δk,G e−ik 4 sinc  ˜ ka 4π  = −2πga ∑ δk,G e−i 2 m˜ sinc  π m˜ lˆ 2  ˜a  G  π  G  since G =  2π a (mx , my , mz ),  lˆ  mα ∈ Z and we have defined m˜ ≡ mx + my + mz . As ex-  pected, only the Fourier components which correspond to reciprocal lattice vectors  164  are non-zero. We write only the non-zero terms by taking k = G π m˜ lˆ 2  BsG = −2πga e−i 2 m˜ sinc π  (A.63)  Note that the average magnetic field is given by the zero Fourier component BsG=0 = −2πgalˆ  (A.64)  The vector potential’s Fourier components (for G = 0) can be found from (4.32) AG =  iG × BsG = −2πiga G2  G ˆ −i π m˜ π × l e 2 sinc m˜ 2 G 2  (A.65)  For G = 0 we note that G2 AG = iG × BsG  (A.66)  Below we will prove that AG = A∗−G . Therefore, A0 is real, so that on the left side we have a real expression. The average magnetic field BsG=0 is real and is given  by (A.64). Therefore, the right side is pure complex. We replace the continuous Fourier transform with a discrete Fourier transform. This can be seen either by replacing the integral by a sum d 3 k → ∑k (2π/a)3 , or  by calculating  A(r) =  1 (2π)3  =  1 (2π)3  Ak eik·r d 3 k  ∑ δk,G AG G 3 (2π)  =  1 (2π)3 a3  =  1 A eiG·r a3 ∑ G  ∑ G  G  165  (A.67) eik·r d 3 k  δ (k − G)AG eik·r d 3 k  From here we also see that since A(r) is real, we must have AG = A∗−G since 1 A∗ e−iG·r a3 ∑ G  A∗ (r) =  (A.68)  G  1 1 A∗−G eiG·r = 3 ∑ AG eiG·r ∑ 3 a a  =  G  G  and indeed (A.65) fulfills this condition. We calculate the Peierls phase factors on an underlying cubic lattice with lattice constant a0 . This is done by computing the line integral of the vector potential from site c to site d. For simplicity, we compute the integral for two nearest neighbors θcd =  2π Φ0  = −ig  Rd Rc  A(r) · dr =  (2π)2 a2 Φ  1  ∑ G2 e 0  1 2π A · nˆ a3 Φ0 ∑ G G  −i π2 m˜  G  sinc  π m˜ 2  Rd  eiG·r dr  (A.69)  Rc  G × lˆ · nˆ eiG·Rc  Rd −Rc  eiG·r dr  0  where the unit vector nˆ points from site c to site d. Look at    (m − mz )dx   y  2π G × lˆ · ndr ˆ = √ (mz − mx )dy a 3   (m − m )dz x y  nˆ = xˆ nˆ = yˆ  (A.70)  nˆ = zˆ  √ recall that lˆ ≡ (1, 1, 1)/ 3. Compute the integral for example for Rc = a0 (cx , cy , cz ), Rd = Rc + a0 xˆ so nˆ = x, ˆ dr = dx and y = z = 0 in the integral G × lˆ · nˆ  Rd −Rc  eiG·r dr =  0  =  166  a0 2πi 2π √ (my − mz ) e a mx x dx 0 a 3 1 my − mz 2πimx a0 a −1 √ e i 3 mx  So that a0 1 my − mz −i π m˜ π (2π)2 g x e 2 sinc m˜ eiG·Rc e2πimx a − 1 θcd = −√ ∑ 2 2 mx 2 3a Φ0 G G  (A.71) which can be rewritten by putting in another sinc function, and cancelling the imaginary part x θcd  ga0 = −√ 3Φ0 2ga √ 0 3Φ0  =  3  2π a 2π a  ∑ G=0 3  ∑ |G|>0  a0 a0 π 1 ˜ i(G·Rc − π2 m+πm x a ) (m − m )sinc m ˜ sinc πm (A.72) i e y z x G2 2 a  π 1 a0 π a0 (my − mz )sinc m˜ sinc πmx sin G · Rc − m˜ + πmx 2 G 2 a 2 a  which can be seen by noting that i  ∑ |G|=0  cG eidG  = i =  ∑  cG [cos(dG ) + i sin(dG )]  (A.73)  |G|=0  ∑ [icG cos(dG ) − cG sin(dG )] = −2 ∑  |G|=0  cG sin(dG )  |G|>0  where cG and dG are odd in G, and the cos term cacelled since it is an odd function summed over an even range, and the sin term is even so it gives twice the value for positive G. Note that by G > 0 we mean that we sum over only a half space, defined for example as mz > 0. Also, by dimensional analysis we see that θ ∑G 1/G3  A.3.2  G2 dG/G3  log G. So we expect slow convergence.  Adding a single monopole  We wish to add a magnetic monopole to the lattice system by calculating the Peierls factors. We start by finding the magnetic field of a single magnetic monopole, using Gauss’ law. The magnetic field is spherically symmetric, so B(r) = B(r)ˆr. Therefore, if we calculate the flux through a spherical Gaussian surface of radius r  167  centered on the monopole, we find B · ds = 4πr2 B(r) = 4πg  (A.74)  so that the magnetic field of a single monopole is B(r) =  g r rˆ = g 3 2 r r  (A.75)  We wish to add a monopole of flux 2π (flux is dimensionless in these units), so that 4πg = 2π and therefore g = 1/2. We now find the vector potential. The magnetic field is spherically symmetric, so we can write A(r) = Aθ ∇ϕ  (A.76)  and find the magnetic field by taking the curl, using a vector identity, noting that the curl of a gradient is zero and substituting the gradient and curl in spherical coordinates B(r) = ∇ × A = ∇ × (Aθ ∇ϕ)  (A.77)  = ∇Aθ × ∇ϕ + Aθ ∇ × ∇ϕ 1 ∂ Aθ = 2 rˆ = B(r)ˆr r sin θ ∂ θ  We solve B(r) =  g 1 ∂ Aθ = 2 r2 sin θ ∂ θ r  (A.78)  To find, after choosing an arbitrary constant Aθ = −g(1 + cos θ )  (A.79)  A(r) = −g(1 + cos θ )∇ϕ  (A.80)  So  168  which can be rewritten using ∇ϕ =  1 1 ϕˆ = (− sin ϕ, cos ϕ, 0) r sin θ r sin θ  (A.81)  g(1 + cos θ ) (sin ϕ, − cos ϕ, 0) r sin θ  (A.82)  so that A(r) =  which is singular for θ = 0. The vector potential can also be rewritten in a form which is sometimes referred to as the ”Dirac Pontential” A(r) =  g r × nˆ g yxˆ − xyˆ = r r − r · nˆ r r−z  (A.83)  where r = (x, y, z), r = |r|, nˆ = (0, 0, 1) so that the Dirac string stretches in the zˆ direction to infinity. The vector potential is singular for x = y = 0 if z ≥ 0, since in  this case r = |z| so r − z = |z| − z = 0 for z ≥ 0.  We want to calculate the Peierls phase factors given by Θcd =  Rd Rc  A · dl  (A.84)  where R are the lattice vectors and the integral must be done for all nearest neighz = 0 for all bors. Note that for d l = dzˆz the integral is zero, since Az = 0, so that θcd  bonds. For d l = dxxˆ and Rd − Rc = ax, ˆ we have y = const ≡ yc and z = const ≡ zc , r=  x2 + y2c + z2c so  Θxcd = gyc  xc +a xc  = g tan−1 + tan−1  gyc xc +a 1 1 dx = dx − r(r − z) zc xc r − zc r xc + a x c − tan−1 yc yc (xc + a)zc yc  (xc + a)2 + y2c + z2c  169  − tan−1  yc  (A.85)  xc zc 2 xc + y2c + z2c  where we have used the identity 1 1 = r(r − z) z  1 1 − r−z r  (A.86)  and the integrals dx r dx r−z  = ln(x + r) z tan−1 y  =  (A.87) x + tan−1 y  xz yr  + ln(x + r)  (A.88)  Similarly, for d l = dyyˆ y θcd = −gxc  yc +a yc  = −g tan−1 + tan−1  gxc dy =− r(r − z) zc yc + a − tan−1 xc  xc  yc +a  dx yc  yc xc  (yc + a)zc 2 xc + (yc + a)2 + z2c  1 1 − r − zc r  − tan−1  xc  (A.89)  yc zc 2 xc + y2c + z2c  x by replacing g → −g and x ↔ y. which can be obtained from θcd  A.3.3  Adding a planar monopole  We wish to add a planar magnetic monopole to the lattice system by calculating the Peierls factors. We assume that the planar magnetic monopole is localized on the two closest layers. We start by finding the magnetic field using Gauss’ law. The magnetic field is spherically symmetric, so B(r) = B(r)ˆr. Therefore, if we calculate the flux through a cylindrical Gaussian surface of radius r and height a (the spacing) centered on the planar monopole, we find B · ds = 2πraB(r) = 4πg  170  (A.90)  so that the magnetic field of a planar monopole is 2g r rˆ = 2g ra a  B(r) =  (A.91)  on the two closest layers, and zero elsewhere. We wish to add a planar monopole of flux 2π (flux is dimensionless in these units), so that 4πg = 2π and therefore g = 1/2. We choose a vector potential A=  2g ϕ zˆ a  (A.92)  on the two closest layers and zero elsewhere, where ϕ is the angular coordinate in cylindrical coordinates. Note that the vector potential has a discontinuity at ϕ = 0 since the angle jumps from 0 to 2π, leading to the conclusion that the Dirac string for this monopole lies along this line. To check, we find the magnetic field from this vector potential by taking the curl in cylindrical coordinates 2g 1 B = ∇×A = a r  rˆ  rϕˆ  zˆ  ∂ ∂r  ∂ ∂ϕ  ∂ ∂z  0  0  ϕ  =  2g rˆ ra  (A.93)  We want to calculate the Peierls phase factors given by Rd  Θcd =  Rc  A · dl  (A.94)  for d l = dxx, ˆ dyyˆ get Θxcd = Θycd = 0 since the vector potential is zero in the x and y directions. In addition, the vector potential in the z direction is zero except on the two closest layers, so that Θzcd = 0 for all bonds except those that cross the plane of the planar monopole. We calculate the non-zero Peierls factors Θzcd =  2g a  zc +a  ϕdz = 2gϕc  (A.95)  zc  Note that for g = 1/2 we find Θzcd = ϕ, so that the Peierls factors are simply equal to the angle in cylindrical coordinates. For example, the four vertical bonds closest 171  to the monopole will have the Peierls factors 0, π/2, π, 3π/2 which are their angles with respect to the line connecting the monopole and the first of these bonds. This result can also be written in cartesian coordinates by noting that x = r cos ϕ and y = r sin ϕ so that yc xc  Θzcd = 2g tan−1  (A.96)  In order to avoid problems arising from using tan−1 , we can calculate eiΘ directly without the use of trigonometric functions iΘzcd  e  A.4  =e  2gϕc i  xc + iyc  =  2g  xc2 + y2c  (A.97)  Derivation of the self consistent Born approximation  In this section we follow [4] sections 5.3-5.4. First we derive Dyson’s equation. We can write the Green’s function (with disorder) as a sum of irreducible diagrams which depend on the self energy Σ and the bare Green’s function G0 (single electron, no disorder) G = G0 + G0 ΣG0 + G0 ΣG0 ΣG0 + . . .  (A.98)  = G0 + G0 Σ(G0 + G0 ΣG0 + . . . ) = G0 + G0 ΣG where G,G0 and Σ depend on k and can be scalars or matrices. In the last step we notice that the expression in the parenthesis is equal to the first line and so we replace it with G. We can now solve this equation for G. We rearrange G0 = G − G0 ΣG0 = (I − G0 Σ)G  (A.99)  and multiply by G−1 0 on the left I = (G−1 0 − Σ)G 172  (A.100)  so we have a closed expression for the Green’s function as a function of the self energy and bare Green’s function −1 G = (G−1 0 − Σ)  (A.101)  By definition, the bare Green’s function G0 = (εF − H0 + i0+ )−1  (A.102)  We substitute in the previous equation, to find G = (εF − H0 − Σ + i0+ )−1  (A.103)  This is the exact one particle Green’s function with disorder. Now our objective is to calculate the self energy. We write it as an infinite sum of irreducible diagrams. In order to evaluate this sum, we make some approximations. The first approximation is to assume a low concentration of impurities nimp ≡ Nimp/V , and then keep only terms that are first order in nimp . Since each scat-  tering event carries a factor of nimp , this limits us to diagrams with one impurity. This means that we are assuming that impurities are far enough away one from another, that electrons scattering off one impurity are not simultaneously scattering off another impurity. The second approximation is the weak scattering limit. We assume that the scattering potential is weak enough, so that at most two electrons scatter off a single impurity simultaneously. This limits the sum of diagrams further, and we are left with just two diagrams: the first with a single electron and a single impurity (a line), and the second with two electrons and a single impurity (a triangle). The first order diagram gives Σ1 = NimpU(k = 0) = Nimp  1 V  U(x)d d x = nimp  U(x)d d x  (A.104)  which simply shifts the energies by a constant, which is proportional to the spatial average of the potential. Therefore, it can be absorbed into the definition of the energies. In some cases (such as our case), the average is zero, and then we don’t 173  need to worry about it at all. The second diagram (the triangle), gives Σ2 = Nimp ∑ U 2 (k − k )G0 (k )  (A.105)  k  which is the self energy in this approximation. In the self consistent Born approximation, we now replace G0 by G in this expression, to find Nimp ∑ U 2 (k − k )G(k )  Σ  (A.106)  k  = Nimp ∑ U 2 (k − k )[εF − H0 (k ) − Σ(k ) + i0+ ]−1 k  where in the last line we have substituted (A.103). For point-like impurities, the interaction is a delta function U(x) = uδ (d) (x), so U(k) =  1 V  uδ (d) (x) eik·x d d x = u  (A.107)  Therefore, the self energy simplifies Σ = Nimp u2 ∑[εF − H0 (k) − Σ(k) + i0+ ]−1  (A.108)  k  It is often convenient to evaluate this sum in the large system limit, where the sum becomes an integral Σ = Nimp  a 2π  d  u2  d d k[εF − H0 (k) − Σ(k) + i0+ ]−1  (A.109)  Note that these sums and integrals are all over the full momentum space. For a periodic system, momentum space can be written as N 3 copies of the BZ, so we can write Σ = Nimp N d u2  ∑ [ε  F  k∈BZ  − H0 (k) − Σ(k) + i0+ ]−1  174  (A.110)  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073186/manifest

Comment

Related Items