UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A study of topological insulators in two and three dimensions Weeks, William Conan 2012

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

Item Metadata


24-ubc_2012_fall_weeks_conan.pdf [ 7.62MB ]
JSON: 24-1.0072924.json
JSON-LD: 24-1.0072924-ld.json
RDF/XML (Pretty): 24-1.0072924-rdf.xml
RDF/JSON: 24-1.0072924-rdf.json
Turtle: 24-1.0072924-turtle.txt
N-Triples: 24-1.0072924-rdf-ntriples.txt
Original Record: 24-1.0072924-source.json
Full Text

Full Text

A Study of Topological Insulators in Two and Three Dimensions by William Conan Weeks  B.Sc., Carleton University, 2004 M.Sc., The University of British Columbia, 2007  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) July 2012 c William Conan Weeks 2012  Abstract In this work we introduce a new model for a topological insulator in both two and three dimensions, and then discuss the possibility of creating a fractional topological insulator in the three dimensional version. We also explore the possibility of engineering a quantum spin Hall phase in Graphene through the application of heavy metal adatoms. We show that the two dimensional model on the Lieb lattice, and its three dimensional counterpart, the so called edge centered cubic or Perovskite lattice, possess non-trivial Z2 invariants, and gapless edge modes, which are the signatures of the topological insulating state. Having established that these lattices can become topological insulators, we then tune several short range hoppings in the model and show that it is possible to flatten the lowest energy bands in each case. After flattening the bands we add in a Hubbard term and then use a mean field decoupling to show that there is a portion of phase space where the system remains nonmagnetic, and then conjecture that the many body ground state in three dimensions could become a fractional topological insulator. For the model on graphene, we start by using density functional theory (DFT) to find a pair of suitable heavy elements that are non-magnetic, have a strong spin orbit coupling and prefer to sit at the center of the hexagonal lattice. We then establish that for adatoms distributed in a periodic configuration, again with DFT, that a gap will open at the Dirac point in the presence of spin orbit coupling. To prove the gap is topologically non-trivial, we show that it is possible to adiabatically connect this model to the original Kane-Mele model, a known topological insulator. Lastly, we show that for adatoms distributed randomly the effect survives.  ii  Preface This thesis is based on a series of publications authored by me, my supervisor Marcel Franz (Chapters 2, 3), and our collaborators at the University of California, Irvine, Jason Alicea, Jun Hu and Ruqian Wu (Chapter 5). Chapter 1- Introduction This chapter contains a brief overview of the field of topological insulators, and is provided to motivate the later chapters. It contains no original work, and is inspired by the several excellent review articles currently available on the subject. Chapter 2- Topological insulators on the Lieb and perovskite lattices This chapter is based on our publication in Physical Review B, 82, 085310 (2010). I had the original idea to explore the lattice models, and was then helped by my supervisor’s expertise to work through all of the details and write the paper. Chapter 3- Flat bands with non-trivial topology in three dimensions This chapter is based on our publication in Physical Review B, 85, 041104(R) (2012). It was based on an original idea by my supervisor, to apply the same techniques seen in a series of recent papers to our model from Chapter 2. I carried out all of the numerical simulations and wrote the initial draft, which was then polished into its final form with help from my supervisor. Chapter 4- Quantum transport This chapter is another review on the quantum transport theory required to calculate various quantities of interest in chapter 5. Again, there were several excellent exisiting reviews of the technique that we have benefitted from reading. Chapter 5- Engineering a quantum spin Hall state in graphene iii  Preface This chapter is based on our publication in Physical Review X, 1, 021001 (2011). It was based on an original idea by Jason Alicea, with individual contributions coming from all authors to arrive at the final product. I was responsible for carrying out the transport calculations involving randomly distributed adatoms and writing that section of the paper.  iv  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  Abstract Preface  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . .  ix  1 Introduction . . . . . . . . . . . . . . . . . . . . . . 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . 1.2 Quantum Hall effect and topological band theory 1.3 The bulk-boundary correspondence . . . . . . . . 1.4 Graphene . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Low energy theory . . . . . . . . . . . . . . 1.4.2 Nanoribbons and edge modes . . . . . . . . 1.5 Haldane model . . . . . . . . . . . . . . . . . . . . 1.6 Z2 topological insulator . . . . . . . . . . . . . . . 1.7 Kane-Mele model . . . . . . . . . . . . . . . . . . 1.8 HgTe quantum wells . . . . . . . . . . . . . . . . . 1.9 3D topological insulators . . . . . . . . . . . . . . 1.9.1 The first 3D topological insulator: Bi1−x Sbx  . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  . . . . . . . . . . . . .  1 1 2 4 4 6 8 9 11 13 14 16 19  2 Topological insulator on the Lieb and perovskite lattices 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Lieb lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Perovskite lattice . . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . .  23 23 23 28 33  v  Table of Contents 3 Flat bands with non-trivial topology in three dimensions 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Lieb lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Perovskite lattice . . . . . . . . . . . . . . . . . . . . . . . 3.5 Magnetic instabilities . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . .  . . . . . .  35 35 37 37 40 42 44  4 Quantum transport . . . . . . . . . . . . . . 4.1 Landauer formalism and Fisher-Lee relation 4.2 Introduction to Green’s functions . . . . . 4.3 Semi infinite leads-Self energy description . 4.4 Local current . . . . . . . . . . . . . . . . . 4.5 Recursive Green’s function technique . . . 4.6 Surface Green’s function . . . . . . . . . .  . . . . . . .  46 46 48 50 54 55 59  5 Engineering a quantum spin Hall state in graphene . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Physics of a single heavy adatom . . . . . . . . . . . . . . . . 5.3 Periodic adatom configurations . . . . . . . . . . . . . . . . . 5.4 Randomly distributed adatoms . . . . . . . . . . . . . . . . . 5.5 Relation to previous work . . . . . . . . . . . . . . . . . . . . 5.6 Discussion of the effective Hamiltonian for graphene with a single adatom . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Density functional theory . . . . . . . . . . . . . . . . . . . . 5.8 Adiabatic continuity of the adatom and Kane-Mele models . 5.9 Graphene-only effective Hamiltonian . . . . . . . . . . . . . . 5.10 Solution by canonical transformation . . . . . . . . . . . . . 5.11 Effective Hamiltonian SOC terms . . . . . . . . . . . . . . . 5.12 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .  62 62 64 67 71 79  . . . . . . . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  . . . . . . .  82 89 90 92 95 97 99  6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103  vi  List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 1.10 1.11 1.12 1.13 1.14 1.15  Graphene lattice and corresponding Brillouin zone Graphene dispersion . . . . . . . . . . . . . . . . . Zigzag and armchair nanoribbon . . . . . . . . . . Bandstructure for zigzag and armchair nanoribbon Haldane model . . . . . . . . . . . . . . . . . . . . Phase diagram for Haldane model . . . . . . . . . Single Kramers pair . . . . . . . . . . . . . . . . . Double Kramers pair . . . . . . . . . . . . . . . . . Kane-Mele dispersion plots . . . . . . . . . . . . . Kane-Mele edge modes . . . . . . . . . . . . . . . . HgTe bandstructure . . . . . . . . . . . . . . . . . HgTe conductance . . . . . . . . . . . . . . . . . . Fermi surface - Weak and strong TI . . . . . . . . Bi1−x Sbx . . . . . . . . . . . . . . . . . . . . . . . ARPES data for Bi2 Se3 . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  . . . . . . . . . . . . . . .  5 7 8 9 10 10 12 12 13 15 16 17 18 20 22  2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8  Lieb lattice . . . . . . . . . . . . . . . . . . Topologically trivial insulators - Lieb . . . . Lieb lattice edge modes . . . . . . . . . . . Lieb lattice wave function . . . . . . . . . . Perovskite lattice . . . . . . . . . . . . . . . Brillouin zone and perovskite bandstructure Topologically trivial insulators - perovskite Perovskite edge modes . . . . . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  . . . . . . . .  24 25 28 29 30 31 32 33  3.1 3.2 3.3 3.4 3.5  Lieb lattice flat band . . . . Flat band edge modes . . . Perovskite lattice flat band Gap evolution perovskite . . Phase diagram for magnetic  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  . . . . .  38 39 40 41 43  . . . . . . . . . . . . . . . . . . . . . . . . . . . . instabilities  . . . . .  . . . . .  vii  List of Figures 3.6  Lieb and perovskite bandstructure in magnetic phase . . . . .  44  4.1 4.2  Recursive Green’s function schematic . . . . . . . . . . . . . . Full recursive Green’s function schematic . . . . . . . . . . .  56 58  5.1 5.2 5.3 5.4 5.5  Adatom on graphene . . . . . . . . . . . . . . . . . . . . . . . DFT and tight binding plots with and without SOC . . . . . Gap size vs adatom concentration . . . . . . . . . . . . . . . . Two terminal conductance schematic . . . . . . . . . . . . . . Nanoribbon bandstructure with SOC and staggered sublattice potential . . . . . . . . . . . . . . . . . . . . . . . . . . . Conductance reference plots for armchair nanoribbon . . . . . Averaged conductance for random adatom distributions at differenct concentrations . . . . . . . . . . . . . . . . . . . . . Current distribution . . . . . . . . . . . . . . . . . . . . . . . Conductance plot for largest system size . . . . . . . . . . . . Hexagonal rashba dispersion plots . . . . . . . . . . . . . . . Hexagonal rashba conductance plots . . . . . . . . . . . . . . Correlated and uncorrelated disorder . . . . . . . . . . . . . . Charge density for 4 × 4 supercell . . . . . . . . . . . . . . . . Adiabatic deformation of energy gap . . . . . . . . . . . . . . Graphene only Hamiltonian on 4 × 4 supercell . . . . . . . . .  63 68 70 73  5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15  74 75 76 77 78 79 80 81 86 91 95  viii  Acknowledgements It has been a great privilege for me to work with my supervisor, Marcel Franz, and I would like to thank him, above all, for giving me the opportunity to be a part of his group over the last five years. His incredible level of skill and knowledge is truly inspiring, and his extraordinary ability to somehow always be available to field questions has proven invaluable over the course of this degree. I am grateful to all of my friends and colleagues, but would especially like to thank Andy Bond, Matt Bouchard, Glen Goodvin, Justin Malecki, Dominic Marchand, Matt Meadows, Wes Miron, Anindya Mukherjee and Gili Rosenberg for keeping life balanced over the years, and my girlfriend Melinda Medina, for giving me a home, and traveling to far away places with me, and for letting me win a few matches of tennis. I would also like to single out Bill Young, for reading several of our papers of his own free will. I would like to thank the members of my supervisory/examining committee, Mona Berciu, Kirk Madison, Mark Van Raamsdonk, Joshua Folk and Alireza Nojeh, for donating their valuable time and helping me with constructive feedback. In addition, I would like to thank Jason Alicea, for having a great idea and giving us all an excuse to collaborate, which was one of the highlights of the degree. My family, scattered in various locations around the globe, have always been supportive and so I would very much like to thank them, and in particular my wonderful father, for paving the way with his own adventures in physics. Lastly, I would like to dedicate this humble thesis to the memory of my twin brother, Ashley Kenneth Weeks.  ix  Chapter 1  Introduction 1.1  Overview  The field of topological insulators has taken off in the last few years, proving to be an exciting area of research in condensed matter physics, while also having ties to fundamental physics due to the possibility of providing a testbed for exotic particles like Majorana fermions, [1] fractionally charged vortices,[2] axions [3, 4] and magnetic monopoles.[5, 6] In order to motivate the topics discussed within, and for general interest, we include a brief history of topological insulators and some of their defining characteristics. There are several excellent review articles available that we have benefitted from reading and modelled this introductory section after, [7–9] and we encourage anyone wishing to go into greater depth to start there. In general, topological insulators are bulk insulators having an energy gap between the valence and conduction bands, but with gapless edge states (2D) or surface states (3D), caused by spin-orbit coupling, which are immune to non-magnetic impurities and geometric perturbations. Similar in spirit to the quantum Hall effect, the novelty of a topological insulator (TI) lies in the fact that it can be characterized by a topological invariant, and is not the result of a spontaneously broken symmetry. However, the invariant can not take on any integer value, as in the quantum Hall case, but is instead a so called Z2 topological quantum number that we label ν, which can be either 0 or 1. In 2D a single Z2 invariant is sufficient for the job, whereas in 3D one needs 4 Z2 invariants (ν0 ; ν1 ν2 ν3 ), which distinguish 16 possible topological phases and fall into the two general classes: weak (WTI) and strong (STI) topological insulators.[10] In 3D if there are an odd number of surface states then we have ν0 = 1 and the system is in the STI phase. If there is an even number of surface states (possibly zero) then ν0 = 0 and we have a trivial insulator or WTI if any of the νi are nonzero. The 2D topological insulator, also known as the quantum spin Hall (QSH) phase, was first predicted to arise in graphene, [11, 12] though later 1  1.2. Quantum Hall effect and topological band theory it was realized that, due to carbon’s small atomic number, spin-orbit coupling is too weak to produce an observable effect under realistic conditions. (First-principles calculations[13–17] on pristine graphene all predict subKelvin gaps for this phase). It was then predicted theoretically,[18] and later confirmed experimentally,[19] that this effect is present in the HgTe/CdTe quantum wells. Fu and Kane extended the concept to 3D with a toy model on the diamond lattice,[10] and also predicted that the alloy Bi1−x Sbx , which has large spin orbit coupling, would be a three dimensional TI.[20] This was later verified experimentally by ARPES measurements which directly observed the topological surface states.[21] Other models with a low energy Dirac structure, namely the 2D kagome and 3D pyrochlore systems, have been recently shown to support topological states. [22–24] Experimentally the most promising system is currently Bi2 Se3 , which is known to have a large bandgap and a single surface Dirac cone. [25] Although not experimentally verified, there is also some theoretical work on topological insulators in cold atoms, which looks promising.[26, 27] On a practical level, these gapless surface states are robust against disorder, since time reversal (T ) symmetry disallows back scattering from non magnetic impurities. When gapped by time-reversal symmetry breaking perturbations, superconducting or excitonic pairing, these surface states give rise to interesting insulating phases with exotic quasiparticle excitations already mentioned above. Electric manipulation of the spin degrees of freedom in topological insulators should also be possible and they can therefore be of value in spintronics applications. Some of the existing topological insulators are also known to be among the best thermoelectric materials, such as the material Bi2 Te3 discussed by Zhang et. al., [28] and thus discovery of new materials or further understanding of existing ones would be of value.  1.2  Quantum Hall effect and topological band theory  Until quite recently, experimental studies of topological phases were confined to the quantum Hall effect, seen in two-dimensional electron systems subjected to strong magnetic fields at low temperatures, where the theoretical understanding of the topological nature of band structures also began. This remarkable effect, discovered by von Klitzing in 1980,[29] gave to incredibly  2  1.2. Quantum Hall effect and topological band theory high precision, quantized values of the Hall conductance σ=ν  e2 . h  (1.1)  Electrons in such a system end up with quantized circular orbits having for energies the famous Landau Levels, En =  n+  1 2  ωc  (1.2)  where ωc is the cyclotron frequency. If N landau levels are completely filled, then there should be an energy gap separating occupied and unoccupied states, resulting in an insulator. However, another key feature of the quantum Hall phase is that metallic edges form at the boundaries of the system, resulting in the observed conductance. Naturally the search then began to understand how such precise quantization could arise independent of the geometry or purity of the system. It was in 1982 in another remarkable paper that Thouless et al,[30] discovered a connection between the Hall conductivity, given by the Kubo formula, and a topological invariant, the so called first Chern number. This Chern number is given by the winding number of the Berry phase of the electron wavefunctions around the Brillouin Zone, and is necessarily an integer, insensitive to details of the band structure itself. More concretely, for a Bloch Hamiltonian H(k) having energy bands En (k), the Chern number of the n-th band can be written as Cn =  1 2π  d2 kF12 (k),  (1.3)  BZ  where F12 (k) is the Berry curvature defined as F12 (k) =  ∂ ∂ A2 (k) − A1 (k), ∂k1 ∂k2  (1.4)  and Aµ (k) = −i unk | ∂k∂µ |unk is the Berry connection where |unk are the  the Bloch wavefunctions of the nth band. The full Chern number is then the sum of Cn over all occupied bands. For a more intuitive picture, the Chern number can also be thought of as an analogue to the genus of a surface, g, which counts the number of holes in the surface. For instance, a sphere has g = 0, while a torus has g = 1.  3  1.3. The bulk-boundary correspondence Just as the genus changes only if a hole closes or opens, the Chern number cannot change unless the energy gap closes.  1.3  The bulk-boundary correspondence  A fundamental consequence of these gapped band structures having non trivial topology is the existence of gapless modes at interfaces where the topological invariant changes.[31–33] As such, the edge states mentioned above, existing at the interface between the integer quantum Hall state and vacuum, turn out to be inextricably linked to the topology of the bulk. Furthermore, the electronic states giving rise to this motion are chiral. That is to say they propagate only in one direction along the edge, a property also related to the underlying topology. These states are insensitive to disorder because there are no states available for backscattering, which underlies the perfectly quantized electronic transport. If we imagine an interface where we smoothly transform between a quantum Hall state (ν = 1) and a trivial insulator (ν = 0), the energy gap has to vanish somewhere along the way, or else it would be impossible for the topological invariant to change. Therefore there must be at least one band comprised of states existing at the boundary, crossing the Fermi energy and connecting the conduction and valence bands of the bulk bandstructure. By changing the Hamiltonian near the interface it is possible to modify the dispersion of the edge states. For example, the edge modes could develop a kink so that they intersect the Fermi energy three times – twice with a positive group velocity and once with a negative group velocity. However, the difference NR − NL between the number of right and left moving modes, can not change, and is determined by the topological structure of the bulk states. This is summarized by the bulk-boundary correspondence: NR − NL = ∆n,  (1.5)  where ∆n is the difference in the Chern number across the interface.  1.4  Graphene  Although interesting for many other reasons, graphene will provide us with a concrete example of the quantum Hall effect in a band theory (section 1.5) and was the first lattice shown to host a topological insulating phase 4  are (a = 1.42 ˚ A): �a1 =  √  3a  �  √ � 1 3 , , 2 2  �a2 =  √  �  √ � √ 1 3 3a − , , �a1 − �a2 = 3a(1, 0a) 2 2  (1)  �√  � � √ � 3 1 3 1 �δ1 = a �δ2 = a (0, 1) , �δ3 = a − ,− ,− 2 2 2 2 1.4. Graphene �√ � � √ � 4π 3 1 4π 3 1 �1 = �2 = K , , K − , , 2 entire 2 2 5,2 so we would like to (section 1.7). It is 3a also the focus of3achapter � √ � � � provide a brief introduction here. 4π 4π 1 3 � Graphene isK made atoms √ of, carbon √ = out 0 , K = arranged , in hexagonal structure 2 2 3 3a as shown in Fig. 1.1(a). The structure is 3not3aa Bravais lattice but can be �√ � � √ � seen as a triangular lattice with 2π 3 1 a basis of two 2πatoms 3per1 unit cell, giving us M1 B(white) = , , M2 = − , , the A (blue) and sublattices. 3a 2 2 3a 2 2  (a)  (2) (3) (4) (5)  (b)  Figure 1.1: Graphene lattice (a) and corresponding Brillouin zone (b). Figure 1: Left: Lattice geometry, with primitive vectors denoted a1,2 . The solid and hollow circles represent the position of the A and B atoms. The shaded parallelogram marks the The lattice vectors can be written as: primitive unit cell. Right: The first BZ of the honeycomb lattice. a √ a √ a1 = ( 3, 3) , a2 = (− 3, 3) , (1.6) 2 2 where a ≈ 1.42 ˚ A is the carbon-carbon distance, while the reciprocal lattice vectors are given by: 1 2π √ 2π √ K1 = ( 3, 1) , K2 = (− 3, 1) . (1.7) 3a 3a Of particular importance for the physics of graphene are the two points K and K , seen in Fig. 1.1(b). Their positions in momentum space are given by: √ 4π 1 3 4π √ ,0 , K= K = √ , . (1.8) 3 3a 3 3a 2 2 The tight-binding Hamiltonian for electrons on the lattice looks as fol-  5  1.4. Graphene lows  a†σ,i bσ,j + h.c.  H = −t  (1.9)  i,j ,σ  where ai,σ (a†i,σ ) annihilates (creates) an electron with spin σ (σ =↑, ↓) on site Ri on sublattice A, and similarly for sublattice B, and t is the nearest neighbor hopping energy (approx. 2.7 eV). Going over to momentum space, we can then write the hamiltonian in the form a†k b†k hk  H= k  where hk =  ak bk  (1.10)  0 f (k) f ∗ (k) 0  (1.11)  and f (k) = 1 + e−ik·a1 + e−ik·a2  (1.12)  Solving for the energy, we then have bands of the form: E± (k) = ±t  3 + 2 cos [k · a1 ] + 2 cos [k · a2 ] + 2 cos [k · (a1 − a2 )] (1.13)  or expanding out the cosine terms E± (k) = ±t  1 + 4 cos  √ 3 kx a cos 2  3 ky a + 4 cos2 2  √  3 kx a 2  (1.14)  consisting of a valence band and conduction band that are particle hole symmetric. They touch at six points, only two of which are distinct, the so called Dirac points, K and K singled out above. The others can be related by the reciprocal lattice basis vectors. And so technically graphene is a metal as there is no gap, but since the density of states vanishes at the K points, it is more accurately described as being a semi-metal or a gapless semiconductor.  1.4.1  Low energy theory  We now explore the low energy regime of the model, where the dispersion is approximately linear, and a connection can be made with the Dirac equation. If we decompose the vector k into one of the vectors K or K plus a small 6  1.4. Graphene  Figure 1.2: Graphene dispersion with magnified Dirac cone. (Figure from Castro Neto et al. [34])  vector q, and then substitute into the Hamiltonian above, to first order in q we find hq =  0 αqx + iqy αqx − iqy 0  (1.15)  Eq± = ± vF |q|  (1.16)  where α = ±1, for K and K respectively. The energy dispersion then reduces to a linear form  where  t3a ≈ 8.7 × 105 m/s (1.17) 2 Comparing this to the relativistic dispersion for particles resulting from the Dirac equation E = ± c2 2 k2 + m2 c4 (1.18) vF =  we see that for massless particles where m = 0, E = ± c|k|  (1.19)  and the dispersion has the same form as above, with the speed of light c taking the place of the Fermi velocity vF .  7  1.4. Graphene  1.4.2  Nac  Nanoribbons and edge modes  We introduce the two standard types of graphene nanoribbon, the zigzag and armchair(Figs. 1.3 (a) and (b) respectively), as the zigzag variety is known studied to support of its own, and the armchair nanoribbons 28 5.2 Band structure studied wit 5.2 Band structure withedge the modes Dirac equation are used exclusively for the numerical simulations in Section 5.4. They are both elongated strips of graphene with a finite width, obtained by cutting a graphene sheet along one of two directions with a difference of 30 degrees y between y 1 2 3 y 1 2 Nac Nzz 1the2 two cuts. 3 (a) (b) B  A  A B  B  B  A  A  B  B  B  A  A  a  Wac+  √  3a 2  x √ Wac+ 3a  A  A  a  a x√  √  x √  3a and armchair edge (b). (Figures Figure zigzag edge 0 a1.3: Schematic forW zz +a W zz +2a 0 (a), Wac+ 23a Wac+ 3a 2 from Dietl [35])  (a)  (b)  0  a  (b  We provide sample bandstructure plots in Figs. 1.4(a) (zigzag) and  Figure ture of graphene nanoribbons. The dotted lines show the bonds5.1: Lattice structure of graphene nanoribbons. The dotted (b) (armchair). From Fig. 1.4(a) we can see the lowest conduction and that have to be cut and the circles the atoms that have to be remo the circles the atoms that have to be removed in order to create highest valence band form two partially flat degenerate bands at zero energy, (a) an armchair nanoribbon, (b) a zigzag nanoribbon. bon, (b) a zigzag nanoribbon.  which turn out to be localized at the edges of the sample. [36] Unlike the topologically protected variety of edge modes encountered so far, these 5.2 Band structure studied with the D zigzag edge theoretically proposed by Fujita et al., [37] are not robust ucture studied with thestates, Dirac equation in the face of perturbations. Furthermore, no current will flow along the Wefor will deriveflat an bands. analytical expression for the band structure of tical expression for the band of armchair and zigzag edgesstructure as the velocity goes to zero perfectly nanoribbons usingthe thedispersion Dirac equation (see Eq. (2.23) in Section 2.2. irac equation (see Eq. (2.23) in following Ref. [28]. AsSection for the 2.2.2), armchair nanoribbon, Fig. 1.4(b), is gapped This expression is onlygapped valid for energies. The resulting band s valid for low energies. Thefor resulting band structure be comthis particular case, butwill alternates between being andlow metallic pared with results of numerical tight-binding calculations. These c merical tight-binding calculations. These calculations consist of depending on the width of the sample in the finite direction. The size of the exact diagonalisation of the tight-binding Hamiltonian and give of the tight-binding Hamiltonian and giveon thethe proper gap also depends widthdispersion ofanthe sample and reaches the zero-gap value relation. We will see that the analytical approach can not capture at the analytical approachofcan not capture effects. graphene for an all infinite ribbon.  ribbon  5.2.1  Armchair ribbon  8  We wave consider the armchair ribbon to be aligned along the y-direct air ribbon to be aligned along the y-direction, so that the function x-direction. To find the dispersion relation for a ribbon, we is restricted in x-direction. To find the dispersion relat replace qx → −i∂x in the Dirac Hamiltonian (Eq. (2.23) in Section e Dirac Hamiltonian (Eq. (2.23) in Section 2.2.2):        ψA 0 i∂x + iqy 0 0 ψA ψA + iqy 0 0  ψ   ψ  3a  i∂ − iq  ψ  0 0 0 0 0 0  1.5. Haldane model For a more in depth look at the physics of nanoribbons we suggest the review articles by Dubois et al., [38] and Wakabayashi et al.. [39] (a)  (b)  3  3  2  2  1 E t 0  1 E t 0  ï1  ï1  ï2  ï2  ï3 0  2π  π  kx  ï3 −π  kx  0  π  Figure 1.4: Sample bandstructure plots for graphene nanoribbons having a zigzag edge (a), and armchair edge (b). Both ribbons are 12 unit cells wide, with a single unit cell corresponding to the red lines seen if figures 1.4(a)and (b).  1.5  Haldane model  Theoretically, another important contribution to topological insulator lore was made by Duncan Haldane,[40] who asked whether it was possible to establish the quantum Hall effect in the absence of a net magnetic field. As a result he discovered that the breaking of T invariance was the essential ingredient to obtain a non-zero Hall conductance, while also providing a lattice model for a Chern insulator to explore, seen in Fig. 1.5. He proposed a tight-binding model on the honeycomb lattice with real nearest-neighbor hopping, identical to the model for graphene above, but with complex valued next-nearest neighbor hopping, responsible for breaking T symmetry, along with an on-site inversion symmetry breaking term. The Hamiltonian is given by c†i cj + t2  HHaldane = t1 i,j  e−iνij φ c†i cj + M i,j  ξi c†i ci .  (1.20)  i  ˆ1 × d ˆ 2 )z = ±1, where dˆ1 and dˆ2 are unit vectors along the with νij = sgn(d two bonds, which constitute the next-nearest neighbor hopping. ξi repre9  1.5. Haldane model  Figure 1.5: Four unit cells of the Haldane model. The arrows on the next nearest neighbor hopping terms indicate the direction for picking up a positive phase. (Figure from Thonhauser and Vanderbilt [41])  sents a staggered on-site potential, which takes the values ±1 depending on whether the i-th site belongs to the A or B sublattice, respectively. The result of the next nearest neighbor hopping is that going around the hexagonal plaquette in the clockwise direction results in the phase eiφ , whereas in the counterclockwise direction one gets (e−iφ ). Just as in the standard quantum Hall system, a non-zero chern number will emerge if the gap is dominated by the T symmetry breaking terms, seen in Fig. 1.6, and chiral edge modes will appear if the system size is finite.  Figure 1.6: Phase diagram for Haldane model. (Figure from Murakami [42])  10  1.6. Z2 topological insulator  1.6  Z2 topological insulator  As Haldane showed, the topologically non trivial states described above can occur only when T symmetry is broken. However, it turns out if the symmetry is restored, a new topological class can emerge, which was discovered in the groundbreaking work by Kane and Mele.[11, 12] They found that by doubling the Haldane model with the inclusion of spin, one could obtain an insulator with a robust pair of helical edge states. To understand the result, discussed below in section 1.7, it is useful to first examine how spin 1/2 particles particles behave under T . T symmetry is represented by an antiunitary operator Θ = exp(iπσy / )K, where σy is the standard Pauli matrix and K is complex conjugation. For spin 1/2 electrons, Θ has the property Θ2 = −1. This leads to an important constraint, known as Kramers’ theorem, that all eigenstates of a T invariant Hamiltonian are at least twofold degenerate. This follows because if a non degenerate state |χ existed then Θ|χ = c|χ for some constant c. This would mean Θ2 |χ = |c|2 |χ , which is not allowed because |c|2 = −1. In the absence of spin orbit interactions, Kramers’ degeneracy is simply the degeneracy between up and down spins which is fulfilled automatically in the absence of a magnetic field. In the presence of spin orbit interactions, however, it has nontrivial consequences. If there are electronic states bound to the edges (or surface in the case of a 3D topological insulator), then Kramers’ theorem requires they be at least twofold degenerate at the T invariant momenta. If there is only one state available for each spin (Fig. 1.7(a)), then opening a gap necessarily violates Kramers’ theorem(Fig. 1.7(b)). If there were two Kramer’s pairs of edge states however(Fig. 1.8(a)), every energy would be fourfold degenerate, and one could imagine transforming the state into a gapped configuration without spoiling the degeneracy at the T invariant momenta (Fig. 1.8(b)). Therefore the only way Kramers’ theorem will allow the band structure to be continuously deformed into a gapped configuration, is if there are an even number of Kramers’ pairs. Thus we have two topologically distinct configurations, and can define a new Z2 topological invariant, ν, having values 0 or 1 corresponding to the even and odd number of Kramer’s pairs, respectively. There are several mathematical formulations of the ❩2 invariant ν [5, 12, 20, 44–49]. However, if the crystal has inversion symmetry there is a shortcut to computing ν [20]. One simply has to calculate the parity 11  k  k  k=0  k=0  (a)  (b)  Figure 2.2: A gapless (a) and a gapped (b) configuration of edge bands in the gap between the bulk bands, with two Kramers’ pairs. Red and blue colours denote opposite spins. 1.6. Z2 topological insulator (a)  (b)  E  E  EF  EF  k  k  k=0  k=0  Figure 1.7: A gapless (a) and a forbidden, gapped, (b) configuration with (b) opposite a12single Kramers (a) pair on the edge. RedChapter and blue spins. 2 denote Topological insulators (Figures from Str¨ om [43]) Figure 2.3: A gapless (a) and a forbidden, gapped, (b) configuration with one (b) blue colours denote opposite (a)Kramers’ pair on the edge. Red and spins. E E  uration, one with an even number ofE Kramers’ pairs, where a gap can E form, and one with an odd number of pairs, where one pair will always be left gapless. Which of these two alternatives that occurs depends on the topological class of the bulk structure. There is thus a Z2 topological ink variant ν, defined modulo 2, whichk tells us if the system is a topologically k=0 k=0 trivial insulator or a topological insulator with conducting edge states. Figure 1.8: A gapless (a) and a gapped (b) configuration of edge states (a) the integer quantum Hall case, the (b) Chern number Comparing with between the bulk bands, with two Kramers pairs. Red and blue denote is n = 0 in these systems, but we now have a new topological invariant opposite spins. (Figures(a) from m [43]) (b) configuration of edge bands Figure 2.2: A gapless andStr¨ a ogapped in the gap between the bulk bands, with two Kramers’ pairs. Red and blue coloursξdenote opposite spins. eigenvalues, m (Λa ) = ±1, at the four special time-reversal invariant points F  F  Λa in the bulk 2D Brillouin zone. The ❩2 invariant is then given by the expression δa = ξm (Λa ), (1.21) E  E  m  where the product is over the Kramers pairs of occupied bands. This has proven very useful for identifying topological insulators from band structure calculations [20, 28, 50, 51]. This is alsoE the method we shall use in Chapter E 2, along with its generalization to three dimensions utilizing the 8 timereversal invariant points in the 3D Brillouin zone. F  F  k  k  k=0  k=0  (a)  (b)  12  Figure 2.3: A gapless (a) and a forbidden, gapped, (b) configuration with one Kramers’ pair on the edge. Red and blue colours denote opposite spins.  1.7. Kane-Mele model  1.7  Kane-Mele model  The Hamiltonian for the original Kane-Mele model on the honeycomb lattice looks as follows, [12]  ij  (1.22)  ij  ij  +iλR  νij c†i σ z cj  c†i cj + iλSO  H=t  ˆ ij )z cj + λv c†i (s × d  ξi c†i ci . i  The first two terms double the Haldane model, as advertised above, where the spin up and spin down components correspond respectively to the values φ = −π/2 and φ = π/2 in Eqn. 1.20, and σ z is a Pauli matrix describing the electron’s spin. The third term is the nearest neighbor Rashba term, which will arise if a perpendicular electric field is applied or upon interaction with a substrate. Lastly, as in Eqn. 1.20, a staggered sublattice potential is included to describe the transition between the QSH phase and the trivial insulator. (a)  (b)  2  2  E t 0  E t 0  ï2  ï2 Γ  M  K  Γ  Γ  M  K  Γ  Figure 1.9: (a) Dispersion resulting from Eqn. 1.22 for λSO = 0.04t . (b) √ Dispersion for λR = 2 3λSO . √ For λR =√0, an energy gap of magnitude |6 3λSO − 2λv | will√arise, and so long as 3 3λSO > λv , the QSH phase survives. For λv > 3 3λSO the gap will be dominated by λv , and the system becomes a trivial√insulator. Although the Rashba term violates σ z conservation, for λR < 2 3λSO the resulting gapped phase will be adiabatically connected to the QSH phase at λR = 0. In Fig. 1.9(a) we show the bulk gap for an infinite system aris13  1.8. HgTe quantum wells ing when λSO = 0.04t, √ and in Fig. 1.9(b), we show how the gap vanishes at the point λR = 2 3λSO . Lastly we solve the lattice model in an armchair strip geometry 10 unit cells wide for λSO = 0.04t (seen in Fig. 1.10), showing clearly the gapless edge states resulting from the bulk boundary correspondence. As the spin-orbit coupling terms in Eq. 1.9 are so widely used, it is also worth mentioning how to connect their tight binding form to the standard spin-orbit coupling terms arising from the non-relativistic limit of the Dirac equation, 1 Hso = (∇V0 × p) · S + V (r), (1.23) 2(mc)2 where V0 is the periodic potential from the underlying lattice, m is the mass of an electron, p is the momentum of the electron, S is the spin operator and V (r) is an external applied potential (e.g. gate voltage). Treating these terms as a perturbation around the points K and K in the brillouin zone (Fig. 1.1(b)) using second-order degenerate state perturbation theory, the effective spin orbit interaction, Hef f = λso σz τz sz + λR (σx τz sy − σy sx )  (1.24)  then follows, where the Pauli matrices, σ, act in the A/B sublattice space, τz describe states at the K and K points and sz represents the electron’s spin. If one expands the tight-binding Hamiltonian (Eq. 1.22) around K and K , it will reduce to Eq. 1.24, and we have the desired connection to Eq. 1.23. A more in depth account of this connection can be found in the following references: [52–54]. In Chapter 2, we apply the Hamiltonian above (Eq. 1.22) to the Lieb lattice, and its 3D counterpart, the perovskite, and are able to show that a topologically non trivial gap emerges, as well as the gapless edge modes for a finite strip of the lattice.  1.8  HgTe quantum wells  Although graphene has been shown to host a QSH insulating phase, as mentioned in the overview carbon is a light element with very weak spin orbit coupling. [13–17] So, a more promising place to look for this physics would be in materials having strong spin orbit interactions, made from heavy elements near the bottom of the periodic table. This was the brilliant approach 14  1.8. HgTe quantum wells  3 2 1 E t 0 ï1 ï2 ï3 −π  kx  0  π  Figure 1.10: Armchair nanoribbon with λSO = 0.04t showing gapless edge modes.  taken by Bernevig, Hughes and Zhang (BHZ) who considered quantum well structures where HgTe is sandwiched between layers of CdTe, and in turn paved the way to the experimental discovery of the QSH phase. We will keep things brief, but an excellent review of the physics of the QSH effect in HgTe/CdTe quantum wells can be found in the article by K¨ onig et al. [55]. In CdTe (Fig. 1.11(a)) the valence band is the so called Γ8 band with J = 3/2, degenerate at k = 0 with another Γ8 band having J = 1/2. Above the Fermi level lies the conductance band, Γ6 . In HgTe (Fig. 1.11(b)) on the other hand, it turns out the spin-orbit interaction is large enough to push the Γ8 bands above the Γ6 band, yielding an inverted band structure. The upper, J = 3/2, Γ8 band becomes the conduction band, and the material becomes a zero gap semiconductor because of the degeneracy at k = 0, with the J = 1/2 band becoming the valence band. BHZ found that when the thickness of the HgTe layer is d < dc = 6.3 nm the 2D electronic states bound to the quantum well have the band order described above for CdTe. For d > dc , however, the 2D bands invert back to their natural state, and it was found that this inversion of the bands as a function of increasing d signals a quantum phase transition between the trivial insulator and the quantum spin Hall insulator. This theory was tested experimentally by K¨onig et al. [19], by measuring the electrical conductance due to the edge states of a HgTe quantum well of thickness d sandwiched between Hg0.3 Cd0.7 T e layers (see Fig. 1.12(a)). In 15  1.9. 3D topological insulators 14  Chapter 2 Topological insulatorsspin Hall effect in HgTe quantum wells 2.4 The quantum  (a)  15  (b)  E  E Γ6  EF  EF  Γ8  Γ8 Γ6 Γ7  Γ7 k k=0  k k=0  Figure (a)structure The band CdTe near k = 0. (b) The inverted Figure 2.4:1.11: The band of CdTestructure near k =Figure 0. of 2.5: The inverted band structure of HgTe near k = 0. band structure of HgTe near k = 0. (Figures from Str¨om [19])  (2 + 1)-dimensional topological insulators, or QSH insulators. thewith J = 1/2, split from the J = 3/2 band by spinis also a Here, Γ8 band characteristic topological invariant previously arises directly theZ2experiment, thediscussed thickness oforbit the HgTe layermeans wasthat varied between d= 5.5 except at coupling, which this bands has a lower energy, from the field theory [26]. = 0, where the two Γ8 bands are degenerate. Above the Fermi level lies nm and d = 12 nm, so that kthe the critical thickness d would be covered. conductance band, called Γ6 . cThis band structure is sketched in fig. Fig. 1.12(b) shows the resistance a the series of samples a enough 2.4. measurements In HgTe on the otherfor hand, spin-orbit interactionas is large 2.4 The quantum spin Hall effect in theHgTe to push Γ8 bands aboveenergy the Γ6 band, thus shifting the valence and function of a gate voltage which tunes the Fermi through the bulk conductance bands one step, as shown in fig. 2.5. The conductance band quantum wells energy gap. Sample I is a narrow with d <anddcsoshowing is nowquantum the upper, Jwell = 3/2, Γ8 band, the valence insuand conductance bands are degenerate k = 0. No gap can For a goodlating review ofbehaviour the physics ofwith the QSH effect inresistance HgTe/CdTe quana large in the at gap. Samples II,now IIIform andwithout IV, breaking and hence also time-reversal symmetry. tum wells (QW:s), including the earliest experiments,this see degeneracy the article by however, are wells the inverted III and both This theoryregime. was testedSamples experimentally by K¨oIV nig et al. [7]. They K¨onig et al. [27]. We will herewider sketch only thein basics. used HgTe with quantum well, where the HgTe well, with thickness d, was The idea of usinga HgTe QW:s as QSH was apresented by the exhibit conductance 2e2systems /h, consistent Landauer-Buttiker framesandwiched between Hg Bernevig, Hughes and Zhang in 2006 [28]. These authors were looking 0.3 Cd0.7 Te layers (see figure 2.6). The idea was which predicts a inquantized conductance e2 /h associated each set band gap, that the QW should behave more like HgTe,with with the inverted for a band work structure of the type depicted fig. 2.3a. As shown by Kane for aof thick layeraand morethickness) like CdTe, with normal and Mele [6], could beThe a candidate for thisboth type bandHgTe strucof graphene edge states. fact that samples (at fixed arethethe theband gap, when HgTe is thinner. In ref. [27], the critical thickness, where the ture and indeed, the degeneracy at the Dirac point gives thethe appropriate same length L = 1µm but different widths = material 0.5µmpushes and the 1µm, indicates SO coupling of thewQW Γ8 band above the Γ6 band, helical edge states, protected by time-reversal symmetry. The problem was calculated to be dc = 6.3 nm. In the experiment, the thickness of thatis that the the transport must be light along the edge. with graphene carbon atoms are too to produce a spinthe HgTe was varied between d = 5.5 nm and d = 12 nm, so that orbit coupling strong enough to give a large enough energy gaplayer for the It was a huge achievement establishing that the HgTe quantum wells are the critical thickness dc should be covered. The fundamental quantum bulk states. indeed topological insulators, but the fact that single of conductance is e2 /h.aSince the group edge of aworldwide QSH system carry two Instead, Bernevig et al. started to look for a heavier material, withG =only channels of dissipationless current, the Hall conductance should be twice the same Kramer’s degeneracy at k = 0, i.e. the Γ point in the Brillouin are capable of carrying out the experiments shows that there is definitely value, σxy = 2e2 /h. A four-terminal conductance measurement zone. Let us have a look at the energy bands near the this Fermi leveli.e. in two room for a less complicated alternative. So this section can to also serve as a was In performed was possible detect. heavy semiconductor materials, namely CdTe and HgTe. CdTe (astoinsee if this 2 What they found was a quantised Landauer conductance GaAs), themotivation valence band isfor the the so called Γ band with J = 3/2. There work8 carried out in Chapter 5 where we have shown, at of 2e /h,  least theoretically, that it is possible to engineer a much simpler quantum spin Hall state by depositing adatoms on the surface of graphene.  1.9  3D topological insulators  As stated in the overview, a 3D topological insulator is characterized by four ❩2 topological invariants (ν0 ; ν1 ν2 ν3 ), which can be most easily understood 16  1.9. 3D topological insulators 16  Chapter 2 Topological insulators  (a)  (b) 2  20  G = 0.01 e /h  cap layer  Barrier  Doping layer Barrier  I  Vg R14,23 / Ω  Insulator  10  6  T = 30 mK  Doping layer  T = 0.03 K 2  G = 2 e /h  10 T = 1.8 K 5  II 0 -1.0  10  -0.5  5  0.0 0.5 (Vg - Vthr) / V  1.0  2  HgTe QW Barrier  III  15  R14,23 / kΩ  10  7  III  10  4  10  3  G = 0.3 e /h 2  G = 2 e /h  IV  Barrier Buffer  -1.0  Substrate  -0.5  0.0  0.5  1.0  1.5  2.0  (Vg- Vthr) / V  Figure 2.6: A simple sketch of a HgTe QW. For a more detailed figure, see ref. [27].  Figure 1.12: (a)Schematic of the layer sequence of the MBE-grown quantum structures experiment (Figure from Str¨om [43]) (b)Experimental as predicted, well for d > dc , while the used QW wasin shown to be an ordinary insulator for two d < dcterminal . The edges of a QSH system carry conductance as toa channels function of a gate voltage that tunes EF of electrons, so the conductance value This predicted conductance was through bulksamples gap. (∼Sample with d < dc , shows insulating behavior, found in samples of size ∼ 1the µm. Larger 20 µm) alsoI, showed the behaviour of a conventional insulator for d < dc and a finite conwhile samples III and IV show quantized transport associated with edge ductance for d > dc , even though this conductance was lower than the predicted value of 2e2 /h found in the from smaller samples. states. (Figure K¨ onig [19]) Importantly, they also found that breaking time-reversal symmetry with even a small magnetic field destroyed the conductivity, again making the QW an ordinary insulator.  by studying the surface states. These are the higher dimensional analogue of the edge states described above. There will be four T invariant points Γ1,2,3,4 in the surface Brillouin zone where surface states, if present, must be Kramers degenerate [Figs. 1.13(a) and (b)]. Away from these special points, the spin-orbit coupling lifts the degeneracy and 2D Dirac points then form in the vicinity (Fig. 1.13(c)). From here one can then determine whether the surface Fermi surface intersects a line joining any pair Γa to Γb an even or an odd number of times. If it is odd, then the surface states are topologically protected. Which of these two alternatives occurs is determined by the four bulk ❩2 invariants. The simplest non trivial 3D topological insulators can be constructed by stacking layers of the 2D quantum spin Hall insulator, and a possible surface Fermi surface for weakly coupled layers stacked along the y direction is sketched in Fig. 1.13(a). This layered state is referred to as a weak topological insulator, and has ν0 = 0, with the other indices (ν1 ν2 ν3 ) serving as Miller type indices describing the orientation of the layers. Unlike the 2D helical edge states of a single layer, T symmetry does not protect the 17  1.9. 3D topological insulators  ky Γ3  ky Γ3  Γ4 kx  Γ1 (a)  EF  Γ4 kx  Γ1  Γ2 (b)  E  Γ2 (c)  Figure 1.13: Fermi circles in the surface Brillouin zone for a weak topological insulator (a) and a strong topological insulator (b). (c) Dirac cone. (Figure from Hasan and Kane [8] )  surface states in this case, and although they must be present for a clean surface, they can be localized in the presence of disorder. ν0 = 1 identifies a distinct phase, called a strong topological insulator, which can not be interpreted as a descendent of the 2D quantum spin Hall insulator. ν0 determines whether an even or an odd number of Kramers points is enclosed by the surface Fermi circle. In a strong topological insulator the surface Fermi circle encloses an odd number of Kramers degenerate Dirac points. The surface states of a strong topological insulator form a unique 2D topological metal [10, 20] that is essentially half an ordinary metal. Unlike an ordinary metal, which has up and down spins at every point on the Fermi surface, the surface states are not spin degenerate. Since T symmetry requires that states at momenta k and −k have opposite spin, the spin must rotate with k around the Fermi surface, as indicated in Fig. 1.13(b). This leads to a non trivial Berry phase acquired by an electron going around the Fermi circle. T symmetry requires that this phase be 0 or π. When an electron circles a Dirac point, its spin rotates by 2π, which leads to a π Berry phase. The Berry phase has important consequences for the behavior in a magnetic field and for the effects of disorder. In particular, in an ordinary 2D electron gas the electrical conductivity decreases with decreasing temperature, reflecting the tendency towards Anderson localization in the presence of disorder [56]. The π Berry phase changes the sign of the weak localization correction to the conductivity leading to weak antilocalization [57]. In  18  1.9. 3D topological insulators fact, the electrons at the surface of a strong topological insulator can not be localized even for strong disorder, as long as the bulk energy gap remains intact [58].  1.9.1  The first 3D topological insulator: Bi1−x Sbx  Although the 3D topological insulating model we introduce in Chapter 2 and explore further in Chapter 3 can not currently be realized in the lab, there exists an impressive body of experimental work on other materials. After first being predicted by Fu and Kane, [10] the first 3D topological insulator to be identified experimentally was the semiconducting alloy Bi1−x Sbx , whose surface bands were mapped in an angle resolved photoemission spectroscopy (ARPES) by Hsieh et al. [21] The charge transport experiments, which were successful in identifying HgTe as a 2D topological insulator [19], are problematic in 3D materials because it is difficult to separate the surface contribution to the conductivity from that of the bulk. Angle resolved photoemission spectroscopy (ARPES), however, which uses a photon to eject an electron from a crystal and then determines the surface or bulk electronic structure from an analysis of the momentum of the emitted electron, is an ideal tool for probing the topological character of the surface states. High-resolution ARPES performed with modulated photon energy allows for a clear isolation of surface states from that of the bulk 3D band-structure because surface states do not disperse along a direction perpendicular to the surface where as the bulk states do. Moreover, unlike in a transport experiment, ARPES carried out in a spin resolution mode can, in addition, measure the distribution of spin orientations on the Fermi surface which can be used to estimate the Berry phase on the surface. Spin sensitivity is critically important for probing the existence of spin-momentum locking on the surface expected as a consequence of bulk topological order. The experiments by Hasan et al. probed both the bulk and surface electronic structure of Bi.09 Sb.91 with ARPES. Fig. 1.14(a) shows the ARPES spectrum, which can be interpreted as a map of the energy of the occupied ¯ to electronic states as a function of momentum along the line connecting Γ ¯ in the projected surface Brillouin zone (Fig. 1.14(b)). Bulk energy bands M associated with the L point are observed that reflect the nearly linear 3D Dirac like dispersion. Kramers’ theorem requires surface states to be doubly degenerate at 19  1.9. 3D topological insulators  EB(eV)  0.1  3D Topo. Insulator (Bi1-xSbx) 1  bulk gap  3 4,5  2  0.0  -0.1  0.0  (a)  0.2  0.4  0.6 -1  G (KP)  - kX (Å )  0.8  1.0  M (KP)  z  K  M  ky  G kx  T X  L  E T  L X  La 4% 7% 8%  8  x=0 x=0.1  6 4  ´ 80  2  LS  Bi  2 1  (c)  r (mW cm)  (b)  x  0 0  100  200  300  T (K)  Figure 1.14: Topological surface states in Bi1−x Sbx : (a) ARPES data on the 111 surface of Bi0.9 Sb0.1 which probes the occupied surface states as a ¯ and function of momentum on the line connecting the T invariant points Γ ¯ in the surface Brillouin zone. (b) Schematic of the 3D Brillouin zone and M its (111) surface projection. (c) Comparing the resistivity of semimetallic pure Bi with the semiconducting alloy. (Figure from Hasan and Kane [8])  20  1.9. 3D topological insulators ¯ and each of the three equivalent M ¯ points. Such a the T invariant points Γ ¯ approximately 15±5meV below EF . Kramers point is indeed observed at M As expected for a system with strong spin orbit interactions, the degeneracy ¯ . The observed surface bands cross the Fermi energy 5 is lifted away from M ¯ and M ¯ indicating that these surface states are topologically times between Γ protected. Accounting for the threefold rotational symmetry and mirror symmetry of the 111 surface, this data shows that the surface Fermi surface ¯ an odd number of times, while it encloses the three equivalent encloses Γ ¯ points an even number of times. This establishes Bi1−x Sbx as a strong M topological insulator, with ν0 = 1. The data is consistent with the predicted (1; 111) topological class. Establishing the first topological insulator in 3D was a remarkable achievement, but, as seen in the Fig. 1.14(a), the surface structure of Bi1−x Sbx is quite complicated, and the band gap rather small, and so a search began for new materials with a larger bandgap and a simpler surface spectrum. This led to the discovery of a second generation of 3D topological insulator materials, with the most promising being Bi2 Se3 (shown in Fig. 1.15), having a single Dirac cone and topological behaviour visible at room temperature. We have just scratched the surface here, but, again, encourage anyone wishing for greater detail to consult the excellent review by Hasan and Kane. [21]  21  1.9. 3D topological insulators  νi={1,000}  Bi2Se3 Γ  M  K  Γ  K  EB (eV)  0.0  M  Bi2Se3  0.4 0.2 (a) -0.2 k 0.0 -1 x (Å ) 0.2  0.0  -1  ky (Å )  0.2  (c)  k y (Å -1 )  0.1  ky  E kx  0.0  -0.1  -0.2 -0.2 -0.1 0.0 0.1 0.2 (b) -1 kx (Å )  (d)  Figure 1.15: (a) ARPESBidata for Bi2 Se3 showing surface electronic states 2Te 3 with a single spin-polarized Dirac cone. (b) The Surface Fermi surface with chiral left-handed spin texture. (c) Surface electronic structure of Bi2 Se3 computed in the local density approximation. (d) Schematic of the spin (b) polarized surface state dispersion in Bi2 X3 (1; 000) topological insulators. (Figure from Hasan and Kane [8])  22  Chapter 2  Topological insulator on the Lieb and perovskite lattices 2.1  Introduction  In this chapter, we present another possible host for a topological insulator, namely the Lieb lattice in 2D and its 3D counterpart the perovskite or edge centered cubic lattice. We show that a simple tight-binding model for electrons on these lattices leads to TI behavior in the presence of spin-orbit coupling. These toy models are novel in that unlike most of the existing models, they have a simple cubic symmetry, and exhibit only a single Dirac node in the spectrum intersecting a degenerate flat band precisely at the Dirac point. We also briefly discuss other perturbations that lead to topologically trivial gapped phases.  2.2  Lieb lattice  We begin with the tight-binding Hamiltonian H0 = −t  c†iα cjα + h.c,  (2.1)  ij α  where c†iα creates an electron of spin α on site i of the so called Lieb lattice, seen in Fig. 2.1(a), and t is the hopping amplitude for nearest neighbour sites. † 0 In momentum space Eq. (2.1) becomes H0 = kσ Ψkσ Hk Ψkσ with Ψkσ = (c1kσ , c2kσ , c3kσ )T and    0 cos(kx ) cos(ky ) . Hk0 = −2t  cos(kx ) 0 0 cos(ky ) 0 0 23  2.2. Lieb lattice (a)  (b)  a2 1  2  a1  3  Figure 2.1: (a) Lieb lattice showing 3-site basis in unit cell, NNN hopping (red dotted lines) and basis vectors a1 and a2 . (b) Tight-binding dispersion for Hk0 . (3)  The spectrum of Hk0 consists of one degenerate flat band Ek = 0 and two dispersive bands (1,2)  Ek  = ±2t  cos2 kx + cos2 ky ,  (2.2)  where the Brillouin zone spans − π2 ≤ kx ≤ π2 , − π2 ≤ ky ≤ π2 . Plotting these bands yields a low energy Dirac like dispersion (seen in Fig. 2.1(b)), with a single Dirac cone in the Brillouin zone. At one third filling then, band 1 will be completely filled, with the degenerate flat band and the upper band being empty, and this state will behave as a gapless band insulator. There are several perturbations to the basic Hamiltonian (2.1) that open up a gap while respecting the translational symmetry of the lattice. The simplest option, recently discussed in the context of ultracold fermionic atoms [59] by varying laser intensities along different directions, is to shift the on site energy of site 1 in Fig. 2.1(a), relative to sites 2 and 3. In this case, one of the dispersing bands remains touching the flat band while the other dispersing band becomes isolated, as seen in Fig. 2.2(a). Another option is to add a dimerization term that staggers the hopping amplitude along both the x ˆ and yˆ directions so that t → t±α. The additional terms, once Fourier transformed, look as follows    0 i sin(kx ) i sin(ky ) . Hk0 = 2α  −i sin(kx ) 0 0 −i sin(ky ) 0 0 24  2.2. Lieb lattice and the flat band becomes isolated from the dispersing bands, which will be gapped symmetrically (seen in Fig. 2.2(b)). One possible way to create such a term would be to include interactions in the model. It has been shown, at the mean field level, that both the square and honeycomb lattices are unstable towards the formation of a dimerization phase when including short range Hubbard terms. [60] Lastly, including the Rashba spin orbit term HR = iλR ij αβ  c†iα (σ × dij )z cjβ ,  (2.3)  where dij is the unit vector along the nearest neighbour bonds connecting site i to its nearest neighbour j and σ is the vector of Pauli spin matrices, results in an isolated flat band while also splitting the spin degeneracy among the valence and conduction bands, which can be seen in Fig. 2.2(c). In this case, the spin mixing results in a Hamiltonian with off diagonal terms, and so the top right block of the full 6 × 6 Hamiltonian reads    0 i sin(kx ) − sin(ky ) . HkR = 2λR  i sin(kx ) 0 0 − sin(ky ) 0 0  with the bottom left block equal to the hermitian conjugate. (a)  (b)  (c)  3 2 1  E 0 t ï1 ï2 ï3 Γ  X  M  Γ  Γ  X  M  Γ  Γ  X  M  Γ  Figure 2.2: Bandstructure examples for topologically trivial insulators. (a) On site energy = −0.5t for site 1 seen in Fig. 2.1(a) (b)Dimerization α = 0.2t (c)Rashba spin orbit interaction, λR = 0.2t In all of the above cases the resulting insulating phases are characterized by conventional broken symmetries and are topologically trivial. The recent paper by Green, Santos and Chamon[61] discusses a similar 25  2.2. Lieb lattice set of energy bands having a flat band directly through the middle of a single Dirac cone, which have been created via staggered flux phases on the kagome lattice. Although their interests lie with the flat band itself, they discuss possible insulating phases that can be introduced in order to isolate the flat band and determine that, for spinless fermions, time reversal symmetry must be broken in order to achieve this. Seeing as the Lieb lattice requires no magnetic field to create the same energy band structure, perhaps it is of some interest then that we can isolate the flat band without breaking time reversal symmetry by way of the dimerization or Rashba terms mentioned above. Now, to see if this model can support a TI phase we include the intrinsic SO interaction term HSO = iλ ij αβ  (d1ij × d2ij ) · σαβ c†iα cjβ ,  (2.4)  where λ is the amplitude for the next-nearest neighbour spin-orbit-induced interaction (shown in Fig. 2.1(a)), the term νij =(d1ij × d2ij )z =±1, where d1ij and d2ij are again the two unit vectors along the nearest neighbour bonds connecting site i to its next-nearest neighbour j with σ the vector of Pauli spin matrices. This will also lead to the formation of a gap at the Dirac point while preserving T and the translational symmetry of H0 . Fourier transforming, we have HkSO    0 0 0 = ±4λ 0 0 i sin (kx ) sin (ky ) , 0 −i sin (kx ) sin (ky ) 0  (2.5)  where the +(−) sign refers to spin up (down) electrons, and the spectrum for the full Hamiltonian, Hk = Hk0 + HkSO , then consists of the same degenerate (3) flat band Ek = 0, and the modified doubly degenerate dispersive bands (1,2)  Ek  = ±2  t2 (cos2 kx + cos2 ky ) + 4λ2 sin2 kx sin2 ky ,  (2.6)  with a gap ∆SO = 4|λ| at the Dirac point. We note that the flat band eliminates the possibility of writing down a low energy Dirac expansion for the system in terms of the Pauli matrices as discussed e.g. in Ref. [22] and other topological insulator models, as it intersects the Dirac cone. However, following Ref. [61], it is possible to 26  2.2. Lieb lattice express the full Hamiltonian above in terms of a set of 3 × 3 matrices that form a spin-1 representation of SU (2). Although not a focus here, this makes it possible, for instance, to calculate the Chern numbers of the bands analytically. To prove that our system is indeed a topological insulator, we now show by an explicit calculation that it possesses a nontrivial Z2 invariant. There are several ways to do this in practice, as stated in Section 1.6, although they are found to be equivalent in the end, so we use the method that is most convenient. According to Fu and Kane[20] when a crystal possesses inversion symmetry the Z2 topological invariant ν is related to the parity eigenvalues ξ2m (Γi ) of the 2m-th occupied energy band at the four T -invariant momenta Γi . Our system is inversion symmetric and so we can use this method to find ν. If we select site 1 of the unit cell as the center of inversion then the parity operator acts as P[ψ1 (r), ψ2 (r), ψ3 (r)] = [ψ1 (−r), ψ2 (−r − a1 ), ψ3 (−r − a2 )]  (2.7)  on the triad of the electron wavefunctions in the unit cell labeled by vector r. In momentum space the parity operator becomes a diagonal 3 × 3 matrix    1 0 0 Pk = 0 e−ia1 ·k 0 , −ia 0 0 e 2 ·k  (2.8)  Γi = π(ˆ xni + yˆmi )/2  (2.9)  and the four T -invariant momenta can be expressed as  with ni , mi = 0, 1. The eigenstates of HΓi can be found analytically in this case and it is then easy to determine the parity eigenvalues of the occupied bands. At 13 filling we find that three ξ’s are positive and one is negative. Which of the four ξ’s is negative depends on the choice of the inversion center but the product Πi ξ(Γi ) = (−1)ν is independent of this choice and determines the non-trivial Z2 invariant ν = 1, confirming that the system is indeed a topological insulator. Similar considerations for 23 filling also yield ν = 1. Next, we solve the model in a strip geometry numerically using exact diagonalization, in order to incorporate edge effects and demonstrate the 27  2.3. Perovskite lattice bulk-boundary correspondence, which states that whenever ν = 1 there will be a pair of topologically protected gapless modes along each edge in the system. For the trivial insulating state with α = 0.2t, we find no edge states (Fig. 2.3(a)), but for λ = 0.2t, we indeed find a pair of spin-filtered gapless states associated with each edge, as seen in Fig. 2.3(b). As there are 3 atoms in the unit cell, the two edges of the slab are not equivalent, and therefore the edge states are not degenerate. These edge states are not entirely localized, although the wavefunction will decay rapidly as can be seen Fig. 2.4, borrowed from the paper by Goldman et al. [84] (a)  (b)  2  2  E t 0  E t 0  −2  −2 0  π  π /2  kx  0  π /2  π  kx  Figure 2.3: (a) Bandstructure for strip of Ny = 10 unit cells with open boundary conditions along the y direction and infinite along x with α = 0.2t. (b) As in (a), but with λ = 0.2t instead.  2.3  Perovskite lattice  The perovskite lattice, depicted in Fig. 2.5, is a straightforward generalization of the Lieb lattice into 3D and can be viewed as a simple cubic lattice with additional sites positioned at the centers of all edges. As with the Lieb lattice, our starting point is the tight binding model given in Eq. (3.1), which in momentum space gives    0 cos kx cos ky cos kz  cos kx 0 0 0  . Hk0 = −2t   cos ky 0 0 0  cos kz 0 0 0  (2.10)  28  5  tions, the Z2 r nσ [14, 15], where nσ = numbers ashe numerical = 1 in agree-  bove by coni.e. a finite by a unique or a Lieb latedges. This g in the range CNP, within e eigenvalues at the edge of 5b, where the lar edge-state h |ψ↓ (x, y)|2 . lical property sociated cur-  (8)  n) �  n) − h.c , (9)  n) �  n) − h.c .  , ψσ (m, n, B))  (10a)  (a)  2  (b)  2.3. Perovskite lattice  | ! (x,y)| !  1 E(ky ) [t]  ant characteris defined on n numbers inwing Ref. [18], the inversion at SOI opens ν = 1, thereum spin-Hall  0  -1  -2  0  ! ky  2!  y [ 0]  x [ 0]  0  Figure 2.4: Edge-states amplitude |ψ↑ (x, y)|2 for the Lieb lattice with 39×39 sites with (a): λ = Energy 0.5t, corresponding energy E = 0.5t which lies within FIG. 5: Panel spectrum Eto=an E(k y ) in the cylinthe gap. (Figure from Goldman et al. [84]) der geometry with tSO = 0.2t. The bulk gap is traversed by  gapless edge-states [i.e. helical states]. Panel (b): Edge-states amplitude |ψ↑ (x, y)|2 for the open Lieb lattice with 39 × 39 (3,4) sites, straight edges then and consists tSO = 0.5t. This localizedflat eigenstate The spectrum of two degenerate bands Ek = 0 and the satisfies Hψ (x, y) = E ψ (x, y) and corresponds to the enλ λ λ two dispersive bands ergy Eλ = 0.5t which lies within the gap. (1,2)  Ek  = ±2t cos2 kx + cos2 ky + cos2 kz .  (2.11)  Again, we have a single Dirac point in the Brillouin zone and at 14 filling the III. LOW-ENERGY FERMIONS: THE system behaves as a gapless band insulator. The spectrum is shown in Fig. QUASIRELATIVISTIC REGIME 2.6(b). As above, we can create a topologically trivial insulating phase by including an anisotropic on site term or a dimerization (Figs. 2.7(a)/(b)). In thisToSection, focus on the the properties of the spin-orbit couopen upwe a topologically non-trivial gap, we consider 1,2 Lieb lattice non-interacting fermions at low energy,space, the Hamilpling Eq.for (2.4). Since the dij now lie in three-dimensional i.e. close to the CNP. We perform a long wavelength tonian does not decouple for the two spin projections asapabove, and instead proximation to the Schr¨ o dinger equation underlying the becomes an 8 × 8 matrix in k-space. The 4 × 4 blocks on the diagonal, have TB Hamiltonian the terms seenwhich in theconsists 2D case, in expressing the spatial  part of the wave function as the product of a fast-varying   0 0 0 part times a slow-varying part. 0Within this approxima 0 SO tion, the H wave function as i sin (kx ) sin (ky ) 0   0 can be written . k,1 = ±4λ   0 −i sin (kx ) sin (ky ) Ψα (Rα0) ∝ eik·Rα0Fα (Rα ) ,  0 0  0 (11) 0  (2.12)  the +(−) refers toαthe spinlattice up (down) The off diagonal blocks wherewhere α ∈ {A,B,H} and R is the siteblock. coordinate. are hermitian conjugates of one another, with Schr¨ the top right block having the We substitute this wave function into the odinger equation and expand the slow-varying part as  Fα (Rα� ±dj ) � Fα (Rα� )±dj ·∇r Fα (r)|r=Rα� +O(|d|2 ) . Collecting all terms, we are left with ˜ = vF Σ · p , H  (12)  where vF = 2�0 t is the Fermi velocity, and p =  29  2.3. Perovskite lattice  a3 4 3 1  a2 2  a1  Figure 2.5: Perovskite lattice showing 4 sites in unit cell along with basis vectors.  form    0 0 0 0  0 0 0 − sin (kx ) sin (kz )  SO . = 4λ  Hk,2  0 0 0 i sin (ky ) sin (kz )  0 sin (kx ) sin (kz ) −i sin (ky ) sin (kz ) 0 (2.13) Unlike the model on the Lieb lattice, the addition of the spin orbit term creates a dispersion in the degenerate flat bands and a gap that depends on the sign of λ, which is no longer symmetric around zero energy. Switching the sign of λ inverts the band structure. We now study the topological classes of these insulating phases. As above, the Z2 topological invariants (ν0 ; ν1 ν2 ν3 ) are easy to evaluate when a crystal possesses inversion symmetry and can be determined from knowledge of the parity eigenvalues ξ2m (Γi ) of the 2m-th occupied energy band at the 8 T -invariant momenta (TRIM) Γi that satisfy Γi = Γi + G. The 8 TRIM in our system can be expressed in terms of primitive reciprocal lattice vectors as Γi=(n1 n2 n3 ) = (n1 b1 + n2 b2 + n3 b3 )/2,  (2.14)  with nj = 0, 1. Then να is determined by the product  30  2.3. Perovskite lattice (a)  (b) kz A 2  L  H  E t 0  Γ  ky ï2  X  M  Γ  kx  A  L  H  Γ  Figure 2.6: (a) High symmetry points in the Brillouin zone. (b) Bandstructure inside bulk along path of high symmetry for Hk0 .  (−1)ν0 =  δn1 n2 n3 ,  (−1)νi=1,2,3 =  δn1 n2 n3 ,  (2.15)  nj=i =0,1;ni =1  nj =0,1  where  N  δi =  ξ2m (Γi ).  (2.16)  m=1  If we select site 1 of the unit cell, Fig. 2.5, as the center of inversion then the parity operator acts as P[ψ1 (r), ψ2 (r), ψ3 (r), ψ4 (r)] = [ψ1 (−r), ψ2 (−r−a1 ), ψ3 (−r−a2 ), ψ4 (−r−a3 )] (2.17) on the four-component electron wave function in the unit cell labeled by vector r. In momentum space and including spin the parity operator becomes a diagonal 8 × 8 matrix    1 0 0 0 0 e−ia1 ·k 0 0  ⊗ 1 0 Pk =  −ia ·k 0 0 e 2 0  0 1 −ia ·k 3 0 0 0 e  (2.18)  It is straightforward to obtain the eigenstates of HΓi and the parity eigenvalues of the occupied bands numerically, then determine the Z2 in31  2.3. Perovskite lattice variants. At quarter filling, we find that δ = −1 at the H point and δ = 1 at the other TRIM, so the spin-orbit phase is a (1; 111) strong topological insulator. An example of the bandstructure in the bulk for λ = −0.2t can be seen in Fig. 2.8(a). As in Ref.[23], we note that it is also possible to obtain a weak topological insulator (WTI) with this model. In the limiting case where the system is tuned so that electrons can only move in decoupled parallel planes, each forming a 2D Lieb lattice, we know from our work above that we have a collection of 2D topological insulators. Such a collection results in a WTI, which will survive even if the interplane coupling is restored. Using the FuKane method[20] we have confirmed that if hopping along the z direction is eliminated, the system at quarter filling becomes a (0; 001) WTI. For non-zero interplane coupling parametrized as t = ζt for the NN hopping and λ = ζλ for the spin-orbit interaction terms, where 0 ≤ ζ ≤ 1, we find that the WTI survives as long as ζ < 2λ/t. At ζ = 2λ/t the two lowest energy bands touch at the point M in the Brillouin zone and there is a sharp transition to the (1; 111) STI phase. We note that the system remains insulating at quarter filling for all other values of ζ. Incorporating edge effects into the system with a slab geometry and plotting the band energies along lines connecting the four surface TRIM (Fig. 2.8(b)), we can also see the bulk energy bands and an odd number of surface states which traverse the gap, a behavior characteristic of STI. (a)  (b)  2  2  E t 0  E t 0  ï2  ï2  Γ  A  L  H  Γ  Γ  A  L  H  Γ  Figure 2.7: (a) Bandstructure inside bulk along path of high symmetry for = −0.5t. (b) As in (a), but with α = 0.2t.  32  2.4. Conclusions (a)  (b)  2  2  E t 0  E t 0  ï2  ï2  Γ  A  L  H  Γ  Γ  X  M  Γ  Figure 2.8: (a) Bulk bandstructure shown along path connecting high symmetry points for λ = −0.2t. (b) Bandstructure for slab geometry along lines connecting four surface TRIM for λ = −0.2t.  2.4  Conclusions  We have established a model for a topological insulator on lattices with simple cubic symmetry. In 2D we have shown that a tight-binding model with spin-orbit coupling on the Lieb lattice supports gapless topologically protected edge modes, and possesses a non-trivial Z2 invariant. In 3D, we have demonstrated that a similar model on the perovskite lattice supports 2D gapless surface states when in the (1;111) STI phase, brought on by spin orbit coupling. We have also identified other gapped phases in these models that are topologically trivial. In addition to the topologically non-trivial band structures, the models considered here exhibit flat bands, which could give rise to interesting strongly correlated electron states even for relatively weak interaction strength.[62] A necessary condition for the non-trivial correlated physics we are interested in is that the single-electron states forming the flat band be delocalized in space (localized states typically lead to formation of a Wigner crystal). Delocalized single-electron states are guaranteed to exist if the flat band possesses a non-zero topological invariant, as happens e.g. in quantum Hall liquids. Unfortunately, this is not the case in models considered in this study. In the Lieb lattice, spin-orbit coupling preserves the flat band and separates it from other bands. The resulting flat band however turns out to be topologically trivial, i.e. it has ν = 0, while the TI behavior derives from the two dispersive bands which have ν = 1. In the perovskite lattice, spin33  2.4. Conclusions orbit coupling produces a significant dispersion in the original flat band and this will limit the importance of correlations in the system. Although we see no general reason that would prevent a flat band from having a non-trivial Z2 invariant, this situation does not seem to occur in commonly used simple models. Can our model be realized in a physical system? There exist many perovskites in nature as well as many layered perovskites composed of weakly coupled 2D planes with Lieb lattice structure. The most prominent example of the latter are the CuO2 planes in high-Tc cuprate superconductors such as YBa2 Cu3 O7 or Bi2 Sr2 CaCu2 O8 . In 3D, our model system is an idealization of the naturally occurring perovskites, most of which also have a heavy central atom in the middle of each cubic cell, such as SrTiO3 or the double perovskite structure Ba2 NaOsO6 . In these real materials the electron behavior near the Fermi level derives from the eg and t2g orbitals of the transition metal element occupying the cubic site while the edge sites are normally oxygens whose p-orbitals are far away from the Fermi level. The resulting band structure for the active orbitals is then significantly more complex than that captured in our simple tight binding model. Nevertheless, our model calculations demonstrate that this class of tight-binding Hamiltonians on lattices with cubic symmetry can support topological phases, both in 2D and 3D. We hope that our work will stimulate detailed band structure calculations of perovskites with heavy transition metal elements in search for new families of topological insulators. In the more exotic realm, it might be possible to artificially engineer the 2D system by modulating a two-dimensional electron gas with a periodic potential having Lieb symmetry, as achieved recently in constructing the ‘artificial graphene’.[63] Another possibility lies with cold Fermionic atoms in optical lattices as discussed in Refs. [26, 27].  34  Chapter 3  Flat bands with non-trivial topology in three dimensions 3.1  Introduction  When single-particle states of bosons or fermions comprise a flat band, the effect of interactions is generically non-perturbative and can lead to new states of quantum matter. For 2D electrons in a strong perpendicular magnetic field the flatness of Landau levels combined with non-zero Chern numbers leads to the spectacular fractional quantum Hall (FQH) effects[64, 65] and other exotic states such as the Wigner crystal[66] and various striped and bubble phases. Recently, models for spinless fermions moving in twodimensional lattices have been constructed that exhibit nearly-flat bands with non-zero Chern numbers, but remarkably in zero external magnetic field.[67–70] It has been proposed, on the basis of heuristic arguments[67, 68] and numerical simulations,[69, 71, 72] that at partial filling repulsive interactions will drive these systems into the zero-field FQH states. Although it is not clear apriori how to experimentally realize the above situation in a physical system, the prospect of engineering systems that can support fractional excitations without involving a magnetic field has led to significant interest recently in flat or nearly flat bands in various lattice models. In this chapter we explore whether it is possible to achieve a flat band with non-trivial topology in a system of electrons in three spatial dimensions. It is well known that the concept of the first Chern number (and the FQH physics) does not generalize to three dimensions. However, recent advances in the theory of time-reversal (T ) invariant band insulators have established a uniquely three-dimensional ❩2 -valued topological invariant[7–9] called ν0 which, as we argue below, can play a similar role. Taking a specific example of electrons moving in the 3D edge-centered cubic (perovskite) lattice we show that a suitable combination of short-range hoppings and spin-orbit coupling (SOC) terms can give rise to a nearly flat band with a non-trivial 35  3.1. Introduction  ❩2 index ν0 = 1. With the chemical potential inside the bandgap such crystal  would be a strong topological insulator with an odd number of topologically protected gapless states associated with all of its surfaces. At partial filling and in the presence of repulsive interactions we argue that the system may become a three-dimensional ‘fractional’ topological insulator.[24, 73–75] Many lattice models are known to support completely flat bands. In 2D these include p-orbitals in the honeycomb lattice,[76] as well as s-orbitals in the kagome, dice and Lieb lattices.[22, 77, 78] 3D examples include the pyrochlore and the perovskite lattice.[50, 77, 78] In their simplest form however the above tight-binding models do not constitute a platform suitable for the study of exotic correlated phases for at least two reasons. First, in all cases the flat bands touch other dispersing bands, usually in a quadratic fashion. The interactions are thus bound to mix states from adjacent bands and this is clearly detrimental to the emergence of correlated phases. Second, even when the Hamiltonian is modified to open a gap, e.g. by adding SOC as in Ref. [78], the resulting flat (or nearly-flat) bands are in all known cases topologically trivial. In such topologically trivial flat bands the eigenstates can be chosen as exponentially localized around lattice sites and, at partial filling, interactions will typically select a crystalline ground state with broken translational symmetry and not an exotic featureless liquid that underlies the FQH effect. As already argued in Refs. [67–69, 78] a natural place to seek such exotic phases is in a flat band that is topologically non-trivial. It is well known that in 2D a non-zero Chern number represents a topological obstruction to the formation of exponentially localized Wannier states.[79] It is precisely this obstruction that tilts the balance in favor of FQH states in partially filled Landau levels (although Wigner crystal phases are still known to arise at very low filling fractions). For ❩2 -odd phases of T -invariant band insulators in 2D and 3D a similar topological obstruction exists[80] if one insists on Wannier states that respect T . This consideration suggests that interacting electrons partially filling a flat band with a non-trivial ❩2 index in a 3D crystal could either (i) spontaneously break T and form a Wigner crystal, or (ii) remain T -invariant and form one of the proposed 3D fractional topological insulators characterized either by spin-charge separation[24] or by fractional magneto-electric effect controlled by fractional values of the axion parameter θ accompanied by ground state degeneracy. [73–75] In the following we lay groundwork for future investigations of these 3D exotic phases by constructing a simple tight-binding model whose spectrum 36  3.2. The model has a flat band characterized by ν0 = 1 and is separated from other bands by a large gap. We remark that constructing a flat band with non-trivial topology in a 3D crystal is a mathematically more constrained problem than in 2D and it is by no means obvious that this can be achieved with shortrange hoppings only.  3.2  The model  As a first step in the program outlined above we study a simple tight-binding model with intrinsic SOC terms on the same lattices seen above in Chapter 2, but with additional short range hoppings on the 2D Lieb lattice, seen in Fig. 3.1a, and the 3D perovskite lattice seen in Fig. 3.3a. The relevant Hamiltonian reads H=−  ij α  tij c†iα cjα − δ  + iλ  (d1ij ij αβ  ×  c†iα ciα i∈corner  d2ij )  · σαβ c†iα cjβ ,  (3.1)  where tij are the hopping amplitudes which we take equal to t1 , t2 and t3 for the first, second and third nearest neighbor sites respectively, and δ is an onsite energy for the corner sites (indicated as filled circles in Figs.3.1a and 3.3a). In Chapter 2, we showed that for nearest neighbor hopping (t1 only) and non-zero λ both Lieb and perovskite lattices became topological insulators possessing a nontrivial ❩2 invariant. In both cases the bands carrying nontrivial invariants were strongly dispersing. From here then, our goal is to flatten out these bands by tuning the hopping strengths t2 , t3 , λ and the onsite energy δ without leaving the topological phase. Also, we want the flat band to be separated from all other bands by a large gap ∆ so that interactions do not mix states from different bands. The relevant figure of merit,[67–69] then, is the ratio of the flat-band width W to the bandgap ∆.  3.3  Lieb lattice  As a warm-up exercise we first consider the valence band in the 2D Lieb lattice. Going over to momentum space and diagonalizing the Hamiltonian  37  3.3. Lieb lattice Hk that follows from Eq. (3.1), now having the additional terms Hkf lat   −δ + c1 0 0 = 0 c1 c2  , 0 c2 c1   (3.2)  where c1 = −2t3 (cos(2kx ) + cos(2ky )) c2 = −2t2 cos (kx ) cos (ky ), the result of the tuning procedure can be seen in Fig. 3.1b. (a)  (b) kx  a2 1  2  a1  3 ky  Figure 3.1: (a) Lieb lattice showing 3-site basis in unit cell and basis vectors a1 and a2 . Second and third neighbor hoppings are indicated by red and green dotted lines, respectively. (b) Tight-binding dispersion for Hk after tuning parameters to achieve maximal figure of merit ∆/W . For clarity the bulk energy bands are plotted over the shifted Brillouin zone 0 ≤ kx ≤ π/a, 0 ≤ ky ≤ π/a where a represents the length of the nearest neighbor bond. The procedure used for the tuning was the simplest available, namely introducing the parameters one by one and modifying the values by hand until the band was as flat as we could achieve. For the following set of values: t1 = 1, t2 = −0.17, t3 = −0.114, λ = 0.303 and δ = −0.725, the ratio of the bandwidth to the bandgap, ∆/W ≈ 40.45. A systematic numerical minimization technique could probably improve upon our result, but comparing to previous results on the kagome, honeycomb, square, checkerboard and ruby lattices, [67–70] this value is already quite respectable. To show that the system has not migrated away from the topological phase during the parameter tuning stage, we first solve the system in a strip 38  3.3. Lieb lattice (a)  (b)  4 1  10  E t  2  ∆ W 0  10  0  ï1  −2  0  π /2  kx  π  10  0.0  0.2  0.4  x  0.6  0.8  1.0  Figure 3.2: (a) Bandstructure for strip of Ny = 10 unit cells with open boundary conditions along the y direction and infinite along x with the same parameters as in Fig. 3.1b. (b) Semi-Log plot showing the evolution of the gap as one moves between a known topological insulating phase and the final set of finely tuned parameters. geometry numerically using the same parameters as above. This verifies that the spin-filtered gapless edge states, seen in Fig. 3.2a, indeed persist. We also consider the Hamiltonian H (x) = (1 − x)H0 + xH,  (3.3)  where H0 includes only t1 and λ and is known to support the topological phase.[78] The result for ∆/W , starting at λ = 0.1 (t1 = 1) then adiabatically tuning x from 0 to 1, for the same final set of parameters in H, can be seen in Fig. 3.2(b). The ratio increases from its initial value ∆/W ≈ 0.17 to its largest value 40.45 at x = 1, without closing the gap. The system is seen to remain in the topological phase. What can one say about the fate of the ground state in the presence of repulsive interactions when the flat band is partially filled? Close to half filling, since the spin up and down electrons are decoupled, one expects the Stoner instability towards the ferromagnetic state. If we parametrize the interaction by the usual on-site Hubbard U then a gap ∼ U should open between the majority and minority bands when W < U < ∆. For large enough U one can focus on the spin-polarized electrons in the majority band. In the presence of sufficiently strong residual interactions (e.g. nearest-neighbor repulsion V ) and at fractional filling, this problem becomes similar to that considered in 39  3.4. Perovskite lattice (a)  (b) 4  a3  E 2 t  4  0  3 1  a2  ï2  2  a1  Γ  A  L  H  Γ  Figure 3.3: (a) Perovskite lattice showing 4 sites in unit cell along with basis vectors. Second and third neighbor hoppings are indicated by red and green dotted lines, respectively. (b) Bandstructure inside bulk along path of high symmetry for Hk with finely tuned parameter set. Refs.[67–69] with spinless fermions. Analogous arguments then suggest the emergence of FQH states under favorable conditions. Another possibility is the formation of a T -invariant fractional topological insulator[81] which can be pictured as two decoupled FQH states for the two spin species. Such a state might be favored if U V.  3.4  Perovskite lattice  We now approach the main objective of this work: constructing flat bands with a non-trivial ❩2 index in a 3D lattice. Our starting point is again the Hamiltonian given in Eq. (3.1) but now on the perovskite lattice, shown in Fig. 3.3(a). In Ref. [78] we showed that for nearest-neighbor hopping and λ < 0 the system becomes a (1; 111) strong topological insulator when the lowest doubly degenerate band is filled. In an analogous fashion to the two dimensional case then, we aim to flatten out the lowest energy band by means of including additional parameters t2 , t3 and δ. So, going over to k-space we  40  3.4. Perovskite lattice  1  10  ∆ W  0  10  ï1  10  0.0  0.2  0.4  x  0.6  0.8  1.0  Figure 3.4: Semi-Log plot showing the evolution of the gap as one moves between a known topological insulating phase and the final set of finely tuned parameters in the perovskite lattice.  have for the 4 × 4 blocks on the diagonal Hkf lat     −δ + c1 0 0 0  0 c1 c2 c3  , =  0 c2 c1 c4  0 c3 c4 c1  (3.4)  where c1 = −2t3 (cos(2kx ) + cos(2ky ) + cos(2kz )), c2 = −2t2 cos (kx ) cos (ky ), c3 = −2t2 cos (kx ) cos (kz ) and c4 = −2t2 cos (ky ) cos (kz ). The bandstructure for the set of parameters that maximize the bandgap to bandwidth ratio is shown in Fig. 3.3b. The following set of values: t1 = 1, t2 = −0.416, t3 = −0.047, λ = −0.331 and δ = 1.318, yielded the ratio ∆/W ≈ 17.56. We remark that the ratio was calculated over the entire Brillouin zone and in Fig. 3.3b we have simply chosen one possible path here to illustrate the result. Including additional parameters in H, such as longer range hoppings, could no doubt further improve the above figure of merit but we do not pursue this here. To ensure that the system remains in the topological phase, we adiabatically tune the parameter x from 0 to 1 in the Hamiltonian H (x) defined in Eq. (3.3). The result of this procedure can be seen in Fig. 3.4. The ratio, again, remains finite across the entire range from its initial value ∆/W ≈ 0.15 to its largest value 17.56 at x = 1. We conclude that the flat-band carries ❩2 index (1;111). Now imagine that the (doubly degenerate) flat band is partially filled and 41  3.5. Magnetic instabilities consider the effect of repulsive interactions. Unlike the 2D case discussed above where the residual U(1) spin symmetry is preserved despite the presence of SOC, in the 3D strong topological insulator the SU(2) spin symmetry is completely broken. In this situation it is not clear what the leading instability might be. In a similar flat-band setting describing Ir-based pyrochlore Y2 Ir2 O7 Pesin and Balents[24] argued for an exotic T -invariant spin-charge separated topological Mott insulator while others found more conventional magnetic phases.[82, 83] Below, we investigate the magnetic instabilities of our model in the simplest case with on-site repulsion and at exact half filling of the flat bands. We show that the system remains non-magnetic when interaction strength is below the critical value which is large compared to the bandwidth W . In this regime, therefore, the true many-body ground state could become a T -invariant fractional topological insulator.  3.5  Magnetic instabilities  We extend the Hamiltonian (3.1) above by including a Hubbard term (1) (1)  Hint = U  (l) (l)  ni↑ ni↓ + U  ni↑ ni↓  i  (3.5)  i,l>1  with the goal of mapping out the magnetic phases as a function of U and U by means of a standard mean-field calculation. The superscript l denotes the basis sites in the lattice, Figs. 3.1(a) and 3.3(a). We use the following decoupling ni↑ ni↓ → ni↑ ni↓ + ni↓ ni↑ − ni↑ ni↓ , (3.6) and definine the magnetization  (l)  (l)  (l)  (3.7)  (l)  (l)  (l)  (3.8)  mi = ni↑ − ni↓ , and occupation number n ¯ i = ni↑ + ni↓ . We can rearrange so that (l)  (l) ni↑  =  (l)  n ¯ i + mi 2  ,  (3.9)  42  3.5. Magnetic instabilities (a)  (b) 20  12  16  8  U/W  U/W  12  4  8 4  0 0  4  U/W  8  0 0  12  1  U/W  2  Figure 3.5: (a) Phase diagram for Lieb lattice (W ≈ 0.034), and (b) perovskite lattice (W ≈ 0.072). Critical values Uc , Uc for the onset of magnetic order reflect the charge distribution on the basis sites. The latter is approximately uniform for the Lieb lattice implying Uc ≈ Uc . On perovskite lattice most of the charge density is on site 1 causing Uc Uc . and  (l)  (l) ni↓  =  (l)  n ¯ i − mi  . 2 Substituting back into Eq. 3.5 then yields the terms  (3.10)  1 2 n ¯ − m2i 2 i i i (3.11) with a similar expansion for U . Furthermore, we focus on uniform magnetic (l) phases, such that mi = m(l) independent of the site index i. To map out the phase diagram, we minimize the ground state energy with respect to m(l) . The result for half filling of the flat bands can be seen in Figures 3.5(a) and 3.5(b), for the Lieb and perovskite lattices respectively. Here m = l m(l) is the total magnetization per unit cell. The key observation to be made here is that both systems remain non-magnetic over a significant range of interaction strengths. U  (1) (1)  ni↑ ni↓ =  U 2  n ¯ i (ni↑ + ni↓ ) + mi (ni↓ − ni↑ ) −  43  3.6. Conclusions (a)  (b)  4 4 2  E 2 t  E t 0  0 ï2  ï2 Γ  X  M  Γ  Γ  A  L  H  Γ  Figure 3.6: (a) Band structure for Lieb lattice with U/W = U /W = 15, m(1) ≈ 0.35 and m(2) = m(3) ≈ 0.32, and (b) Perovskite lattice with U/W = 5, U /W = 25, m(1) ≈ 0.62, m(2) = m(3) ≈ 0.11 and m(4) ≈ 0.05.  3.6  Conclusions  We have established that it is possible to engineer nearly flat bands characterized by non-trivial ❩2 topological invariants in both the 2D Lieb lattice and its three dimensional counterpart, the edge centered cubic (perovskite) lattice. With our 2D example we have simply added to the growing number of lattices capable of producing this behavior, whereas the 3D result is to the best of our knowledge completely novel. What makes these findings potentially broadly relevant is that they provide a concrete framework for future studies addressing the nature of the many-body ground state that occurs in the presence of repulsive interactions and at partial filling. Unlike the 2D case where the Laughlin liquid[65] furnishes a well established paradigm for the topologically ordered correlated ground state, in 3D such a template is presently lacking. Therefore, understanding the fate of the ground state of electrons in the topologically non-trivial 3D flat band in the presence of strong interactions beyond the simple mean-field analysis is clearly a problem outside the scope of this study. In this work we made a first step in this direction by demonstrating, within the mean-field approximation, that the ground state remains non-magnetic over a significant portion of the U -U phase diagram (Fig 3.5b). Since mean-field calculations tend to overestimate the stability of ordered phases the actual non-magnetic region will likely be somewhat larger and can in principle host a T -invariant fractional topological insulator. Our result in 3D thus provides a concrete 44  3.6. Conclusions foundation for addressing these interesting questions using the suite of analytical approaches or numerical many-body techniques and we hope that it might inspire further studies. Finally we remark that there exist many perovskites in nature. Our theoretical results demonstrate that a near-ideal situation (in terms of strong interaction physics) can arise in a very simple model in this family of lattices. Also alluded to previously[78] was the possibility of engineering a relevant 2D system by modulating a two-dimensional electron gas with a periodic potential having Lieb lattice symmetry, in an analogous fashion to the ‘artificial graphene’ created recently.[63] Another possibility lies with cold fermionic or bosonic atoms in optical lattices as discussed in Refs. [[26, 27]], with a detailed study on the Lieb lattice provided in Ref. [84].  45  Chapter 4  Quantum transport Before beginning the final chapter on our efforts to engineer a quantum spin Hall state in graphene, we take a detour to describe the techniques required to simulate the transport phenomena encountered in Section 5.4.  4.1  Landauer formalism and Fisher-Lee relation  The theoretical tools commonly used for dealing with transport in mesoscopic devices were developed by Landauer-B¨ uttiker, and we aim here to briefly introduce the technique. For a more in depth review we turn the interested reader towards Ferry and Goodnick or Datta,[85, 86] but below we try to include enough detail to keep this thesis self contained. Throughout this section, and the rest of the chapter, we have also benefitted greatly from the efforts by Georgo Metalidis and Petra Dietl, [35, 87] to distill the treatments found in the textbooks and follow their approach quite closely. The standard treatment is to have the central mesoscopic device hooked up to electron reservoirs via leads, and the current through the device is then directly related to the probability that an electron has been transmitted through it. The leads are considered to be free from defects and translationally invariant in the direction of propagation, although their finite width will yield quantized modes in the transverse direction. If these channels are labelled by an index m the result for the conductance looks as follows G=  2e2 h  Tm (E).  (4.1)  m  If we had a uniform electron gas travelling at speed v with n electrons per unit length, then we would have a current equal to env. In our case, the transverse modes from the leads will each have a dispersion relation Em (k) and if the modes are occupied according to the Fermi distribution function  46  4.1. Landauer formalism and Fisher-Lee relation f (E), then we can write for the current I1→2 =  2e L  m,k  Tm (E)v(k)f (E − ν1 )  (4.2)  where v(k) is the group velocity, ν1 is the chemical potential in lead 1 and we have used the constant electron density 1/L for a single mode propagating along a wire of length L. The overall current can then be found by subtracting the contribution from the second lead, with chemical potential ν2 , giving 2e Tm (E)v(k) [f (E − ν1 ) − f (E − ν2 )] (4.3) I= L m,k  Using the rule for converting a discrete sum in k-space to an integral  k  →  L 2π  dk  (4.4)  and the expression for the group velocity v(k) =  1 dE h dk  (4.5)  allows us to rewrite the above as an integral over the energy I=  2e L  m  dE Tm (E) [f (E − ν1 ) − f (E − ν2 )] .  (4.6)  For low temperatures and voltages, eV = ν1 − ν2 , and we are in the linear response regime allowing us to employ the Taylor expansion f (E − ν1 ) − f (E − ν2 ) ≈ δ(E − EF )(ν1 − ν2 ).  (4.7)  giving us the desired expression G=  2e2 I = V h  Tm (E)  (4.8)  m  So to calculate the transport properties of our device, we will need these transmission amplitudes. However, instead of calculating the scattering amplitude directly, it is common practice to find the Green’s function, which can then be used to find the transmission coefficient. The so called Fisher-Lee 47  4.2. Introduction to Green’s functions relations accomplish this task,[88] which we introduce in the next section.  4.2  Introduction to Green’s functions  For a system described by Hamiltonian H, the time-independent Schr¨odinger equation reads [E − H] Ψ (r) = 0, (4.9) and the Green’s function operator can then be defined as [E − H] G(E) = ✶,  (4.10)  in this case giving the response of the system to a delta function perturbation. So long as E is not an eigenvalue of H the solution is given by G(E) = [E − H]−1 .  (4.11)  To remedy the fact that the expression is not well defined for all energies, we introduce the so called retarded and advanced Green’s functions G± (E) = lim [E ± iη − H]−1 ,  (4.12)  η→0+  and as they are related by hermitian conjugation, we shall use G for the retarded function and drop the ± notation for the remainder of the chapter. We shall also require the spectral function, which is defined as A = i G − G† .  (4.13)  However, we shall need to use it in one of its other guises, and thus introduce the expansion of the Green’s function in terms of the eigenbasis: G(E) =  1 = E + iη − H  k  |k k| E + iη −  (4.14) k  where the |k ’s are all eigenvectors of H with the corresponding eigenvalues  48  4.2. Introduction to Green’s functions k.  Expanding the spectral function in the eigenbasis yields 1 1 − E + iη − H E − iη − H 1 1 |k k| − E + iη − k E − iη −  A(E) = i = i  k  = k  |k k|  2η (E − k )2 + η 2  (4.15) (4.16) k  (4.17)  and letting η go to zero then gives A = 2π k  δ (E −  k ) |k  k|,  (4.18)  the desired expression. In order to motivate the effort below where we go on to calculate the Green’s function, we also introduce the Fisher-Lee relation for the transmission coefficient (4.19) T12 = T r Γ1 G12 Γ2 G†12 where  Γ1(2) = i Σ1(2) − Σ†1(2)  (4.20)  and the so called self-energies Σ1(2) can be thought of as effective Hamiltonians describing the interaction between the device and leads, with G12 a matrix of Green’s function elements connecting the first column of the central device to the last. The expression, 4.19, actually turns out to be intractable if H is of infinite dimension, which will occur if we take the semi-infinte leads into consideration. However, it is often possible to find the self energy analytically for clean leads in a simple geometry, leading to a finite Hamiltonian, or alternatively using a numerical iteration procedure. This is the method we shall use and introduce later in section 4.6 after first deriving an expression for the self energies.  49  4.3. Semi infinite leads-Self energy description  4.3  Semi infinite leads-Self energy description  We start by writing full Hamiltonian of our system, consisting of two leads and the central device, in the following block form   HL †  H= VLC 0  VLC HC VRC  0    † , VRC HR  (4.21)  where HL , HC , and HR are Hamiltonians for the left lead, the central device, and the right lead, and the off-diagonal terms, VLC and VRC contain the coupling terms between the device and the leads. The Hamiltonian is, of course, hermitian, and so the following relations are true † , VCL = VLC  † . VCR = VRC  The Green’s function can be subdivided in a can then be written as   − HL −VLC 0 GL †    −V † G − HC −VRC CL LC 0 0 −VRC − HR  (4.22)  similar fashion, and Eq. (4.10) GLC GC GRC   0 GCR  = ✶, GR  (4.23)  where we introduce = (E+iη)✶. We need to solve for GC , and so expanding the matrix equation above we find the series of equations (4.23) ( − HL ) GLC − VLC GC = 0,  (4.24)  −VRC GC + ( − HR ) GRC = 0.  (4.26)  † † −VLC GLC + ( − HC ) GC − VRC GRC = ✶,  (4.25)  From the first and the third equations we have GLC = ( − HL )−1 VLC GC , −1  GRC = ( − HR )  VRC GC ,  (4.27) (4.28)  and substituting into the second equation we then find ( − HC − Σ) GC = ✶,  (4.29)  50  4.3. Semi infinite leads-Self energy description where we have introduced the self energy, Σ = Σ1 + Σ2 , having the form † † Σ = VLC ( − HL )−1 VLC + VRC ( − HR )−1 VRC , Σ1  (4.30)  Σ2  And so for a central device coupled to two leads the Green’s function can be found by solving the following equation GC (E) = [(E + iη)✶ − HC − Σ]−1 ,  (4.31)  where the effects of the leads are included through the self-energy. It is also worth noting that the self-energy, (4.30), is determined solely by the coupling Hamiltonians and the retarded Green’s functions of the isolated leads Gi = ( − Hi )−1 (i = L, R). And so the self-energy is independent of the state of the central device itself. As was noted earlier, although we have managed to successfully account for the leads through the self-energy term, the solution still involves the inversion of an infinite matrix, and so some further step is required to find a tractable alternative. However, because we are limiting ourselves to nearestneighbour hopping from lead to device, the matrices VLC and VRC in Eqn. 4.30 will only have nonzero elements between the sites on the surface of each lead and their neighbouring sites in the central device. Therefore only the 11 surface Green’s functions, G11 L and GR are required, which we solve for in section 4.6. Lastly, to show that the self energy terms described above are equivalent to the ones found in the Fisher-Lee relation, Eqn. 4.19, we now aim to calculate the current through the system in a complementary fashion to section 4.1. We assume, as before, that the leads are in equilibrium, and there is small voltage V applied between them. It is usually necessary to use the nonequilibrium Green function method, as utilized below for the local current, but as this is a non-interacting system, we can get by with a simplified approach, following Paulsson. [89] So, we start by writing the Schr¨odinger equation, Eqn. 4.9, in the following block form   HL  V† LC 0  VLC HC VRC  0 † VRC HR       ΨL ΨL   ΨC  = E  ΨC  , ΨR ΨR  (4.32)  51  4.3. Semi infinite leads-Self energy description where ΨL , ΨC , and ΨR are vector wave functions of the left lead, the central device, and the right lead respectively. We then proceed by assuming we have an incoming wave in the left lead, Ψ0L , that is known, and use the ansatz ΨL = Ψ0L + Ψ1L , where the reflected component Ψ1L , as well as ΨC , and ΨR , appear only as a result of the interactions between the subsystems. The other necessary condition is that we find a retarded solution. With these conditions, we solve for equation (4.32) yielding † Ψ1L = 1 + G0L VLC GC VLC Ψ0L ,  (4.33)  † ΨR = G0R VRC GC VLC Ψ0L ,  (4.34)  ΨC =  † Ψ0L . GC VLC  (4.35)  where GC is the full Green’s function of the central device including the contribution from the self energies of the leads. Assuming we are in steady-state, the probability to find an electron in the device, |ψC |2 , will be conserved, and therefore †  ∂ψC ∂|ψC |2 † ∂ψC = ψC + ψC . (4.36) ∂t ∂t ∂t Using this condition, and the time-dependent Schrodinger equation, 0=  i  ∂ψ = Hψ, ∂t  (4.37)  we find  ∂ψC i † † ψR , (4.38) ψL + HC ψC + VRC =− VLC ∂t with a similar expression for the hermitian conjugate. Substituting back into Eqn. 4.36 then results in an expression of the form  where ji=L,R =  i  j = jL + jR  (4.39)  † Ψ†i ViC ΨC − Ψ†C ViC Ψi ,  (4.40)  which we interpret as probability currents coming from the left and right leads. The electrical current then follows from multiplying through by the charge e. 52  4.3. Semi infinite leads-Self energy description To find the current in terms of the Green’s functions we’ll now substitute in our expressions for the wave functions above (Eqn. 4.33-4.35)). Starting with the right lead, IR→L =  ie  =  ie  =−  † Ψ†R VRC ΨC − Ψ†C VRC ΨR  e  † † 0† † 0 0 Ψ0† L VLC GC VRC GR − GR VRC GC VLC ΨL † † 0 Ψ0† L VLC GC ΓR GC VLC ΨL .  (4.41)  The full current of all possible left eigenstates is given by I=  e  jλ = λ  † † 0 Ψ0† Lλ VLC GC ΓR GC VLC ΨLλ fL (Eλ ),  (4.42)  λ  the distribution function fL (Eλ ) describes the population of the left states, the distribution function of the right lead is absent here, because we consider only the current from the left to the right. The same current is given by the Landauer formula through the transmission function T¯(E) I=  e h  ∞  T (E)fL (E)dE.  (4.43)  −∞  If one compares these two expressions for the current, the transmission function at some energy is obtained as T (E) = 2π λ  † † 0 δ(E − Eλ ) Ψ0† Lλ VLC GC ΓR GC VLC ΨLλ  = 2π λ  δ  δ(E − Eλ ) Ψ0† Lλ VLC Ψδ  † Ψ†δ G†C ΓR GC VLC  = δ  = Tr ΓL G†C ΓR GC  2π λ  † Ψ†δ G†C ΓR GC VLC Ψ0Lλ  δ(E − Eλ )Ψ0Lλ Ψ0† Lλ  VLC Ψδ (4.44)  which is equivalent to the Fisher-Lee relation. As advertised, to evaluate the sum in brackets we used the eigenfunction expansion (4.18) for the left contact.  53  4.4. Local current  4.4  Local current  The expression for the local current relies on the non equilibrium Green’s function (NEGF) formalism, which is significantly more involved than the simple result used in the simulations. To go in depth here would take us quite far off track, so we will attempt to give a brief derivation following Jiang et al. [90] A more complete treatment of the NEGF formalism can be found in the textbooks by Ferry and Goodnick and Datta, [85, 86] as well as the more recent textbook by Di Ventra. [91] So, for a model on the lattice, the local current originating from site i to j is given as the expectation value of the time rate of change of the local charge at that point, which is then expanded using the Heisenberg equation of motion: Ji→j = e N˙ i =  ie  H,  Niα  (4.45)  α  =− where  e h  j  α,β  < Hiα,jβ G< iα,jβ (t, t) − Hjα,iβ Giα,jβ (t, t)  † G< iα,jβ (t, t ) = i ciα (t, t )cjβ (t, t )  (4.46)  is the Keldysh Green’s function and α, β represent the state of the system. Taking the Fourier transform of the above, so that we go from time to energy representation, yields Ji→j = −  2e h  ∞ α,β  −∞  Re Hiα,jβ G< iα,jβ (E) dE  (4.47)  For small applied voltages and in the zero temperature limit, which is the regime we are working in, applying the Keldysh equation G< = Gr (iΓL fL + iΓR fR ) Ga  (4.48)  54  4.5. Recursive Green’s function technique allows us to split up the integral above such that Ji→j =  2e h +  eVR α,β 2e2  h  Im [Gr (ΓL + ΓR ) Ga ] dE  −∞  (4.49) Hiα,jβ Gnjβiα (EF )](VL  Im[ αβ  − VR )  where VL , VR are the voltages at the left and right leads and α, β denote the state indices. Gn (EF ) = Gr ΓL Ga is electron correlation function with the line width function Γα (EF ) = i[Σrα (EF ) − (Σrα (EF ))+ ] , the Green function Gr (EF ) = [EF − Hcen − ΣrL (EF ) − ΣrR (EF )]−1 and the Hamiltonian in the central region denoted as Hcen . Because of time reversal symmetry, the first part of the equation above vanishes and thus we are left with the expression Ji→j =  4.5  2e2 Im[ h  α,β  Hiα,jβ Gnjβ,iα (EF )](VL − VR )  (4.50)  Recursive Green’s function technique  It is possible to solve Eqn. 4.31 directly, once the effects of the leads have been made finite, but this is computationally expensive and would limit the lattice sizes one can study, while also being inefficient as only part of this information is needed to find the conductance via Eqn. 4.19. So in lieu of finding the entire Green’s function in one step, we introduce the recursive Green’s function technique, where the lattice is divided into smaller sections which can be computed exactly and then ’glued’ together using Dyson’s equation, G = g + gV G  (4.51)  Here, g is the Green’s function of the isolated columns, which are known, V contains the hopping terms between columns and can be thought of as the perturbation, and G is the full Green’s function of the connected system we are interested in. Starting at one end of the sample then, it is possible to proceed through the entire lattice by adding one column at a time, arriving at the full Green’s function for the system once the opposite end has been attached. To begin, 55  inverting a C × C matrix. However, as the system size (i.e., the number of lattice sites) increases, this is a rather costly operation. Fortunately, as explained next, there is a more efficient way.  4.3.2  Recursive Green’s function technique  Instead of evaluating the total Green’s function at once, the lattice of the device is divided into several sections. Presuming that the Green’s function of each of these sections is 4.5. Recursive Green’s function technique known, they are ”glued together” by using Dyson’s equation (see, e.g. Ref. [20]) like Dietl, we illustrate theGtechnique = g + gV using G. the simplest possible case, a (4.19) system consisting of just two sections (illustrated in Fig. 4.1) which we would likeg to connect using Dyson’s as above.sections, Once it clear the In this equation is the Green’s function of equation the disconnected V isdescribes hoppinghow between the sections while G is we thecan Green’ connected system, to connect these two pieces, thenfunction move onoftothe finding the full i.e. the set Green’s functionfor wethe want to find (see Fig.function 4.2). Lettechnique us give anapplied example of equations recursive Green’s to of anhow to attach two isolated section with the So, aidinofthe Dyson’s equation. The is depicted arbitrary number of sections. basis of eigenstates forexample the columns, in Fig. 4.2. we start by projecting Dyson’s equation between the first and last columns, 1 and N of the system  gn1  gnn  gn+1,n+1  GN1  V 1  n  n+1  N  1  N  Figure 4.2: Two isolated sections are attached the help sections of Dyson’s in order Figure 4.1: Schematic illustrating how with two separate areequation joined toto obtain the Green’s function for the connected system. gether using Dyson’s equation to give the full connected Green’s function. (Figure from Dietl [35]) Suppose the propagators gn1 , gnn and gn+1,n+1 describing the isolated systems are known. The Green’s function for the coupled system one would like to obtain is GN 1 . The notation gn,m(i, j) stands for the propagator between site i of column n and site j of column m of the disconnected system: GN 1 = N |G|1  j) = ni|g|mj = Ngn,m |g|1(i,+ N |g|α .α|V |β β|G|1 |α ,|β  (4.20)  (4.52)  = N |g|n + 1 n + 1|V |n n|G|1  = gN,n+1 Vn+1,n Gn,1  Complete sets of lattice site states have been inserted, and we have used the fact that V has only non-zero elements between columns n and n + 1 because we are only considering nearest neighbour hopping. Since g is only available for the disconnected system, we have also set gN 1 and gN n to zero. Implied above is the splitting of the two separate pieces from 1 to n and then from n + 1 to N. To proceed, we can repeat the above procedure to find the unknown Green’s function Gn,1 and then again to find the closed set of equations: 56  4.5. Recursive Green’s function technique  Gn1 = gn1 + gnn Vn,n+1 Gn+1,1  (4.53)  Gn+1,1 = gn+1,n+1 Vn+1,n Gn1 ,  (4.54)  Substituting the second into the first and re-arranging we then have: Gn1 = [1 − gnn Vn,n+1 gn+1,n+1 Vn+1,n ]−1 gn1  (4.55)  which we can then substitute back into Eqn. 4.52 to get an expression for GN 1 in terms of the known Green’s functions for the isolated columns as desired. With this in hand, we can now proceed to derive the full blown technique for use on an arbitrarily large system in two dimensions. For convenience, we consider a rectangular tight binding lattice with M rows and N columns, but this will be perfectly compatible with the required honeycomb lattice used for our simulations, where each column will be constructed from graphene unit cells instead of just the individual sites. The effect of the leads attached to the device is accounted for through their self energies, which will yield a finite Hamiltonian, as shown in section 4.6, and is applicable so long as the leads are attached only to the left and right edges. If they were attached centrally, then one could end up with hopping terms that were not nearest neighbour and thus the current technique would break down. The first step is to divide the system into isolated columns, and for each of these from i = 1, 2...., N we calculate green’s function using direct inversion:   giiisol = (E + iη) ✶ − i|HD +  p=L,R  −1  Σp |i   (4.56)  where for the first and last columns, we have included the self-energy of the leads. We then start adding the columns one by one from the left to build the complete device, as shown in Fig. 4.2. If we imagine that the first n columns have been attached and we have the L L L Green’s functions GL n1 and Gnn , the Green’s functions Gn+1,1 and Gn+1,n+1 can then be derived by projecting Dyson’s equation (Eqn. 4.51) between the appropriate columns as shown above. The results are as follows:  57  4.3 Lattice Green’s function method  25  4.5. Recursive Green’s function technique isol g11  isol gnn  GL2,1  isol gNN  + + + +  + + +  + +  +  GLn+1,1  GLn1  GN1  Figure 4.2: Schematic illustrating the full recursive Green’s function technique deviceGreen’s is stitched together from left to rightofuntil the propFigure 4.3:where The the recursive function technique consists dividing the device agator between the first and last columns is obtained. (Figure from Dietl in single columns. After calculating the Green’s function of the isolated columns the [35])functions are added together from left to right. In the end the propagator Green’s between the first and the last column is obtained. When the last column N is connected, we obtain the Green’s function GLN 1 = GN 1 that we were looking for (see Fig. 4.3). From these Green’s functions the transmission −1 L extracted by the amplitudes canGbe expression (see Section = ✶ − g isolfollowing Vn+1,n GL Vn,n+1 g isol 4.2): n+1,n+1  GL n+1,1  n+1,n+1  =  n,n  † TRL = Tr[Γ GL Vn+1,n GNL RG 1 ΓL GN 1 ]. n+1,n+1 n,1  n+1,n+1  (4.26)  isol Starting from the far left side with n = 1 and GL 1,1 = g1,1 , it is then possible to proceed through the entire sample in this manner. After the last column has been connected we have GL N 1 = GN 1 , and our work is done. These set of steps give us all that is required to describe transport phenomena within the Landauer-B¨ uttiker formalism. To deal with describing the current across a sample, we’ll require an extension on this technique, which we now address. Above, we have found only GN 1 , which was all we needed for the overall conductance. However, for local properties we’ll require Gn1 for any n, as can be seen from Eqn. 4.50 below. The procedure starts from the right side of the device, stitching the  58  4.6. Surface Green’s function Green’s functions for the isolated columns together with Dyson’s equation R as was done above. At each step, GR nn is given in terms of Gn+1,n+1 as follows: isol R GR nn = ✶ − gnn Vn,n+1 Gnn Vn+1,n  −1  isol gnn  (4.57)  isol To start off, we set GR N N = GN N , and then apply the above for all n = N − 1, N − 2, ..., 1. From here, we can find the final set of Green’s functions by stitching together the GL and GR in pairs, attaching a strip L of connected columns 1 to n with known Green’s functions GL n1 and Gnn to the strip of columns N + 1 to N with known Green’s functions GR n+1,n+1 . Again, this is done for all n = 1, ...., N leading to the expression R Gn1 = ✶ − GL nn Vn,n+1 Gn+1,n+1 Vn+1,n  4.6  −1  GL n1  (4.58)  Surface Green’s function  As discussed in section 4.3, to incorporate the effect of the leads we need to find the surface Green’s functions. Finding the analytic solution is possible in most cases, but we will rely on a numeric technique to achieve the same end. For an isolated lead, the total retarded Green’s function looks the same as Eqn. 4.10 g (E) = [(E + iη) ✶ − Hl ]−1 (4.59)  The fact that these leads are perfect conductors in the Landauer-B¨ uttiker formalism, and thus translationally invariant, means we can write the Hamiltonian in block diagonal form      (E + iη) ✶ − Hl =       d −A 0 0 ··· −B D −A 0 · · ·   0 −B D −A · · ·   .. ..  .  . 0 −B D  .. .. .. .. . . . .  (4.60)  where d, D, A and B are all 2M × 2M matrices, with M being the width of the lead and the factor of 2 coming from spin. The matrix D contains the hoppings for a single isolated column, with A and B providing the hopping between neighbouring columns. The matrix d has been singled out because it contains the hoppings for the surface of the lead, but it is identical to D. 59  4.6. Surface Green’s function If the total Green’s function is then partitioned in a similar way, consisting of 2M × 2M submatrices, we can re-write the above as follows:          d −A 0 0 ··· −B D −A 0 · · · 0 −B D −A · · · .. .. . . 0 −B D .. .. .. .. . . . .          ·        g11 g21 g31 .. .  g12 g13 · · · g22 g23 g24 · · · g32 g33 g34 · · · . g42 g43 g44 . . .. .. .. .. . . . .        = ✶ (4.61)     If we then focus on the first column of the matrix product, we end up with an infinite set of equations for g11 , the surface Green’s function: dg11 = 1 + Ag21 , Dgp1 = Bgp−1,1 + Agp+1,1 ,  (4.62) ∀p ≥ 2.  (4.63)  g2r,1 = D−1 Bg2r−1,1 + D−1 Ag2r+1,1  (4.64)  Thus, it appears that in order to find g11 we would need to know all matrices gp1 with p ≥ 2. However, using the second equation, we can express the Green’s functions gp1 with even indices p = 2r (r = 1, 2, ...) as a function of the Green’s functions with odd indices:  Inserting this expression back into Eqn. 4.62, we can then write g11 as a function of only the gp1 with p odd: d − AD−1 B g11 = 1 + AD−1 A g31  (4.65)  D − AD−1 B − BD−1 A g2r+1,1 = BD−1 B g2r−1,1 + AD−1 A g2r+3,1 (4.66) Comparing these equations with the original expressions in Eqn. 4.62, they can be seen to have the same form if we define the following set of  60  4.6. Surface Green’s function renormalized matrices: d = d − AD−1 B,  D = D − AD−1 B − BD−1 A, A = AD−1 A,  B = BD  −1  (4.67)  B,  gr,1 = g2r−1,1 ,  r = 2, 3...  And so it is possible to iterate this procedure again and again, which will eventually allow us to arrive at an approximation for the surface Green’s function as the effective interaction A becomes vanishingly small. The equation for g11 after n iterations looks as follows: dn g11 = 1 + An g2n +1,1  (4.68)  and once An vanishes for sufficiently large n, we have the approximation for the surface Green’s function: g11 ≈ d−1 n  (4.69)  61  Chapter 5  Engineering a quantum spin Hall state in graphene 5.1  Introduction  As mentioned in the introduction, while landmark experiments indeed observed this phase in HgTe[19], experimental activity on the QSH effect has remained limited by notorious practical difficulties associated with this material. Given the comparative ease with which graphene can be fabricated, it would be highly desirable to artificially enhance graphene’s spin-orbit coupling strength to elevate the gap protecting the QSH state to observable levels. A practical means of doing so would pave the way to a new generation of QSH experiments accessible to a broad spectrum of researchers worldwide. Previous work has explored the possibility that strong interactions[60, 92, 93] or radiation[94] might drive graphene into a robust topological phase. Our goal in this chapter is to explore a practical new route to this end. Specifically, we will ask whether graphene can ‘inherit’ strong spin-orbit coupling from a dilute concentration of heavy adatoms (whose innate spinorbit coupling strength can be on the eV scale) deposited randomly into the honeycomb lattice. The basic principle underlying our proposal can be understood by considering processes in which an electron from graphene tunnels onto an adatom— whereupon it ‘feels’ enormous spin-orbit coupling—and then returns to the graphene sheet. Such processes in effect locally enhance graphene’s spinorbit strength manyfold. A viable proposal must, however, address numerous competing effects. For instance, adatoms often form local magnetic moments[95] which potentially spoil the time-reversal symmetry protecting the QSH effect. Moreover, adatoms generically mediate both intrinsic and Rashba spin-orbit coupling. The latter is believed to be detrimental to the QSH phase[11], and previous work has indeed established that certain kinds 62  5.1. Introduction  (a)  B  H T  (b)  (c) 4  3  5  2 6  1  c  Figure 5.1: (a) Depending on the element, adatoms favor either the highsymmetry ‘bridge’ (B), ‘hollow’ (H), or ‘top’ (T) position in the graphene sheet. (b) Detailed view of an H position adatom, which is best suited for inducing the intrinsic spin-orbit coupling necessary for stabilizing the topological phase. The desired spin-orbit terms mediated by the adatom are illustrated in (c). Red and yellow bonds represent the induced secondneighbor imaginary hopping, whose sign is indicated by the arrows for spin up electrons. For spin down electrons the arrows are reversed.  63  5.2. Physics of a single heavy adatom of adatoms do generate substantial Rashba coupling in graphene that typically overwhelms the intrinsic contribution [96–98]. (As Ref. [98] showed, however, magnetic adatoms inducing strong Rashba coupling can induce an interesting ‘quantum anomalous Hall’ state in graphene.) The adatoms may also favor competing, ordinary insulating states depending on their precise locations in the lattice. And finally, since spin-orbit coupling is generated non-uniformly in graphene, the stabilization of a QSH phase even in an otherwise ideal situation is unclear a priori. After an extensive search using tight-binding and first-principles analyses, we have found that two elements, indium and thallium, are capable of stabilizing a robust QSH state in graphene. Neither element forms a magnetic moment, and although they do generate significant Rashba coupling, for symmetry reasons this remarkably does not suppress the QSH state. We find that gaps many orders of magnitude larger than that predicted in pure graphene can form even with coverages of only a few percent; for example, at 6% coverage indium yields a gap on the order of 100K, while for thallium the gap approaches room temperature. These predictions revive graphene as a viable QSH candidate, and can be verified by probing the gap and associated edge states using spectroscopic and conductance measurements.  5.2  Physics of a single heavy adatom  To set the stage, let us first briefly review the Kane-Mele model[11] describing pure, undoped graphene with spin-orbit coupling. The Hamiltonian can be expressed as HKM = Ht + Hso , where Ht describes the usual nearestneighbor hopping and Hso encodes intrinsic spin-orbit coupling. In terms of operators c†rα that add electrons with spin α to site r of the honeycomb lattice and Pauli matrices sx,y,z that act on the spin indices, Ht and Hso explicitly read Ht = −t  (c†r cr + h.c.)  (5.1)  rr  (iνrr c†r sz cr + h.c.).  Hso = λso  (5.2)  rr  Here and below, spin indices are implicitly summed whenever suppressed. In Eq. (5.2) the sum runs over second-nearest-neighbor lattice sites, and νrr are signs that equal +1 if an electron hops in the direction of the arrows 64  5.2. Physics of a single heavy adatom in Fig. 5.1(c) and −1 otherwise. Thus Hso describes ‘chiral’ spin-dependent second-neighbor electron hopping. When λso = 0, the band structure exhibits the familiar gapless Dirac cones centered on momenta ±Q, resulting in semimetallic behavior. Turning on λso = 0 generates an energy gap[11] √ ∆ = 6 3|λso | at the Dirac points, transforming the system into a (very fragile[13–17]) QSH insulator. Importantly, if mirror symmetry with respect to the graphene plane is broken, then Rashba coupling—which involves spin flips and thus breaks the U(1) spin symmetry enjoyed by H—will also be present[11]. Rashba coupling competes with the intrinsic spin-orbit term in pure graphene, and beyond a critical value closes the gap and destroys the QSH state. If heavy adatoms are to stabilize a more robust QSH phase in graphene, then at a minimum they should be non-magnetic (to preserve T ) and modify the physics near the Dirac points primarily by inducing large intrinsic spinorbit coupling. The latter criterion leads us to focus on elements favoring the ‘hollow’ (H) position in the graphene sheet indicated in Fig. 5.1(a). Compared to the ‘bridge’ (B) and ‘top’ (T) positions, adatoms in the H position can most effectively mediate the spin-dependent second-neighbor hoppings present in the Kane-Mele model, while simultaneously avoiding larger competing effects such as local sublattice symmetry breaking generated in the T case. Since H-position adatoms generically reside on one side of the graphene sheet, they will clearly mediate Rashba spin-orbit coupling as well, leading to a potentially delicate competition. If the adatoms’ outer-shell electrons derive from p-orbitals, however, the induced intrinsic spin-orbit terms always dominate over the induced Rashba interactions. One can establish this key result by studying graphene with a single adatom of this type. To model this setup, we employ operators dmα for the adatom states, with m = 0, ±1 and α =↑, ↓ labeling the orbital and spin angular momentum quantum numbers. The coupling of these orbitals to graphene is conveniently expressed in terms of the following operators: 1 Cmα = √ 6  6  π  e−i 3 m(j−1) crj α ,  (5.3)  j=1  where the sum runs over the six sites surrounding the adatom shown in Fig. 5.1(b). Guided by symmetry (see section 5.6 for details), we consider the  65  5.2. Physics of a single heavy adatom following minimal Hamiltonian for this single-adatom problem: H = Hg + Ha + Hc  (5.4)  6  Hg = Ht − δµ Ha = m=0,±1  + Hc =  √  c†rj crj  (5.5)  j=1  † |m| dm dm  + Λso (d†1 sz d1 − d†−1 sz d−1 )  2Λso (d†0 s− d−1 + d†0 s+ d1 + h.c.)  † − (t|m| Cm dm m=0,±1  + h.c.),  (5.6) (5.7)  with s± = (sx ± isy )/2. Here Hg represents the nearest-neighbor hopping model of graphene supplemented by a chemical potential δµ for the six sites surrounding the adatom. Physically, a non-zero δµ leads to an excess electron density at those sites, screening any net charge from the adatom. Crystal field effects and spin-orbit coupling split the adatom p-orbitals through the |m| and Λso , Λso terms in Ha . Finally, Hc allows electrons to tunnel between the adatom and its neighboring carbon sites. All couplings in H are real except t1 , which is pure imaginary. Note also that we have ignored the exceedingly weak spin-orbit terms that couple the electrons in graphene directly, as well as Coulomb interactions for the adatom. Exclusion of the latter terms is justified for the non-magnetic adatoms we will consider below. Recalling that intrinsic spin-orbit coupling in graphene conserves the S z spin component while Rashba coupling does not, one can readily understand how the adatom mediates both kinds of interactions by viewing the tunnelings in Hc perturbatively. For example, when an electron hops onto the m = 0 orbital via t0 , flips its spin through the Λso coupling, and then hops back to the honeycomb lattice via t1 , spin-orbit coupling of Rashba type is locally generated in graphene. In sharp contrast to the situation for the Kane-Mele model, however, Rashba terms mediated in this fashion are irrelevant for the low-energy physics. This crucial feature can be understood by Fourier transforming C0α ; remarkably, the components with Dirac-point momenta ±Q vanish identically. More physically, electrons near the Dirac points interfere destructively when hopping onto the m = 0 adatom orbital. We stress that this argument holds only for p-orbitals. If the relevant adatom orbitals carry higher angular momentum, spin-flip processes which do affect the Dirac points will generically appear. The induced Rashba terms may 66  5.3. Periodic adatom configurations still be subdominant, though whether this is the case depends on details of the Hamiltonian unlike the situation for p-orbitals. In contrast, S z -conserving events whereby an electron hops both onto and then off of an m = ±1 orbital via t1 locally mediate the desired intrinsic spin-orbit interactions. No obstruction exists for tunneling onto the m = ±1 orbitals via t1 , so these processes can effectively modify the physics near the Dirac points. It is important to note, however, that t1 hopping mediates additional couplings as well, so the dominance of these spin-orbit terms remains unclear. Nevertheless, given the high symmetry preserved by the H-site adatom, none of the conventional broken-symmetry gapped phases of graphene—such as charge-density wave or ‘Kekule’ orders—are obviously favored here, so it is reasonable to expect these additional terms to play a relatively minor role. We explicitly confirm this intuition in the multiadatom situation, to which we now turn.  5.3  Periodic adatom configurations  As a first step in understanding the multi-adatom case, we examine a periodic system with one adatom residing in a large N × N supercell. This situation allows us to utilize density functional theory (DFT) to ascertain suitable heavy elements and obtain a quantitative understanding of their effects on graphene. (See section 5.7 for computational details.) To ensure large spin-orbit coupling, we focused on elements in rows five and six of the periodic table, including In, Sn, Sb, Te, I, La, Hf, Pt, Au, Hg, Tl, Pb and Bi. For each element, we calculated the total energy in the three adsorption geometries shown in Fig. 5.1(a) along with the adatom’s spin moment. Our calculations reveal that two elements—indium (atomic number Z = 49) and thallium (Z = 81)—satisfy our criteria of both favoring the high-symmetry H position and being non-magnetic. Furthermore, both elements exhibit partially filled p-shells, ensuring that the Rashba coupling they mediate in graphene is benign at the Dirac points. It is instructive to first examine the electronic properties of indium on graphene in a 4 × 4 supercell, without spin-orbit coupling. As Fig. 5.2(a) illustrates, the Dirac cones characteristic of pure graphene indeed remain massless—despite the reduced translation symmetry, conventional gapped phases are not stabilized here, consistent with the intuition developed in the single-adatom case above. Indium does, however, electron-dope graphene and shifts the Fermi level EF to 0.95 eV above the Dirac points. From the 67  5.3. Periodic adatom configurations  pz  1  px+py  EF 7 meV  0 -1  7 meV  Energy (eV)  2  (c) tight-binding, with SOC (In)  (b) DFT, with SOC (In)  (a) DFT, no SOC (In)  3  -2 -3  M  Γ  3  Γ 0  2  4  Γ  M  K  Γ 0  2  4 Γ  M  K  Γ 0  4  8  (f) tight-binding, with SOC (Tl)  (e) DFT, with SOC (Tl)  2 1  EF  0 -1  21 meV  21 meV  Energy (eV)  K  (d) DFT, no SOC (Tl)  -2 -3  Γ  M  K  Γ 0  3  6  Γ  M  K  Γ 0  3  6  Γ  M  K  Γ 0  3  6  Figure 5.2: All data correspond to one adatom in a 4 × 4 supercell, with the upper row corresponding to indium and the lower row corresponding to thallium. The left panels in (a) and (d) correspond to the band structure and LDOS computed using DFT without spin-orbit coupling. The horizontal dashed lines indicate the Fermi level (EF ), which shifts due to electrondoping from the adatoms. Insets zoom in on the band structure near the K point within an energy range −35 to 35meV, showing that without spin-orbit interactions neither indium nor thallium open a gap at the Dirac points. The central panels in (b) and (e) are the corresponding DFT results including spin-orbit coupling. Remarkably, in the indium case a gap of 7meV opens at the Dirac point, while with thallium the gap is larger still at 21meV. Finally, the right panels in (c) and (f) were obtained using the tight-binding model described in the text.  68  5.3. Periodic adatom configurations adatom’s local density of states (LDOS) displayed in Fig. 5.2(a), one can see that indium’s 5p orbitals lie almost entirely above EF , implying that the 5p electron in neutral indium nearly completely transfers to graphene. (The charge of an indium adatom is +0.8e from the Bader charge division scheme.) Note that the relatively diffuse pz LDOS indicates that this orbital hybridizes more strongly with graphene compared to the px,y orbitals. Replacing indium with thallium, again without spin-orbit coupling, leads to the band structure and LDOS shown in Fig. 5.2(d). Clearly the electronic structure is modified very little by this substitution; importantly, the Dirac cones remain massless with thallium as well. Thus any gap opening at the Dirac points must originate from spin-orbit coupling. Fig. 5.2(b) displays the band structure and LDOS for spin-orbitcoupled indium on graphene. Note the sizable spin-orbit splitting in the LDOS for the px,y orbitals. More remarkably, a gap ∆so ≈ 7meV now appears at the Dirac points, which already exceeds the spin-orbit-induced gap in pure graphene[13–17] by several orders of magnitude. The analogous results for thallium—whose atomic mass is nearly twice that of indium—are still more striking. As Fig. 5.2(e) illustrates, p-orbital splittings of order 1eV are now evident in the LDOS, and a gap ∆so ≈ 21meV opens at the Dirac points. We emphasize that these results apply for adatom coverages of only 6.25%. To explore the dependence of ∆so on the adatom coverage, we additionally investigated systems with one adatom in 5 × 5, 7 × 7 and 10 × 10 supercells. (For N × N supercells with N a multiple of 3, the Dirac points reside at zero momentum and can thus hybridize and gap out even without spin-orbit coupling. We thus ignore such geometries.) The values of ∆so along with the Fermi level EF computed for the coverages we studied appear in Fig. 5.3. Quite naturally, the gap decreases as one reduces the coverage, as does the Fermi level since fewer electrons are donated to graphene at lower adatom concentrations. It is worth highlighting that the gap decreases sublinearly and that a sizable ∆so ≈ 10meV remains even with a mere 2.04% thallium concentration. Since spin-orbit coupling clearly underlies the gap onset, it is tempting to associate the insulating regime opened by the adatoms with a QSH state. To verify this expectation, we appeal to tight-binding simulations of the 4×4 supercell geometry. We model the system by a Hamiltonian H4×4 consisting of uniform hopping for graphene (including weak second- and third-neighbor tunneling) supplemented by the additional local terms in Eq. (5.4) at each adatom. Figure 5.2(c) displays the tight-binding band structure and adatom 69  1.0  25  0.8  20  0.6  15  0.4  10  0.2  5  0.0  0  1  2  3  4  5  6  7  ∆so (meV)  EF (eV)  5.3. Periodic adatom configurations  0  Adatom Coverage (%) Figure 5.3: Spin-orbit-induced gap ∆so opened at the Dirac points and Fermi level EF (measured relative to the center of the gap) for different indium and thallium adatom coverages. The open and filled symbols represent data for indium and thallium, respectively.  LDOS for indium with the following parameters: first, second, and thirdneighbor graphene hoppings t = 2.82eV, t = 0.22eV, and t = 0.2eV; δµ = 0.5eV; 0 = 2.5eV, 1 = 1.8eV; t0 = 2eV, it1 = 0.95eV; and Λso = Λso = 0.1eV. The corresponding data for thallium appear in Fig. 5.2(f); parameters are the same as for indium, except for 1 = 1.9eV and Λso = Λso = 0.31eV. The good agreement with DFT provides strong evidence that the underlying physics is well-captured by our tight-binding model. With this in hand, one can then demonstrate the topological nature of the insulating regime by showing that the Dirac point gap remains non-zero upon adiabatically deforming the Hamiltonian to the Kane-Mele model. This can be achieved by defining a new Hamiltonian H(λ) = (1−λ)H4×4 +λHKM that interpolates between our adatom model at λ = 0 and the Kane-Mele model at λ = 1. For either indium or thallium parameters, the Dirac point gap never closes as λ varies from 0 to 1 (see section 5.8 for details). Thus either element can indeed stabilize a robust QSH state in graphene, at least when periodically arranged.  70  5.4. Randomly distributed adatoms  5.4  Randomly distributed adatoms  In an experiment adatoms will likely occupy random H-site locations in a graphene sheet that is of course not pristine. It is therefore important to understand the impact of such randomness and disorder on the stability of the topological phase that was found above in pure graphene with periodic adatoms. To address general features of the problem in a way that reduces its computational complexity, we will describe the disordered, random-adatom system using the following graphene-only Hamiltonian, H = Ht + I  δHI −  δµr c†r cr .  (5.8)  r  Here δµr represents a random on-site potential arising, e.g., from the substrate (and not the adatoms). In the second term, I labels the random plaquettes occupied by adatoms and δHI  = −δµ  c†r cr + λso r∈I  +iλR r,r ∈I  (iνrr c†r sz cr + h.c.) rr  ∈I  c†r (s × drr )z cr  (5.9)  with drr = r − r . The first term in δHI describes the chemical potential that screens charge from the adatoms, while the last two terms capture the local intrinsic (λso ) and Rashba (λR ) spin-orbit couplings induced by electrons hopping from graphene to an adatom and back, as discussed in Sec. II. Formally, these spin-orbit terms can be derived by integrating out perturbatively the adatom degrees of freedom in Hamiltonian (5.5-5.7), which we have done explicity in section 5.9 below. It is important to notice that unlike the conventional nearest-neighbor Rashba term considered in Ref. [11] for pristine graphene, the adatom-generated Rashba term connects all sites in the hexagon. Such a ‘hexagon Rashba term’ has the property of vanishing at the Dirac points as already argued in Sec. II. The adatom also induces other symmetry-allowed terms, such as further-neighbor spin-independent hoppings, which we disregard because they are either weak or do not lead to qualitative changes in the results reported below. To study the robustness of the QSH phase in the presence of random adatoms and disorder, we first provide a Born approximation for Graphene with random SOC and potential disorder. We start by considering the low71  5.4. Randomly distributed adatoms energy theory following Kane and Mele, [11] with Hamiltonian H = H0 + H  (5.10)  where d2 rψ † (r) (σx τz ∂x + σy ∂y ) ψ(r)  H0 = −i v  (5.11)  representing pure graphene, and H =  d2 rψ † (r) i  [u0 + λτz σz sz ] δ(r − Ri )ψ(r)  (5.12)  ˆ U  representing the potential and Dresselhaus interactions (we have ignored the Rashba term due to it vanishing at the Dirac point). Here, σ, τ and s are Pauli matrices in sublattice, valley and spin spaces, respectively, the Ri give the random positions of the adatoms, S0 is the area of a unit cell, needed for correct dimensionality and u0 , λ are couplings with dimension of energy. In momentum space (taking = 1), H0 = S  d2 k † ψ v (σx τz kx + σy ky ) ψk , (2π)2 k  (5.13)  ˆk h  and H =S  d2 q (2π)2  e−iRi ·q  S0  i  d2 k † ˆ ψ U ψk , (2π)2 k+q  (5.14)  S = N S0 being the area of the system where N is the number of unit cells. d2 k Using the equivalence S1 k ↔ (2π) 2 we may then write ˆ k ψk + ψk† h  H= k  † ˆq ψk ψk+q U  pq q  (5.15)  k  −iRi ·q . Aside from the matrix structure of h ˆk ˆq = S0 and pq = with U ie S ˆq , this Hamiltonian has a standard form for treating random disorder in and U metals (see Eqn. 5.2.2 in Doniach and Sondheimer [99]). We can therefore apply standard formulas to evaluate the disorder contribution to the self energy Σ(k, ω) in the first and second Born approximation. Then  G(k, ω) = G0 (k, ω)−1 − Σ(k, ω)  −1  .  (5.16) 72  5.4. Randomly distributed adatoms To first order ˆ (q = 0) Σ(1) (k, ω) = NI U ˆ = nI U  (5.17)  = nI [u0 + λτz σz sz ] , where nI = NI S0 /S is the impurity concentration. So, looking at the individual terms, nI U0 simply shifts the Fermi energy of Graphene away from the neutral point, which can be compensated by gating, and the Dresselhaus term, as hoped, opens up a gap ∆SO = nI λ resulting in a 2D topological phase. We now verify the above claims, while also probing for the signature gapless edge states by simulating the classic two-terminal conductance measurement in the geometry of Fig. 5.4.  LEAD 1  LEAD 2  W  L Figure 5.4: Configuration used to determine the two-terminal conductance. Semi-infinite clean graphene leads are connected to the sample, which has adatoms distributed randomly across it. The length L of the system is set by the number of columns (indicated by dotted lines). The width W refers to the number of sites along either edge of the figure, corresponding to the number of unit cells of an armchair nanoribbon required to create an individual column. Within the Landauer-B¨ uttiker formalism [100, 101] we employ the recursive Green’s function method [102] to evaluate the conductance G of a length-L, width-W graphene strip as we vary the Fermi energy. The details of this technique are supplied above in Chapter 4. As a benchmark we first 73  5.4. Randomly distributed adatoms provide a simulation for the conductance of nanoribbon including spin orbit coupling throughout its entirety, as well as one having a topologically trivial insulating gap arising from a staggered sublattice potential. In this way we can provide bandstructure plots, seen in Figs. 5.5a and b, to be compared directly with the respective conductance plots, Figs. 5.6a and b. (a)  (b)  3  3  2  2  1 E t 0  1 E t 0  ï1  ï1  ï2  ï2  ï3 −π  kx  0  π  ï3 −π  kx  0  π  Figure 5.5: Bandstructure plots for graphene nanoribbons having intrinsic spin orbit coupling, λ = 0.05t (a), and a staggered sublattice potential ν = 0.25t (b). As one sweeps across the chemical potential, for each new band crossed in Figs. 5.5(a) and (b), the conductance increases by exactly 2e2 /h, consistent with the Landauer theory. And for the ballistic edge states across the gap in Fig. 5.5(a) , we see a steady 2e2 /h, which is the conductance signature of a topological insulator. Having confidence in the numerical technique for a uniform system, we now consider the simplest case for the random-adatom setup with δµr = 0 (corresponding to an otherwise clean graphene sheet), δµ = 0, and λR = 0. Figure 5.7 illustrates the conductance for this case at several adatom concentrations ni , and model parameters indicated in the caption. A 2e2 /h plateau with width proportional to ni clearly emerges, strongly suggesting the onset of a bulk mobility gap and quantized ballistic conduction by edge states. This picture is corroborated by Fig. 5.8, which displays the current distribution for EF = 0.15eV across a smaller system size (chosen for clarity). These results establish that in principle adatoms need not be periodic to stabilize a QSH phase, though the parameters employed were unrealistic. To make contact with our DFT results, we also considered a system with 74  5.4. Randomly distributed adatoms (a)  (b)  20  18  18  16  16  14 12  12  G(e2/h)  G(e2/h)  14 10 8  10 8 6  6 4  4  2  2  0 ï1.0  ï0.5  0.0  E/t  0.5  1.0  0 ï1.0  ï0.5  0.0  0.5  1.0  E/t  Figure 5.6: Conductance plots corresponding to Figs. 5.5a and b above. λso = 0.02t and δµ = 0.1t, which would yield a similar bulk mobility gap to that seen in our thallium simulations if the adatom coverages were the same. While finite-size effects prevent us from studying the low-adatom coverages considered in Sec.5.3, Fig.5.9 shows that for ni = 0.15 a robust conductance plateau indeed persists for these more realistic parameter values. (At much lower coverages the width of the edge states approaches the system sizes we are able to simulate.) We have also studied the effects of the hexagon Rashba term and residual on-site potential disorder on the stability of the QSH phase, and for the remainder of this section we will consider random adatom arrangements, modeled using the graphene-only Hamiltonian, Eq. (5.8), with δµ = 0 for simplicity. Even a relatively large λR ∼ 2λso has a weak effect on the width of the conductance plateau, as expected given our analysis in Sec. II. With uncorrelated on-site disorder δµr , the topological phase is also remarkably robust—δµr can fluctuate on the scale of t while degrading the plateau only marginally. The topological phase is more sensitive, however, to correlated disorder which is likely more relevant experimentally. In this case the conductance plateau survives only when δµr varies on a scale smaller roughly than the mobility gap for the clean case. To observe the QSH phase experimentally, the chemical potential should therefore fluctuate on scales smaller than the ∼ 10meV gaps we predicted here. While this might be difficult to achieve with graphene on standard SiO2 substrate, recent experiments using hexagonal boron nitride substrates[103] show disorder energy scales as low as 15K, which should be sufficient to observe the QSH  75  5.4. Randomly distributed adatoms  2  G(e /h)  20 18  ni=0.1  16  ni=0.2  14  ni=0.3  12 10 8 6 4 2 0  ï0.6  ï0.4  ï0.2  0.0  0.2  0.4  0.6  EF(eV)  Figure 5.7: Conductance G as a function of the Fermi energy EF , averaged over 40 independent random adatom realizations at different concentrations ni = 0.1, 0.2, 0.3 for a system of size W = 80 and L = 40 with λso = 0.1t. (Error bars reflect range of conductance obtained over 40 random adatom distributions)  gap. So, we first consider the influence of the hexagon Rashba term defined in Eq. (5.8). While it is true that its effects identically vanish at the Dirac points, it nevertheless causes some splitting of the bands in their vicinity. As a result these Rashba terms can reduce somewhat the size of the mobility gap observed with either periodic or random adatoms, though to a far lesser extent than the traditional Rashba coupling. The splitting will depend on the ratio of λso to λR , and, more importantly, the absolute strength of λR compared to t. For purposes of comparison, we first solve for a system having uniform coverage so that we can find the dispersion and see how the splitting behaves for two different strengths of λR , keeping the ratio λso to λR the same. We start with λso = 0.04t and λR = 0.08t, their ratio chosen to be similar to that estimated in Eq. (5.50) below, and then halve the respective values for the second plot. The results can be seen in Figs. 5.10 (a) and (b). In Figures 5.11(a) and (b) we show the results for a set of parameters chosen to illustrate the effect in the random case while avoiding the complications due to the finite size of the system. As in Fig. 5.10(a), the mobility gap has clearly been somewhat reduced in Fig. 5.11(a), but a ro76  5.4. Randomly distributed adatoms  40 35 30 25  y  20 15 10 5 0  5  10  15  20  x  Figure 5.8: Current distribution across a sample of size W = 40 and L = 20 at ni = 0.2, λso = 0.1t and EF = 0.15eV . (The color and arrow length represent the magnitude of the current)  bust plateau remains. We remark that the standard Rashba term with only nearest neighbor coupling (introduced in Chapter 1), at the same coupling strength, would completely destroy the plateau. For the smaller values of λR and λso , seen in Fig. 5.11(b) the effect of hexagonal Rashba coupling becomes even less pronounced, although the finite size effects prevent us from obtaining reliable results at low adatom coverage in this regime. It is also worth emphasizing that one should not view λR as reducing the gap compared to those found in our DFT and tight-binding results displayed in Fig.5.2. Those simulations incorporated both the intrinsic and hexagon Rashba spin-orbit couplings (among other couplings neglected for simplicity here), and the results already reflect the influence of the latter. To study residual disorder unrelated to the adatoms (e.g., arising from impurities in the substrate), two models were used, both of which can be captured by an on-site potential term of the form specified in Eq. (5.8). In the first model, we incorporate an uncorrelated random chemical potential in the range [−∆, ∆] on every lattice site. In the second, we employ a longer-range, correlated disorder potential of the form[104] Nimp  δµri = j=1  Vj exp −  |ri − rj |2 , 2d2  (5.18) 77  5.4. Randomly distributed adatoms 20 18 16  G(e2/h)  14 12 10 8 6 4 2 0  ï0.3  ï0.2  ï0.1  0.0  0.1  EF(eV)  Figure 5.9: Conductance for the largest simulated system size using realistic parameters for thallium adatoms (λso = 0.02t and δµ = 0.1t) estimated from DFT data. Here W = 200, L = 100, t = 2.7eV and the coverage is ni = 0.15.  where Nimp is the number of impurity sites, Vj is also uniformly random over the range [−∆, ∆], and d is the radius of the impurity potential. As in Ref. [104], we characterize Nimp by the ratio ρ = Nimp /N , with N being the total number of carbon atoms in the sample. While we have not averaged over multiple disorder realizations as in Fig. 5.7, with only a single disorder realization one can still see the general trend towards instability of the QSH phase as ∆ increases. Our simulations show that this state, despite already being stabilized by a disordered arrangement of adatoms, in fact exhibits remarkable resilience against additional uncorrelated substrate-induced disorder in the graphene sheet. As Fig. 5.12(a) illustrates, in this case a very high value of ∆ ≈ t can be applied before the first signs of disorder affecting the topological phase appear, and even here a considerable conductance plateau survives. The topological phase is, however, more sensitive to long-range correlated disorder, which is likely more relevant for experiment when charged impurities from the substrate present the main disordering mechanism. Figure 5.12(b) illustrates that here the deterioration of the plateau is more rapid. At ∆ = 0.2t ≈ 0.54eV, with ρ = 0.07 and d = 5, there is still a remnant of the gapless edge states, but for stronger disorder the topological phase is destroyed. Notice that this 78  5.5. Relation to previous work (a)  (b)  2  2  E 0 t  E 0 t  ï2  ï2  Γ  M  K  Γ  Γ  M  K  Γ  Figure 5.10: Dispersion plots for uniform system with (a) hexagonal Rashba coupling λR = 0.08t and intrinsic spin-orbit strength is λso = 0.04t, (b) λR = 0.04t and λso = 0.02t. Inset shows dispersion enlarged at the Dirac point for each case. value of ∆ is comparable to the width of the plateau in the ‘clean’ case where the only disorder source is the randomly distributed adatoms.  5.5  Relation to previous work  Since numerous recent studies have explored the possibility of employing adatoms to control graphene’s electronic properties, here we briefly comment on the distinctions between our proposal and related earlier works. Castro Neto and Guinea[96] previously argued that graphene’s spin-orbit strength can be greatly enhanced by adatom deposition, though the physical mechanism and nature of induced spin-orbit coupling is quite different from what we considered here. The enhancement discussed in Ref. [96] stemmed from structural changes in the honeycomb lattice induced by a T-position adatom (such as hydrogen whose innate spin-orbit coupling is weak). We explicitly neglected such contributions, as our DFT simulations demonstrated that the adatoms we considered only weakly affect the carbon-carbon bonds. Perhaps even more importantly, the effect was discussed solely in the context of enhanced Rashba spin-orbit interactions, which do not stabilize a QSH state in graphene. More recently, Abdelouahed et al.[97] examined the influence of heavy Au adatoms on graphene’s electronic structure using first-principles calcu79  5.5. Relation to previous work (a)  (b)  20  16  18  hR=0  16  hR=0.08t  12  G(e2/h)  G(e2/h)  hR=0.04t  12  14 10 8 6  10 8 6 4  4  2  2 0  hR=0  14  ï0.3 ï0.2 ï0.1  0.0  EF(eV)  0.1  0.2  0.3  0  ï0.2  ï0.1  0.0  0.1  0.2  EF(eV)  Figure 5.11: Conductance G as a function of the Fermi energy EF for a system including (a) λso = 0.04t and λR = 0.08t, for a system of size W = 140 and L = 70 with adatom coverage ni = 0.2(b) λso = 0.02t and λR = 0.04t, for a system of size W = 160 and L = 80 with adatom coverage ni = 0.25 In both cases the nearest-neighbor hopping strength t = 2.7eV. lations. These authors considered dense coverages, with one T-position adatom residing above each carbon site on graphene’s A sublattice, and found substantial enhancement of graphene’s Rashba coupling to values of order 10meV. Again, this type of induced spin-orbit coupling is not suitable for driving a QSH state, and in any case this geometry is very likely to produce only a topologically trivial band insulator due to explicit breaking of sublattice symmetry by the adatoms. (Lighter Ni adatoms were also studied, but were found to modify graphene’s spin-orbit coupling only very weakly.) Closest in spirit to our study is the very interesting work of Qiao et al.[98], who considered the possibility that adatoms can stabilize another topological phase in graphene—a quantum anomalous Hall state. Unlike the QSH phase, this state breaks time-reversal symmetry and exhibits chiral edge modes leading to a non-trivial Hall conductivity as in the integer quantum Hall effect. Thus the requirements for stabilizing a quantum anomalous Hall phase in graphene are quite different. Qiao et al. showed that this can be accomplished by a combination of strong Rashba spin-orbit coupling and an exchange-induced Zeeman field. They further demonstrated via firstprinciples calculations that 6.25% coverage of periodically arranged iron adatoms (which form magnetic moments) provide both ingredients and generate a quantum anomalous Hall phase protected by a 5.5meV gap. Note 80  5.5. Relation to previous work (a)  (b)  14  14 6=0 6=0.25t 6=0.5t 6=t  12  10  G(e2/h)  G(e2/h)  10 8 6  8 6  4  4  2  2  0  ï0.6  ï0.4  ï0.2  0.0  EF(eV)  0.2  6=0 6=0.05t 6=0.1t 6=0.2t  12  0.4  0.6  0  ï0.6  ï0.4  ï0.2  0.0  0.2  0.4  0.6  EF(eV)  Figure 5.12: (a) Uncorrelated disorder randomly distributed over the range [−∆, ∆] and (b) Correlated disorder as defined in Eq. (5.18) with ρ = 0.07 and d = 5. In all cases the system is of size W = 80 and L = 40 with adatom coverage ni = 0.2, intrinsic spin orbit coupling strength λso = 0.1t and nearest-neighbor hopping strength t = 2.7eV. that while iron favors the H site like indium and thallium, its outer-shell electrons derive from d rather than p orbitals. Thus our argument for the irrelevance of Rashba spin-orbit coupling at the Dirac points does not apply to iron adatoms, consistent with the numerical findings of Ref. [98]. Although the gap produced by iron is smaller than that produced by indium and especially thallium at the same coverage, a promising feature of their proposal is that iron adatoms do not alter the Fermi level. We also wish to point out a connection of our work to ‘topological Anderson insulators’—systems which are non-topological in the clean limit but acquire non-trivial topology when disordered. Such states were first predicted in Ref. [105] and explored further in several subsequent studies[90, 106– 108]. The robust QSH state that can emerge upon randomly depositing dilute adatoms onto pure graphene (which by itself is a non-topological semimetal) can be viewed as an example of a topological Anderson insulator. Note, however, that this categorization is somewhat loose, as this phase in fact smoothly connects to the topological phase of the pure KaneMele model[11], which has no disorder (see also Ref. [108]). Nevertheless, this viewpoint is useful in that it highlights the remarkable fact that disorder, with suitable properties, can provide a practical avenue to generating topological phases. 81  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom  5.6  Discussion of the effective Hamiltonian for graphene with a single adatom  In this section we will consider a sheet of graphene located in the (x, y) plane, with a single adatom sitting at the H position (recall Fig.5.1). An important feature of such an adatom is the high symmetry that it preserves compared to the cases where the adatom resides at the B or T positions. In particular, the following symmetries remain here: π/3 rotation (Rπ/3 ) and x-reflection (Rx ) about the adatom site, as well as time-reversal (T ). Note that since the adatom is displaced vertically out of the (x, y) plane, mirror symmetry with respect to the graphene plane (Rz ) is broken. It will nevertheless be useful below to understand how operators transform under Rz , since this will give us insight as to which couplings mediate intrinsic versus Rashba spin-orbit coupling. We will continue to assume that p-orbitals form the active states for the adatom. We employ operators d†mα which add electrons with spin α and orbital angular momentum quantum number m = 0, ±1 to the adatom, and operators c†rα which add electrons with spin α to site r of graphene’s honeycomb lattice. These operators transform under Rπ/3 according to Rπ/3 : Rπ/3 :  π z  cr → ei 6 s cr  (5.19)  i π6 sz i mπ 3  dm → e  e  dm  (5.20)  where r corresponds to r rotated by π/3. Note that here and below, both cr and dm are considered to be two-component objects, with the upper and lower components corresponding to spin up and down, respectively. Similarly, under Rx we have Rx : Rx :  π x  c(x,y) → ei 2 s c(−x,y) i π2 sx  dm → e  d−m .  (5.21) (5.22)  Time-reversal, which is antiunitary, sends T :  T :  π y  cr → ei 2 s cr  m i π2 sy  dm → (−1) e  (5.23) d−m .  (5.24)  Note that the (−1)m factor above arises because one conventionally defines spherical harmonics such that [Y m ]∗ = (−1)m Y −m . Finally, the operators 82  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom transform under Rz as Rz :  π z  cr → −ei 2 s cr  (5.25)  1−m i π2 sz  Rz :  dm → (−1)  e  dm  (5.26)  In the last equation, the factor of (−1)1−m reflects properties of spherical harmonics under z → −z. Our aim now is to construct a minimal tight-binding model for graphene with a single H-position adatom. We write the Hamiltonian as H = Hg + Ha + Hc , where Hg involves only the carbon degrees of freedom, Ha involves only the adatom degrees of freedom, and Hc couples the two. Consider first Ha . The most general quadratic form for the adatom Hamiltonian can be written 1  uαm,m d†m sα dm ,  Ha =  (5.27)  α=0,x,y,z m,m =−1  where s0 is the identity, sx,y,z are Pauli matrices, and uαm,m = [uαm ,m ]∗ due to Hermicity. Using the transformation rules obtained above, one can show that symmetry strongly constrains the allowed uαm,m terms. For example, in the α = 0 sector, Rπ/3 implies that u0m,m must vanish except when m = m , while Rx further requires that u0m,m = u0−m,−m . (Time-reversal constrains these couplings no further.) The α = x, y, z terms can be similarly constrained, leading to the expression Ha = m=0,±1  +  √  † |m| dm dm  + Λso (d†1 sz d1 − d†−1 sz d−1 )  2Λso (d†0 s− d−1 + d†0 s+ d1 + h.c.).  (5.28)  Let us elaborate on the physics here in greater detail. If the adatom were in isolation, continuous rotation symmetry would further restrict |m| to be independent of m and require that Λso = Λso . In this case the Λso , Λso terms above simply encode an isotropic L · S spin-orbit interaction that split the orbitals into a lower j = 1/2 doublet and a j = 3/2 quadruplet that is higher in energy by 3Λso . Crystal field effects arising from the graphene environment lead to m-dependence of |m| and allow Λso to differ from Λso . Note also that we have assumed that the adatom’s s-orbitals (and all other inner levels) are far in energy from the Dirac points and can be safely neglected.  83  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom This is indeed justified by our supercell density functional theory (DFT) calculations for both indium and thallium, where the 5s and 6s electrons contribute significantly to the density of states only at several eV below the Dirac points. We have further ignored four-fermion interaction terms in Ha . This, too, is justified for indium and thallium, since both elements donate their outer p electron to graphene and form non-magnetic configurations. In a perturbative tight-binding picture, electrons hopping onto the adatom from graphene therefore lead to a singly occupied p orbital, for which Coulomb repulsion is expected to be unimportant. Finally, we note that Ha is even under Rz symmetry even though we did not enforce this. In other words, the fact that the adatom is displaced away from the graphene sheet leads to no additional terms in the Hamiltonian Ha compared to the (fictitious) case for an H-position adatom residing directly in the (x, y) plane. This fact is crucial for the irrelevance of induced Rashba terms for graphene, as we discuss further below. For the graphene-only part of the Hamiltonian, we posit that the usual nearest-neighbor hopping model Ht is modified primarily through a change in chemical potential for the six sites adjacent to the adatom: 6  Hg = Ht − δµ  c†rj crj .  (5.29)  j=1  Clearly Hg is invariant under not only Rπ/3 , Rx , and T , but also Rz despite this symmetry being broken by the adatom. (Terms which violate Rz necessarily involve spin-orbit coupling, which we have justifiably neglected because it is exceedingly weak.) Since the adatom generically modifies the bond lengths between carbon atoms in its vicinity, the hopping amplitudes will be modulated as well near the adatom. According to our DFT calculations, however, perturbations to graphene’s lattice structure are quite small (the carbon bond lengths in the immediate vicinity of either indium or thallium change by less than 1%). We thus expect this effect to be minor and have neglected it for simplicity. We have also neglected in Hg any spin-orbit terms which directly couple the carbon atoms, which is orders of magnitude weaker than the spin-orbit interaction indirectly mediated by the adatom. In contrast, the induced chemical potential δµ is not a weak effect— we find that our tight-binding simulations best reproduce the DFT band structure when δµ is on the order of 1eV. Physically, δµ appears because 84  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom the indium or thallium adatoms donate their outer p-orbital electron to the graphene sheet, leaving behind a net positive charge. Electrons in the graphene sheet thus prefer to conglomerate in its vicinity in order to screen the positively charged adatom. We note that of course there is no fundamental reason why the induced chemical potential should be confined to only the nearest six sites to the adatom. Our DFT simulations demonstrate, however, that the induced charge modulation occurs very locally around the adatom, so that this is expected to be a good approximation. As an illustration, consider the geometry with one thallium adatom in a 4 × 4 supercell. Figure 5.13 displays the induced charge density ∆ρ ≡ ρ(graphene + thallium) − ρ(graphene) − ρ(thallium) relative to the charge densities obtained when graphene and thallium decouple completely. Here the yellow/blue regions correspond to isosurfaces with ∆ρ = ±0.01eV/˚ A3 . Clearly the charge modulation occurs predominantly within the first two ‘rings’ of carbon sites surrounding the adatom. Finally, let us discuss the Hamiltonian Hc that allows electrons to tunnel between the adatom and graphene sheet. For simplicity, we only allow for tunneling events which couple the adatom to its six nearest carbon atoms, and assume (because of carbon’s weak intrinsic spin-orbit coupling) that such processes are spin-independent. It is useful to introduce the linear combinations 6 π 1 (5.30) e−i 3 m(j−1) crj α Cmα = √ 6 j=1 that, like dmα , carry angular momentum quantum number m under Rπ/3 . It follows from the symmetry transformation rules for cr obtained above that Cmα transforms as follows: Rπ/3 : Rx : T :  Rz :  π  (5.31)  C−m  (5.32)  π z  Cm → e i 6 s e i 3 m Cm i π2 sx  Cm → e  i π2 sy  Cm → e  Cm → −e  C−m  i π2 sz  Cm  (5.33) (5.34)  Using the first three symmetries along with the corresponding transformation rules for dm , it is straightforward to show that the spin-independent  85  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom  (a)  (b)  Figure 5.13: Adatom-induced charge modulation. (a) Top and (b) side views of the charge density induced by a thallium adatom in a 4 × 4 supercell. Yellow/blue surfaces correspond to a positive/negative induced charge density (relative to the case where graphene and thallium decouple completely). The charge modulations relax rather rapidly away from the adatom and occur mostly within the two innermost ‘rings’ of carbon sites near the adatom.  86  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom hoppings of interest must take the form Hc = −  † (t|m| Cm dm + h.c.),  (5.35)  m=0,±1  with t0 real and t1 pure imaginary. Under Rz symmetry, the t0 hopping is even while the t1 hopping terms are odd. This is easy to understand physically by considering the fictitious case where the adatom resides exactly in the graphene plane so that Rz symmetry is present. In this case the t1 terms reflect hopping between the adatom’s px,y orbitals and the carbon pz orbitals, which are orthogonal in this pathological limit. That is, tunneling between graphene’s pz orbitals and the adatom’s m = ±1 orbitals arises only when Rz is lifted. In contrast, the adatom’s pz orbitals overlap nontrivially with carbon’s even when Rz is preserved. Remarkably, the only terms in H that do not respect Rz are precisely these t1 couplings. To see what kind of spin-orbit processes the adatom can mediate, let us now consider the coupling between graphene and the adatom from a perturbative perspective. An electron from the graphene sheet can tunnel onto the adatom via either t0 or t1 . Due to the adatom’s strong spin-orbit coupling, the electron need not return to the graphene sheet via the same tunneling process by which it arrived. For example, an electron can tunnel via the t0 hopping, flip its spin, and then return to the graphene sheet via t1 . Consequently, at second order in perturbation theory, the adatom will mediate effective couplings between the carbon atoms that are proportional to t20 , t21 , and t0 t1 . The electron’s S z spin component is conserved in the former two cases, but not the third. Furthermore, since the t0 coupling in Hc is even under Rz while the t1 coupling is odd, it is only the induced terms for graphene proportional to t0 t1 that violate Rz . It follows that such spin-flip processes mediate Rashba spin-orbit coupling, while the t21 events mediate intrinsic spin-orbit interactions. (The t20 processes do not generate spin-orbit coupling.) It is worth noting as an aside that our ability to classify such processes in this simple way relies on the fact that all terms in Ha are even under Rz . If Ha also involved terms that violated Rz , then any second-order tunneling event would generically be capable of mediating Rashba spin-orbit coupling. To generate a quantum spin Hall (QSH) state in graphene, the physics near the Dirac points must be dominated by intrinsic rather than Rashba spin-orbit interactions. An important feature of our proposal is that Rashba coupling mediated by an adatom with a partially filled p shell is always 87  5.6. Discussion of the effective Hamiltonian for graphene with a single adatom subdominant to the induced intrinsic spin-orbit terms. This follows because electrons with momentum near the Dirac points hop very ineffectively onto the adatom’s m = 0 orbital. Previously we discussed how this effect arises because electrons destructively interfere when tunneling onto this orbital. Here we provide a complementary perspective using the language of lowenergy continuum Dirac fields. Focusing on physics near the Dirac points, one can express the lattice fermion operators for the A and B sublattices of graphene as follows, cr∈A,α ∼ eiQ·r ψR,Aα + e−iQ·r ψL,Aα iQ·r  cr∈B,α ∼ e  −iQ·r  ψR,Bα + e  ψL,Bα .  (5.36) (5.37)  Here ψR/Lα are slowly varying, two-component fields that describe lowenergy excitations in the vicinity of the Dirac points at momenta ±Q = 4π x ˆ (a is the separation between neighboring carbon atoms). Let us ± 3√ 3a use this decomposition to rewrite the operator Cm=0 that appears in the t0 tunneling. Neglecting pieces involving derivatives of ψR/Lα , we obtain  +  1 √ eiQ·r2 + eiQ·r4 + eiQ·r6 ψR,A 6 eiQ·r1 + eiQ·r3 + eiQ·r5 ψR,B  +  e−iQ·r2 + e−iQ·r4 + e−iQ·r6 ψL,A  +  e−iQ·r1 + e−iQ·r3 + e−iQ·r5 ψL,B .  Cm=0 ∼  (5.38)  Remarkably, all terms in parenthesis above vanish identically. The leading terms in the low-energy expansion of Cm=0 thus involve derivatives and become, consequently, increasingly unimportant as one approaches the Dirac points. This implies that the Rashba coupling mediated by the adatom—which necessarily involves tunneling via t0 —is highly ineffective at influencing the physics at the Dirac points. In contrast, Cm=±1 admits a non-trivial expansion in terms of ψR/Lα even when terms involving derivatives are dropped. It follows that electrons with momenta near the Dirac points can effectively tunnel onto and off of the adatom via t1 , and the intrinsic spin-orbit coupling mediated by these hoppings thus always dominates over the induced Rashba terms. This is not necessarily the case for adatoms with partially filled d or f shells, because there are then additional orbitals onto which electrons 88  5.7. Density functional theory from graphene can tunnel. For instance, here electrons can tunnel onto the adatom’s m = ±2 orbitals via a tunneling t2 ; such processes, unlike t0 tunneling, are not suppressed at the Dirac points. Effective couplings between carbon atoms proportional to t1 t2 are odd under Rz , involve spin flips, and correspond to Rashba spin-orbit interactions that may significantly affect the low-energy physics at the Dirac points. The importance of such terms compared to the induced intrinsic spin-orbit interactions then depends on details of the Hamiltonian. If, say, the additional orbitals hybridize very weakly with graphene or reside far in energy from the Dirac points due to crystal field effects, then d- or f -shell adatoms may still be suitable for stabilizing a QSH state in graphene.  5.7  Density functional theory  All DFT calculations were carried out with the Vienna ab-initio simulation package (VASP) [109, 110], at the level of the generalized-gradient approximation (GGA) [111]. We used the projector augmented wave (PAW) method for the description of the core-valence interaction [112, 113]. A vacuum space of 15˚ A was employed in all calculations, along with a 15 × 15 k-grid mesh[114] for integrals in the two-dimensional Brillouin zone. The energy cutoff of the plane wave expansion was set to 500eV. For pure graphene, the optimized lattice constant is 2.468˚ A, consistent ˚ with the experimental value of 2.46A. When addressing the preferred location (H, B, or T) of different elements on graphene, a 4 × 4 supercell was employed with one adatom per cell. Recall that we have two minimal criteria for stabilizing a QSH state in graphene: the adatoms should (i) prefer to occupy the H position and (ii) not form a net spin moment. Among the elements we searched in this work, In, I, La, Hf, Au, Hg and Tl satisfy the first criterion, while the remaining elements prefer either the T or B site. However, I, La, Hf and Au form magnetic moments and thus fail the second criterion, while the s and d bands of Hg reside far below graphene’s Fermi level, implying a negligible interaction between the two. The only two remaining elements—indium and thallium—satisfy both criteria and enter the H site non-magnetically. Including spin-orbit coupling, the binding energies Eb = E(adatom + graphene) − E(graphene) − E(atom)  (5.39)  89  5.8. Adiabatic continuity of the adatom and Kane-Mele models of indium and thallium on graphene at the H site are -0.525 and -0.133 eV, respectively. The total energy per site when indium resides at the T and B sites is higher by 87 and 80meV, respectively, than when indium resides at the H site. This is qualitatively consistent with the results of Ref. [95] where spin-orbit coupling was neglected. The conclusion is the same for thallium, though the preference for the H position is weaker; here the T and B sites are higher in energy by 32 and 22meV. At low temperatures compared to these energy differences indium and thallium are expected to occupy only the H position, with indium-carbon and thallium-carbon bond lengths of 2.83 and 2.90˚ A as determined by GGA calculations. We also calculated atomic structures of indium and thallium on graphene within the local density approximation (LDA). The indium-carbon and thalliumcarbon bonds are shortened by only 0.09 and 0.06˚ A, respectively, compared with the GGA results. It is thus reasonable to conclude that both GGA and LDA can describe the structural and electronic properties well.  5.8  Adiabatic continuity of the adatom and Kane-Mele models  We now turn to a periodic system with one indium or thallium adatom in a 4 × 4 supercell. Our first-principles calculations demonstrated that, with either element, the Dirac points remain massless in the non-relativistic limit but acquire a substantial gap upon including spin-orbit coupling. We further showed that the band structure and local density of states at the adatom site could be well-accounted for by an effective tight-binding model H4×4 . This model included first-, second-, and third-neighbor hopping amongst carbon atoms in the graphene sheet, as well as the additional couplings discussed in the previous section around each adatom. Note that inclusion of these further-neighbor tunnelings is necessary for obtaining reasonable quantitative agreement with first-principles band structure at energies away from the Dirac points, even for pure graphene[115]. Since the adatoms’ local density of states is concentrated at energies greater than 1eV above the Dirac points (see Fig. 5.2), one needs these additional hoppings to accurately model their hybridization with the graphene bands. (We stress that this is important for a quantitative description only. None of the qualitative features discussed in this chapter depend on these fine details. For instance, with only nearest-neighbor hopping a gap still opens at the Dirac points, 90  5.8. Adiabatic continuity of the adatom and Kane-Mele models  Energy gap [meV]  20 15 10 5 0 0  0.2  0.4  λ  0.6  0.8  1  Figure 5.14: Energy gap at the Dirac point as the Hamiltonian adiabatically deforms from the periodic adatom Hamiltonian H4×4 at λ = 0 to the Kane-Mele model at λ = 1. The solid blue and dashed red lines correspond to H4×4 evaluated with parameters appropriate for indium and thallium, respectively. In both cases the gap remains finite as λ varies from 0 to 1, indicating that the Kane-Mele model and adatom Hamiltonian are adiabatically connected and thus support the same quantum spin Hall phase.  though the local density of states for the adatoms’ pz states is difficult to accurately fit in this case.) One advantage of this effective tight-binding model is that it enables one to readily verify that the insulating regime generated by the spin-orbitinduced gap indeed corresponds to a QSH state. This can be perhaps most easily proven by demonstrating that H4×4 and the Kane-Mele model HKM [see Eqs. 5.1 and 5.2] can be smoothly deformed into one another without closing the Dirac point gap. To this end, consider the Hamiltonian H(λ) = (1 − λ)H4×4 + λHKM that interpolates between these models as λ varies from 0 to 1. For HKM (which includes only nearest-neighbor hopping and uniform intrinsic spin-orbit coupling for graphene), we choose t = 2.7eV and λso = −0.0015eV. As Fig. 5.14 illustrates, the Dirac point gap computed for H(λ) indeed remains finite for all λ between 0 to 1 with either indium (solid blue line) or thallium (dashed red line) parameters input into H4×4 . The QSH state known to be supported by HKM and the insulating state stabilized by either type of adatom are, consequently, the same topological phase of matter.  91  5.9. Graphene-only effective Hamiltonian  5.9  Graphene-only effective Hamiltonian  In this section we outline the derivation of the graphene-only Hamiltonian, Eqs. (5.8-5.9), from the tight binding model in Eqs. (5.4-5.7) describing the graphene sheet with H-position adatoms. We assume, for the purposes of this section only, that one adatom resides at each hexagonal plaquette but disallow any direct hopping of electrons between adatom p orbitals. Under these assumptions we may exploit the translation symmetry of the problem and write the Hamiltonian (4) in momentum space as H = k ψk† Hk ψk with A/B B T ψk = (cA remove electrons from the A/B k , ck ; d1k , d0k , d−1k ) , where ck graphene sublattices while dmk remove electrons from the adatom orbitals. As before, spin is treated implicitly. The Hamiltonian matrix in this basis reads hg T Hk = , (5.40) T † ha where hg =  0 τk , τk∗ 0  2  τk = −t  eik·δj  (5.41)  j=0  is the standard nearest-neighbor tight-binding Hamiltonian for graphene 2πj with δj = (sin 2πj 3 , − cos 3 ) unit vectors pointing from a site on sublattice A to its three neighbors on the B sublattice. Also, √   − z s 0 + Λ s 2Λ 1 so so √ √ ha =  2Λso s+ √ 0 (5.42) 2Λso s−  + z 2Λso s − Λ s 0 1 so  describes the adatom p orbitals. Finally T =  V1 V0 V−1 , ∗ V−1 V0∗ V1∗  2  Vm = t|m|  ei  2πj m 3  e−ik·δj  (5.43)  j=0  represents transitions between graphene and the adatom orbitals as described by Hc in Eq. (5.7). We now wish to ‘integrate out’ the adatom degrees of freedom and thus determine their effect on the electrons in graphene. The simplest way to ˜= accomplish this goal is to first perform a unitary transformation H → H ˜ g and h ˜ a in the transformed Hamiltonian e−S HeS with S chosen such that h ˜ are decoupled, i.e. T = 0. In this basis the adatom degrees of freedom can 92  5.9. Graphene-only effective Hamiltonian be integrated out trivially. Following the steps outlined in Ref. [15] (shown explicitly in the section 5.10 below), we find, to second order in T , ˜ g = hg − 1 T h−1 T † + hg T h−2 T † + h.c. + O(T 4 ). h a a 2  (5.44)  This result is equivalent to treating Hamiltonian (5.40) perturbatively to second order in T . We are interested in the effect of adatoms on the low-energy fermonic modes in graphene occuring in the vicinity of the Dirac points ±Q. In this limit, clearly, the second term in the brackets becomes negligible (since hg vanishes as k → ±Q), so we therefore focus on the first term. Although it is possible to find an exact inverse of ha the result is cumbersome and tends to obscure the simple physics underlying the formation of the spinorbit-induced gap in the system. To avoid this complication we make an additional assumption that spin-orbit splitting of the adatom orbitals is small compared to the crystal field effects, |Λso |, |Λso | | m |, and write ha = h0 + h1 with h0 = diag( 1 , 0 , 1 ). To first order in Λso and Λso we then −1 −1 have h−1 h−1 a 0 − h0 h1 h0 and ˜g h  −1 −1 † † hg − T h−1 0 T + T h0 h1 h0 T .  (5.45)  The adatom-induced correction to hg is seen to have a simple structure. The first term is spin-independent and mediates on-site potential terms and hoppings between sites in the hexagon surrounding an adatom. The second, spin-dependent term in turn produces the intrinsic and Rashba spin-orbit couplings. Various terms resulting from Eq. (5.45) can be explicitly worked out in a somewhat tedious but straightforward calculation, seen in section 5.11 below. The spin-orbit terms take the form δhSOC = g with  R1 + D R2 , † R2 R1 − D  (5.46)  2  |t1 |2 D = Λso √ 2 sz sin k · (δj − δj+1 ), 3 1 j=0  (5.47)  93  5.9. Graphene-only effective Hamiltonian R1 = Λso and  √  2t0 |t1 | 3 0 1  j,l  √ 2t0 |t1 | R2 = −iΛso 3 0 1  [s × (δl − δj )]z sin k · (δj − δl ),  j,l  [s × (δl + δj )]z e−ik·(δj +δl ) .  If we now recall that δ0 + δ1 + δ2 = 0 and identify √ 2t0 |t1 | |t1 |2 λR = Λso λso = Λso √ 2 , 3 0 1 2 3 1  (5.48)  (5.49)  (5.50)  we see that D and R1,2 correspond respectively to the intrinsic and hexagon Rashba spin-orbit couplings in Eq. (5.9). One can readily verify √ that R1,2 vanish at k = ±Q while D approaches a nonzero value of 3 3λso , as expected on the basis of arguments presented above. Using the tight-binding parameters obtained in Sec.5.3 we may estimate λso 23meV and λR = 58meV for thallium (for indium both are about factor of 3 smaller). At 6% coverage, assuming uniform averaging of√ λso over all plaquettes, this would imply a spin-orbit-induced gap of 0.06 × 6 3λso 14meV, somewhat smaller that the 21meV gap predicted by DFT at the same coverage. We attribute this discrepancy to the neglect of higher-order terms in Eq. (5.45); since t0,1 are of similar magnitude as 0,1 this expansion cannot be expected to yield quantitatively accurate values of the coupling constants. Nevertheless the procedure outlined above is useful in that it illustrates how various symmetry-allowed spin-orbit terms emerge in the effective graphene-only model. It also confirms that although the Rashba coupling exceeds the intrinsic coupling by more than a factor of two, its special form, dictated by symmetry, forces it to vanish at the Dirac point and renders λR irrelevant for the low-energy physics. Using the same methods one can also estimate various spin-independent terms induced by the adatoms. These include the on-site potential and additional hopping terms mentioned above. Since the latter are already present in pristine graphene and since adatom corrections are generally small compared to these bare values, we expect their effect on the low-energy physics to be minimal and do not discuss them further. To make a connection with the band structure plots provided in Fig. 5.2, we now apply the graphene only Hamiltonian to a 4 × 4 supercell, replacing the single adatam with a single plaquette carrying an imaginary next-nearest 94  5.10. Solution by canonical transformation  2  Energy (eV)  1 0 ï1 ï2 ï3 ï4  Γ  M  K  Γ  Figure 5.15: Energy gap for graphene only Hamiltonian applied to a 4 × 4 supercell with a single plaquette carrying imaginary next-nearest neighbor hopping. Here, λso 23meV resulting in a gap of ∼ 15meV. neighbor hopping with strength λso . For simplicity, as we only wish to verify that a gap persists even for dilute coverage, we have ignored the Rashba and onsite terms. As expected, a gap opens (Fig. 5.15), with a magnitude in agreement with the rough estimates made above.  5.10  Solution by canonical transformation  Expanding using the Hausdorff-Baker transform, we have for our canonical transform ˜ = e−S HeS = H − [S, H] + 1 [S, [S, H]] + .... H 2!  (5.51)  where we require ˜= H  ˜g 0 h ˜a 0 h  (5.52)  95  5.10. Solution by canonical transformation to leading order in T . So if we take S=  0 M M† 0  (5.53)  ensuring that eS is unitary, calculating the leading order term yields H − [S, H] =  hg T T † ha  MT† M ha † −M hg −M † T  −  +  −T M † hg M −ha M † T † M  (5.54)  And since the off-diagonal elements must be zero, we then have T = M ha − hg M,  (5.55)  T † = ha M † − M † hg ,  (5.56)  which we can use to isolate for M. First though we need to find the contribution to second order in T coming from 21 [S, [S, H]]. To proceed, we re-write Eqn. 5.54 in the form H − [S, H] =  hg − M T † − T M † 0 † 0 ha + M T + T † M  MT† + TM† T =H− † † T −M T − T † M  (5.57)  so that we have an expansion for [S, H], and neglecting the diagonal terms as they are already second order in T we find 1 1 0 T [S, [S, H]] = S, T† 0 2 2 1 = 2  MT† + TM† 0 † 0 −M T − T † M  (5.58)  Combining with the above then yields ˜ g = hg − 1 M T † + T M † h 2  (5.59)  96  5.11. Effective Hamiltonian SOC terms Finally, we can go ahead and solve for M iteratively using Eqn. 5.55 above, yielding −1 M = T h−1 a + hg M ha  (5.60)  −2 2 = T h−1 a + hg T ha + O(T )  which then gives us Eqn. 5.44 after substituting back into Eqn. 5.59.  5.11  Effective Hamiltonian SOC terms  To calculate the SOC terms in Eqn. 5.45 above we first introduce the pam to reduce clutter. Expanding, we then have rameter ξ = Vm −1 † T h−1 0 h1 h0 T  √    z 2Λso s− √ 0 ξ1 ξ−1 √Λso s +  2Λso s ξ0  2Λso s−   ξ0∗ √ 0 + ∗ z ξ−1 ξ1∗ 0 2Λso s −Λso s  =  ξ1 ξ0 ξ−1 ∗ ξ0∗ ξ1∗ ξ−1  =  R1 + D R2 † R1 − D R2  (5.61)  where R1 =  √  ∗ 2Λso s+ (ξ0 ξ1∗ + ξ0∗ ξ−1 ) + s− ξ0∗ ξ1 + ξ0 ξ−1  ,  D = Λso sz |ξ1 |2 − |ξ−1 |2 , and R2 =  √  2Λso s+ (ξ0 ξ−1 + ξ0 ξ−1 ) + s− (ξ0 ξ1 + ξ0 ξ1 )  (5.62) (5.63) (5.64)  97  5.11. Effective Hamiltonian SOC terms Working on the individual terms, we have for R1 iΛ t0 |t1 | + R1 = so s 6 0 1 =  Λso t0 |t1 | 6 0 1  Λ t0 |t1 | = so 6 0 1  2 j,l=0  −e−i  2πl 3  e−ik·(δj −δl ) + e−i  2 j,l=0  sin k · (δj − δl ) s+ e−i  2πl 3  2πl 3  eik·(δj −δl ) + h.c  + s− ei  2πl 3  2 j,l=0  sin k · (δj − δl ) −sx δly + sy δlx  (5.65)  2  Λ t0 |t1 | (s × δl )z sin k · (δj − δl ) = so 6 0 1 j,l=0 √ 2t0 |t1 | [s × (δl − δj )]z sin k · (δj − δl ) = Λso 3 0 1 j,l  matching Eqn. 5.48 above. For the off diagonal terms iΛ t0 |t1 | R2 = so 6 0 1  2  s+ e−i  2πl 3  e−ik·(δj +δl ) + s− ei  2πl 3  e−ik·(δj +δl )  j,l=0 2  Λso t0 |t1 | (s × δl )z e−ik·(δj +δl ) 6 0 1 j,l=0 √ 2t0 |t1 | = −iΛso [s × (δl + δj )]z e−ik·(δj +δl ) 3 0 1 =  (5.66)  j,l  98  5.12. Discussion again, matching Eqn. 5.49 as desired. Lastly, for the Dresselhaus term 2  2π 2π |t1 |2 D = Λso √ 2 sz e−i 3 (l−j) − ei 3 (l−j) eik·(δj −δl ) 3 1 j,l=0  2  |t1 |2 2π = Λso √ 2 sz sin (l − j) eik·(δj −δl ) − e−ik·(δj −δl ) 3 3 1 j,l=0 2  |t1 |2 2π = Λso √ 2 sz sin (l − j) sin k · (δj − δl ) 3 3 1 j,l=0 |t1 |2 = 2Λso √ 2 sz 3 1  sin l>j  (5.67)  2π (l − j) sin k · (δj − δl ) 3  2  |t1 |2 sin k · (δj − δj+1 ) = 2Λso √ 2 sz 3 1 j=0 yielding Eqn. 5.47, as desired.  5.12  Discussion  In this chapter we have demonstrated that dilute heavy adatoms—indium and thallium in particular—are capable of stabilizing a robust QSH state in graphene, with a band gap exceeding that of pure graphene by many orders of magnitude. Aside from fabrication ease, another virtue graphene possesses as a potential QSH system is the large number of experimental probes available in this material because the electrons are not buried in a heterostructure. Angle resolved photoemission spectroscopy (ARPES), for instance, actually benefits from electron doping by the adatoms, and thus could be employed to detect the onset of the Dirac point gap at arbitrary coverages. Scanning tunneling microscopy (STM) would be similarly wellsuited for probing the bulk gap, the spatial structure of the LDOS around an adatom, and even the topologically protected gapless edge modes. Detecting the edge states via transport will be more challenging, since this would require repositioning the Fermi level inside of the bulk gap. Nevertheless, the technology for doping graphene by the requisite amount does exist (see, e.g., Refs. [116–120]). Although the most commonly used ionliquid techniques [116–118] work best near room temperature, recently developed approaches employing ferroelectric substrates[119] and solid poly99  5.12. Discussion mer electrolytes[120] show comparable results and remain applicable down to cryogenic temperatures. The requisite doping may also be achievable by introducing additional high-electronegativity adatoms that absorb electrons. A quantum spin Hall phase in an appropriately doped graphene sample can be detected through the classic two-terminal conductance measurement akin to the historic observation performed in HgTe quantum wells[19]. A more complete characterization of the topological phase can be obtained through non-local edge state[121] and magnetotransport measurements which provide unambiguous evidence for topologically protected helical edge states [122, 123]. Fulfilling Kane and Mele’s original vision by artificially turning graphene into a strong spin-orbit system could have remarkable technological implications. Among these, the possibility of employing the protected edge states for realizing topological superconductivity, Majorana fermions, and related phenomena such as non-Abelian statistics and ‘fractional Josephson effects’ are particularly tantalizing[124]. Interestingly, one of the primary requirements for observing these phenomena—generating a superconducting proximity effect in graphene—has already been achieved by numerous groups[125–128]. These exciting possibilities, along with potential spintronics and low-dissipation devices, provide strong motivation for pursuing the line of attack we introduced for stabilizing a robust QSH state in graphene.  100  Chapter 6  Conclusion 6.1  Summary  In this work we established a new model for a topological insulator in both two and three dimensions, and using these models as a starting point we then went on to show that it was possible to flatten out the lowest energy bands in each case and explored the possibility of establishing a fractional topological insulator in the three dimensional version. For the final project we explored the possibility of engineering a quantum spin Hall phase in graphene via adatom deposition. We summarize the key results from each chapter below. In Chapter 2 we employed a simple tight binding model for electrons on the Lieb and perovskite lattices, and were able to show that in the presence of spin-orbit coupling they both entered the topological insulating phase. We calculated the Z2 invariants from the parity eigenvalues of the Hamiltonian at the time reversal invariant momenta for both lattice models, and provided band structure plots showing the existence of gapless edge states in each case. In Chapter 3 we added further range hoppings to the topological insulator models established in Chapter 2, and then maximized the relevant figure of merit, the ratio of the flat-band width W to the bandgap ∆, by fine tuning this set of hoppings. By adding interactions in the form of a mean-field decoupled Hubbard term, we then went on to show that a region of phase space remained that was non magnetic, which allowed us to speculate on the possibility of establishing a fractional topological insulating phase. In Chapter 5 we began by using DFT to find a suitable heavy adatom having large spin-orbit coupling that preferred to sit in the H position, above the center of the hexagonal plaquette. Indium and thallium both emerged as promising candidates, and further DFT studies, utilizing a periodic super cell containing an individual adatom, established that in the presence of spin orbit coupling a gap opened up at the Dirac point in the spectrum. For a modest coverage of ∼ 6%, a gap of 7meV was found for indium and 101  6.2. Future work 21meV for thallium, both orders of magnitude larger than what one can expect from pristine graphene. By fitting the DFT results with a tightbinding model, we were then able to show that the adatom model could be adiabatically connected to the Kane-Mele model, a known topological insulator, without closing the gap. We verified the effect could survive in a more realistic situation by performing transport calculations using the classic Landauer-B¨ uttiker formalism where the adatoms were randomly distributed, and completed the study by exploring the effects of disorder.  6.2  Future work  The topological insulator models developed in chapter 2 are currently toy models, in the sense that there are no known materials, to the best of our knowledge, that have the requisite lattice geometry and spin orbit coupling. So the hope is that our work will stimulate detailed band structure calculations of perovskites with heavy transition metal elements in search for new families of topological insulators. The prospects for finding naturally occurring materials that can support the flat bands discussed in chapter 3 are less likely still, due to the fine tuning required among the various short range hoppings. However, our result in 3D does provide a concrete foundation for future studies that may be able to answer whether it is possible to create a fractional topological insulating phase using the suite of available analytical approaches or numerical many-body techniques. Our adatom work looks more promising from an experimental verification standpoint, as the project has received some attention from the Max Planck-UBC center for quantum materials. Studies are currently underway by the Folk group here at UBC to probe the effects of depositing indium on graphene using low temperature electronic transport techniques and there are also plans for Kern’s group at MPI to explore the same setup using photoemission.  102  Bibliography [1] L. Fu and C. L. Kane. Superconducting proximity effect and majorana fermions at the surface of a topological insulator. Phys. Rev. Lett., 100: 096407, 2008. [2] B. Seradjeh, J.E. Moore, and M. Franz. Exciton condensation and charge fractionalization in a topological insulator film. Phys. Rev. Lett., 103:066402, 2009. [3] X. L. Qi, R. Li, J. Zang, and S.-C. Zhang. Inducing a magnetic monopole with topological surface states. Science, 323:1184, 2009. [4] A.M. Essin, J.E. Moore, and D. Vanderbilt. Magnetoelectric polarizability and axion electrodynamics in crystalline insulators. Phys. Rev. Lett., 102:146805, 2009. [5] X. L. Qi, T. L. Hughes, and S. C. Zhang. Topological field theory of time-reversal invariant insulators. Phys. Rev. B, 78:195424, 2008. [6] G. Rosenberg and M. Franz. Witten effect in a crystalline topological insulator. Phys. Rev. B, 82:035105, 2010. [7] J. E. Moore. The birth of topological insulators. Nature, 464:194, 2010. [8] M. Z. Hasan and C. L. Kane. Colloquium: Topological insulators. Rev. Mod. Phys., 82:3045, 2010. [9] X.-L. Qi and S.-C. Zhang. Topological insulators and superconductors. Rev. Mod. Phys., 83:1057, 2011. [10] L. Fu, C. L. Kane, and E. J. Mele. Topological insulators in three dimensions. Phys. Rev. Lett., 98:106803, 2007. [11] C. L. Kane and E. J. Mele. Quantum spin hall effect in graphene. Phys. Rev. Lett., 95:226801, 2005. 103  Chapter 6. Bibliography [12] C. L. Kane and E. J. Mele. Z2 topological order and the quantum spin hall effect. Phys. Rev. Lett., 95:146802, 2005. [13] D. Huertas-Hernando, F. Guinea, and A. Brataas. Spin-orbit coupling in curved graphene, fullerenes, nanotubes, and nanotube caps. Phys. Rev. B, 74:155426, 2006. [14] H. Min, J. E. Hill, N. A. Sinitsyn, B. R.Sahu, L. Kleinman, and A. H. MacDonald. Intrinsic and rashba spin-orbit interactions in graphene sheets. Phys. Rev. B, 74:165310, 2006. [15] Y. Yao, F. Ye, X. L. Qi, S.-C. Zhang, and Z. Fang. Spin-orbit gap of graphene: First-principles calculations. Phys. Rev. B, 75:041401, 2007. [16] J. C. Boettger and S. B. Trickey. First-principles calculation of the spin-orbit splitting in graphene. Phys. Rev. B, 75:121402, 2007. [17] 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, 2009. [18] B.A. Bernevig, T.L. Hughes, and S.-C. Zhang. Quantum spin hall effect and topological phase transition in hgte quantum wells. Science, 314:1757, 2006. [19] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang. Quantum spin hall insulator state in HgTe quantum wells. Science, 318:766, 2007. [20] L. Fu and C. L. Kane. Topological insulators with inversion symmetry. Phys. Rev. B, 76:045302, 2007. [21] 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:970, 2008. [22] H. M. Guo and M. Franz. Topological insulator on the kagome lattice. Phys. Rev. B, 80:113102, 2009. [23] H. M. Guo and M. Franz. Three-dimensional topological insulators on the pyrochlore lattice. Phys. Rev. Lett., 103:206805, 2009. 104  Chapter 6. Bibliography [24] D.A. Pesin and L. Balents. Mott physics and band topology in materials with strong spinorbit interaction. Nature Physics, 6:376, 2010. [25] Y. Xia, L. Wray, D. Qian, D. Hsieh, A. Pal, H. Lin, A. Bansil, D. Grauer, Y.S. Hor, R.J. Cava, and M.Z. Hasan. Observation of a large-gap topological-insulator class with a single dirac cone on the surface. Nature Physics, 5:398, 2009. [26] T. D. Stanescu, V. Galitski, J. Y. Vaishnav, C. W. Clark, and S. Das Sarma. Topological insulators and metals in atomic optical lattices. Phys. Rev. A, 79:053639, 2009. [27] S. L. Zhu, H. Fu, C. J. Wu, S. C. Zhang, and L. M. Duan. Spin hall effects for cold atoms in a light-induced gauge potential. Phys. Rev. Lett., 97:240401, 2006. [28] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang. Topological insulators in bi2se3, bi2te3 and sb2te3 with a single dirac cone on the surface. Nature Physics, 5:438, 2009. [29] K. v. Klitzing, G. Dorda, and M. Pepper. New method for highaccuracy determination of the fine-structure constant based on quantized hall resistance. Phys. Rev. Lett., 45:494, 1980. [30] 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, 1982. [31] D. B. Kaplan. A method for simulating chiral fermions on the lattice. Physics letters B, 288:342, 1992. [32] E. Fradkin and E. Dagotto. Physical realization of the parity anomaly in condensed matter physics. Phys. Rev. Lett., 57:2967, 1986. [33] B. A. Volkov and O. A. Pankratov. Jounal of Experimental and Theoretical Physics, 42:178, 1985. [34] 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, 2009. [35] P. Dietl. Numerical studies of electronic transport through graphene nanoribbons with disorder. Master’s thesis, Karlsruhe, 2008. 105  Chapter 6. Bibliography [36] C. Tao, L. Jiao, O. V. Yazyev, Y.-C. Chen, J. Feng, X. Zhang, R. B. Capaz, J. M. Tour, A. Zettl, S. G. Louie, H. Dai, and M. F. Crommie. Spatially resolving edge states of chiral graphene nanoribbons. Nature Physics, 7:616, 2011. [37] M. Fujita, K. Wakabayashi, K. Nakada, and K. Kusakabe. Peculiar localized state at zigzag graphite edge. Journal of the Physical Society of Japan, 65:1920, 1996. [38] S. M.-M. Dubois, Z. Zanolli, X. Declerck, and J.-C. Charlier. Electronic properties and quantum transport in graphene-based nanostructures. Eur. Phys. J. B, 74:235111, 2009. [39] K. Wakabayashi, K.-I. Sasaki, T. Nakanishi, and T. Enoki. Electronic states of graphene nanoribbons and analytic solutions. Science and Technology of Advanced Materials, 11:054504, 2010. [40] 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, 1988. [41] T. Thonhauser and D. Vanderbilt. Insulator Chern-insulator transition in the Haldane model. Phys. Rev. B, 74:235111, 2006. [42] S. Murakami. Quantum spin hall phases. Prog. Theor. Phys. Supplement, 176:279, 2008. [43] A. Strom. Interactions on the edge of a quantum spin Hall insulator. PhD thesis, Gothenburg, 2011. [44] L. Fu and C.L. Kane. Time reversal polarization and a z2 adiabatic spin pump. Phys. Rev. B, 74:195312, 2006. [45] J. E. Moore and L. Balents. Topological invariants of time-reversalinvariant band structures. Phys. Rev. B, 75:121306, 2007. [46] T. Fukui and Y. Hatsugai. Quantum spin hall effect in three dimensional materials: Lattice computation of z2 topological invariants and its application to bi and sb. Journal of the Physical Society of Japan, 76:053702, 2007.  106  Chapter 6. Bibliography [47] T. Fukui, T. Fujiwara, and Y. Hatsugai. Topological meaning of z2 numbers in time reversal invariant systems. Journal of the Physical Society of Japan, 77:123705, 2008. [48] R. Roy. Z2 classification of quantum spin hall systems: An approach using time-reversal invariance. Phys. Rev. B, 79:195321, 2009. [49] Z. Wang, X.-L. Qi, and S.-C. Zhang. Equivalent topological invariants of topological insulators. New J. Phys., 12:065007, 2010. [50] H. M. Guo and M. Franz. Three-dimensional topological insulators on the pyrochlore lattice. Phys. Rev. Lett., 103:206805, 2009. [51] J. C. Y Teo, L. Fu, and C. L. Kane. Surface states and topological invariants in three-dimensional topological insulators: Application to bi1?xsbx. Phys. Rev. B, 78:045426, 2008. [52] H. Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, L. Kleinman, and A. H. MacDonald. Intrinsic and Rashba spin-orbit interactions in graphene sheets. Phys. Rev. B, 74:165310, 2006. [53] S. Konschuh, M. Gmitra, and J. Fabian. Tight-binding theory of the spin-orbit coupling in graphene. Phys. Rev. B, 82:245412, 2010. [54] D. Huertas-Hernando, F. Guinea, and A. Brataas. Spin-orbit copuling in curved graphene, fullerenes, nanotubes and nanotube caps. Phys. Rev. B, 74:155426, 2006. [55] M. Konig, H. Buhmann, L. W. Molenkamp, T. Hughes, C.-X. Liu, X.-L. Qi, and S.-C. Zhang. The quantum spin Hall effect: Theory and experiment. Journal of the Physical Society of Japan, 77:031007, 2008. [56] P. A. Lee and T. V. Ramakrishnan. Disordered electronic systems. Rev. Mod. Phys., 57:287, 1985. [57] H. Suzuura and T. Ando. Crossover from symplectic to orthogonal class in a two-dimensional honeycomb lattice. Phys. Rev. Lett., 89: 266603, 2002. [58] K. Nomura, M. Koshino, and S. Ryu. Topological delocalization of two-dimensional massless dirac fermions. Phys. Rev. Lett., 99:146806, 2007. 107  Chapter 6. Bibliography [59] R. Shen, L. B. Shao, B. Wang, and D. Y. Xing. Single dirac cone with a flat band touching on line-centered-square optical lattices. Phys. Rev. B, 81:041410, 2010. [60] C. Weeks and M. Franz. Interaction-driven instabilities of a dirac semimetal. Phys. Rev. B, 81:085105, 2010. [61] D. Green, L. Santos, and C. Chamon. Isolated flat bands and spin-1 conical bands in two-dimensional lattices. Phys. Rev. B, 82:075104, 2008. [62] D.L. Bergman, C. Wu, and L. Balents. Band touching from real-space topology in frustrated hopping models. Phys. Rev. B, 78:125104, 2008. [63] M. Gibertini, A. Singha, V. Pellegrini, M. Polini, G. Vignale, A. Pinczuk, L. N. Pfeiffer, and K. W. West. Engineering artificial graphene in a two-dimensional electron gas. Phys. Rev. B, 79:241406, 2009. [64] D. C. Tsui, H. L. Stormer, and A. C. Gossard. Two-dimensional magnetotransport in the extreme quantum limit. Phys. Rev. Lett., 48: 1559, 1982. [65] R. B. Laughlin. Anomalous quantum hall effect: An incompressible quantum fluid with fractionally charged excitations. Phys. Rev. Lett., 50:1395, 1983. [66] P. K. Lam and S. M. Girvin. Liquid-solid transition and the fractional quantum-hall effect. Phys. Rev. B, 30:473, 1984. [67] E. Tang, J.-W. Mei, and X.-G. Wen. High-temperature fractional quantum hall states. Phys. Rev. Lett., 106:236802, 2011. [68] K. Sun, Z. Gu, H. Katsura, and S. Das Sarma. Nearly flatbands with nontrivial topology. Phys. Rev. Lett., 106:236803, 2011. [69] T. Neupert, L. Santos, C. Chamon, and C. Mudry. Fractional quantum hall states at zero magnetic field. Phys. Rev. Lett., 106:236804, 2011. [70] X. Hu, M. Kargarian, and G. Fiete. Topological insulators and fractional quantum hall effect on the ruby lattice. Phys. Rev. B, 84:155116, 2011. 108  Chapter 6. Bibliography [71] D.N. Sheng, Z.-C. Gu, K. Sun, and L. Sheng. Fractional quantum hall effect in the absence of landau levels. Nature Communications, 2:389, 2011. [72] N. Regnault and B.A. Bernevig. Fractional chern insulator. Phys. Rev. X, 1:021014, 2011. [73] J. Maciejko, X.-L. Qi, A. Karch, and S.-C. Zhang. Fractional topological insulators in three dimensions. Phys. Rev. Lett., 105:246809, 2010. [74] B. Swingle, M. Barkeshli, J. McGreevy, and T. Senthil. Correlated topological insulators and the fractional magnetoelectric effect. Phys. Rev. B, 83:195139, 2011. [75] X.-L. Qi. Generic wave-function description of fractional quantum anomalous hall states and fractional topological insulators. Phys. Rev. Lett., 107:126803, 2011. [76] C. Wu, D. L. Bergman, L. Balents, and S. Das Sarma. Flat bands and wigner crystallization in the honeycomb optical lattice. Phys. Rev. Lett., 99:070401, 2007. [77] D. L. Bergman, C. Wu, and L. Balents. Band touching from real-space topology in frustrated hopping models. Phys. Rev. B, 78:125104, 2008. [78] C. Weeks and M. Franz. Topological insulators on the lieb and perovskite lattices. Phys. Rev. B, 82:085310, 2010. [79] D.J. Thouless. Wannier functions for magnetic sub-bands. J. Phys. C, 17, 1984. [80] A.A. Soluyanov and D. Vanderbilt. Wannier representation of z2 topological insulators. Phys. Rev. B, 83:035108, 2011. [81] M. Levin and A. Stern. Fractional topological insulators. Phys. Rev. Lett., 103:196803, 2009. [82] X. Wan, A.M. Turner, A. Vishwanath, and S.Y. Savrasov. Topological semimetal and fermi-arc surface states in the electronic structure of pyrochlore iridates. Phys. Rev. B, 83:205101, 2011.  109  Chapter 6. Bibliography [83] W. Witczak-Krempa and Y.-B. Kim. Topological and magnetic phases of interacting electrons in the pyrochlore iridates. Phys. Rev. B, 85: 045124, 2011. [84] N. Goldman, D.F. Urban, and D. Bercioux. Topological phases for fermionic cold atoms on the lieb lattice. Phys. Rev. A, 83:063601, 2011. [85] D. K. Ferry and S. M. Goodnick. Transport in Nanostructures. Cambridge, 1998. [86] S. Datta. Electronic transport in mesoscopic sysetms. Cambridge, 1995. [87] G. Metalidis. Electronic transport in mesoscopic sysetms. PhD thesis, Halle, 2007. [88] D. S. Fisher and P. A. Lee. Relation between conductivity and transmission matrix. Phys. Rev. B, 23:6851, 1981. [89] M. Paulsson. Non equilibrium Green’s functions for dummies: Introduction to the one particle NEGF equations. arXiv:1010.4679, 2002. [90] H. Jiang, L. Wang, Q. f. Sun, and X. C. Xie. Numerical study of the topological anderson insulator in hgte/cdte quantum wells. Phys. Rev. B, 80:165316, 2009. [91] M. Di Ventra. Electronical transport in nanoscale systems. Cambridge, 2008. [92] S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang. Topological mott insulators. Phys. Rev. Lett., 100:156401, 2008. [93] T. Pereg-Barnea and G. Refael. Inducing topological order in a honeycomb lattice. Phys. Rev. B, 85:075127, 2010. [94] J.-I. Inoue and A. Tanaka. Photoinduced transition between conventional and topological insulators in two-dimensional electronic systems. Phys. Rev. Lett., 105:017401, 2010. [95] K. T. Chan, J. B. Neaton, and M. L. Cohen. First-principles study of metal adatom adsorption on graphene. Phys. Rev. B, 77:235430, 2008. 110  Chapter 6. Bibliography [96] A. H. Castro Neto and F. Guinea. Impurity-induced spin-orbit coupling in graphene. Phys. Rev. Lett., 103:026804, 2009. [97] S. Abdelouahed, A. Ernst, J. Henk, I. V. Maznichenko, and I. Mertig. Spin-split electronic states in graphene: Effects due to lattice deformation, rashba effect, and adatoms by first principles. Phys. Rev. B, 82:125424, 2010. [98] Z. Qiao, S. A Yang, W. Feng, W.-K. Tse, J. Ding, Y. Yao, J. Wang, and Q. Niu. Quantum anomalous hall effect in graphene from rashba and exchange effects. Phys. Rev. B, 82:161414, 2010. [99] S. Doniach and E. H. Sondheimer. Green’s functions for solid state physicists. Imperial College Press, 1998. [100] R. Landauer. Electrical resistance of disordered one-dimensional lattices. Philos. Mag., 21:863, 1970. [101] M. Buttiker. Absence of backscattering in the quantum hall effect in multiprobe conductors. Phys. Rev. B, 38:9375, 1988. [102] G. Metalidis and P. Bruno. Greens function technique for studying electron flow in two-dimensional mesoscopic samples. Phys. Rev. B, 72:235304, 2005. [103] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L K. L. Shepard, and J. Hone. Boron nitride substrates for high-quality graphene electronics. Nature Nanotechnology, 5:722, 2010. [104] S. Yuan, H. De Raedt, and M. I. Katsnelson. Electronic transport in disordered bilayer and trilayer graphene. Phys. Rev. B, 82:235409, 2010. [105] J. Li, R.-L. Chu, J. K. Jain, and S.-Q. Shen. Topological anderson insulator. Phys. Rev. Lett., 102:136806, 2009. [106] C. W. Groth, M. Wimmer, A. R. Akhmerov, J. Tworzydlo, and C.W. J. Beenakker. Theory of the topological anderson insulator. Phys. Rev. Lett., 103:196805, 2009.  111  Chapter 6. Bibliography [107] H.-M. Guo, G. Rosenberg, G. Refael, and M. Franz. Topological anderson insulator in three dimensions. Phys. Rev. Lett., 105:216601, 2010. [108] E. Prodan. Three-dimensional phase diagram of disordered HgTe/CdTe quantum spin-hall wells. Phys. Rev. B, 83:195119, 2011. [109] G. Kresse and J. Furthmuller. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci., 6:15, 1996. [110] G. Kresse and J. Furthmuller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54:11169, 1996. [111] J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 77:3865, 1996. [112] P. E. Blochl. Projector augmented-wave method. Phys. Rev. B, 50: 17953, 1994. [113] G. Kresse and D. Joubert. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B, 59:1758, 1999. [114] H. J. Monkhorst and J. D. Pack. Special points for brillouin-zone integrations. Phys. Rev. B, 13:5188, 1976. [115] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordej´on. Tight-binding description of graphene. Phys. Rev. B, 66:035412, 2002. [116] A. Das, S. Pisana, B. Chakraborty, S. Piscanec, S. K. Saha, U. V. Waghmare, K. S. Novoselov, H. R Krishnamurthy, A. K. Geim, A. C. Ferrari, and A. K. Sood. Monitoring dopants by raman scattering in an electrochemically top-gated graphene transistor. Nature Nanotechnology, 3:210, 2008. [117] J. Yan, T. Villarson, E. A. Henriksen, P. Kim, and A. Pinczuk. Optical phonon mixing in bilayer graphene with a broken inversion symmetry. Phys. Rev. B, 80:241417, 2009. [118] J. T. Ye, M. F. Craciun, M. Koshino, S. Russo, S. Inoue, H. T. Yuan, H. Shimotani, A. F. Morpurgo, and Y. Iwasa. Accessing the transport 112  Chapter 6. Bibliography properties of graphene and its multi-layers at high carrier density. arXiv:1010.4679, 2010. [119] Y. Zheng, G.-X. Ni, S. Bae, C.-X. Cong, O. Kahya, C.-T. Toh, H. R. ¨ Kim, D. Im, T. Yu, J. H. Ahn, B. H Hong, and B. Ozyilmaz. Waferscale graphene/ferroelectric hybrid devices for low-voltage electronics. EuroPhys. Lett., 93:17002, 2011. [120] D. K. Efetov and P. Kim. Controlling electron-phonon interactions in graphene at ultra high carrier densities. Phys. Rev. Lett., 105:256805, 2010. [121] A. Roth, C. Brune, H. Buhmann, L. W. Molenkamp, J. Maciejko, X.L. Qi, and S.-C. Zhang. Nonlocal transport in the quantum spin hall state. Science, 325:294, 2009. [122] G. Tkachov and E. M. Hankiewicz. Ballistic quantum spin hall state and enhanced edge backscattering in strong magnetic fields. Phys. Rev. Lett., 104:166803, 2010. [123] J. Maciejko, X.-L. Qi, and S.-C. Zhang. Magnetoconductance of the quantum spin hall state. Phys. Rev. B, 82:155310, 2010. [124] L. Fu and C. L. Kane. Josephson current and noise at a superconductor/quantum-spin-hall-insulator/superconductor junction. Phys. Rev. B, 79:161408, 2009. [125] H. B. Heersche, P. Jarillo-Herrero, J. B. Oostinga, L. M. K. Vandersypen, and A. F. Morpurgo. Bipolar supercurrent in graphene. Nature, 446:56, 2007. [126] X. Du, I. Skachko, and E. Y. Andrei. Josephson current and multiple andreev reflections in graphene sns junctions. Phys. Rev. B, 77:184507, 2008. [127] C. Ojeda-Aristizabal, M. Ferrier, S. Gu´eron, and H. Bouchiat. Tuning the proximity effect in a superconductor-graphene-superconductor junction. Phys. Rev. B, 79:165436, 2009. [128] D. Jeong J.-H. Choi, G.-H. Lee, S. Jo, Y.-J. Doh, and H.-J. Lee. Observation of supercurrent in pbin-graphene-pbin josephson junction. Phys. Rev. B, 83:094503, 2011. 113  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items