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 θ 6= 0 is predicted to bind a fractional electric charge. We con- duct 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/2pih)~B ·~E, and that θ = pi 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 disor- der 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 topo- logical ‘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 experi- ments. ii Preface This thesis is based almost exclusively on extensive notes written by myself through- out 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 super- visor 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 re- search 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 topo- logical 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 super- visor’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 Bi2Se3 ” 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 ap- proximation 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 Re- view 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 Introduction to Topological Insulators . . . . . . . . . . . . . . . . 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 A brief history of topological insulators . . . . . . . . . . . . . . 1 2 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 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Our cubic model . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Explicit calculation of θ . . . . . . . . . . . . . . . . . . . . . . 36 4 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 5 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 6 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 Bi2Se3 . . . . . . . . . . . . . . . . 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Table 7.1 Lattice constants for Bi2Se3 . . . . . . . . . . . . . . . . . . . 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 Bi2Se3 . . . . . . . . . . . . . . . 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 topologi- cal 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 4.4 3D charge density plot for the layer just below the planer monopole for weak disorder . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 4.5 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 . . . . . . . . . . . . . . . . . . . 89 Figure 6.4 Quasi 1D system with a monopole and anti-monopole with pe- riodic 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 temperature113 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 coordi- nate, for different temperatures . . . . . . . . . . . . . . . . . 122 Figure 7.5 Impurity spin expectation value as a function of temperature for the discretized Bi2Se3 lattice . . . . . . . . . . . . . . . . 123 Figure 7.6 Critical temperature as a function of the chemical potential µ for the discretized Bi2Se3 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 Per- ovskite 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 tempera- ture 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 Grad- uate Entrance Scholarship, the International Partial Tuition Scholarship, the PhD Tuition Fee Award, the University Graduate Fellowship, the Pacific Century Grad- uate 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 atten- tion and research in the condensed matter community. The advance has been rapid, on both the theoretical and experimental fronts, and the excitement in the commu- nity 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 pub- lished 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 dimen- sions (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 sym- metry and a ferromagnet breaks the continuous rotation symmetry of the vacuum. In contrast, a topological phase includes all states which can be deformed contin- 1 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 perpendicu- lar 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-Nightingale- Nijs) can be calculated by a surface integral of the Berry flux over the Brillouin zone (BZ) n= 1 2pi∑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 Ek EF 0/a−pi Conduction Band Valence Band Quantum Hall State n=1 Insulator n=0 (a) (b) /a−pi 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 BInsulating State Quantum Hall State E k 0 E k EG (a) (b) (c) (d) (e) (f) /a−pi/a−pi 0 /a−pi/a−pi hωc 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 spin- orbit 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 con- 4 QH QSH spinless 1D chain spinful 1D chain 2=1+1 4=2+2 ba 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 6=−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 pi or −pi . The difference between these two paths is 2pi , which for spin 5 1/2 gives a relative minus sign between them. Therefore, when adding the con- tributions 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 topolog- ically 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. ba Figure 1.4: (a) Destructive interference on a lens with an antireflective coat- ing. 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 anti- clockwise, giving a spin rotation of ±pi , netting a difference of 2pi . When spin 1/2 particles rotate by 2pi 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 pres- ence 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. How- ever, 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 ap- parent that the native spin orbit coupling is too weak to support the topological order [21–25]. The search continued, focusing on heavier elements which typi- cally 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 topo- logical 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] ∆Laxion = θ ( e2 2pih ) ~B ·~E (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 θ → θ +2pi [31]. However, if TRI is present, only θ = 0,pi are allowed, since T takes ~E→ ~E and ~B→−~B. This led to 7 the realization that θ is in fact the Z2 invariant, θ = ν0pi (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θ/2pi . In this context, topological insulators present a platform for testing fundamental physical paradigms. As in two dimensions (see Eq. 1.2), the three- dimensional topological invariant can be calculated by an integral over the BZ of a Berry connection[29, 30] θ = 1 4pi ∫ BZ ε i jkTr [ Ai∂ jAk+ 2i 3 AiA jAk ] d3k (1.6) where A αβi = −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 cumber- some 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 =∏ m,i ξ2m(~Γi) (1.7) 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 topologi- cally classify our model. The first model shown to be a TI in 3D was a toy model on the diamond lat- tice, followed soon after by a prediction and observation that the semiconducting alloy Bi1−xSbx 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 can- didate for experiments, which led to the discovery of the so called “second gen- eration” 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 Bi2Se3 , with a relatively large band gap of ∼ 0.3eV (equivalent to about 3600K) [34, 35]. The spectrum of Bi2Se3 and other TIs has been studied using ARPES [35–37] (see Fig. 1.6) and scanning tunnelling microscopy (STM) [38–40], show- ing 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 Bi2Se3 , 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 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 that indeed time reversal invariance is fulfilled. 1 (a) (b) (a) (b) (a) (b) 0 20 40 60 80 0.0 0.5 1.0 1.5 2.0 2.5 T ￿S F M ￿(µ B )   surf1 surf2 z bulk T surfc T bulkc 0 20 40 60 80 ï2 ï1 0 1 2 T ￿S A F ￿(µ B )   surf1 surf2 z avg bulk T surfc T bulkc (a) (b) 0 20 40 60 80 100 120 ï1 0 1 2 T ￿S ￿(µ B )   avg2 avg4 bulk surf T surfc T bulkc 1 2 3 4 5 6 7 8 9 10 ï1 0 1 2 ￿S ￿(µ B )   0 40 60 70 80 90 100 120 Figure 1.5: (a) Energy spectrum of a massless Dirac fermion. The bottom (red) cone is the valence band, fully occupied at half filling. The top (green) cone is the conduction band, which is assumed to be vacant. (b) Massive Dirac fermion. This figure shows how the opening of a gap tends to lower the free energy. 9 NATURE PHYSICS DOI: 10.1038/NPHYS1274 LETTERS L FZ Γ Γ Γ Γ Γ M M M3 M2 M1 K M K K SS SS Low High EF0.1 0 ¬0.1 ¬0.2 ¬0.3 ¬0.4 ¬0.15 0 0.15 ky (Å ¬1) ¬0.15 0 0.15 Theoretical calculations kx (Å ¬1) ky (Å ¬1) E B  ( eV ) E B  ( eV ) 0.1 0 ¬0.1 ¬0.2 ¬0.3 ¬0.4 In te ns ity  ( a. u. ) ¬0.12 0 0.12 kz kx ky ky kx (111) 0.5 0 ¬0.5 ¬1.0 1.0 En er gy  ( eV ) Γ MK Bulk No SOC SOC No SOC SOC Surface a d e f b c 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 Bi2Se3(111). Electron dispersion data measured with an incident photon energy of 22 eV near the 0-point along the 0–M (a) and 0–K (b) momentum-space cuts. c, The momentum distribution curves corresponding to a suggest that two surface bands converge into a single Dirac point at 0. The V-shaped pure SS band pair observed in a–c is nearly isotropic in the momentum plane, forming a Dirac cone in the energy–kx–ky space (where kx and ky are in the 0–K and 0–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 Bi2Se3 and the two-dimensional BZ of the projected (111) surface. e, The surface Fermi surface (FS) of the two-dimensional SSs along the K–0–M momentum-space cut is a single ring centred at 0 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 0 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 0 and M; this ensures the existence of a pi 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 0 and M when SOC is included (red dotted lines). experiment. The ‘V’ bands cross EF at 0.09Å−1 along 0–M and at 0.10Å−1 along 0–K, and have nearly equal band velocities, approx- imately 5×105 ms−1, along the two directions. A continuum-like manifold of states—a filled U-shaped feature—is observed inside the V-shaped band pair. All of these experimentally observed features can be identified, to first order, by a direct one-to-one comparison with the LDA band calculations. Figure 1f shows the theoretically calculated (see the Methods section) (111)-surface electronic structure of bulk Bi2Se3 along the K–0–M k-space cut. The calculated band structure with and without SOC are overlaid together for comparison. The bulk band projection continuum on the (111) surface is represented by the shaded areas, blue with SOC and green without SOC. In the bulk, time-reversal symmetry demands E(k,↑) = E(−k,↓) whereas space inversion symmetry demands E(k,↑) = E(−k,↑). Therefore, all the bulk bands are doubly degenerate. However, because space inversion symmetry is broken at the terminated surface in the experiment, SSs are generally spin-split on the surface by spin–orbit interactions except NATURE PHYSICS | VOL 5 | JUNE 2009 | www.nature.com/naturephysics 399 © 2009 Macmillan Publishers Limited.  All rights reserved.  Figure 1.6: ARPES measurements of electronic dispersion on the [111] sur- face of Bi2Se3 . The incident photon energy was 22 eV, and data was collectd near the Γ-point along the Γ-M (a) and Γ-K (b) lines in momen- tum space. (c) The two surface lines meet at a Dirac point at Γ. Figure and caption from Ref. [35]. Breaking TRI, for example by adding magnetic dopants, is expected to open a gap in the spectrum of the surface states, provided that the magnetic moments are ordered perpendicular to the surface. The resulting spectrum then resembles that of a ‘massive’ Dirac fermion (Fig. 1.5b). There is considerable interest in having a system with an odd number of massive Dirac fermions, since it is predicted to exhibit many interesting topological phenomena, including the alf quantum Hall effect on the urface (e2/2h Hall conductance) [30], the image magnetic onopole (an electric charge adjacen to a TI results in the field of a magn tic monopole emb dded n the TI) [44, 45], and a Kerr/Faraday angle quantization in u it of the fine structur constant [46, 47]. A tunable gap ould lso allow th control of the surface transport, and could in a dition lead to unique practical applications associated with purely electric control of the surface magnetization [48, 49]. A signature of the massive Dirac fermion has been observed recently using ARPES in magnetically doped Bi2Se3 , [50, 51] although the interesting effects associated with it have yet to be seen in a laboratory. A surprising feature of these experiments 10 NATURE PHYSICS DOI: 10.1038/NPHYS1270 ARTICLES 0.15 0.10 0.05 0 ¬0.05 ¬0.10 ¬0.15 En er gy  ( eV ) En er gy  ( eV ) K Γ M 8 6 4 2 0 ¬2 ¬4 ¬6 ¬8 En er gy  ( eV ) En er gy  ( eV ) K Γ M 0.4 0.3 0.2 0 ¬0.1 ¬0.2 ¬0.3 0.1 K Γ M 8 6 4 2 0 ¬2 ¬4 ¬6 8 6 4 2 0 ¬2 ¬4 ¬6 0.2 0 ¬0.1 ¬0.2 ¬0.3 0.1 8 6 4 2 0 ¬2 ¬4 K Γ M Sb2Se3 Bi2Se3 Sb2Te 3 Bi2Te 3 a c d b 0.4 0.3 0.2 0 ¬0.1 ¬0.2 0.1 Figure 4 | Surface states. a–d, Energy and momentum dependence of the LDOS for Sb2Se3 (a), Sb2Te3 (b), Bi2Se3 (c) and Bi2Te3 (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 0 point as red lines dispersing in the bulk gap for Sb2Te3, Bi2Se3 and Bi2Te3. No surface state exists for Sb2Se3. be carried out on the other three materials, from which we see that Sb2Te3 and Bi2Te3 are qualitatively the same as Bi2Se3, whereas the SOCof Sb2Te3 is not strong enough to induce such an inversion. Topological surface states The existence of topological surface states is one of the most important properties of the topological insulators. To see the topological features of the four systems explicitly, we calculate the surface states of these four systems on the basis of an ab initio calculation. First we construct the maximally localized Wannier function (MLWF) from the ab initio calculation using the method developed by Marzari and co-workers21,22. We divide the semi- infinite system into a surface slab with finite thickness and the remaining part as the bulk. The MLWF hopping parameters for the bulk part can be constructed from the bulk ab initio calculation, and the ones for the surface slab can be constructed from the ab initio calculation of the slab, in which the surface correction to the lattice constants and band structure have been considered self-consistently and the chemical potential is determined by the charge neutrality condition.With these bulk and surfaceMLWFhopping parameters, we use an iterative method23,24 to obtain the surface Green’s function of the semi-infinite system. The imaginary part of the surface Green’s function is the local density of states (LDOS), from which we can obtain the dispersion of the surface states. The surface LDOSon the [111] surface for all four systems is shown in Fig. 4. For Sb2Te3, Bi2Se3 andBi2Te3, one can clearly see the topological surface states that form a single Dirac cone at the 0 point. In comparison, Sb2Se3 has no surface state and is a topologically trivial insulator. Thus, the surface-state calculation agrees well with the bulk parity analysis, and confirms conclusively the topologically non-trivial nature of the three materials. For Bi2Se3, the Fermi velocity of the topological surface states is vF' 5.0×105 ms−1, which is similar to that of the other two materials. Low-energy effective model As the topological nature is determined by the physics near the 0 point, it is possible to write down a simple effective Hamiltonian to characterize the low-energy long-wavelength properties of the system. Starting from the four low-lying states |P1+z ,↑ (↓)〉 and |P2−z ,↑ (↓)〉 at the 0 point, such a Hamiltonian can be constructed by the theory of invariants25 for the finite wave vector k. On the basis of the symmetries of the system, the generic form of the 4× 4 effective Hamiltonian can be written down up to the order of O(k2), and the tunable parameters in the Hamiltonian can be obtained by fitting the band structure of our ab initio calculation. The important symmetries of the system are time-reversal symmetry T , inversion symmetry I and three-fold rotation symmetry C3 along the z axis. In the basis of (|P1+z ,↑〉, |P2−z ,↑〉, |P1+z ,↓〉, |P2−z ,↓〉), the representation of the symmetry operations is given by T =K · iσ y ⊗ I2×2, I = I2×2⊗ τ3 andC3= exp(i(pi/3)σ z⊗I2×2), whereK is the complex conjugation operator, σ x,y,z and τ x,y,z denote the Pauli matrices in the spin and orbital space, respectively. By requiring these three symmetries and keeping only the terms up to quadratic order in k, we obtain the following generic form of the effective Hamiltonian: H (k) = 0(k)I4×4+ M(k) A1kz 0 A2k−A1kz −M(k) A2k− 00 A2k+ M(k) −A1kz A2k+ 0 −A1kz −M(k)  + o(k2) (1) with k± = kx ± iky , 0(k)= C+D1k2z +D2k2⊥ andM(k)=M −B1 k2z − B2k2⊥. By fitting the energy spectrum of the effective Hamiltonian with that of the ab initio calculation, the parameters in the effective model can be determined. For Bi2Se3, our fitting leads to M = 0.28 eV, A1 = 2.2 eVÅ, A2 = 4.1 eVÅ, B1 = 10 eVÅ2, B2= 56.6 eVÅ2, C =−0.0068 eV, D1= 1.3 eVÅ2, D2= 19.6 eVÅ2. Except for the identity term 0(k), the Hamiltonian (1) is nothing but the 3D Dirac model with uniaxial anisotropy along the z-direction and k-dependent mass terms. From the fact NATURE PHYSICS | VOL 5 | JUNE 2009 | www.nature.com/naturephysics 441 © 2009 Macmillan Publishers Limited.  All rights reserved.  Figure 1.7: Local density of states (warmer colours are higher values) as it depends on omentum nd energy for Sb2Se3 (a), Sb2T 3 (b), Bi2Se3 (c), Bi2Te3 (d) on the [111] surface. Red areas correspond to bulk energy bands and blue regions correspond to bulk energy gaps. The red lines near the Γ point show the surface states clearly for Sb2Te3 , Bi2Se3 , Bi2Te3 , but they do not exist for Sb2Se3 . Figure from Ref. [34]. 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. We discuss and analyze a possible solution to this puzzle in Chapter 7. Another peculiar effect that occurs in a strong topological insulator is the ’worm- hole’ effect. It is well known that 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 STI, 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 insula- tor. We establish this ‘wormhole’ effect by both general qualitative considerations and by numerical calculations within a minimal lattice model in Chapter 5 Topological insulators survive under non-magnetic disorder, as long as it is not too strong compared with the band gap. Once the disorder is of the order of 11 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 simula- tions [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 introduc- ing 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 exam- ples 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 =∑ i, j c†iMi jc j (2.1) 13 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† =∑ i, j c†jM ∗ i jci =∑ i, j c†iM ∗ jic j (2.2) 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 hop- ping 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 inver- sion. 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 opera- tors for each sublattice, and rewrite the real space representation using the sublat- tice 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 = 1√ N∑i ei~k·~Ric†i (2.3) c†i = 1√ N∑~k e−i~k·~Ric†k (2.4) 14 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 =∑ i, j c†iMi jc j = ∑ ~k,~k′ c†kck′ 1 N∑i, j Mi j ei ~k ′ ·~R j−i~k·~Ri (2.5) = ∑ ~k,~k′ c†kck′∑ ~δ Mi j ei ~k·~δ ( 1 N∑j ei(~k ′−~k)·~R j ) = ∑ ~k c†kck∑ ~δ Mi j ei ~k·~δ =∑ ~k c†kH(~k)ck where ~δ ≡ ~R j−~Ri, and in the last line we have defined the Fourier transform of M as H(~k)≡∑ ~δ Mi j ei ~k·~δ (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 ei ~k·~δ ei~K·(~R j−~Ri) =∑ ~δ Mi j ei ~k·~δ = H(~k) (2.7) since for any lattice vector ~R and reciprocal lattice vector ~K we have ei~K·~R = 1. Us- ing the definition (2.6) and the Hermitianity of M, we prove that H(~k) is Hermitian too H†(~k) =∑ ~δ M∗i j e −i~k·~δ =∑ ~δ M∗ji e i~k·~δ =∑ ~δ Mi j ei ~k·~δ = H(~k) (2.8) In (2.5) we used the fact that δ~k~k′ = 1 N∑i ei(~k−~k ′ )·~Ri (2.9) 15 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) c†k = 1√ N∑i ei~k·~Ric†i =∑ ~k′ c† k′ ( 1 N∑i ei(~k−~k ′ )·~Ri ) (2.10) On the other hand, we can write c†k =∑ ~k′ c† k′ δ~k~k′ (2.11) By comparing the last two equations we get (2.9). Similarly, by substituting (2.3) into (2.4) we get δi j = 1 N∑ ~k ei~k·(~R j−~Ri) (2.12) 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 Ψ†kH(~k)Ψk (2.13) 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 E ik and N wavefunctions ϕ i k. 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 = ∑ i, j,σ ,σ ′ c†iσMi jσσ ′ c jσ ′ (2.14) 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 = ∑ ~k,σ ,σ ′ c†kσH(~k)ckσ ′ (2.15) where H(~k)≡∑ ~δ Mi jσσ ′ e i~k·~δ (2.16) Once again, we write the Hamiltonian as H = ∑~kΨ † kH(~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 = ε∑ i c†i ci− t ∑ 〈i, j〉 c†i c j (2.17) 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 = 2pi Φ0 ∫ ~R j ~Ri ~A · ~dl (2.18) and H0 = ε∑ i c†i ci− t ∑ 〈i, j〉 eiθi jc†i c j (2.19) The hopping term is gauge invariant, as can be seen by making a gauge transfor- mation ci → eiχici (2.20) θi j → θi j+χi−χ j So that eiθi jc†i c j→ eiθi j+iχi−iχ jc†i e−iχic j eiχ j = eiθi jc†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 = ε∑ i,σ c†iσciσ − t ∑ 〈i, j〉,σ eiθi jc†iσc jσ (2.22) 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 correspond- ing to each value of sz. For example, for spin 12 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,σ ≡ ε∑ i c†iσciσ − t ∑ 〈i, j〉 eiθi jc†iσc jσ (2.23) Spin-orbit terms are common in tight-binding models of topological insula- tors. 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α ∑ 〈〈i, j〉〉 eiθik eiθk j ( c†i↑ c†i↓ )T ~σ · (d̂ik× d̂k j) ( c j↑ c j↓ ) (2.24) where 〈〈i, j〉〉 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 ma- trices 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 un- changed P≡ τx⊗ I (2.25) For example, we define the operator c†k,σ that creates an electron with momentum ~k and spin σ on sublattice A, and similarly d†k,σ creates an electron on sublattice B. Now define a vector Ψk ≡ (ck↑,ck↓,dk↑,dk↓), and check that P indeed swaps sublattices PΨk =  1 0 0 1 1 0 0 1   ck↑ ck↓ dk↑ dk↓ =  dk↑ dk↓ ck↑ ck↓  (2.26) Note that P2 = I, and therefore P= P−1. Now define a time-reversal operator Θ≡ i(I⊗σy)K (2.27) 19 where K is the complex conjugation operator, and the inverse time-reversal opera- tor Θ−1 ≡ iK(I⊗σy) (2.28) So that ΘΘ−1 = [i(I⊗σy)K] [iK(I⊗σy)] =−i2(I⊗σ2y ) = I⊗ I (2.29) 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 † iMi jc j, we define a Bloch Hamiltonian H(~k), as in (2.6) H(~k)≡∑ i, j Mi j ei ~k·(~Ri−~R j) (2.32) 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 =∑ i, j ΘMi j ei ~k·(~Ri−~R j)Θ−1 =∑ i, j Mi j e−i ~k·(~Ri−~R j) = H(−~k) (2.33) 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 (2.34) 20 with i = (nxnynz) 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 δi (2.35) where δi =∏ m ξ2m(~Γi) (2.36) The index i = n1n2n3 (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 6=i=0,1 δn1n2n3 (2.37) We use this method in Chapter 3 to topologically classify our model. 21 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 = λ ∑ µ=x,y,z i ( iΨ†i σµΨi+µ +h.c. ) (3.1) with Ψi ≡ (ci↑,ci↓). This can also be written as HSO = λ ∑ µ=x,y,z i,σ ,σ ′ (ic†i,σσµ,σσ ′ci+µ,σ +h.c.) (3.2) where λ is the strength of the spin-orbit term, σµ are the Pauli matrices, c†iσ cre- 22 ates an electron on site i with spin σ , and site i+ µ is taken to mean the nearest neighbor in direction µ̂ . Note that the hopping amplitude is iλ for bonds pointing in directions x̂, ŷ, ẑ, and−iλ for−x̂,−ŷ,−ẑ. Therefore, the sum of all phase factors (±i) on each face is zero. We rewrite the Hamiltonian in momentum space HSO = iλ∑ ~k ( c†k↑ c†k↓ )T [ σx(eikxa− e−ikxa)+σy(eikya− e−ikya)+σz(eikza− e−ikza) ](ck↑ ck↓ ) = −2λ∑ ~k ( c†k↑ c†k↓ )T [σx sin(kxa)+σy sin(kya)+σz sin(kza)] ( ck↑ ck↓ ) (3.3) = −2λ∑ ~k ( c†k↑ c†k↓ )T [ ∑ µ σµ sin(kµa) ]( 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Ψ†kHSO(~k)Ψk, and HSO(~k) =−2λ∑ µ σµ sin(kµa) =−2λ ( sin(kza) sin(kxa)− isin(kya) sin(kxa)+ isin(kya) −sin(kza) ) (3.4) This is a sum of anti-commuting matrices, so that the energies can be written im- mediately εk =±2λ √ sin2(kxa)+ sin2(kya)+ sin2(kza) (3.5) The two energy bands are degenerate at Dirac points located at ~K = pi a (nx,ny,nz) (3.6) where nµ ∈ Z. The reciprocal lattice primitive vectors are ~bµ = (2pi/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 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 each site. If we define creation operators c†iσ and d † iσ on the two orbitals, and Ψ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 = ε∑ i,σ (c†iσdiσ +h.c.)− t ∑ 〈i, j〉,σ (c†iσd jσ +h.c.) (3.8) In momentum space Ht =∑ k,σ {c†kσdkσ (ε−2t[cos(kxa)+ cos(kya)+ cos(kza)])+h.c.} (3.9) so that Ht(~k) = mk(τx⊗ I) (3.10) 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 =± √ ε2k +m 2 k (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 E~K0+~q =± √ (2λa)2(q2x+q2y+q2z )+(ε−6t)2 (3.14) We see that near ~K0 the spectrum is gapped with mass m~K0 = ε − 6t. Near one of the Dirac points ~K, the spectrum will have a mass m~K = ε−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(kxa)+ cos(kya)+ cos(kza)] d2 = 0 d3 = −2λ sin(kxa) d4 = −2λ sin(kya) d5 = −2λ sin(kza) In order to classify this insulator, we must find δi = −sgn(d1(~k =~Γi)) for each of the 8 TRIM, defined by ~Γi = 12(nx~bx+ny~by+nz~bz) with i = (nxnynz) and nµ ∈ {0,1}. For a cubic lattice of spacing a, the reciprocal lattice vectors are ~bµ = (2pi/a)µ̂ , so that we find ~Γi=(nxnynz) = pi a (nx,ny,nz) (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µpi) (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 nx ny nz ∑µ(−1)nµ d1(~Γi) δi δ<i δi δ>i δi 1 0 0 0 3 ε−6t + + + + − 2 0 0 1 1 ε−2t + + + − − 3 0 1 0 1 ε−2t + + + − − 4 1 0 0 1 ε−2t + + + − − 5 0 1 1 -1 ε+2t + + − − − 6 1 0 1 -1 ε+2t + + − − − 7 1 1 0 -1 ε+2t + + − − − 8 1 1 1 -3 ε+6t + − − − − Table 3.1: Time reversal invariant momenta and the parity eigenvalues for our model (−1)ν1 = ∏ nx=1 ny,nz=0,1 δi = δ4δ6δ7δ8 (3.19) (−1)ν2 = ∏ ny=1 nx,nz=0,1 δi = δ3δ5δ7δ8 (−1)ν3 = ∏ nz=1 nx,ny=0,1 δi = δ2δ5δ6δ8 (−1)ν0 = 8 ∏ i=1 δi So that (ν0;ν1ν2ν3) =  (0;000) ε <−6t orε > 6t (1;111) −6t < ε <−2t (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 !1;000"!0;111"!0;000"!1;111"!0;111" Ε"2tΕ"#2t Ε"6tΕ"#6t 2 t Ε 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 la- bels 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 λ 6= 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 topo- logical 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 semi- 28 metal on this plane. This is the domain wall. We define m(z) ≡ ε(z)− 6t and choose a solitonic profile with asymptotic values lim z→∞m(z) = +2t (3.21) lim z→−∞m(z) = −2t and m(0) = 0. By looking at (3.11), we see that in the orbital subspace, only two Pauli matri- ces 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 −pi/4 around τx, which will therefore leave τx unchanged. If we define U ≡ e−i(pi/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 Ψk→U†Ψk = 1√ 2  1 0 i 0 0 1 0 i i 0 1 0 0 i 0 1   ck↑ ck↓ dk↑ dk↓ = 1√2  ck↑+ idk↑ ck↓+ idk↓ ick↑+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 hap- pen 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 representa- 29 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, andHq 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 0 m̃(z)−∂z )( a(z) b(z) ) = 0 (3.29) so that (m̃(z)+∂z)a(z) = 0 (3.30) (m̃(z)−∂z)b(z) = 0 30 By rearranging and integrating both sides, we find a(z) = e− ∫ z 0 m̃(z ′)dz′ (3.31) b(z) = e+ ∫ z 0 m̃(z ′)dz′ However, when the argument of the exponent is positive, the resulting wave func- tion 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′ (3.32) Similarly for the lower block of hz, ϕ−(z) = ( 0 1 ) e− ∫ z 0 m̃(z ′)dz′ (3.33) So that the two wavefunctions of the full matrix hz are ϕ+(z) = N  1 0 0 0  e− ∫ z 0 m̃(z ′)dz′ (3.34) ϕ−(z) = N  0 0 0 1  e− ∫ z 0 m̃(z ′)dz′ whereN 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 haveHq= 0, so thatH = 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 ofH . 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 exponen- tial (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 onlyHq. P≡  1 0 0 0 0 0 0 1  (3.35) So that the projected Hamiltonian is h± = P†HqP=−vF ( 1 0 0 0 0 0 0 1 ) −iqx−qy −iqx+qy iqx+qy iqx−qy   1 0 0 0 0 0 0 1  = vF ( 0 qy+ iqx qy− iqx 0 ) = vF(qyσx−qxσy) (3.36) We project the operator basis vector too Ψ̃q = P†Ψq = ( 1 0 0 0 0 0 0 1 ) 1√ 2  cq↑+ idq↑ cq↓+ idq↓ icq↑+dq↑ icq↓+dq↓ = 1√2 ( cq↑+ idq↑ icq↓+dq↓ ) (3.37) So that the full projected Hamiltonian, which we interpret as the effective Hamil- tonian for the edge states, can be written as HB = ∑ ~q Ψ̃†qh±Ψ̃q (3.38) = vF 2 ∑ ~q ( c†q↑− id†q↑ −ic†q↓+d†q↓ )T ( 0 qy+ iqx qy− iqx 0 )( cq↑+ idq↑ icq↓+dq↓ ) 32 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 ε± =±vF √ q2x+q2y (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 e−iα ) (3.41) φ− = 1√ 2 ( −eiα 1 ) We construct the unitary matrix U that has φ± as its columns U = 1√ 2 ( 1 −eiα e−iα 1 ) (3.42) which is guaranteed to diagonalize h± U†h±U = ε 2 ( 1 eiα −e−iα 1 )( 0 eiα e−iα 0 )( 1 −eiα e−iα 1 ) = ε ( 1 0 0 −1 ) (3.43) and the operator basis vectors γq =U†Ψ̃q = 1 2 ( 1 eiα −e−iα 1 )( cq↑+ idq↑ icq↓+dq↓ ) (3.44) = 1 2  ei α2 [e−i α2 (cq↑+ idq↑)+ ei α2 (icq↓+dq↓)] e−i α 2 [ −e−i α2 (cq↑+ idq↑)+ ei α2 (icq↓+dq↓) ] 33 If we define γk ≡ ( γk+ γk− ) then HB = ∑ q,i=± εiγ†qiγqi (3.45) 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+} = 1 4 [ {cq↑,c†q↑}+{cq↓,c†q↓}+{dq↑,d†q↑}+{dq↓,d†q↓} ] = 1 {γq−,γ†q−} = 1 4 [ {cq↑,c†q↑}+{cq↓,c†q↓}+{dq↑,d†q↑}+{dq↓,d†q↓} ] = 1 (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) = 1√ 2 [ ϕ+(z)+ e−iαϕ−(z) ] =N  1 0 0 e−iα  e− ∫ z 0 m̃(z ′)dz′ (3.47) Ψ−(~q,z) = 1√ 2 [−eiαϕ+(z)+ϕ−(z)]=N  −eiα 0 0 1  e− ∫ z 0 m̃(z ′)dz′ (where the normalization constant has been redefinedN →√2N ) To find the energies, we note that the Hamiltonian is a sum of anticommuting matricesH = hz+Hq. Therefore, it is convenient to square it, soH 2 = h2z+H 2 q . 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 √ q2x+q2y (3.48) 34 Note that these bound states are written in the block off-diagonal basis given by (3.23) Ψk = 1√ 2  ck↑+ idk↑ ck↓+ idk↓ 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′ = e− 2t vF |z| (3.51) 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 − 2tvF ∫ z 0 tanh ( z′ ξ ) dz′ = e − 2tvF ξ lncosh(z ′) ∣∣∣z/ξ 0 = cosh ( z ξ )− 2tvF ξ (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∓ 2t vF z (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 en- ergy, 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) = m~K(τ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 = pi a (nx,ny,nz) (3.58) 36 We define the momentum space non-abelian gauge field, in this case a set of three matrices A j of dimensions 4×4 A µνj =−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 4pi ∫ BZ d3kε i jkTr [ Ai∂ jAk+ 2 3 iAiA jAk ] (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 matricesA j and substitute into this equa- tion. For each Dirac point ~K, we find θK = ∫ Λ 0 d3q 2mKv3F (m2K +q2v2F)2 = tan−1 ( vFΛ mK ) − mKvFΛ m2K + 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 = pi 2 sign(mK) (3.62) and θ =∑ K θK = pi 2∑K sign(mK) (3.63) where the result is understood to be mod 2pi , so that 0 ≤ θ < 2pi . 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 = pi 2∑K sign [ ε−2t∑ µ (−1)nµ ] (3.64) We now vary ε and find that θ = 0 |ε|< 2t or |ε|> 6tpi 2t < |ε|< 6t (3.65) which verifies our earlier findings (3.20) if we identify ν0 = θ/pi . 38 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/2pih)~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 fundamen- tal 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 θ 6= 0 is predicted to bind a (generally fractional) electric charge −e(θ/2pi + 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 θ = pi emerges naturally in the descrip- tion 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 in- sulator 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 electro- magnetic response of a STI is characterized by the axion term ∆Laxion = θ ( e2 2pih ) ~B ·~E (4.1) with θ = pi , 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 mag- netism, θ can acquire an arbitrary value. Fluctuations in the magnetic order param- eter 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 mag- netic flux Φ0 = hc/e) immersed in an axion medium must carry electric charge 40 −e(θ/2pi+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 sec- tion 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 ob- stacles. 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 2pih ) ~B ·~E (4.2) We recall Maxwell’s equations with magnetic monopoles (and currents) and θ = 0 (see page 273 in Jackson [69], for example) ~∇ ·~E = ρe (4.3) ~∇ ·~B= ρm ~∇×~E =−∂ ~B ∂ t − ~Jm ~∇×~B= ∂ ~E ∂ t + ~Je 41 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− e4pi2 ~∇θ ·~B (4.4) ~∇ ·~B= 0 ~∇×~E =−∂ ~B ∂ t ~∇×~B= ∂ ~E ∂ t + ~Je+ e 4pi2 ( ∂θ ∂ t ~B+~∇θ ×~E ) (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 κ =α/4pi2 = e/4pi2). We infer that Maxwell’s equations in the presence of a non-zero θ and magnetic charges and currents should be ~∇ ·~E = ρe− e4pi2 ~∇θ ·~B (4.5) ~∇ ·~B= ρm ~∇×~E =−∂ ~B ∂ t − ~Jm ~∇×~B= ∂ ~E ∂ t + ~Je+ e 4pi2 ( ∂θ ∂ t ~B+~∇θ ×~E ) 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 θ → θ +2pin of the axion actionSaxion, implying that θ can be chosen from the interval [0,2pi). 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 transforma- tion 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= 4pigδ (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 ∂ t + e 4pi2 ∂θ ∂ t ~∇ ·~B= 0 (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∆θ pi (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 ∆θ = pi , thus one expects a unit magnetic monopole (g= 1/2) to bind fractional charge Q=−e ( 1 2 +n ) (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− e4pi2 ~∇θ ·~B (4.9) ~∇ ·~B= ρm ~∇×~E =−∂ ~B ∂ t − ~Jm ~∇× ~H = ∂ ~D ∂ t + ~Je+ e 4pi2 ( ∂θ ∂ t ~B+~∇θ ×~E ) As before, we take the divergence of the last equation, finding ~∇ · ∂ ~D ∂ t + e 4pi2 ∂θ ∂ t ~∇ ·~B= 0 (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/20piε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 = 4pigδ (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) = pi[Θ(R− r)−Θ(ε− r)] (4.13) where Θ is the Heaviside function. We note that ~∇θ(r) = dθ dr r̂ = pi [ dΘ(R− r) r − dΘ(ε− r) r ] (4.14) = −pi[δ (r−R)−δ (r− ε)]r̂ 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 ·d~s= ∫ ρm d3r (4.15) to find the magnetic field ~B(~r) = g r2 r̂ (4.16) So that the first equation of (4.5) gives the effective charge density ~∇ ·~E =− e 4pi2 ~∇θ ·~B= ge [ δ (r−R) 4piR2 − δ (r− ε) 4piε2 ] ≡ ρe f (r) (4.17) 45 By integrating over a narrow region around each sphere we find the charge on the outer and inner surfaces Qo = ∫ R+δ R−δ ρe f (r)4pir2dr = ge (4.18) Qi = ∫ ε+δ ε−δ ρe f (r)4pir2dr =−ge (4.19) Therefore, for g = 12 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 = 4pigδ (3)(~r) (4.20) This equation is fulfilled by a Coulombic magnetic field ~Bg = g ~r r3 = g r2 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 6= 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 = gr2 r̂−4pigθ(z)δ (x)δ (y)ẑ (4.23) The singular part ~Bs can be thought of as an infinitely long infinitesimal solenoid leading flux lines in the −ẑ direction. It is often referred to as a “Dirac string”. 46 Equation (4.22) is now fulfilled since ~∇ ·~Bs =−4pigδ (x)δ (y)dθ(z)dz =−4pigδ (3)(~r) (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 Φ= ∫ ~B ·d~s = g [∫ r̂ r2 ·d~s−4pi ∫ θ(z)δ (x)δ (y)ẑ ·d~s ] (4.25) = g [ 1 R2 ∫ ds−4pi ∫ θ(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 r2 r̂ =−g~∇ ( 1 r ) (4.26) since in general, for a function of the radial coordinate only, we have ~∇ f (r) = d f dr r̂ (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 (4.29) 47 On the other hand ~∇×~B = ~∇× (~∇×~A) (4.30) = ~∇(~∇ ·~A)−∇2~A=−∇2~A 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 6= 0) ~Ak = i ~k k2 ×~Bsk (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 situ- ation 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 = (2pi/Φ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 = eh̄/2mec is the Bohr magneton and ~S denotes the electron spin. For free electrons g is close to 2 but in solids the effec- tive g can be substantially larger. The Zeeman coupling thus leads to an additional 48 term in the Hamiltonian, HZ =−gµB 12∑j ~B j · ( Ψ†j~σΨ 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/4pir2)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 14(4× 203) = 8,000. In order to calculate the charge density at half filling we require knowledge of all the occupied eigenstates. ï10ï8 ï6 ï4 ï2 0 2 4 6 8 10ï10ï8 ï6ï4ï2 024 6810 x y ln (δ ρ + c) 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 lat- tice with 203 sites with a unit monopole at its center, with parameters t = λ , ε = 4t, leading to a bulk gap ∆= 4t. 49 0 2 4 6 8 10 12 14 160.0 0.1 0.2 0.3 0.4 0.5 r δQ (r )   g=0 g=2 g=10 Figure 4.2: The excess charge δQ(r) (in units of e) for different Zeeman cou- pling 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 topologi- cal 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 conduc- tion bands. Thus, the valence band shows a net deficit of half a state in the vicinity 50 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90ï3.4 ï3.2 ï3.0 ï2.8 ï2.6 ï2.4 ï2.2 ï2.0 ï1.8 ï1.6 log r lo g |δQ (r ) − δQ 0 (r )|   g=2 g=6 g=10 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 pi 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,pi 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 struc- ture with a roughly exponential profile that is controlled by the properties of the microscopic Hamiltonian and is thus non-universal. At the intermediate length- scales 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 in- serted 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 cou- pling 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 sug- gested by exploiting emergent instead of fundamental monopoles. A classic exam- ple of such an emergent behavior in a crystalline solid is the 2007 theoretical pre- diction [81] and the subsequent experimental observation [82–84] of monopoles in frustrated magnetic systems called ‘spin ice’, realized in certain magnetic py- rochlore compounds such as Dy2Ti2O7 or Ho2Ti2O7. 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 cre- ates 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 conden- sate, 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 Φ=Φ0eiχ can contain vortices – point-like topological defects with the phase χ winding by ±2pi 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. ï7 ï5 ï3 ï1 1 3 5 7ï7ï5 ï3ï1 13 57 ï0.01 0.00 0.01 0.02 0.03 0.04 x y 2 − ρ 1 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 fash- ion. 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 cal- culation 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 ï5 ï3 ï1 1 3 5 7ï7 ï5 ï3 ï1 1 3 5 7 x y   ï0.01 0.00 0.01 0.02 0.03 0.04 (a) ï7 ï5 ï3 ï1 1 3 5 7ï7 ï5 ï3 ï1 1 3 5 7 x y   ï0.06 ï0.04 ï0.02 0.00 0.02 0.04 0.06 (b) ï7 ï5 ï3 ï1 1 3 5 7ï7 ï5 ï3 ï1 1 3 5 7 x y   ï0.01 0.00 0.01 0.02 0.03 0.04 (c) ï7 ï5 ï3 ï1 1 3 5 7ï7 ï5 ï3 ï1 1 3 5 7 x y   ï0.06 ï0.04 ï0.02 0.00 0.02 0.04 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 2 4 6 8 10 120.0 0.1 0.2 0.3 0.4 0.5 r δQ (r )   µ=0.00 µ=0.02 µ=0.05 µ=0.20 4 5 60.494 0.496 0.498 0.500 Figure 4.6: The excess charge δQ(r) (in units of e) for different values of dis- order strength µ . The inset shows a close up of the saturation. At this scale a small deviation from the expected asymptotic value 1/2 that in- creases 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 spec- tral gap. model it here by adding a term HD =∑ j µ+j Ψ † jΨ j+∑ j µ−j Ψ † jτzΨ j, (4.34) 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/2pir)δ (z)r̂ (in cylindrical coordinates) and the 56 vector potential can be chosen as ~A= (Φ0/2pi)ϕδ (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 in- plane 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 polyacety- lene [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 frac- tional 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. Al- though at the microscopic level θ is no longer quantized in the presence of time reversal symmetry (TRS)-breaking perturbations, we expect its effective value, rel- evant to the physics at long lengthscales, to remain pinned at pi 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. Quali- tatively, this suggests that the Witten effect may survive inclusion of a moderate concentration of magnetic dopants with randomly oriented moments. At low tem- peratures moments may order ferromagnetically [89]. In this case both TRS and the inversion symmetry are broken (the latter due to the random position of mag- netic 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 be- lieve 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 exper- imental challenge is reduced to designing a suitable method for the detection of fractional charge bound to the monopole. The fractional charge of elementary ex- citations 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 pro- vide 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 insula- tor. 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 sur- face described by Eq. (5.1). This argument was devised to establish the fractional charge of quasiparticles in fractional quantum Hall liquid (FQHL) [93] and in- volves 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 Fara- day’s law ~∇×~E =−(1/c)(∂~B/∂ t). This induces a Hall current on the STI surface ~j = σxy(~E× ẑ) which brings electric charge δQ= σxy Φ0 c = ( n+ 1 2 ) e (5.2) to the solenoid. Since the flux tube carrying a full flux quantum Φ0 can be re- moved 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 ex- citations 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 quasi- particles. 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 EFM STI \\ a) b) 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 in- dicated. b) Flux tube threading a cylindrical hole in a STI. Arrows illustrate the helical spin state for upward moving electrons (for down- movers 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 η = 12 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 H = 1 2 v [ h̄~∇ · n̂+ n̂ · (~p×~σ)+(~p×~σ) · n̂] (5.3) where v is the Dirac velocity, ~p = −ih̄~∇ is the momentum operator and ~σ = (σ1,σ2,σ3) is the vector of Pauli spin matrices. The magnetic flux is included by replacing ~p with ~pi = ~p− (e/c)~A, where ~A= ηΦ0(ẑ×~r)/2pir2 is the vector po- tential. For a cylindrical inner surface n̂=−(cosϕ,sinϕ,0), the Hamiltonian (5.3) becomes, in cylindrical coordinates and taking v= h̄= 1, Hk =− 12R + ( 1 R(i∂ϕ +η) −ike−iϕ ikeiϕ − 1R(i∂ϕ +η) ) . (5.4) We assumed a plane-wave solution eikz along the cylinder axis and replaced−i∂z→ k. The eigenstates ofHk are of the form Ψk(ϕ) = ( fk eiϕgk ) eiϕl (5.5) with l integer. The spinor Ψ̃k = ( fk,gk)T is an eigenstate of H̃kl = σ2k−σ3(l+ 1 2 −η)/R with an energy eigenvalue Ekl =±vh̄ √ 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 ∆= 2vh̄ R ∣∣∣∣12 −η ∣∣∣∣ . (5.7) 62 When η = 12 , i.e. at half flux quantum, the l = 0 mode becomes gapless, E η=1/2 k0 = ±vh̄|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 pi 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]. Be- low, 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 inter- section of the flux tube with the magnetized surface as a function of η = Φ/Φ0. For small η we observe δQ = 12eη , consistent with the fractional Hall conductiv- ity σxy = e2/2h expected on the basis of Eq. (5.1). At η = 12 a charge e/2 travels along the flux tube and combines with the negative charge that has built up on the opposite surface. For η > 12 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 indi- cated 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 η = 12 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. Sweep- ing 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 ob- servable 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/piB)1/2 for N oscillations. Second, R must be sufficiently small so that the maximum spectral gap Eq. (5.7) is large compared to kBT , or otherwise the oscillations in the conductance will be washed out by thermal broadening. This gives R. h̄v/ √ 2kBT . Taking typical val- ues B = 10T, N = 10, v = 5× 105m/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 mag- netic 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 Bm∗∆/eh̄' 6T , by substituting the bulk gap ∆' 0.3eV and the effective mass m∗ ' 0.14me, for Bi2Se3 [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 suffi- ciently 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 phe- nomenology 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, con- ductance 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. Sb2Se3, 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 - non- inverted 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 impuri- ties. The self energy is given by (see derivation in Appendix A.4) Σ(~k) = Nu2 ( a 2pi )2 ∫ d2k [ε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 u2 = 1 NU0 ∫ U0 2 −U02 u2du= 2 NU0 [ u3 3 ]U0 2 0 = U20 12N (6.3) and the infinitesimal imaginary part is required in order to make the integral con- verge. 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. [εF−H0(~k)−Σ(~k)]−1 = ( ε̃F−∑ i α̃iσi )−1( ε̃F +∑ i α̃iσi )−1( ε̃F +∑ i α̃iσi ) = ε̃F +∑i α̃iσi ε̃2F −∑i α̃2i (6.5) We substitute this expression into (6.2), to find Σ= U20 12 ( a 2pi )2 ∫ BZ d2k ε̃F +∑i α̃iσi ε̃2F −∑i α̃2i (6.6) which we can separate into four equations Σ0 = A ∫ BZ d2k εF−Σ0− γk2 D(~k) (6.7) Σ1 = A ∫ BZ d2k αkx+Σ1 D(~k) Σ2 = A ∫ BZ d2k −αky+Σ2 D(~k) Σ3 = A ∫ BZ d2k m+Σ3+βk2 D(~k) where we have defined the prefactor A ≡ U2012 ( a 2pi )2 and the denominator D(~k) ≡ ε̃2F −∑i α̃2i . 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−ARe [∫ BZ d2k ε̃F− γk2 D(~k) ] (6.8) m̄ = m+ARe [∫ BZ d2k m̄+βk2 D(~k) ] (6.9) 71 We write out the denominator D(~k)≡ ε̃2F −∑ i α̃2i = (ε̄F− γk2)2− (αkx+Σ1)2− (αky−Σ2)2− (m̄+βk2)2 → (εF− γk2)2−α2k2− (m+βk2)2 (6.10) = (ε2F −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, d2k = 2pikdk, and impose a cutoff Λ, to find ε̄F = εF− 2piAγγ2−β 2 Re ∫ Λ 0 dk εF γ k− k3 ε2F−m2 γ2−β 2 − α2+2(εFγ+mβ ) γ2+β 2 k 2+ k4  (6.11) The integral we need to solve is of the form ∫ Λ 0 dk Ak− k3 B+Ck2+ k4 '−1 4 ln ∣∣∣∣Λ4B ∣∣∣∣ (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 Λ= pi/a, to find ε̄F = εF−2pi [ U20 12 ( a 2pi )2] γ β 2− γ2 1 4 ln ∣∣∣∣ γ2−β 2ε2F −m2 (pi a )4∣∣∣∣ (6.13) = εF− 12 U20 a 2 48pi γ β 2− γ2 ln ∣∣∣∣β 2− γ2ε2F −m2 (pi a )4∣∣∣∣ and similarly m̄ = m− 1 2 U20 a 2 48pi β β 2− γ2 ln ∣∣∣∣β 2− γ2ε2F −m2 (pi a )4∣∣∣∣ (6.14) 72 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 = U20 a 2 48pi β β 2− γ2 ln ∣∣∣∣β 2− γ2ε2F −m2 (pi a )4∣∣∣∣ (6.15) which can be solved to give |ε2F −m2|= (pi a )4 |β 2− γ2|e− 48pim a2 β2−γ2 β 1 U20 (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−U 2 0 a 2 48pi 1 γ±β ln ∣∣∣∣β 2− γ2ε2F −m2 (pi 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 so that m̄= 0 which is the condition for the critical line. space and time. The quadratic term !p2 ¼ "@2!r2 in H, acting on the decaying state / e"x=", adds a negative correction #m to the topological mass. The renormalized mass !m ¼ mþ #m can therefore have the opposite sign as the bare mass m. Topological mass renormalization by disorder, and the resulting change in the phase diagram, has previously been studied without the terms quadratic in momentum [11]. The sign of !m and m then remains the same and the TAI phase cannot appear. We extract the renormalized topological mass !m, as well as the renormalized chemical potential !$, from the self- energy " of the disorder-averaged effective medium. To make contact with the computer simulations [2,3], we discretize H on a square lattice (lattice constant a) and take a random on-site disorder potential U, uniformly distributed in the interval ("U0=2, U0=2). We denote by H0ðkÞ the lattice Hamiltonian of the clean quantum well in momentum representation [5,7]. The self-energy, defined by ðEF "H0 " "Þ"1 ¼ hðEF "HÞ"1i; (2) with h& & &i the disorder average, is a 2' 2matrix which we decompose into Pauli matrices: " ¼ "0%0 þ "x%x þ "y%y þ "z%z. The renormalized topological mass and chemical potential are then given by !m ¼ mþ lim k!0 Re"z; !$ ¼ EF " limk!0 Re"0: (3) The phase boundary of the topological insulator is at !m ¼ 0, while the Fermi level enters the (negative) band gap when j !$j ¼ " !m. In the self-consistent Born approximation, " is given by the integral equation [12] " ¼ 112U20ða=2&Þ2 Z BZ dk½EF þ i0þ "H0ðkÞ " ")"1: (4) (The integral is over the first Brillouin zone.) The self- energy is independent of momentum and diagonal (so there is no renormalization of the parameters ', !, (). By calculating !m and !$ as a function of EF and U0 we obtain the two curves A and B in Fig. 1. We have also derived an approximate solution in closed form [13], !m ¼ m" U 2 0a 2 48&@2 !!2 " (2 ln!!!!!!!!!2 " (2E2F "m2 " &@ a # 4 !!!!!!!!; (5a) !$ ¼ EF " U 2 0a 2 48&@2 (!2 " (2 ln!!!!!!!!!2 " (2E2F "m2 " &@ a # 4 !!!!!!!!; (5b) showing that the correction #m ¼ !m"m to the topologi- cal mass by disorder is negative—provided !> (. For !< ( the clean HgTe quantum well would be a semi- metal, lacking a gap in the entire Brillouin zone. Neither the TAI phase nor the QSH phase would then appear. In HgTe the parameter ! is indeed larger than (, but not by much [10]. Equation (5) implies that the lower branch of curve B (defined by !$ ¼ !m< 0) is then fixed at EF * m> 0 independent of U0. This explains the puzzling ab- sence of the TAI phase in the valence band (EF < 0), observed in the computer simulations [2,3]. To quantitatively test the phase diagram resulting from the effective medium theory, we performed computer simulations similar to those reported in Refs. [2,3]. The conductance G is calculated from the lattice Hamiltonian [5,7] in a strip geometry, using the method of recursive Green functions. The strip consists of a rectangular disor- dered region (width W, length L), connected to semi- infinite, heavily doped, clean leads [14]. Theory and simu- lation are compared in Figs. 1 and 2. Figure 1 shows the phase diagram. The weak-disorder boundary of the TAI phase observed in the simulations is described quite well by the self-consistent Born approxi- mation (curve B)—without any adjustable parameter. Curve B limits the region where (A) the renormalized topological mass !m is negative and (B) the renormalized chemical potential !$ lies inside the band gap: j !$j<" !m. Condition (A) is needed for the existence of edge states with a quantized conductance. Condition (B) is not needed for an infinite system, because then Anderson localization suppresses conductance via bulk states as well as coupling of edge states at opposite edges. In the relatively small FIG. 1 (color online). Computer simulation of a HgTe quan- tum well [10], showing the average conductance hGi as a function of disorder strength U0 (logarithmic scale) and Fermi energy EF, in a disordered strip of width W ¼ 100a and length L ¼ 400a. The TAI phase is indicated. Curves A and B are the phase boundaries resulting from the effective medium theory. Curve A separates regions with positive and negative renormal- ized topological mass !m, while curve Bmarks the crossing of the renormalized chemical potential !$ with the band edge (j !$j ¼ " !m). Both curves have been calculated without any adjustable parameter. The phase boundary of the TAI at strong disorder is outside of the regime of validity of the effective medium theory. PRL 103, 196805 (2009) P HY S I CA L R EV I EW LE T T E R S week ending 6 NOVEMBER 2009 196805-2 Figure 6.1: Average conductance 〈G〉 as a function of disorder strength U0 and the Fermi energy EF in a 2D strip of dimensions 400a× 100a. Curves A and B are the phase boundaries from the self consistent Born approximation. Figure from Ref. [3]. 6.4 TAI in 3D 6.4.1 General treatment After considering the 2D case, it is natural to ask whether the TAI phase can be found in 3D also. We consider a generalization of the Hamiltonian in the previous su section to three dimensions. The Hamiltonian we will use was presented in Section 3.2 and is the same one we used previously to study the Witten effect and the Wormhole effects, with an additional constant term to shift the chemical 74 potential, and a hopping term Hγ = 6γ∑ i,σ c†iσciσ − γ ∑ 〈i, j〉,σ c†iσc jσ (6.18) 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Ψ † kH0(~k)Ψk, with H0(~k) = −2γ∑ i cos(kia)−2λ∑ i sin(kia)τ3⊗σi+mkτ1⊗σ0 (6.19) = 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(kia), di ≡−2λ sin(kia), d4 = 2γ[3−∑i cos(kia)]. 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 = ( ∑ µ dµγµ )2 = ∑ µ>ν dµdν{γµ ,γν}=∑ µ d2µ (6.20) So that the energies are E0(~k) = d4± √ ∑ µ d2µ = d4± √ m2k+ ε 2 k (6.21) where εk =±2λ √ sin2(kxa)+ sin2(kya)+ sin2(kza) (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 = (pi,pi,pi), 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 2pi )3 ∫ d3k [ε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 u2 = 1 U0 ∫ U0 2 −U02 u2du= 2 U0 [ u3 3 ]U0 2 0 = U20 12 (6.25) and the infinitesimal imaginary part is required in order to make the integral con- verge. 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(kia)], 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. [εF−H0(~k)−Σ(~k)]−1 = ( ε̃F−∑ µ d̃µγµ )−1( ε̃F +∑ µ d̃µγµ )−1( ε̃F +∑ µ d̃µγµ ) = ε̃F +∑µ d̃µγµ ε̃2F −∑µ d̃2µ (6.27) We substitute this expression into (6.24), to find Σ= U20 12 ( a 2pi )3 ∫ BZ d3k ε̃F +∑µ d̃µγµ ε̃2F −∑µ d̃2µ (6.28) which we can separate into five equations Σ4 = A ∫ BZ d3k εF−Σ4−d4 D(~k) (6.29) Σ0 = A ∫ BZ d3k mk+Σ0 D(~k) Σi = A ∫ BZ d3k di+Σi D(~k) where we have defined the prefactor A ≡ U2012 ( a 2pi )3 and the denominator D(~k) ≡ ε̃2F −∑µ d̃2µ . 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−ARe [∫ BZ d3k ε̄F− γck D(~k) ] (6.30) and m̄ = m+ReΣ0 = m+ARe [∫ BZ d3k m̄+ tck D(~k) ] (6.31) 77 where we have defined ck ≡ 2 [ 3−∑ i cos(kia) ] (6.32) and we have assumed that the smallest gap is at the Γ point, so that m= ε−6t and then we rewrite mk = ε−2t∑ i cos(kia) = m+ tck (6.33) We consider first the three equations for Σi. In the zero order Born approximation we have Σi = A ∫ BZ d3k di D(~k)|Σi=0 (6.34) Write out the denominator D(~k)≡ ε̃2F −∑ µ d̃2µ = (εF−d4−Σ4)2−∑ µ (dµ +Σµ)2 (6.35) = (εF−d4−Σ4)2− (mk+Σ0)2−∑ i (−2λ sin(kia)+Σi)2 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(kia) 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 D(~k) = [ ε̄F−2γ ( 3−∑ i cos(kia) )]2 − [ m̄+2t ( 3−∑ i cos(kia) )]2 − (2λ )2∑ i sin2(kia) = (ε̄F− γck)2− (m̄+ tck)2+λ 2sk (6.36) 78 where we have defined sk ≡ −4∑ i sin2(kia) (6.37) 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| (6.38) Σ4 = U20 12N3 ∑k∈BZ εF + i0+− γck−Σ4 D(~k)| εF→εF+i0+ Σ0 = U20 12N3 ∑k∈BZ m+ tck+Σ0 D(~k)| εF→εF+i0+ For the second and third equations, we notice that we are summing over expres- sions 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 calcu- late 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 = 2pi niN , ni = 1, . . . ,N and i= x,y,z and N = 50 (for exam- ple) 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 = (U20 /12N3)I0, Σ̃4 = (U20 /12N3)I4, ε̃F = Re Σ̃4±|m+Re Σ̃0| (e) If ∣∣∣ ε̃F−εFεF ∣∣∣ < tol, ∣∣∣Re Σ̃0−ReΣ0ReΣ0 ∣∣∣ < tol and ∣∣∣ Im Σ̃0−ImΣ0ImΣ0 ∣∣∣ < tol (and the same for Σ4), we have reached self consistency - save these values (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 val- ues 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 parame- ters. 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 (6.39) U20 12N3 ∑k∈BZ εF + i0+− γck−Σ4 D(~k)| εF→εF+i0+ −Σ4 = 0 U20 12N3 ∑k∈BZ m+ tck+Σ0 D(~k)| εF→εF+i0+ −Σ4 = 0 80 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+λ 2sk (6.41) = [ε̄F + m̄+(t− γ)ck] [ε̄F− m̄− (t+ γ)ck]+λ 2sk This expression is more convenient to work with when looking for the phase bound- aries, 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 ∫ BZ d3k ε̄F− γck D(~k) (6.42) m̄ = m+A ∫ BZ d3k m̄+ tck D(~k) and we must remember to take only the real part of the integrals. Add and subtract m̄+ ε̄F = m+ εF +A ∫ BZ d3k (m̄− ε̄F)+(t+ γ)ck D(~k) (6.43) m̄− ε̄F = m− εF +A ∫ BZ d3k (m̄+ ε̄F)+(t− γ)ck 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 ∫ BZ d3k 2m̄+(t+ γ)ck D+(~k) (6.44) 2m̄ = m− εF +A ∫ BZ d3k (t− γ)ck D+(~k) 81 and D+(~k) =−(t− γ)ck [2m̄+(t+ γ)ck]+λ 2sk (6.45) Now define the integrals I+1 ≡ ∫ BZ d3k 1 D+(~k) (6.46) I+2 ≡ ∫ BZ d3k ck D+(~k) Rewrite the two equations in terms of these integrals 0 = m+ εF +A [ 2m̄I+1 +(t+ γ)I + 2 ] (6.47) 2m̄ = m− εF +A(t− γ)I+2 Solve for A and εF (in terms of A) to find A = m̄−m m̄I+1 + tI + 2 (6.48) εF = −m−A [ 2m̄I+1 +(t+ γ)I + 2 ] Similarly, for the second case, ε̄F− m̄= 0 we have 2m̄ = m+ εF +A ∫ BZ d3k (t+ γ)ck D−(~k) (6.49) 0 = m− εF +A ∫ BZ d3k 2m̄+(t− γ)ck D−(~k) We continue as for the previous case, by defining the integrals I−1 and I − 2 , and solving for A. The final result can be written as A± = m̄−m m̄I±1 + tI ± 2 (6.50) ε±F = ∓m∓A [ 2m̄I±1 +(t± γ)I±2 ] 82 and the denominator (we write both cases together now) D±(~k) =− [(t∓ γ)ck] [2m̄+(t± γ)ck]+λ 2sk (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 =∑ µ d2µ (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−∑ µ d2µ = 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 I±1 and I ± 2 numerically (b) Calculate U±0 = √ 12A±(2pi/a)3 and ε±F The advantages of this method are the elimination of the iterative process, prob- ing 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 ∫ BZ d3k εF− γck D(~k)|Σ=0 (6.54) Σ0 = A ∫ BZ d3k m+ tck D(~k)|Σ=0 and the denominator D(~k)|Σ=0 = (εF− γck)2− (m+ tck)2+λ 2sk (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 ck = 2 [ 3−∑ i cos(kia) ] ' 2 [ 3−3+ 1 2 (ak)2 ] = (ak)2 (6.56) sk = −4∑ i sin2(kia)'−4(ak)2 Substitute in the self energy equations Σ4 = A ∫ Λ d3k εF− γ(ak)2 D(~k)|Σ=0 (6.57) Σ0 = A ∫ Λ d3k 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 d3k = 4pik2dk, so that Σ4 = 4piA ∫ Λ 0 dk [εF− γ(ak)2]k2 D(~k)|Σ=0 (6.58) Σ0 = 4piA ∫ Λ 0 dk [m+ t(ak)2]k2 D(~k)|Σ=0 We rearrange the integrals Σ4 = 4piA a2 γ γ2− t2 ∫ Λ 0 dk εF γa2 k 2− k4 ε2F−m2 (γ2−t2)a4 − 2(εFγ+mt)−v2 (γ2−t2)a2 k 2+ k4 (6.59) Σ0 = 4piA a2 t γ2− t2 ∫ Λ 0 dk m ta2 k 2+ k4 ε2F−m2 (γ2−t2)a4 − 2(εFγ+mt)−v2 (γ2−t2)a2 k 2+ k4 The integrals we need to solve are of the form ∫ Λ 0 dk Ak2± k4 B+Ck2+ k4 ' ±Λ (6.60) 85 Where we have kept only the part that diverges linearly for Λ→ ∞. The additional terms are constant in the limit Λ→ ∞. We set Λ= pi/a, substitute the result of the integral and the prefactor A to find Σ4 =−4pia2 U20 12 ( a 2pi )3 γ γ2− t2 pi a = U20 24pi γ t2− γ2 (6.61) and similarly Σ0 = − U 2 0 24pi t t2− γ2 (6.62) Therefore, the renormalized Fermi energy and mass, in the zero order Born ap- proximation are given by ε̄F ≡ εF−ReΣ4 = εF− U 2 0 24pi γ t2− γ2 (6.63) m̄ ≡ m+ReΣ0 = m− U 2 0 24pi t t2− γ2 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 nega- tive (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− U 2 0 24pi t t2− γ2 = 0 (6.64) 86 or U0 = √ 24pim t2− γ2 t ≡Uc (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− U 2 0 24pi 1 γ± 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 con- tains an unusual ‘axion’ term ∼ θ~E ·~B with θ = pi . According to Witten [32] a magnetic monopole inserted into a medium with non-zero θ binds electric charge −e(θ/2pi+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 0 20 40 60 80 100 120 140 160 180 200 ï10 0 10 20 30 40 50 60 U0 E F 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 approx- imation. 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 be- low 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 U0 (meV) 50 100 150 200 250 300 G  (e 2 /h ) 0 1 2 3 4 Bulk+Surfaces Bulk 0 5 10 15 20 25 30 35 40 5 10 15 20 y z x c) d) a) STAI b) ST AI 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 sin- gle 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(ĒF)| = |Re(m̄)|. We use parameters ε = 145meV, t = 24meV, λ = 20meV, and γ = 16meV, corresponding to m = 1meV. The cal- culation of G is outside the scope of this dissertation and is included in Ref. [98]. . with the monopole and anti-monopole low lying states, therefore doubling the lo- calization 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 al- tering 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 ex- pected 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 lo- calized 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 1 25 50 75 100 ï0.50 ï0.25 0.00 0.25 0.50 z δ Q   4997 4999 5000 5003 1 25 50 75 100ï0.20 ï0.10 0.00 0.10 0.20 z δ ρ   δρ δρ̄ a) b) c) Figure 6.4: A prism-shaped sample with monopole (+), anti-monopole (−), and field lines indicated. The asymmetric field distribution allows us to use PBC in all directions and thus avoid difficulties associated with the surface states. nowhere to go. In addition, in order to maximize the localization length probed, 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 on a strong disorder to get a TAI. We measured the induced fractional charge in a configuration containing a monopole and an anti-monopole depicted in Fig. 6.4. Our results presented in Fig. 6.5 and Fig. 6.6 clearly indicate fractional charge ±e/2 bound to the monopole, confirming the expected value of θ = pi . This result lends additional support to our identification of STAI as a genuinely 3D topological phase characterized by the bulk axion term. 90 1 25 50 75 100ï0.20 ï0.10 0.00 0.10 0.20 z δ ρ   δρ δρ̄ 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. 1 25 50 75 100 ï0.50 ï0.25 0.00 0.25 0.50 z δ Q   4997 4999 5000 5003 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 half- filling). The steps of magnitude ±e/2 at the monopole/anti-monopole locations show the expected localized fractional charge due to the Wit- ten 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 useU0 = 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 sur- face 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 Bi2Se3 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 Bi2Se3 [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 order- ing? 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 re- covering 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 Surface magnetic ordering in topological insulators 7.3.1 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ẑ then we get an additional term in the Hamiltonian H =H0+ JM σz 2 , (7.2) 93 where J is the exchange coupling strength. Since H is a sum of anti-commuting matrices we can write the spectrum down immediately Ek =± √ v2k2+(JM/2)2, (7.3) where k2 = k2x + k 2 y . We see that a gap of size JM has opened up. However, even though from a the- oretical standpoint this proposition seems promising, experimentally it has proven very difficult to fabricate a sample with the requisite properties. Two key chal- lenges need to be overcome: first, to observe most of the interesting surface phe- nomena, one requires the surface to remain insulating; however most ferromagnets in nature are metallic. Second, for a ferromagnet in a thin-film geometry, the mag- netization 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 ex- pected since they do not break TRI. Conversely, doping with magnetic impurities, for example Bi2−xFexSe3 , resulted in a spectral gap that increased with the concen- tration of magnetic dopants x, with a gap of 60meV for x = 0.25 (the bulk gap for Bi2Se3 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 (T 3Dc > T 2D c ), 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 T bulkc < T surf c (for another example of this, see Ref. [102, 103]). Therefore, there is a regime T bulkc < T < T surf c 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 sur- face 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 T surfc 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 clos- ing 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 topolog- ical 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 Bi2Se3 Several papers have considered a Hamiltonian which is a so called “discretized” model of Bi2Se3 . We understand the discretization here as projecting the Bi2Se3 lattice on a tetragonal (elongated cubic) lattice. The tetragonal model’s low energy ap- proximation agrees with the low energy approximation of the Bi2Se3 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 Bi2Se3 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 Bi2Se3 . The first is that the kz term of that Hamiltonian can be written (in our notation) as A1kz(σ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 Bi2Se3 we compare the low energy ap- proximation of our model, with the low energy approximation from the literature [34] H(~k) = ε0(~k)+  M(~k) A1kz 0 A2k− A1kz −M(~k) A2k− 0 0 A2k+ M(~k) −A1kz A2k+ 0 −A2kz −M(~k) +O(k2) (7.4) with k± = kx± iky, ε0(~k) =C+D1k2z +D2k2⊥ and M(~k) =M−B1k2z −B2k2⊥. The 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) H(~k) = ε0(~k)+  M(~k) Axyk− 0 Azkz Axyk+ −M(~k) −Azkz 0 0 −Azkz M(~k) Axyk+ Azkz 0 Axyk− −M(~k) +O(k2) (7.5) which can be obtained from the previous matrix by swapping rows and columns (2↔ 4), and ε0(~k) =C+Dxyk2⊥+Dzk2z and M(~k) =M0+Bxyk2⊥+Bzk2z . However, 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 ap- pears 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 pi 4 σ2σx (7.6) So that our model gives H(~k) = d4−µ+  d0 −d1+ id2 0 −d3 −d1− id2 −d0 d3 0 0 d3 d0 −d1− id2 −d3 0 −d1+ id2 −d0  (7.7) We now look at the low energy approximation around the Γ point (0,0,0) sin(kiai) ' kiai (7.8) cos(kiai) ' 1− (kiai) 2 2 97 so that d0 = ε−2∑ i ti cos(kiai)' ε−2∑ i ti+∑ i ti(kiai)2 (7.9) d4 = γ0−2∑ i γi cos(kiai)' γ0−2∑ i γi+∑ i γi(kiai)2 di = −2λi sin(kiai)'−2λikiai =−viki 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 Azkz = −d3 Axyk− = −d1+ id2 which yield simple equations for the parameters of our model. The first equation gives γ0 = C+2 Dz c2 +4 Dxy a2 (7.11) γz = Dz c2 γ⊥ = Dxy a2 The second equation gives ε = M0+2 Bz c2 +4 Bxy a2 (7.12) tz = Bz c2 t⊥ = Bxy a2 98 and the third and fourth equations give λz = Az 2c (7.13) λ⊥ = Axy 2a The parameters were obtained in [34] by fitting to an ab initio calculation of the spectrum of Bi2Se3 and are given by [89] as M0 =−0.28eV, Axy = 2.2eVÅ, Az = 4.1eVÅ, Bxy = 10eVÅ2, Bz = 56.6eVÅ2, C =−0.0068eV, Dxy = 1.3eVÅ2, Dz = 19.6eVÅ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 Bi2Se3 . 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 Bi2Se3 (obtained experimentally) are labelled a′ and c′. The two equations are a2c = √ 3 2 a′2c′ (7.14) a2 c = 3 2 √ 3a′2 c′ and the solutions a = √ 3 2 a′ (7.15) c = c′√ 3 For the lattice constants of Bi2Se3 we use a′ = 4.1387Å (7.16) c′ = 28.629Å 99 x a′ (Å) c′ (Å) 0 4.1387 28.629 0.0023 4.1381 28.627 0.0044 4.1365 28.625 0 4.1396 28.636 0.03 4.106 28.562 Table 7.1: Lattice constants for Bi2Se3 So that we have a = 5.0689Å (7.17) c = 16.5290Å 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) = v~k ·~σ (7.18) where v is the Fermi velocity, which is assumed the same for the x and y directions (as is reasonable for Bi2Se3 ),~k = (kx,ky), and ~σ = (σy,−σx) are Pauli matrices in the spin subspace. 100 On substitution into (7.46) we find HMFe (~k) = He(~k)+ JM σz 2 −µ = v(kxσy− kyσx)+ JM2 σz−µ (7.19) so that the energies are Ee = α √ v2k2+ M̃2 (7.20) where k2 = k2x + k 2 y , 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 (HMFe 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 = α ∂ ∂M √ v2k2+ M̃2 (7.22) = α MJ2 4 ( v2k2+ M̃2 )− 12 = MJ2 4 α uk where we have defined uk ≡ √ v2k2+ M̃2, so that Ee = αuk−µ Substitute into (7.64) to find m = 1 NJ∑kα f (Ee) ∂Ee(~k,α) ∂M (7.23) = JM 4N ∑kα α uk f (αuk−µ) = −JM 4N ∑k 1 uk [1− f (uk+µ))− f (uk−µ)] In the limit µ = 0, which is the case for half filling at T = 0, we find m=−JM 4N ∑k 1−2 f (uk) uk =−JM 4N ∑k tanh(βuk2 ) uk (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 2S(S+1) M βJx (7.25) Which we can now substitute into the expression for m (in the limit µ = 0 for all T ), to find − 3 2S(S+1) M βJx =−JM 4N ∑k tanh(βuk2 ) uk (7.26) After dividing by M and an additional factor on both sides, we find an equation for Tc 1 = S(S+1) 6 βcJ2x 1 N∑k tanh(12βcv|k|) v|k| (7.27) where we have set M = 0 in the energies uk ≡ √ v2k2+ M̃2 ' v|k|, since this term would just give higher order terms in M (we expanded both sides to first order in M). We integrate 1 N∑k tanh(12βcv|k|) v|k| = 2pi ( a 2pi )2 ∫ Λ 0 kdk tanh(12βcvk) vk (7.28) = ( a 2pi )2 4pi βcv2 ∫ 1 2βvΛ 0 dz tanhz = ( a 2pi )2 4pi βcv2 lncosh ( 1 2 βcvΛ ) = ( a 2pi )2 2pi βcv2 βvΛ In the last step we assume that the cutoff is large βcvΛ 1 so that lncoshx ' x. Substituting back into the gap equation, we find T sur fc ' pi S(S+1) 3kB ( Λa 2pi )2 J2x Λv (7.29) 102 For S = 5/2, J = 0.5eV, x = 0.05, v = 2λ⊥a and a cutoff Λa/(2pi) ' 1/5 we find T surfc ' 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, T surfc ' pi S(S+1) 3kB ( Λa 2pi )2 J2x (Λv)2 (Λv−|µ|). (7.30) As expected, T surfc is seen to decrease away from half filling. To complete the argument we would now like to give a similar simple estimate for T bulkc . Unfortunately, the bulk critical temperature is not easy to estimate, since unlike the topologically protected surface state, whose physics is simple and uni- versal, the bulk of a TI can be complicated and T bulkc will generally depend on the details of the band structure and other factors. For Bi2Se3 with 5% concentration of Cr dopants, Ref. [89] estimates T bulkc ' 70K using first-principles numerical cal- culations. Comparing with our rough estimate for T surfc given above we see a clear indication that a T bulkc < T surf c 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−xMnxSe3 (and similar ma- terials), 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 Bi2Te3 [107] (see FIG 11 in that paper). Therefore, we add only an on-site electron-impurity spin-spin interaction term H =∑ ~k Ψ†k [ He(~k)−µ ] Ψk+ J∑ I ~SI ·~sI (7.31) 103 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 ∑ I ~SI ·~sI = ∑ i ηi~Si ·~si (7.32) = ∑ i (ηi~Si− ~M) · (~si−~m)+∑ i ~M ·~si+∑ i ηi~Si ·~m−∑ i ~M ·~m where ηi ≡ 1 i ∈ I0 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 (7.36) 104 The quantity of interest is the average electron spin, given by 〈~s〉= ∑i〈~si〉 2N = ~m 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∑ I ~SI ·~sI ' J ~M ·∑ i ~si+ J~m ·∑ I ~SI−NJ ~M ·~m (7.38) 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 HMF =∑ ~k Ψ†k [ He(~k)+ J ~M · ~σ 2 −µ ] Ψk+ J~m ·∑ I ~SI−NJ ~M ·~m (7.39) where we have used ∑ i ~si =∑ i Ψ†i ~σ 2 Ψi =∑ ~k Ψ†k ~σ 2 Ψk (7.40) Conside the case in which the magnetization is in the z direction, so ~M =Mẑ and ~m= mẑ. We have HMF = ∑ ~k Ψ†k [ He(~k)+ JM σz 2 −µ ] Ψk+ Jm∑ I ~SzI −NJMm (7.41) = HMFe +H MF I −NJMm 105 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(kiai), di ≡−2λi sin(kiai), d4 = γ0−2∑i γi cos(kiai). 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 Bi2Te3 and Bi2Se3 , 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 λI (7.44) 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 HMFe =∑ ~k Ψ†k [ He(~k)+ JM σz 2 −µ ] Ψk =∑ ~k Ψ†kH MF e (~k)Ψk (7.45) and we can write explicitly HMFe (~k) = He(~k)+ JM σz 2 −µ = d4+∑ µ dµγµ + JM σz 2 −µ (7.46) 106 and the energies are Ee = d4+ γ √ ∑ µ d2µ + M̃2+2δM̃ √ d20 +d 2 3 −µ (7.47) = d4+ γ √ d21 +d 2 2 + [√ d20 +d 2 3 +δM̃ ]2 −µ 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± √ ∑ µ d2µ −µ (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 6= 0 we find that the degeneracy is broken by a splitting of 2|M̃|= |JM| 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 β lnZ (7.52) where Z is the many particle partition function, β = 1/(kBT ) and kB is the Boltz- 107 ï4 ï3 ï2 ï1 0 1 2 3 4 5 E Γ 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 =∏ kδγ Zkδγ =∏ kδγ ∑ ni=0,1 e−βniEe(δ ,γ,~k) =∏ kδγ ( 1+ e−βEe(δ ,γ,~k) ) (7.53) So that the electron free energy is Fe =− 1β lnZe =− 1 β ∑kδγ ln ( 1+ e−βEe(δ ,γ,~k) ) (7.54) The partition function for the NI impurities is ZI = ( s ∑ λ=−s e−βJmλ )NI (7.55) so that the impurity free energy is FI =− 1β lnZI =− NI β ln s ∑ λ=−s e−βJmλ (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 ∂m = ∂FI ∂m −NJM = 0 (7.57) So M = 1 NJ ∂FI ∂m = 1 NJ ∂ ∂m [ −NI β ln S ∑ λ=−S e−βJmλ ] (7.58) = −2x ∂ ∂α [ ln S ∑ λ=−S e−αλ ] 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 ∑ λ=−S e−αλ = e−αS ( 1+ eα + e2α + · · ·+ e2αS) (7.59) = e−αS 1− e(2S+1)α 1− eα = e(S+ 1 2 )α − e−α(S+ 12 ) e α 2 − e− α2 = sinh[(S+ 12)α] sinh α2 Substitute into the magnetization, to find M = −2x ∂ ∂α ln [ sinh[(S+ 12)α] sinh α2 ] (7.60) = −x [ (2S+1)coth [ (2S+1) α 2 ] − coth α 2 ] =−2SxBS(αS) where BS(y) is the Brillouin function BS(y) = 2S+1 2S coth 2S+1 2S y− 1 2S coth 1 2S y (7.61) 109 So that our final equation for M is M =−2SxBS(βJmS) (7.62) In the zero temperature limit kBT → 0 we expect the impurity magnetic moments to be aligned and at their extremum value ±S, so that the total magnetization is ±NIS/N = ±2Sx. Indeed, in this limit we get BS → 1 (assuming m > 0) so that M ' −2Sx. In the limit of infinite temperature kBT → ∞ we expect the impurity 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 ∂M = ∂Fe ∂M −NJm= 0 (7.63) So m= 1 NJ ∂Fe ∂M = 1 NJ ∂ ∂M [ − 1 β ∑kδγ ln ( 1+ e−βEe(~k,δ ,γ) )] (7.64) = 1 NJ ∑kδγ e−βEe(~k,δ ,γ) 1+ e−βEe(~k,δ ,γ) ∂Ee(~k,δ ,γ) ∂M = 1 NJ ∑kδγ f (Ee) ∂Ee(~k,δ ,γ) ∂M where f (E) = (1+ eβE)−1 is the Fermi-Dirac distribution. We take the derivative of the energy ∂Ee ∂M = γ ∂ ∂M √ d21 +d 2 2 + (√ d20 +d 2 3 +δM̃ )2 (7.65) = δγJ 2 (√ d20 +d 2 3 +δM̃ )[ d21 +d 2 2 + (√ d20 +d 2 3 +δM̃ )2]−1/2 = δγJ 2 εkδ ukδ 110 where we have defined εkδ = √ d20 +d 2 3 +δM̃ (7.66) ukδ = √ d21 +d 2 2 + ε2kδ So that Ee = γukδ − µ̃ (7.67) and µ̃ ≡ µ−d4. We substitute into the magnetization m to find m = 1 2N ∑kδγ δγ εkδ ukδ f (Ee) (7.68) = 1 2N∑kδ δ εkδ ukδ [ f (ukδ − µ̃)− f (−ukδ − µ̃)] = − 1 2N∑kδ δ εkδ ukδ [1− f (ukδ − µ̃)− f (ukδ + µ̃)] = 1 2N∑k { εk− uk− [1− f (uk−+ µ̃)− f (uk−− µ̃)]− εk+uk+ [1− f (uk++ µ̃)− f (uk+− µ̃)] } where we have used f (−E) = 1 1+ e−βE = 1− 1 1+ eβE = 1− f (E) (7.69) In the limit µ = 0 and γ0 = γi = 0 we find µ̃ = 0 so that m = − 1 2N∑kδ δ εkδ ukδ [1−2 f (ukδ )] (7.70) = − 1 2N∑kδ δ εkδ ukδ tanh ( βukδ 2 ) = 1 2N∑k [ εk− uk− tanh ( βuk− 2 ) − εk+ uk+ tanh ( βuk+ 2 )] 111 where we have used 1−2 f (x) = 1− 2 1+ eβε = eβx−1 eβx+1 = tanh ( βx 2 ) (7.71) The chemical potential µ can be determined by summing over the average occu- pation 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 appropri- ate to discuss the density n= 1 N ∑kδγ f (Ee) = 1 N ∑kδγ f (γukδ − µ̃) (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=∑ ~k,i 〈Ψi|σz2 |Ψi〉=∑ ~k,i (U† σz 2 U)ii f (Ei) (7.73) 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 0 10 20 30 40 50 60 70 80 900 0.5 1 1.5 2 2.5 T (K) 〈S 〉(µ B ) 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 com- pute 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 Bi2Se3 low energy model (which should be correct near the Γ point), we can ex- pect 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 transi- tion. 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 0 10 20 30 40 50 60 70 80 90 1000 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 T (K) M / x (µ B )   4 8 12 16 20 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 2N ∑kδγ δγ εkδ ukδ f (Ee) (7.74) = − 1 2N∑kδ δ εkδ ukδ Now consider a low energy approximation around the Γ point εkδ = √ d20 +d 2 3 +δM̃ ' √ m20+ v2zk2z −δJSx (7.75) ukδ = √ d21 +d 2 2 + ε2kδ ' √ v2⊥k 2 ⊥+ ε 2 kδ 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 Nz ∑ kz,δ δ [ 1 2N⊥ ∑ k⊥ εkδ ukδ ] (7.76) 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 =√ v2k2+ ε2k , so that msur f (T = 0)' JSx2N ∑k 1√ v2k2+(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 2N⊥ ∑ k⊥ εkδ ukδ ' 1 2N⊥ ∑ k⊥ εkδ√ v2⊥k 2 ⊥+ ε 2 kδ (7.78) ' Ω 2 ( a 2pi )2 ∫ Λ⊥ 0 2pik⊥ dk⊥√ k2⊥+Ω2 ' piΩ ( a 2pi )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) = − pi v⊥ ( a 2pi )2 Λ⊥∑ δ δ [ 1 Nz ∑ kz εk ] (7.79) 115 Once again, consider just the expression in the square brackets 1 Nz ∑ kz εkδ ' c 2pi ∫ Λz 0 dkz (√ m20+(vzkz)2−δJSx ) (7.80) = c 2pi vz ∫ Λz 0 dkz √( m0 vz )2 + k2z −δJSxΛz  The integral gives ∫ Λz 0 dkz √( m0 vz )2 + k2z ' 1 2 Λ2z (7.81) so that we find 1 Nz ∑ kz εkδ ' c 2pi [ 1 2 vzΛ2z −δJSxΛz ] (7.82) So m(T = 0) ' − pi v⊥ ( a 2pi )2 Λ⊥∑ δ δ c 2pi [ 1 2 vzΛ2z −δJSxΛz ] (7.83) ' 2pi ( JSx Λ⊥v⊥ )( Λ⊥a 2pi )2(Λzc 2pi ) This result has a quadratic dependence on the cutoff, as expected. Perhaps surpris- ingly, 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 1. We do this by expanding both sides of the equation for m to first order in M̃ m = 1 2N ∑kδγ δγ εkδ ukδ f (Ee) (7.84) We start with the left hand side. Consider the limit of the Brillouin function for 116 small arguments. Since cothx' 1 x + x 3 +O(x3) (7.85) we can approximate the Brillouin function to first order BS(y) = 2S+1 2S coth 2S+1 2S y− 1 2S coth 1 2S y' S+1 3S y (7.86) So that M =−2SxBS(βJmS)'−2S(S+1)3 βJmx (7.87) and therefore m'− 3 S(S+1) M̃ βJ2x (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) ∂ ∂M̃ [ εkδ ukδ f ] = ∂εkδ ∂M̃ f ukδ + ∂ f ∂M̃ εkδ ukδ − ∂ukδ ∂M̃ εkδ u2kδ f (7.89) We take the required derivatives, assuming M̃ > 0 ∂εkδ ∂M̃ = √ d20 +d 2 3 +δ (7.90) and ∂ukδ ∂M̃ = ∂ukδ ∂εkδ εkδ ∂M̃ = εkδ ukδ (√ d20 +d 2 3 +δ ) (7.91) and finally ∂ f ∂M̃ = ∂ f ∂Ee ∂Ee ∂M̃ =−γβ f 2 eβEe εkδ ukδ (√ d20 +d 2 3 +δ ) (7.92) 117 We substitute the last three equations into the previous equation, to find ∂ ∂M̃ [ εkδ ukδ f ] = (√ d20 +d 2 3 +δ ) f ukδ [ 1− γβ f eβEe ε 2 kδ ukδ − ε 2 kδ 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 = (7.94) −βcJ2x S(S+1)6 1 N ∑kδγ δγ (√ d20 +d 2 3 +δ ) f ukδ [ 1− γβc f eβEe ε2kδ ukδ − ε 2 kδ u2kδ ]∣∣∣∣∣ M̃=0 We can take the “surface limit” to check this result. This involves taking all hop- ping parameters along the z direction to zero, taking µ̃ = 0, and taking the low energy approximation. This gives d21 + d 2 2 = v 2k2, d0 = d3 = d4 = 0. So δ = 1, εk = M̃, uk = √ v2k2+ M̃2 and therefore the energies are Ee = γ √ v2k2+ M̃2. On substitution into the previous equation, we find 1 = (7.95) −βcJ2x S(S+1)6 1 N∑kγ γ f√ v2k2+ M̃2 [ 1− γβc f eβEe M̃ 2√ v2k2+ M̃2 − M̃ 2 v2k2+ M̃2 ]∣∣∣∣∣ M̃=0 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 = −βcJ2xS(S+1)6 1 N∑kγ γ f (γv|k|) v|k| (7.96) = βcJ2x S(S+1) 6 1 N∑k tanh (1 2βcv|k| ) v|k| where we have used ∑ γ γ f (γv|k|) = f (v|k|)− f (−v|k|) = 2 f (v|k|)−1 =− tanh(1 2 βv|k|) (7.97) The last equation is identical to the equation derived for the surface case (7.27). 118 7.7 Adding the surfaces Consider now a system with open boundary conditions in the z direction and peri- odic 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 con- ditions 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,µ λµΨ†jτzσµΨ j+µ +h.c. (7.98) Hcd = ε∑ j Ψ†jτxΨ j− ∑ <i, j> tµΨ†i τxΨ j+h.c. Hγ = γ0∑ j Ψ†jΨ j− ∑ <i, j> γµΨ†iΨ j+h.c. 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(kxa)+σy sin(kya)+ iλzT̃τzσz (7.99) + [ε−2t⊥(cos(kxa)+ cos(kya))]τx− tzR̃τx + [γ0−2γ⊥(cos(kxa)+ cos(kya)]− γzR̃ 119 where He =∑ ~k⊥ ϕ†(~k⊥)He(~k⊥)ϕ(~k⊥) (7.100) 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⊥ Ψ(~k⊥)ei ~k⊥·~R j (7.101) (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∑ j ∑ I∈z ~SI ·~sI '∑ j [ J ~M j ·∑ i∈ j ~si+ J~m j ·∑ I∈ j ~SI−N⊥J ~M j ·~m j ] (7.102) where ~M j and ~m j are the magnetization of the impurities and electrons (respec- tively) 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 HMF = (7.103) ∑ ~k⊥ ϕ†(k⊥) [ He(~k⊥)+ J ~M⊗ ~σ 2 −µ ] ϕ(k⊥)+ J ∑ j,I∈ j ~m j ·~SI−N⊥J∑ j ~M j ·~m 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ẑ 120 and ~m= mẑ (M and m are both vectors of length Lz). We have HMF = ∑ ~k⊥ ϕ†(k⊥) [ He(~k⊥)+ JM⊗ σz2 −µ ] ϕ(k⊥)+ J ∑ j,I∈ j m j ·~SI−N⊥J∑ j M j ·m j = HMFe +H MF I −N⊥J∑ j M j ·m j (7.104) We now diagonalize the full electron Hamiltonian He(~k⊥)+JM⊗ σz2 −µ , and find a matrixU 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 =∑ ~k⊥ 〈Ψi|σz2 |Ψi〉 (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 =∑ ~k⊥ 〈ϕi|σz2 |ϕi〉=∑ ~k⊥ 〈ϕ̃i|U†σz2 U |ϕ̃i〉= 4i ∑ j=4i−3 ( U† σz 2 U ) j j 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 T bulkc ' 73K, and the sur- face remains FM ordered up to T surfc ' 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 2 4 6 8 10 12 14 16 18 20 ï1.0 ï0.5 0.0 0.5 1.0 1.5 2.0 2.5 z 〈S 〉(µ B )   0 40 60 80 90 100 120 Figure 7.4: Impurity spin expectation value as a function of the z co- ordinate, 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 Bi2Se3 dispersion close to the Γ point based on first principles calculations [34, 42, 43]: γ0 = 0.3391eV,γ⊥ = 0.0506eV,γz = 0.0717eV,ε = 1.6912eV, t⊥ = 0.3892eV, tz = 0.2072eV,λ⊥ = 0.2170eV and λz = 0.1240eV. of if they are in the bulk or surface, are fully polarized. As we increase the tem- perature, the magnetization of the bulk drops faster than the magnetization of the surface. Eventually we cross T bulkc 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 T surfc , 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 mag- netization as the temperature approached T bulkc 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 → T bulkc , and the proximity of the ordered surfaces, and found T bulkc from the bulk calculation with periodic boundary conditions, the re- sults 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 0 20 40 60 80 100 120 0.0 0.5 1.0 1.5 2.0 2.5 T 〈S 〉(µ B )   surf avg bulk T surfcT bulkc Figure 7.5: Impurity spin expectation value as a function of temperature for the discretized Bi2Se3 lattice. We plot the order parameter on the sur- faces (“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 Bi2Se3 (namely, the fluc- tuations in magnetization), we decided to try another TI model, the so called Per- ovskite lattice [111]. This lattice can be described as a simple cubic lattice with ad- ditional sites on the center of all edges. The nearest neighbor tight binding Hamil- 123 ï0.10 ï0.05 0.00 0.05 0.10 0.15 0.20 0.25 70 75 80 85 90 95 100 µ T c   surf surfa bulk Figure 7.6: Critical temperature as a function of the chemical potential µ for the discretized Bi2Se3 lattice. We plot the surface critical temper- ature (squares, blue), and the approximate surface temperature from Eq. (7.30) (solid line, black). In addition, we plot the bulk critical tem- perature as obtained from an average over the bulk sites (circles, red). The error bars were obtained by fitting the magnetization below the crit- ical temperature to the known behaviour for a second order phase tran- sition, 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/2pi = 0.215. The vertical dashed lines mark the bulk gap. tonian is H(~k) =−2t  0 cos(kxa) cos(kya) cos(kza) cos(kxa) 0 0 0 cos(kya) 0 0 0 cos(kza) 0 0 0  (7.108) 124 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 HSO(~k) = 4λ σx⊗  0 0 0 0 0 0 0 0 0 0 0 i 0 0 −i 0 sin(kya)sin(kza) + (7.109) σy⊗  0 0 0 0 0 0 0 i 0 0 0 0 0 −i 0 0 sin(kxa)sin(kza)+ σz⊗  0 0 0 0 0 0 i 0 0 −i 0 0 0 0 0 0 sin(kxa)sin(kya)  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 =∑ ~k, j 〈Ψ j|σz2 ⊗Oi|Ψ j〉= 1 N/4∑ ~k, j (U† σz 2 U) j j f (E j) (7.110) 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(βJmiS) (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 direc- tion. This phase persists until T bulkc ' 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 0 10 20 30 40 50 60 ï2 ï1 0 1 2 T (K) 〈S 〉(µ B )   1 2 3 4 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 T (K) 〈S 〉(µ B )   AF FM 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 H(~k) = ILz⊗−2t  0 cos(kxa) cos(kya) 0 cos(kxa) 0 0 0 cos(kya) 0 0 0 0 0 0 0 − tR̃⊗ I2 (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 R̃=  Ru+Rd Rd 0 0 0 Ru Ru+Rd Rd 0 0 0 Ru Ru+Rd . . . 0 . . . . . . . . . . . . Rd 0 0 0 Ru Ru+Rd  (7.113) where Ri ju = δi1δ j4, Ri jd = δi4δ j1 and 0 is understood to be a zero matrix of dimen- sion 4. The spin orbit hopping is given by HSO(~k) = 2λ [ σx⊗ T̃x sin(kya)+σy⊗ T̃y sin(kxa)sin(kza) (7.114) + σz⊗  0 0 0 0 0 0 i 0 0 −i 0 0 0 0 0 0 sin(kxa)sin(kya)  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) T̃x =  Tu+Td −Td 0 0 0 −Tu Tu+Td −Td 0 0 0 −Tu Tu+Td . . . 0 . . . . . . . . . . . . −Td 0 0 0 −Tu Tu+Td  (7.115) where T i ju = δi2δ j4, T i jd = δi4δ j2 and 0 is understood to be a zero matrix of dimen- 128 sion 4. Similarly, the matrix T̃y codes the hopping for sites 3 and 4 (site 3 is on the y axis) T̃y =  Tu+Td −Td 0 0 0 −Tu Tu+Td −Td 0 0 0 −Tu Tu+Td . . . 0 . . . . . . . . . . . . −Td 0 0 0 −Tu Tu+Td  (7.116) where T i ju = δi3δ j4, T i jd = δi4δ j3 and 0 is understood to be a zero matrix of dimen- sion 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 SFMi = (Si,1+Si,2+Si,3+Si,4)/4 and an AF order parameter SAFi = (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 ' 4.5K, where it becomes AF ordered. The bulk remains AF ordered until T bulkc ' 60.5K, and the surface remains AF or- dered until T surfc ' 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 0 20 40 60 80 0.0 0.5 1.0 1.5 2.0 2.5 T 〈S F M 〉(µ B )   surf1 surf2 z bulk T surfc T bulkc 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.25eV, t = 1eV, λ = 0.2eV, 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 or- dered. Our conclusions are based on general arguments that involve only universal properties of the topologically protected surface states and on numerical calcula- tions 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 cou- pling. 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 0 20 40 60 80 ï2 ï1 0 1 2 T 〈S A F 〉(µ B )   surf1 surf2 z avg bulk T surfc T bulkc 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 rea- sonable results for T bulkc 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 T bulkc and T surfc . On the other hand, our mean field calculation did not include electron- electron interactions, which would tend to stabilize the long range magnetic order and thereby strengthen the ordered phases. It has been shown recently that electron- electron 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 Bi2Se3 . We note that the critical temperature for the bulk magnetic ordering of the mag- netically doped topological insulator Bi2−xMnxTe3 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 T bulkc ' 70K for Cr-doped Bi2Se3 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 ar- guments 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−xMnxTe3 as a guideline, one may surmise that T bulkc for Fe and Mn doped Bi2Se3 lies in a similar range of temperatures. It is then entirely pos- sible 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 sur- face of a topological insulator [115, 116], and a recent experimental study on Bi2Se3 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 Bi2Se3 , 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 insu- lators (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 θ 6= 0 is predicted to bind a fractional electric charge−eθ/2pi . 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 gen- eral qualitative considerations and by numerical calculations within a minimal lat- tice 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 Bi2Se3 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 experimen- tal 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 Sb2Se3 , 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önig, 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 2pi 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θ/2pi . 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 Bi2Se3 , Bi2Te3 and Sb2Te3 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, Bi2Te3 . 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 Bi2Se3 . 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 CuxBi2Se3. 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án, 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 12 . 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 Ho2Ti2O7. 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 Dy2Ti2O7. 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 Bi2Se3 : 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−xMnxSe3 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 Bi2Se3 and Sb2Se3 . Journal of Magnetism and Magnetic Materials, 304(1):e164 – e166, 2006. → pages 100 [106] T. Jungwirth, Jairo Sinova, J. Mašek, J. Kučera, 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−xMnxTe3. 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 Bi2Se3 . 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 Žitko. 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 Bi2Se3 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 x = z+ z∗ 2 (A.1) y = z− z∗ 2i Therefore ∂x ∂ z∗ = 1 2 (A.2) ∂y ∂ z∗ = − 1 2i 147 If we have a function f (x,y), we can use the chain rule to write d f dz = ∂ f ∂x ∂x ∂ z + ∂ f ∂y ∂y ∂ z (A.3) = 1 2 ( ∂ f ∂x + 1 i ∂ f ∂y ) = 1 2 ( ∂ f ∂x − i∂ f ∂y ) so d dz = 1 2 ( ∂ ∂x − i ∂ ∂y ) (A.4) Similarly, d dz∗ = 1 2 ( ∂ ∂x + i ∂ ∂y ) (A.5) This can also be seen by taking the complex conjugate of the previous equation. Using these two derivatives, we write d2 dzdz∗ = 1 4 ( ∂ ∂x − i ∂ ∂y )( ∂ ∂x + i ∂ ∂y ) (A.6) = 1 4 ( ∂ 2 ∂x2 + ∂ 2 ∂y2 ) = 1 4 ∇2 So that the Laplacian ∇2 = 4 d2 dzdz∗ (A.7) A.1.2 Polar coordinates We want to transform from (z,z∗) to (r,θ), where z = r eiθ (A.8) z∗ = r e−iθ 148 and r = √ zz∗ (A.9) θ = 1 2i ln ( z z∗ ) Using the chain rule d dz = ∂ r ∂ z ∂ ∂ r + ∂θ ∂ z ∂ ∂θ (A.10) = 1 2 e−iθ ( ∂ ∂ r − i r ∂ ∂θ ) and similarly, d dz∗ = 1 2 eiθ ( ∂ ∂ r + i r ∂ ∂θ ) (A.11) Using (A.4) and (A.5) we find the useful result ∂ ∂x − i ∂ ∂y = e−iθ ( ∂ ∂ r − i r ∂ ∂θ ) (A.12) ∂ ∂x + i ∂ ∂y = eiθ ( ∂ ∂ r + i r ∂ ∂θ ) This can also be obtained by going directly from Cartesian coordinates to polar coordinates x = r cosθ (A.13) y = r sinθ so r = √ x2+ y2 (A.14) θ = tan−1 (y x ) 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= tanb, and we can write d tan−1 a da = db da = ( da db )−1 = ( d tanb db )−1 = 1 1+ tan2 b = 1 1+a2 (A.15) Now calculate the derivatives we will need for the chain rule ∂ r ∂x = x r = cosθ (A.16) ∂ r ∂y = y r = sinθ and ∂θ ∂x = ∂ ∂x tan−1 (y x ) =− y r2 =−sinθ r (A.17) ∂θ ∂y = ∂ ∂y tan−1 (y x ) = x r2 = cosθ r Using the chain rule, we find ∂ ∂x = ∂ r ∂x ∂ ∂ r + ∂θ ∂x ∂ ∂θ (A.18) = cosθ ∂ ∂ r − sinθ r ∂ ∂θ and similarly ∂ ∂y = sinθ ∂ ∂ r + cosθ r ∂ ∂θ By combining the last two equations and rewriting using eiθ = cosθ + isinθ we find ∂ ∂x − i ∂ ∂y = e−iθ ( ∂ ∂ r − i r ∂ ∂θ ) (A.19) ∂ ∂x + i ∂ ∂y = eiθ ( ∂ ∂ r + i r ∂ ∂θ ) which is identical to (A.12). 150 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 operatorU (i.e. [H,U ] = 0) that generates the above rotation. By block diagonalizing the Hamiltonian, we break the problem down into sev- eral smaller diagonalization problems. Since direct diagonalization is O(N3), for large matrices this has the potential of drastically reducing the time and mem- ory 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 re- quires 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 diag- onalization 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 num- bering within each quadrant so that a rotation of pi/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  (A.20) and the Hamiltonian will have the form H =  h11 h12 h13 h14 h21 h22 h23 h24 h31 h32 h33 h34 h41 h42 h43 h44  (A.21) where hi j are matrices of side N30 , 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 =  h g f g† g† h g f f g† h 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 =  h a+b f a−b a−b h a+b f f a−b h a+b a+b f a−b h = I4⊗h+Ma⊗a+Ma⊗b+M f ⊗ f (A.24) where Ma =  0 1 0 1 1 0 1 0 0 1 0 1 1 0 1 0 = (I2+σ1)⊗σ1 (A.25) Mb =  0 1 0 −1 −1 0 1 0 0 −1 0 1 1 0 −1 0 = i(I2−σ1)⊗σ2 M f =  1 0 0 1 1 0 0 1 = σ1⊗ I2 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 pi/2, by using u= e−ipi/4σ2⊗ I2 to get M̃a = (I2+σ3)⊗σ1 = ( 2σ1 0 0 0 ) (A.26) M̃b = i(I2−σ3)⊗σ2 = ( 0 0 0 2iσ2 ) M̃ f = σ3⊗ I2 = ( 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 U1 = 1√ 2 ( −1 1 1 1 ) (A.27) and 2iσ2 is diagonalized by U2 = 1√ 2 ( i −i 1 1 ) (A.28) 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,( U†1 0 0 U†2 )( I 0 0 −I )( U1 0 0 U2 ) = ( U†1U1 0 0 −U†2U2 ) = ( I 0 0 −I ) (A.29) Putting the two rotations together we find U = u ( U1 0 0 U2 ) = ( e−i pi 4 σ2⊗ I2 ) 1√ 2  −1 1 0 0 1 1 0 0 0 0 i −i 0 0 1 1  (A.30) 154 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 R̃Ψ1 =Ψ2 (A.32) R̃2Ψ1 =Ψ3 R̃3Ψ1 =Ψ4 and R̃4 = I4. Therefore, we can write Ψ†1gΨ2 =Ψ † 1g(R̃Ψ1) (A.33) Ψ†1 fΨ3 =Ψ † 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.,Ψ†1gΨ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 4N30 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 N30 wavefunctions, and we label the four matrices which have these wavefunctions as their columns as Ψ̃i 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- tion matrix Ψ̃ (thereby rotating each of the wavefunctions separately), so that the matrix Ψ contains the 4N30 wavefunctions of H. Ψ=UΨ̃ (A.35) This can be seen by taking HΨ= EΨ, inserting the identityUU† in between H and Ψ, and multiplying both sides byU† on the left, so that (U†HU)(U†Ψ) =E(U†Ψ), and we identify Ψ̃=U†Ψ, so that Ψ=UΨ̃. 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 R=  0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0  (A.36) which rotates the system by pi/2, thereby taking quadrant one to quadrant two (for example), and in general R  Ψ1 Ψ2 Ψ3 Ψ4 =  Ψ2 Ψ3 Ψ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 di- agonal 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 oper- ator 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 pi/2, we use the rotation operator Rs = e−i pi 4 σ3 = 1√ 2 (I− iσz) = 1√ 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 Rs = ( 1 0 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 H = 4 ∑ i=0 σi⊗Hi (A.40) where σ0 ≡ I2, the identity matrix of rank two. Using this decomposition, we impose the above commutation condition, to find ∑ i [Rs⊗R,σi⊗Hi] = 0 (A.41) 157 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 =  hi ai+bi 0 ai−bi ai−bi hi ai+bi 0 0 ai−bi hi ai+bi ai+bi 0 ai−bi hi = I4⊗hi+Ma⊗ai+Mb⊗bi (A.42) 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, R= 1 2 [I2⊗ (σ1+ iσ2)+σ1⊗ (σ1− iσ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 RH1R† = 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 =  h1 g2 0 −g†1 g†2 h2 g1 0 0 g†1 −h1 −g2 −g1 0 −g†2 −h2  (A.45) H2 =  h2 g1 0 g † 2 g†1 −h1 −g2 0 0 −g†2 −h2 −g1 g2 0 −g†1 h1  (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 σi⊗hi (A.47) 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↓↓ ) (A.49) 159 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 = ( h↑↑+h↓↓ 2 ) (A.51) h1 = ( h↑↓+h↓↑ 2 ) h2 = ( h↓↑−h↑↓ 2i ) h3 = ( h↑↑−h↓↓ 2 ) 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 previ- ous 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 (A.52) 160 and form the matrix that has its eigenvectors as its columns U = 1 2  0 −1 0 i 0 −i 0 1 0 1 0 −1 0 −1 0 1 0 −1 0 −i 0 i 0 1 0 1 0 1 0 1 0 1 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  (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 prob- lem (and any problem whereU 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 ma- trix. 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 U†HU =  B1 B2 B3 B4  (A.54) 161 where B1 = ( h0−h3+2i(b0−b3) ih1−h2+ ã1+ ã∗2+ b̃∗1− b̃2 −ih1−h2+ ã∗1+ ã2− b̃1+ b̃∗2 h0+h3−2(a0+a3) ) B2 = ( h0−h3+2(a0−a3) ih1−h2− ã1− ã∗2+ b̃∗1− b̃2 −ih1−h2− ã∗1− ã2− b̃1+ b̃∗2 h0+h3+2i(b0+b3) ) B3 = ( h0−h3−2(a0−a3) ih1−h2+ ã1+ ã∗2− b̃∗1+ b̃2 −ih1−h2+ ã∗1+ ã2+ b̃1− b̃∗2 h0+h3−2i(b0+b3) ) B4 = ( h0−h3−2i(b0−b3) ih1−h2− ã1− ã∗2− b̃∗1+ b̃2 −ih1−h2− ã∗1− ã2+ b̃1− b̃∗2 h0+h3+2(a0+a3) ) (A.55) where b̃i = (1+ i)bi, ãi = (1+ i)ai, we have decomposed gi = ai+bi into Hermi- tian 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 4N30 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 N30 wavefunctions, and we label the four matrices which have these wavefunctions as their columns as Ψ̃i (i = 1, . . . ,4), and construct the matrix of wavefunctions of H̃ Ψ̃=  Ψ̃1 Ψ̃2 Ψ̃3 Ψ̃4  (A.56) To find the wavefunctions of the original Hamiltonian H, we rotate the wavefunc- tion matrix Ψ̃ (thereby rotating each of the wavefunctions separately), so that the matrix Ψ contains the 4N30 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 simul- taneously 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) =−4pigδ (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) =∑ ~R ~Bs,0(~r−~R) (A.59) 163 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−i ~k·~rd3r = ∫ ∑ ~R ~Bs,0(~r−~R)e−i~k·~rd3r (A.60) = ∑ ~R e−i~k·~R ∫ ~Bs,0(~r)e−i ~k·~rd3r = ∑ ~R e−i~k·~RF [ ~Bs,0(~r) ] =∑ ~G δ~k,~GF [ ~Bs,0(~r) ] 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−i ~k·~rd3r (A.61) = −4pigl̂ ∫ δ (x− y)δ (y− z) [ θ(z)−θ ( z− a 2 )] e−i~k·~rd3r = −4pigl̂ ∫ ∞ −∞ [ θ(z)−θ ( z− a 2 )] e−ik̃zdz = −4pigl̂ ∫ a/2 0 e−ik̃zdz=−2pigae−ik̃ a4 sinc ( k̃a 4 ) l̂ where k̃ ≡ kx+ ky+ kz and sinc(x)≡ sin(x)/x, so that ~Bsk = ∑ ~G δ~k,~GF [ ~Bs,0(~r) ] (A.62) = −2piga∑ ~G δ~k,~G e −ik̃ a4 sinc ( k̃a 4pi ) l̂ = −2piga∑ ~G δ~k,~G e −i pi2 m̃sinc (pi 2 m̃ ) l̂ since ~G= 2pia (mx,my,mz), 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 ~BsG = −2pigae−i pi2 m̃sinc (pi 2 m̃ ) l̂ (A.63) Note that the average magnetic field is given by the zero Fourier component ~BsG=0 =−2pigal̂ (A.64) The vector potential’s Fourier components (for ~G 6= 0) can be found from (4.32) ~AG = i~G G2 ×~BsG =−2piiga ( ~G G2 × l̂ ) e−i pi 2 m̃sinc (pi 2 m̃ ) (A.65) For ~G= 0 we note that G2~AG = i~G×~BsG (A.66) Below we will prove that ~A~G = ~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 ∫ d3k→ ∑k(2pi/a)3, or by calculating ~A(~r) = 1 (2pi)3 ∫ ~A~k e i~k·~r d3k (A.67) = 1 (2pi)3 ∫ ( ∑ ~G δ~k,~G~A~G ) ei~k·~r d3k = 1 (2pi)3 (2pi)3 a3 ∑ ~G ∫ δ (~k− ~G)~A~G ei ~k·~r d3k = 1 a3∑ ~G ~A~G e i~G·~r 165 From here we also see that since ~A(~r) is real, we must have ~A~G = ~A ∗ −~G since ~A∗(~r) = 1 a3∑ ~G ~A∗~G e −i~G·~r (A.68) = 1 a3∑ ~G ~A∗−~G e i~G·~r = 1 a3∑ ~G ~A~G e i~G·~r 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 = 2pi Φ0 ∫ ~Rd ~Rc ~A(~r) ·d~r = 1 a3 2pi Φ0∑~G ~A~G · n̂ ∫ ~Rd ~Rc ei~G·~rdr (A.69) = −ig(2pi) 2 a2Φ0 ∑~G 1 G2 e−i pi 2 m̃sinc (pi 2 m̃ )[( ~G× l̂ ) · n̂ ] ei~G·~Rc ∫ ~Rd−~Rc 0 ei~G·~rdr where the unit vector n̂ points from site c to site d. Look at ( ~G× l̂ ) · n̂dr = 2pi a √ 3  (my−mz)dx n̂= x̂ (mz−mx)dy n̂= ŷ (mx−my)dz n̂= ẑ (A.70) recall that l̂≡ (1,1,1)/√3. Compute the integral for example for ~Rc= a0(cx,cy,cz), ~Rd = ~Rc+a0x̂ so n̂= x̂, dr = dx and y= z= 0 in the integral [( ~G× l̂ ) · n̂ ]∫ ~Rd−~Rc 0 ei~G·~rdr = 2pi a √ 3 (my−mz) ∫ a0 0 e 2pii a mxxdx = 1 i √ 3 my−mz mx ( e2piimx a0 a −1 ) 166 So that θ xcd = − (2pi)2g√ 3a2Φ0 ∑ ~G 1 G2 my−mz mx e−i pi 2 m̃sinc (pi 2 m̃ ) ei~G·~Rc ( e2piimx a0 a −1 ) (A.71) which can be rewritten by putting in another sinc function, and cancelling the imag- inary part θ xcd = − ga0√ 3Φ0 ( 2pi a )3 ∑ ~G6=0 1 G2 (my−mz)sinc (pi 2 m̃ ) sinc ( pimx a0 a ) iei(~G·~Rc− pi 2 m̃+pimx a0 a ) (A.72) = 2ga0√ 3Φ0 ( 2pi a )3 ∑ |~G|>0 1 G2 (my−mz)sinc (pi 2 m̃ ) sinc ( pimx a0 a ) sin ( ~G ·~Rc− pi2 m̃+pimx a0 a ) which can be seen by noting that i ∑ |~G|6=0 cG eidG = i ∑ |~G|6=0 cG[cos(dG)+ isin(dG)] (A.73) = ∑ |~G|6=0 [icG cos(dG)− cG sin(dG)] =−2 ∑ |~G|>0 cG sin(dG) 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 ' ∫ G2dG/G3 ' logG. So we expect slow convergence. A.3.2 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 ·d~s= 4pir2B(r) = 4pig (A.74) so that the magnetic field of a single monopole is ~B(~r) = g r2 r̂ = g ~r r3 (A.75) We wish to add a monopole of flux 2pi (flux is dimensionless in these units), so that 4pig= 2pi 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 r2 sinθ ∂Aθ ∂θ r̂ = B(r)r̂ We solve B(r) = 1 r2 sinθ ∂Aθ ∂θ = g r2 (A.78) To find, after choosing an arbitrary constant Aθ =−g(1+ cosθ) (A.79) So ~A(~r) =−g(1+ cosθ)~∇ϕ (A.80) 168 which can be rewritten using ~∇ϕ = 1 r sinθ ϕ̂ = 1 r sinθ (−sinϕ,cosϕ,0) (A.81) so that ~A(~r) = g(1+ cosθ) r sinθ (sinϕ,−cosϕ,0) (A.82) 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 ~r× n̂ r−~r · n̂ = g r yx̂− xŷ r− z (A.83) where ~r = (x,y,z), r = |~r|, n̂ = (0,0,1) so that the Dirac string stretches in the ẑ 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 ·d~l (A.84) where ~R are the lattice vectors and the integral must be done for all nearest neigh- bors. Note that for d~l = dzẑ the integral is zero, since Az = 0, so that θ zcd = 0 for all 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 dx r(r− z) = gyc zc ∫ xc+a xc dx ( 1 r− zc − 1 r ) (A.85) = g [ tan−1 ( xc+a yc ) − tan−1 ( xc yc ) + tan−1 ( (xc+a)zc yc √ (xc+a)2+ y2c+ z2c ) − tan−1 ( xczc yc √ x2c+ y2c+ z2c )] 169 where we have used the identity 1 r(r− z) = 1 z ( 1 r− z − 1 r ) (A.86) and the integrals∫ dx r = ln(x+ r) (A.87)∫ dx r− z = z y [ tan−1 ( x y ) + tan−1 ( xz yr )] + ln(x+ r) (A.88) Similarly, for d~l = dyŷ θ ycd = −gxc ∫ yc+a yc dy r(r− z) =− gxc zc ∫ yc+a yc dx ( 1 r− zc − 1 r ) (A.89) = −g [ tan−1 ( yc+a xc ) − tan−1 ( yc xc ) + tan−1 ( (yc+a)zc xc √ x2c+(yc+a)2+ z2c ) − tan−1 ( yczc xc √ x2c+ y2c+ z2c )] which can be obtained from θ xcd by replacing g→−g and x↔ y. 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 ·d~s= 2piraB(r) = 4pig (A.90) 170 so that the magnetic field of a planar monopole is ~B(~r) = 2g ra r̂ = 2g ~r a (A.91) on the two closest layers, and zero elsewhere. We wish to add a planar monopole of flux 2pi (flux is dimensionless in these units), so that 4pig = 2pi and therefore g= 1/2. We choose a vector potential ~A= 2g 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 2pi , 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 ~B= ~∇×~A= 2g a 1 r ∣∣∣∣∣∣∣ r̂ rϕ̂ ẑ ∂ ∂ r ∂ ∂ϕ ∂ ∂ z 0 0 ϕ ∣∣∣∣∣∣∣= 2g ra r̂ (A.93) We want to calculate the Peierls phase factors given by Θcd = ∫ ~Rd ~Rc ~A ·d~l (A.94) for d~l = dxx̂,dyŷ get Θxcd =Θ y cd = 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 zc ϕdz= 2gϕc (A.95) 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,pi/2,pi,3pi/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 Θzcd = 2g tan −1 ( yc xc ) (A.96) In order to avoid problems arising from using tan−1, we can calculate eiΘ directly without the use of trigonometric functions eiΘ z cd = e2gϕci = ( xc+ iyc√ x2c+ y2c )2g (A.97) A.4 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−10 on the left I = (G−10 −Σ)G (A.100) 172 so we have a closed expression for the Green’s function as a function of the self energy and bare Green’s function G= (G−10 −Σ)−1 (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 approxi- mations. 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 an- other, 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)ddx ) = nimp ∫ U(~x)ddx (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∑ ~k′ U2(~k−~k′)G0(~k′) (A.105) which is the self energy in this approximation. In the self consistent Born approxi- mation, we now replace G0 by G in this expression, to find Σ ' Nimp∑ ~k′ U2(~k−~k′)G(~k′) (A.106) = Nimp∑ ~k′ U2(~k−~k′)[εF−H0(~k′)−Σ(~k′)+ i0+]−1 where in the last line we have substituted (A.103). For point-like impurities, the interaction is a delta functionU(~x) = uδ (d)(~x), so U(~k) = 1 V ∫ uδ (d)(~x)ei~k·~xddx= u (A.107) Therefore, the self energy simplifies Σ = Nimpu2∑ ~k [εF−H0(~k)−Σ(~k)+ i0+]−1 (A.108) It is often convenient to evaluate this sum in the large system limit, where the sum becomes an integral Σ = Nimp ( a 2pi )d u2 ∫ ddk[ε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 N3 copies of the BZ, so we can write Σ = NimpNdu2 ∑ ~k∈BZ [εF−H0(~k)−Σ(~k)+ i0+]−1 (A.110) 174

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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0073186/manifest

Comment

Related Items