Keywords: Scanning tunneling microscopy

First and foremost, I would like to thank my supervisor Cristina. Thank you for giving me the chance to do this thesis, for supporting me from the start even though we hadn't worked together before. Then also for the three years that followed, for all the projects you've included me in, and for believeing in me and telling me that I could have a

1

Résumé

Le travail présenté dans ce manuscrit se concentre sur une technique particulière, que nous appellerons souvent la technique de l'impureté. Cette technique est née au carrefour de deux domaines différents de la physique de la matière condensée : l'étude des états de bord et de surface et l'étude des effets d'impuretés.

Ces dernières années, un regain d'intérêt a eu lieu pour les phénomènes de bord et de surface en raison de l'avènement de la matière condensée topologique. Ce domaine moderne, né dans les années 80 avec la découverte notable de l'effet Hall quantique (Klitzing et al., 1980), se focalise sur les phases de la matière (et les phénomènes uniques qui y sont associés) qui ne sont pas expliquées par la théorie de Landau des transitions de phase. La théorie de Landau est basée sur la brisure de symétrie : des phases distinctes sont supposées correspondre à des symétries différentes, et donc une transition de phase est toujours associée à un changement de symétrie. Cette théorie peut par exemple expliquer la différence entre un liquide et un solide, un cristal en particulier. Un liquide est une phase désordonnée avec une symétrie de translation continue tandis qu'un cristal est une phase ordonnée qui présente une symétrie de translation discrète (voir Fig. 1.1) : dans un liquide les atomes sont distribués aléatoirement et donc le liquide reste le même si on applique une translation arbitraire aux atomes dans l'espace, alors que dans un cristal les atomes sont organisés en un réseau qui ne reste inchangé que si on applique une translation spécifique (correspondant à un entier multiplié par la constante de réseau). En revanche, les phases topologiques ne diffèrent pas toujours par leurs symétries mais sont plutôt différenciées par la valeur d'un paramètre spécial appelé invariant topologique, ou nombre quantique topologique. 1Un certain nombre de matériaux topologiques présentent des états spéciaux à leur frontière, qu'il s'agisse de coins, de bords, de charnières ou de surfaces. Par exemple, les isolants topologiques possèdent des états de bord métalliques (Kane and Mele, 2005), et les supraconducteurs topologiques présentent des états liés de Majorana. 2 En raison de la nature topologique de ces matériaux, leurs états de bord sont particulièrement résistants aux perturbations et au désordre : ils sont « protégés » par la topologie. La robustesse des états de bord topologiques a été une motivation importante pour l'étude des matériaux topologiques, car ils peuvent avoir des applications très étendues : les états métalliques pourraient transporter du courant sans dissipation, les états protégés en général sont de bons candidats pour les qubits... Même si les matériaux topologiques sont ce qui vient le plus souvent à l'esprit lorsqu'on pense aux "états de bord" de nos jours, d'autres systèmes possèdent également des états qui sont liés à leurs bords, ou dont l'existence est directement liée à la présence de bords. Par exemple, le graphène, l'un des matériaux 2D les plus connus, héberge des états localisés pour certains types de bords. Il peut être particulièrement intéressant d'étudier ces états, par exemple pour voir comment ils se transforment sous l'application de certains champs, ou en fonction des interactions en jeu (Lado et al., 2015).

L'autre domaine principal qui se rapporte à cette thèse est l'étude des effets d'impuretés. Les impuretés ont toujours joué un rôle important dans la matière condensée, car il est tout simplement impossible d'étudier un matériau (expérimentalement) sans qu'il présente aucun défaut. Bien sûr au fil des années, de nouveaux procédés de fabrication nous ont permis de fabriquer des matériaux de plus en plus purs, mais les impuretés restent importantes et peuvent grandement affecter les résultats d'une expérience. Une autre raison fondamentale pour laquelle les impuretés sont étudiées, c'est parce que beaucoup de phénomènes physiques intéressants ne se produisent en réalité qu'en présence de défauts. Un exemple est l'effet Kondo : la diffusion d'électrons de conduction dans un métal due à des impuretés magnétiques peut affecter la résistivité électrique à basse température, ce qui dans certains cas peut produire des fermions lourds, c'est-à-dire des électrons de masse effective élevée (Hewson and Kondo, 2009).

La méthode que nous présentons ici, la "méthode des impuretés" ou "technique des impuretés", vise à apporter une nouvelle façon d'étudier les phénomènes de bord, en utilisant les impuretés. La plupart du temps, la physique des bords est étudiée théoriquement en simulant un système fini avec les bords souhaitées. Dans la méthode des impuretés, nous considérons plutôt un système infini et "créons" une frontière avec des impuretés. Un exemple est donné dans la Fig. 1.2 pour un système 2D. Dans ce cas, nous ajoutons une ligne d'impuretés qui, si elle est suffisamment forte, agira comme un mur et "coupera" le système en deux. On se retrouve ainsi avec un bord effectif de chaque côté de la ligne d'impureté, où des états de bord peuvent se former si le système de base les supporte.

Mathématiquement, nous utilisons des modèles tight-binding pour décrire les systèmes infinis, et utilisons le formalisme des fonctions de Green et des matrices T, couramment utilisés pour étudier la diffusion, pour tenir compte des impuretés et calculer leur effet. La méthode a l'avantage d'être plus analytique que numérique : le calcul suit une série d'équations, et la seule partie qui puisse réellement être considérée comme numérique consiste à évaluer quelques intégrales. Cette dernière partie peut même se faire de manière "légèrement" numérique en approchant simplement les intégrales par des sommes de Riemann ou avec la règle trapézoïdale. 3 Pour des systèmes suffisamment simples, des résultats purement analytiques peuvent même être obtenus, comme nous l'avons montré par exemple dans Kaladzhyan et al. 2020. Il existe bien sûr d'autres méthodes pour étudier les états de bord : diagonalisation exacte des hamiltoniens tight-binding pour des systèmes finis, méthodes itératives pour calculer les fonctions de Green aux limites (voir Lopez Sancho et al. 1985et Peng et al. 2017 pour une application récente), utiliser la correspondance bulk-bords aux limites (Delplace et al., 2011)... Le but de ce travail n'est pas de comparer la technique des impuretés à ces autres méthodes4 , mais plutôt d'explorer les possibilités qu'elle peut offrir : que peut-on en obtenir ? dans quelle mesure ? pouvons-nous l'utiliser pour trouver de nouveaux résultats? Afin de répondre à ces questions, nous avons appliqué la méthode des impuretés à un certain nombre de systèmes, qui correspondent plus ou moins aux différents chapitres de ce manuscrit.

Nous commencerons au chapitre 2 en introduisant quelques caractéristiques de base des états de bord à travers deux exemples bien connus : le graphène et le modèle SSH. Dans ce chapitre d'introduction, nous n'utiliserons pas la technique des impuretés, mais effectuerons seulement quelques calculs simples en utilisant la diagonalisation exacte. L'objectif de ce chapitre est double : présenter les états de bord dans quelques systèmes simples, et introduire des concepts qui seront utilisés dans les chapitres suivants, à savoir le modèle de base du graphène et quelques notions de topologie.

Ensuite, au chapitre 3, nous donnons les détails de la technique des impuretés. Pour cela, nous devons introduire les différentes fonctions de Green pertinentes, ainsi que le formalisme de la matrice T. Ce chapitre contient la plupart des "maths", et nous y ferons souvent référence dans les autres chapitres afin de ne pas réécrire les mêmes équations pour chaque système.

Dans le chapitre 4 nous présentons notre première application de la technique des impuretés. Nous nous concentrons sur deux types de systèmes topologiques : le modèle Kane-Mele, un isolant topologique 2D, et des modèles de semi-métaux de Weyl, qui vivent en 3D. Pour chaque modèle, nous montrons que nous pouvons récupérer les états de bord caractéristiques, avec les propriétés attendues, en comparant la fonction spectrale obtenue par la technique des impuretés avec le spectre d'énergie calculé en utilisant la diagonalisation exacte. Nous montrons ainsi que la technique peut fonctionner correctement aussi bien en 2D qu'en 3D, le cas 1D ayant été prouvé dans de précédents travaux (Kaladzhyan and Bena, 2019).

Nous passons ensuite au chapitre 5 à un système non topologique : la jonction SNS. Ce système est composé de trois parties reliées entre elles : un supraconducteur, un métal normal et un autre supraconducteur. Cette jonction implique donc deux frontières, qui déterminent l'étendue de la région normale. En raison de la diffusion aux frontières, des états spéciaux apparaissent dans la région normale, appelés états liés d'Andreev. Nous montrons que nous pouvons récupérer ces états grâce à notre technique, cependant le processus est un peu différent de celui utilisé dans le chapitre précédent car il implique de nouvelles difficultés : il faut prendre en compte deux frontières, pas une seule, et ces frontières doivent permettre la diffusion particule-trou, un ingrédient principal dans la formation des états liés d'Andreev. Nous constatons que la technique des impuretés ne fonctionne que pour un certain type de jonction, mais cela reste un résultat intéressant car ce régime particulier n'est généralement étudié que par des méthodes numériques plus "lourdes". Dans le chapitre 6 nous introduisons un nouveau type d'application pour la technique des impuretés : l'étude de l'interférence des quasi-particules. Ce phénomène est dû à la diffusion des états de surface par les impuretés ; en quelque sorte on revient ainsi à l'autre origine de la technique : l'investigation des effets d'impureté. Pour ce faire, nous passons par deux étapes : d'abord nous appliquons la technique de l'impureté pour créer une surface et obtenir les états de surface, puis, en utilisant le même formalisme mathématique, nous ajoutons à cette surface une impureté singulière et observons ses effets. Nous prenons comme exemples les modèles semimétalliques de Weyl utilisés dans le chapitre 4, et comparons nos résultats avec ceux obtenus par les auteurs des modèles via des techniques plus conventionnelles. Nous constatons que nos résultats ne sont pas en accord avec ceux des auteurs mentionnés et soutenons que la méthode des impuretés devrait donner les "meilleurs" résultats, en partie parce qu'elle ne reproduit pas les problèmes aujourd'hui bien connus des techniques conventionnelles.

Nous terminons au chapitre 7 par une étude du graphène multicouche et du graphite, en mettant l'accent sur les effets de l'empilement et du trigonal warping. L'idée d'étudier ces systèmes est née d'une collaboration avec des équipes expérimentales, qui ont examiné les effets des impuretés trouvées dans leurs échantillons de graphite. Ce fut donc l'occasion d'exploiter tout ce que nous avions appris des précédentes applications de la technique des impuretés, en complétant par des comparaisons avec des résultats expérimentaux récents, et en concluant convenablement cette thèse.

Dans cette thèse, nous nous sommes intéressés à une nouvelle méthode pour étudier les phénomènes de bord et de surface. Tout d'abord, nous avons montré qu'elle recrée correctement les états de bord connus, tels que ceux caractéristiques des isolants topologiques et des semi-métaux de Weyl ; puis nous avons vu qu'elle peut modéliser des effets de bord plus compliqués, recréant par exemple des états liés d'Andreev dans certaines jonctions SNS ; et enfin, nous avons constaté que la technique constitue en fait une bonne base pour l'étude de l'interférence des quasi-particules, à la fois dans l'espace des moments et dans l'espace réel, en montrant des exemples sur des modèles de semi-métaux de Weyl et sur des systèmes de graphène et de graphite.

A travers l'étude de tous ces systèmes, nous avons observé plusieurs avantages à l'utilisation de cette technique. D'une part, elle pourrait être considérée comme semi-analytique, dans la mesure où elle est exacte jusqu'au calcul de certaines intégrales. Le codage de la technique est également assez simple : à partir d'une petite matrice hamiltonienne, il suffit de suivre les équations successives pour calculer les fonctions de Green et les matrices T, et la fonction spectrale résultante ou la densité d'états locale par exemple. La méthode fonctionne dans toutes les dimensions physiques (1D, 2D et 3D), et avec des systèmes topologiques et non topologiques. Sa plus grande force réside sans doute dans sa construction : elle permet de modéliser des bords en s'affranchissant de la nécessité de considérer des systèmes finis.

Pendant la rédaction de ce manuscrit, d'autres applications de la technique de l'impureté ont été poursuivies. Des résultats entièrement analytiques ont été obtenus pour les interférences de quasi-particules du graphène bicouche à basse énergie (Kaladzhyan et al., 2021b), ce qui permet en particulier de mieux comprendre les effets du trigonal warping. La méthode a également été utilisée pour étudier les chaînes de Shiba et leur topologie [START_REF] Sedlmayr | New tools to determine the topological character of shiba chains[END_REF], où elle peut ignorer les problèmes d'effet de taille finie, connus pour se produire dans les calculs tight-binding. Parmi les applications futures possibles, l'interférence de quasi-particules semble être un choix prometteur, car elle peut prendre en compte des contributions importantes qui sont supprimées lors de l'utilisation des approximations habituelles de JDOS et SSP.

The work presented in this manuscript is centered around a certain technique, which we will often refer to as the impurity technique. This technique arose at the crossroads between two different domains in condensed matter physics: the study of edge and surface states, and the investigation of impurity effects.

In recent years, there has been a lot of regained interest in edge and surface phenomena due to the rise of topological condensed matter. This modern field, born in the 80s with the notable discovery of the quantum Hall effect (Klitzing et al., 1980), focuses on phases of matter (and their associated unique phenomena) which are not explained by the Landau theory of phase transitions. The Landau theory is based on symmetry breaking: distinct phases are assumed to correspond to different symmetries, and so a phase transition is always associated with a change in symmetry. This theory can for example explain the difference between a liquid and a solid, a crystal in particular. A liquid is a disordered phase with a continuous translation symmetry while a crystal is an ordered phase which exhibits a discrete translation symmetry (see Fig. 1.1): in a liquid the atoms are randomly distributed and so the liquid remains the same if we apply an arbitrary translation to the atoms in space, whereas in a crystal the atoms are organized into a lattice which remains unchanged only if we apply a specific translation (corresponding to an integer times the lattice constant). In contrast, topological phases don't always differ by their symmetries but are instead differentiated by the value of a special parameter called the topological invariant, or topological quantum number. 1A number of topological materials show special states at their boundary, be it corners, edges, hinges or surfaces. For example, topological insulators possess metallic edge states (Kane and Mele, 2005), and topological superconductors exhibit Majorana bound states. 2 Because of the topological nature of these materials, their edge states are particularly resistant to perturbations and disorder: they are "protected" by the topology. The robustness of topological edge states has been an important motivation to the study of topological materials, as they can have extensive applications: metallic states could carry current without dissipation, protected states in general are great candidates for qubits... Even though topological materials are what often come to mind nowadays when thinking "edge states", other systems also possess states which are bound to their edges, or which existence is directly linked to the presence of edges. For example graphene, one of the most well-known 2D materials, hosts localized states for certain types of edges. It can be particularly interesting to study these states, for instance to see how they transform under the application of certain fields, or depending on the interactions at play (Lado et al., 2015).

The other main field which relates to this thesis, is the study of impurity effects. Impurities have always played an important role in condensed matter, because it is simply impossible to study a material (experimentally) without it presenting any defects. Of course over the years, new manufacturing processes have allowed us to make purer and purer materials, but impurities remain important and can greatly affect the results of an experiment. Another fundamental reason why impurities are studied, is because a lot of interesting physics actually only occur in the presence of defects. One example is the Kondo effect: scattering of conduction electrons in a metal due to magnetic impurities can affect electrical resistivity at low temperature, which in some cases can produce heavy fermions, electrons with a high effective mass (Hewson and Kondo, 2009).

The method we present here, the "impurity method" or "impurity technique", aims to bring a new way of studying edge phenomena, using impurities. Most of the time, edge physics is studied theoretically by simulating a finite system with the desired edges. In the impurity method, we instead consider an infinite system, and "create" a boundary with impurities. An example is given in Fig. 1.2 for a 2D system. In this case we add a line of impurities which, if strong enough, will act as a wall and "cut" the system in two. We thus end up with an effective edge on each side of the impurity line, where edge states can form if the base system supports them.

Mathematically, we use tight-binding models to describe the infinite systems, and use the formalism of Green's functions and T-matrices, commonly used to study scattering, to account for the impurities and compute their effect. The method has the advantage of being more analytical than numerical: the computation follows a set of equations, and the only part which could really be considered numerical consists in evaluating some integrals. This last part can even be done in a "light" numerical way by simply approximating the integrals with Riemann sums or with the Figure 1.2 Representation of the impurity technique on a 2D system: to an infinite system (shown in grey) we add a line of impurities (shown in red), effectively cutting the system in two. Each subsystem thus possesses an effective edge, where edge states can form (in blue). trapezoidal rule. 3 For not too complicated systems, purely analytical results can even be obtained, as we showed for example in Kaladzhyan et al. 2020. Of course, other methods to study edge states exist: exact diagonalization of tight-binding Hamiltonians for finite systems, iterative methods to compute boundary Green's functions (Lopez Sancho et al. 1985 andPeng et al. 2017 for a recent application), using the bulk-boundary correspondence (Delplace et al., 2011)... The aim of this work is not to compare the impurity technique with these other methods4 , but rather to explore the possibilities that it offers: what can we obtain from it? To what extent? Can we use it to find new results?

As a way to answer those questions, we applied the impurity method to a number of systems, which more or less correspond to the different chapters of this manuscript.

We will start in Chapter 2 by introducing some basic features of edge states through two wellknown examples: graphene and the SSH model. In this introductory chapter, we will not use the impurity technique, but only carry out some simple calculations using exact diagonalization. The aim of this chapter is two-fold: to present edge states in some simple systems, and to introduce concepts which will be used in the following chapters, namely the basic graphene model and some notions of topology.

Then in Chapter 3 we give the details of the impurity technique. For this, we must introduce the different relevant Green's functions, as well as the T-matrix formalism. This chapter contains most of the "maths", and we will often refer to it in the other chapters so as not to rewrite the same equations for each system.

In Chapter 4 we present our first application of the impurity technique. We focus on two kinds of topological systems: the Kane-Mele model, a 2D topological insulator, and Weyl semimetal models, which live in 3D. For each model we show that we can recover the characteristic edge states, with the expected properties, by comparing the spectral function obtained through the impurity technique with the energy spectrum computed using exact diagonalization. We thus show that the technique can function correctly in both 2D and 3D, the 1D case having been proven in a previous work (Kaladzhyan and Bena, 2019).

We then switch in Chapter 5 to a non-topological system: the SNS junction. This system is made of three parts, joined together: a superconductor, a normal metal and another superconductor. This junction thus involves two boundaries, which determine the extent of the normal region.

Due to scattering at the boundaries, special states appear in the normal region, which are called Andreev bound states. We show that we can recover these states using our technique, however the process is a little different from the one used in the previous chapter as it involves new difficulties: we must take into account two boundaries, not just one, and these boundaries must allow particle-hole scattering, a main ingredient in the formation of Andreev bound states. We find that the impurity technique only works for a certain kind of junction, but this is still an interesting result as that particular regime is usually only studied through more "heavy" numerical methods.

In Chapter 6 we introduce a new type of application for the impurity technique: the study of quasiparticle interference. This phenomenon is due to scattering of surface states by impurities; in a way we thus come back to the other origin of the technique: the investigation of impurity effects.

To do this we go through two steps: first we apply the impurity technique to create a surface and obtain the surface states, then, using the same mathematical formalism, we add a singular impurity to that surface and observe its effects. We take as examples the Weyl semimetal models used in Chapter 4, and compare our results with those obtained by the models' authors through more conventional techniques. We find that our results do not agree with those of the mentioned authors, and argue that the impurity method should give the "better" results, in part because it does not reproduce well-known issues from the conventional techniques.

We end in Chapter 7 with a study of multi-layer graphene and graphite, with an emphasis on the effects of stacking and trigonal warping. The idea of investigating these systems originates from a collaboration with experimental teams, who examined the effects of impurities found in their graphite samples. It was thus an opportunity to use all we had learned from the previous applications of the impurity technique, completing with comparisons to recent experimental results, and bringing a fitting conclusion to this thesis. In this thesis, we focus on systems in which the presence of edges causes new states to appear. One can think of edge states, i.e. states which are localized at an edge, with a probability distribution decreasing exponentially when moving away from the edge. But delocalized states can also be linked to the presence of edges, e.g. the Andreev bound states in Superconductor-Normal-Superconductor (SNS) junctions which are discussed in Chapter 5.

Edge states can be separated into two categories: topological states, and non-topological states. Indeed a lot of edge states examples can be found in topological systems, and most of the systems presented in this manuscript are topological. For this reason, in this chapter we give two simple examples of edge states, one from a non-topological system, graphene (in its most basic form), and one from a topological system, the Su-Schrieffer-Heeger (SSH) chain. Through these examples we want to display the main characteristics one can look for when studying edge states, and with the SSH chain also briefly introduce the concept of topology in condensed matter.

| Simple graphene: a non-topological system with edge states

Graphene is one of the most, if not the most, well-known 2D system in condensed matter physics.

It is a single sheet of carbon atoms arranged in a honeycomb lattice, as shown in Fig. 2.1.

The honeycomb lattice has a unit cell of two sites, usually labeled A and B. The tight-binding Hamiltonian with only nearest-neighbor hopping is written as:

H = t i j c † i,A c j,B + h.c. (2.1)
where c † i,α creates an electron at site i of the corresponding type α = A or B. This Hamiltonian can be Fourier transformed to obtain the Bloch Hamiltonian for bulk graphene:

H(k) = t       1 + 2 cos       √ 3 2 k y a       e -i 3 2 k x a       c † k,A c k,B + h.c. (2.2)
where a is the distance between two nearest neighbors, and

c † k,α = 1 √ N i e -ik•R i c † i,α (2.3)
with N the number of unit cells and R i the position of the cell containing the site i.

The Brillouin Zone (BZ) of graphene is a hexagon which corners are called K points (see Fig. 2.2). Out of the six K points, only two are independent, the other four can be obtained by translations involving the reciprocal lattice vectors. We'll call these two independent points K and K .

In Fig. 2.2 we show the band structure of graphene along a line in the BZ which crosses K and K . As we can see, the spectrum around the K points is linear. Such a dispersion is characteristic of massless Dirac fermions, studied in high-energy physics, which is why the K points of graphene are also called Dirac points. Let us now change the boundary condition in one direction from periodic to open. We now have an infinite ribbon of graphene, with two edges. Depending on the direction and the position at which the ribbon ends, we can obtain three different types of edges: zigzag, bearded or armchair (see Fig. 2.3). In the convention we have chosen for the lattice orientation, taking open boundary conditions (OBC) in x will result in either zigzag or bearded edges, while taking OBC in y will result in armchair edges. The Hamiltonian for a graphene ribbon can be obtained by taking the partial Fourier transform of Eq. (2.1) in the infinite direction. For a ribbon infinite along x we thus have:

-π -2π 3 0 2π 3 π k x a -π -4π 3 √ 3 -2π 3 √ 3 0 2π 3 √ 3 4π 3 √ 3 π k y a K K 0.1 0.2 K K -1.0 -0.5 0.0 0.5 1.0 E/t
H(k x ) = t Y c † k x ,Y,A c k x ,Y,B + e -i 3 2 k x a c † k x ,Y,A c k x ,Y- √ 3 2 a,B + c † k x ,Y,A c k x ,Y+ √ 3 2 a,B + h.c., (2.4)
where Y designates the cell position in the y direction, and for a ribbon infinite along y:

H(k y ) = t X       c † X,k y ,A c X,k y ,B + 2 cos       √ 3 2 k y a       c † X,k y ,A c X-3 2 a,k y ,B       + h.c., (2.5)
where X designates the cell position in the x direction. Using these two Hamiltonians, we can obtain the spectrum for a graphene ribbon with any of the three types of edges. The spectra for the four possible edge configurations are given in Fig. 2.4: a ribbon with two zigzag edges, two bearded edges, one zigzag and one bearded edge, or two armchair edges.

As we can see, flat bands appear in the spectrum at zero energy when we have either zigzag or bearded edges (or both), but not with armchair edges. These flat bands actually correspond to edge states, i.e. states localized on the edge. From Fig. 2.4a and 2.4b we observe that zigzag edge states exist for k y a ∈ I with I = -4π This shows that edge states can form in graphene, depending on the type of edge. The presence of zero energy states is often reminiscent of topological systems, however, the existence of topologically protected edge states does not depend on the edge considered. For this reason, graphene is categorized as a non-topological system. Graphene can however be turned into a topological system: one example consists in adding spin-orbit interactions to its Hamiltonian, creating a topological insulator refered to as the Kane-Mele model, which we will present in Chapter 4.

-π -π √ 3 0 π √ 3 π k y a -3 -2 -1 0 1 2 3 E/t (a) Two zigzag edges -π -π √ 3 0 π √ 3 π k y a -3 -2 -1 0 1 2 3 E/t (b) Two bearded edges -π -π √ 3 0 π √ 3 π k y a -3 -2 -1 0 1 2 3 E/t (c) One zigzag edge, one bearded edge -π -2π 3 0 2π 3 π k x a -3 -2 -1 0 1 2 3 E/t (d) Two armchair edges

| A brief introduction to topology: the Su-Schrieffer-Heeger chain

We now introduce the Su-Schrieffer-Heeger (SSH) chain (Su et al., 1979), a simple 1D system, in order to illustrate the main characteristics of topological systems.

t 1 t 2 Figure 2
.5 An SSH chain of 10 sites. The unit cell is made of two neighboring sites, with hopping t 1 (t 2 ) between sites belonging to the same (different) unit cell.

We give a picture of the SSH chain in Fig. 2.5. The system is made from unit cells of two sites, and two different nearest neighbor hoppings exist: t 1 which links two sites from a same unit cell, and t 2 for sites belonging to different unit cells.

Let us first consider the two extreme cases where one of the two hoppings is equal to zero. When t 2 = 0, only hoppings inside the unit cells remain. As we can see in Fig. 2.6, the sites pair

-t 1 t 1 E Figure 2
.6 SSH chain and its energy spectrum for t 2 = 0.

together and the eigenvalues are ±t 1 . Indeed in this case we have a block diagonal Hamiltonian, where each block is the 2 × 2 matrix:

0 t 1 t 1 0 , (2.6)
which eigenvalues are ±t 1 .

The t 1 = 0 case is similar, except that this time the two sites at the ends of the chain do not have partner sites to pair up with and become isolated (see Fig. 2.7). We thus obtain the additional eigenvalue 0 in the spectrum. In this case, the presence of zero energy states is evident and doesn't seem particularly interesting. Let us now turn the hopping t 1 back on. Energy spectra for t 2 = 1 and different values of t 1 are shown in Fig. 2.8. We now have a lot more eigenvalues than in the previous extreme cases, but the important thing is: zero-energy states are present, even though t 1 0.

In order to investigate these peculiar states, we plot their distribution along the chain in Fig. 2.9. As we can see, the zero-energy states are exponentially localized at the ends of the chain. We also see that they become less localized when t 1 increases, but we'll come back to that later.

These plots show that, just like in the case of t 1 = 0, we obtain two zero-energy states localized at the ends of the chain for different non-zero values of t 1 . However, these states do not appear for all values of t 1 : in Fig. 2.10 we show instances where the eigenvalue 0 is absent.

The reader may have noticed the difference between the values taken for t 1 in Fig. 2.8 and 2.10. Indeed, in the first plots we have taken positive values of t 1 which are inferior to t 2 , while in the latter plots t 1 is superior to t 2 . The presence or absence of zero-energy states in these instances can be explained by considering that the system is in different phases when 0 ≤ t 1 < t 2 and when 0 ≤ t 2 < t 1 . However, the phase transition at t 1 = t 2 cannot be explained by the Landau theory because no symmetry was broken in order to go from one phase to the other. We call these types of phases "topological phases".

According to the theory of topological order, the system is in the same phase when t 1 = 0 and when 0 < t 1 < t 2 , which explains why we see the zero-energy edge states in Fig. 2.8. For t 1 > t 2 , an important feature of the system changes: the value of the topological invariant. We will not explain here how to compute the topological invariant, as it was not studied in this thesis, but it can be computed from bulk attributes. This is important, as it means that this invariant is a bulk property, that does not depend on the boundary conditions of the system.

The topological invariant is a number that takes quantized values. In the case of the SSH model, it can only take integer values. That is usually the case in topological systems, although the invariant can also be a rational number or a set of numbers. When it is zero, it means that the system is in a so-called topologically trivial phase. These phases do not show topological effects: the system can appear as a simple insulator for example. The SSH chain is in a topologically trivial phase when 0 ≤ t 2 < t 1 . When the invariant has a non-zero value, the system is in a topologically non-trivial phase: topological effects appear, such as the zero-energy edge states of the SSH chain. The study of topological edge states has been very popular in condensed matter physics lately, as they are by nature robust to perturbations. Indeed these states are protected by certain symmetries of the system, and as long as perturbations respect those symmetries, they will not affect the topological states.

At the topological phase transition, in the case of the SSH chain we saw that this happens at t 1 = t 2 , the topological invariant is not defined. Indeed since it is quantized, with different values for different phases, it cannot change smoothly at the transition. This translates in the bulk spectrum by a closed gap. At a topological phase transition, the bulk spectrum cannot be gapped, otherwise the invariant cannot change. We show this in the case of the SSH chain in Fig. 2.11 where we plot the bulk spectrum across the phase transition: the spectrum is gapped, but closes at t 1 = t 2 .

We now move on to the core of this thesis, introducing the impurity technique and its corresponding mathematical formalism. -2

-1 0 1 2 E (a) t 1 = 1.2 -2 -1 0 1 2 E (b) t 1 = 1.5 -3 -2 -1 0 1 2 3 E (c) t 1 = 1.8 Figure 2
.10 Spectrum of an SSH chain of 100 sites for t 2 = 1 and other values of t 1 . In this chapter, we introduce the main mathematical tools that are used in the impurity method, namely the Green's functions and the T-matrix formalism. This introduction is not very general, and most of the equations that follow are only correct under certain hypotheses, the main one being the consideration of a quadratic Hamiltonian. For a more general or detailed introduction on Green's functions or T-matrices, as this chapter is only meant as a support to understand the constituents of the impurity technique, and does not include much in the way of proofs, we invite the reader to look to the many textbooks or lectures on the subject. Some references are given in the next sections.

0 π 2 π 3π 2 2π k -1.0 -0.5 0.0 0.5 1.0 E (a) t 1 = 0 0 π 2 π 3π 2 2π k -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 E (b) t 1 = 0.5 0 π 2 π 3π 2 2π k -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 E (c) t 1 = 0.8 0 π 2 π 3π 2 2π k -2 -1 0 1 2 E (d) t 1 = 1 0 π 2 π 3π 2 2π k -2 -1 0 1 2 E (e) t 1 = 1.2 0 π 2 π 3π 2 2π k -2 -1 0 1 2 E (f) t 1 = 1.5 Figure 2.

| Green's functions

When studying quantum systems, the Green's function can be a very helpful mathematical tool. Through it, one can obtain properties of a system without having to determine its eigenstates. Depending on the quantities of interest, this can help save a lot of time. For example, one can derive the density of states (DOS) of electrons for a given energy value, and thus recreate the energy spectrum. The Green's function can also be used to compute topological invariants for topological systems.

Here we introduce the Green's functions that were of interest to us, as well as the spectral function and its relation to the electron DOS. Most of the Green's functions' properties are given here without proof. More details can be found in books or courses on Green's functions [START_REF] Rickayzen | Green's Functions and Condensed Matter[END_REF][START_REF] Mahan | Many-Particle Physics[END_REF].

In what follows we use the Heisenberg representation where the time dependence of an operator A is defined by: A(t) = e iHt Ae -iHt (3.1) and its expectation value by:

A = 1 Z tr e -βH A (3.2)
with H the Hamiltonian in the grand-canonical ensemble, β = 1/k B T the thermodynamic beta and Z = tr e -βH the partition function.

We only studied systems with a quadratic Hamiltonian, i.e. where the electrons are considered to be free particles. For this reason, the notions developed here are less general and we consider Hamiltonians of the form:

H = αβ c † α h αβ c β (3.3)
where c † α creates a particle in the state α.

Here we present two types of Green's functions: the retarded Green's function and the Matsubara Green's function.

The retarded Green's function

We define the retarded Green's function as: As we only treat the case of free particles, we need only focus on the single-particle retarded Green's function:

G R t, t = -iθ(t -t ) A(t), B(t ) ζ (3.
g ab R (t, t ) = -iθ(t -t ) [ ψa (t), ψ † b (t )] ζ (3.5)
where ψ † a creates a particle in the state a. Physically, the single-particle retarded Green's function describes the behaviour of an extra particle added to the system. By developing the expectation value in Eq. (3.5), we get two terms ψa (t) ψ † b (t ) and ψ † b (t ) ψa (t) . The first one is the probability amplitude for finding a particle in the state a at time t given that a particle is created in the state b at time t . This term therefore describes the way a particle goes from state b to state a between the times t and t. Likewise, the second term describes a "hole" going from state a to state b between the times t and t .

One can also show that if the Hamiltonian does not depend explicitly on time, which is the case for the systems studied here, the retarded Green's function depends only on tt . Hence in what follows we write the single-particle retarded Green's function as g ab R (tt ) = g ab R (t, t ).

The retarded Green's function can give us the physical quantities we are interested in, such as the spectral function and the DOS. However it is not straight-forward to compute. This is why we use Matsubara Green's functions. The Matsubara Green's function is defined on imaginary times.

It is related to the retarded Green's function and, as we will see in the following section, has a simple expression for the case of free fermions.

The Matsubara Green's function

By analytical continuation with Eq. (3.1), we define the imaginary time evolution as: A(τ) = e Hτ Ae -Hτ .

(3.6)

The Matsubara Green's function is defined as:

G M (τ, τ ) = -T (A(τ)B(τ )) (3.7)
where T is the time-ordered operator, defined by:

T (A(τ)B(τ )) = θ(τ -τ )A(τ)B(τ ) + ζθ(τ -τ)B(τ )A(τ). (3.8)
Once again, since we only consider quadratic Hamiltonians, we only use single-particle Matsubara Green's functions, that we define as:

g ab M (τ, τ ) = -T ( ψa (τ) ψ † b (τ )) .
(3.9) Also, since the Hamiltonian does not depend explicitely on time, we have

g ab M (τ -τ ) = g ab M (τ, τ ).
One can show that:

g ab M (τ + β) = ζg ab M (τ). (3.10)
This means that we can write the Matsubara Green's function as a Fourier series:

g ab M (τ) = 1 β ω n e -iω n τ g ab M (iω n ) (3.11)
where the ω n are the Matsubara frequencies, given by e iω n β = ζ.

One can also show that the Fourier transform g ab R (ω) of the retarded Green's functions is related to g ab M (iω n ) by the analytical continuation:

g ab M (iω n → ω + i0 + ) = g ab R (ω). (3.12)
This direct relationship is very useful in the case of a quadratic Hamiltonian. Indeed one can show that for a general quadratic Hamiltonian (3.3), the Matsubara Green's function is a matrix of the same size as h and has the simple expression:

g M (iω n ) = 1 iω n -h . (3.13)
This means that if h is a small matrix, as is the case for infinite systems with only a few degrees of freedom, the Matsubara Green's function is easy to compute.

The spectral function and DOS

We define the spectral function as:

A(ω) = - 1 π Im {tr [G R (ω)]} (3.14)
where G R (ω) is the Fourier transform of the retarded Green's function. To understand its meaning, let us take the case of a quadratic Hamiltonian. By analytic continuation of Eq. (3.13) we have:

G R (ω) = 1 ω + i0 + -h . (3.15)
By diagonalizing h, we can define a Green's function for a specific state of eigenvalue E:

G R (E, ω) = 1 ω + i0 + -E (3.16)
which in turn gives us:

A(E, ω) = δ(ω -E).
(3.17)

In this case the spectral function is a delta function, which tells us that an excitation can only happen by adding a particle in a state of energy ω = E. This means that for a non-interacting system, the spectral function behaves like a density of states: it only takes non-zero values for values of ω corresponding to the energies available to the particles. When interactions are considered, the spectral function is no longer a delta function, but one can show that it is still similar to the density of states at a given energy.

| The T-matrix formalism

In this section, we introduce the T-matrix formalism. A T-matrix is a tool that can be used together with Green's functions to study the effects of disorder in electronic systems.

Adding an impurity to a system will cause electrons to scatter. But this scattering event can be very complex to describe. In some cases we may not be interested in the scattering itself, but only in its effect. The T-matrix is a way to do just that: it creates a link between the system without scattering, which we call the unperturbed system, and the system with scattering, which we call the perturbed system. These two systems, perturbed and unperturbed, are represented mathematically by Green's functions, while the scattering event is represented by the T-matrix.

In the brief introduction that follows we use a few assumptions, which are valid for all the systems we worked on, and simplify the equations as well as the understanding of the formalism (for a more general theory of T-matrices see for example [START_REF] Mahan | Many-Particle Physics[END_REF][START_REF] Bruus | Many-Body Quantum Theory in Condensed Matter Physics -An Introduction[END_REF].

First we assume that the unperturbed system is described by a quadratic Hamiltonian and is translationally invariant. Typically, we consider a Bloch Hamiltonian H(k), written as a small matrix e.g. 2 × 2 or 4 × 4.

Second we assume that the impurity is localized and is represented by a scattering potential resembling a delta function in space. For example, the impurity can be point-like, represented by a potential: V(r) = Uδ(rr ).

As mentioned above, we first consider an unperturbed system which is translationally invariant, described by a Bloch Hamiltonian H(k). We define the associated unperturbed Matsubara Green's function as:

G 0 (k, iω n ) = 1 iω n -H(k) . (3.18)
As explained in Sec. 3.1.2, we use the Matsubara Green's function because it has a simple expression for the systems we consider.

We now add an impurity to the system. For a localized impurity, the impurity scattering can be modeled using an impurity potential series. That is, the Green's function G of the perturbed system can be expanded in terms of the Green's function G 0 of the unperturbed system and the impurity potential V:

G = G 0 + G 0 VG 0 + G 0 VG 0 VG 0 + ... (3.19) = G 0 + G 0 (V + VG 0 V + ...)G 0 . (3.20)
For now, we omit the parameters k, iω n for simplicity. If we then define T as the infinite sum (V + VG 0 V + ...), we get a more compact expression for G:

G = G 0 + G 0 TG 0 . (3.21)
This T-matrix is the sum of the contributions corresponding to all-order impurity scattering processes. Again though, this method can only be used if the impurity is localized. For extended impurities, if the impurity scattering is weak, the Born approximation can be used instead. In that case only the first order term in the impurity potential series is considered, i.e. we simply have T = V.

Adding again the parameters k, iω n , the equation we use to compute the perturbed Green's function is:

G(k 1 , k 2 , iω n ) = G 0 (k 1 , iω n )δ(k 1 -k 2 ) + G 0 (k 1 , iω n )T (k 1 , k 2 , iω n )G 0 (k 2 , iω n ). (3.22)
Since the unperturbed system is translationally invariant, G 0 depends only on one momentum. However, adding the impurity breaks this invariance, which is why G depends on two distinct momenta k 1 and k 2 .

We can compute G 0 for a given Hamiltonian using Eq. (3.18), so the only thing needed to obtain the perturbed Green's function G from Eq. (3.22) is an expression for the T-matrix. By definition we can write T as:

T (k 1 , k 2 , iω n ) = V(k 1 , k 2 ) + k' V(k 1 , k')G 0 (k', iω n )V(k', k 2 ) + ... (3.23) = V(k 1 , k 2 ) + k' V(k 1 , k')G 0 (k', iω n )        V(k', k 2 ) + k" V(k', k")G 0 (k", iω n )V(k", k 2 ) + ...        (3.24) = V(k 1 , k 2 ) + k' V(k 1 , k')G 0 (k', iω n )T (k', k 2 , iω n ) (3.25)
which means that we need to solve Eq. (3.25) to obtain an expression for the T-matrix.

We consider an impurity localized around r = 0, represented by a delta function impurity potential: V(r) = Uδ(r). By Fourier transforming V, one can see that V is constant in momentum space, and simply equal to U. Then the solution to Eq. (3.25) is:

T (iω n ) = I -U dk V BZ G 0 (k, iω n ) -1 U (3.26)
where V BZ is the volume of the first Brillouin Zone. This result can be verified by plugging Eq. (3.26) into Eq. (3.25), taking V(k, k') = U.

| Using the T-matrix formalism to study edge states

In order to study edge states, or other phenomena linked to the presence of edges, we used this T-matrix formalism in a very specific way. The idea is to take an infinite system, and to create an edge using an impurity. Illustrations are given in Fig. 3.1 for 1D, 2D and 3D. Typically, when adding a localized impurity to a system, some states will form around the impurity. What we observed is that, when the impurity potential strength U goes to infinity, i.e. is much larger than the energy scales of the system, these localized states transform into the edge states one would obtain on "normal" edges. This technique thus presents a new way to obtain and study edge states. In each case, the infinite system is subjected to an impurity, which effectively splits it into two semi-infinite subsystems. If the system possesses edge states, these are recreated on each side of the impurity. These edge states are represented in blue in (a) and (b), or exist on the surface represented in red in (c).

This method offers several advantages compared to the usual ways of obtaining edge states. For one, we consider the system to be infinite. Usually, to study a system with edges, one needs to make it finite in at least one direction. In that case, k is no longer a good quantum number along that direction, and the Hamiltonian matrix becomes much bigger, with a size proportional to N × N with N the number of sites in the finite direction. Such a big matrix makes it impossible to obtain any analytical result, and also makes numerical caluculations slower, mostly when considering big systems with a lot of sites. With our technique, it is possible to obtain some analytical results, such as an expression for the Green's function, the energy dispersion of the edge states or even the edge states themselves. This depends of course on the complexity of the Hamiltonian: the difficulty lies in particular in computing the integral in Eq. (3.26). Numerically, there is a clear advantage since the computation time no longer depends on the size of the system. It is also easier to implement compared to the case of a big Hamiltonian, where one must remember to take into account the open boundary conditions.

Another point worth mentioning is that even the numerical results are not too numerical. The technique is exact up to the integration in Eq. (3.26), which is usually impossible to compute exactly1 .

Point, line and surface Green's functions

Before moving on to the application of this technique, we go over the equations which can be used for the different cases of 1D, 2D and 3D systems. Here we introduce the point, line or surface Green's functions and their corresponding spectral function.

In 1D, the Hamiltonian only depends on one momentum direction k. The unperturbed Green's function can be computed using Eq. (3.18), the T-matrix using Eq. (3.26) with V BZ = L BZ the length of the Brillouin Zone, and then we obtain the Green's function using Eq. (3.22) where T does not depend on k.

Since we are interested in what happens next to the impurity, we must take the Green's function at a corresponding point in space. Typically, for a tight-binding model where we add an impurity potential on the site r = 0, we take the Green's function at the nearest neighbor site i.e. r = a or r = -a, with a the distance between two nearest neighbors. This means we do a Fourier transform of G:

G p (r = ±a, iω n ) = dk 1 dk 2 L 2 BZ G(k 1 , k 2 , iω n )e ±ik 1 a e ∓ik 2 a (3.27)
where G p labels the real space Green's function at a given point. Using G p we can in turn compute the spectral function at r = ±a:

A(r = ±a, E) = - 1 π Im tr G p (r = ±a, E) (3.28)
where we write A and G p as a function of E instead of iω n , such that iω n = E + iδ, δ → 0 + . The positive infinitesimal shift of energy, +iδ, corresponds to an inverse quasiparticle lifetime and is generally related to the width of the energy levels. We often use E as a parameter for a more evident correspondence with energy.

If the system exhibits edge states, this quantity is non zero at energy values E corresponding to the edge states energy. In contrast, if we evaluate A(r, E) with r 0, i.e. far from the impurity, it has zero value at the same energies. Unless of course bulk states exist at those energies, making the edge states indistinguishable from the bulk states in the spectrum. But that is not the case for the systems we considered in this thesis.

In 2D the situation is a bit different: instead of adding a point-like impurity, we add a line-like one. This means that translational invariance is only broken in one direction. For example, for an impurity along the y direction, at x = 0, the impurity potential takes the form V(r) = Uδ(x).

Eq. (3.26) is then modified, such that we only integrate over the direction in which the translational invariance is broken:

T (k 1y , k 2y , iω n ) = δ(k 1y -k 2y ) I -U dk x L BZ,x G 0 (k x , k 1y , iω n ) -1 U (3.29)
where L BZ,x is the length of the BZ along the x direction. This time, the T-matrix depends on k y , but contains the delta term δ(k 1yk 2y ). This term characterizes the conservation of momentum along y, since translational invariance is not broken by the impurity in that direction.

In a similar way to the 1D case, we can compute the real space Green's function G l on a line close to the impurity:

G l (x = ±a, k y , iω n ) = dk 1x dk 2x L 2 BZ,x G(k 1x , k y , k 2x , k y , iω n )e ±ik 1x a e ∓ik 2x a . (3.30)
Again, since momentum is conserved along y, G must be taken at momenta k 1 , k 2 with the same k y component. The spectral function on the x = ±a line is then written as:

A(x = ±a, k y , E) = - 1 π Im tr G l (x = ±a, k y , E) . (3.31)
The 3D case is similar to the 2D one, with translational invariance broken in only one direction and conserved in the two others. If we again write the impurity potential as V(r) = Uδ(x), we have for the T-matrix, the real space surface Green's function and the surface spectral function:

T (k 1y , k 1z , k 2y , k 2z , iω n ) = δ(k 1y -k 2y )δ(k 1z -k 2z ) I -U dk x L BZ,x G 0 (k x , k 1y , k 1z , iω n ) -1 U (3.32) G s (x = ±a, k y , k z , iω n ) = dk 1x dk 2x L 2 BZ,x G(k 1x , k y , k z , k 2x , k y , k z , iω n )e ±ik 1x a e ∓ik 2x a . (3.33) A(x = ±a, k y , k z , E) = - 1 π Im tr G s (x = ±a, k y , k z , E) . (3.34)

Numerical computation of the point, line or surface Green's function

Due to the delta function in the formula for the perturbed Green's function G (3.22), it is not convenient to compute it numerically. Since we are only interested in the resulting point, line or surface Green's function, we decompose its integral so as to remove the delta function. For instance, for the point Green's function we have:

G p (r = ±a, iω n ) = dk 1 dk 2 L 2 BZ [G 0 (k 1 , iω n )δ(k 1 -k 2 ) + G 0 (k 1 , iω n )T (k 1 , k 2 , iω n )G 0 (k 2 , iω n )] e ±ik 1 a e ∓ik 2 a = dk 1 L BZ G 0 (k 1 , iω n ) + dk 1 dk 2 L 2 BZ G 0 (k 1 , iω n )T (k 1 , k 2 , iω n )G 0 (k 2 , iω n )e ±ik 1 a e ∓ik 2 a .
For a delta function impurity localized around r = 0, we can use the T-matrix formula of Eq. (3.26), which does not depend on k 1 and k 2 , thus further simplifying the equations:

G p (r = ±a, iω n ) = dk 1 L BZ G 0 (k 1 , iω n ) + dk 1 dk 2 L 2 BZ G 0 (k 1 , iω n )T (iω n )G 0 (k 2 , iω n )e ±ik 1 a e ∓ik 2 a = dk 1 L BZ G 0 (k 1 , iω n ) + dk 1 L BZ G 0 (k 1 , iω n )e ±ik 1 a T (iω n ) dk 2 L BZ G 0 (k 2 , iω n )e ∓ik 2 a = B(iω n ) + ∆(iω n )
where B(iω n ) corresponds to the "bulk" contribution to G p , involving only the unperturbed bulk system, while ∆(iω n ) embodies the impurity's effect as it corresponds to the perturbed part of the Green's function.

In the same way, the line and surface Green's functions can be split into two similar contributions B and ∆.

Average correction to the spectral function

In some cases, such as the study of the Kane-Mele model in Sec. 4.1, it can be more convenient to compute the average correction to the spectral function. This quantity is actually equal to the average correction to the number of electronic states at given momentum and energy. It therefore allows us to visualize the impurity-induced states more easily, as it takes positive values for those states only. We demonstrate this here in the 2D case but the computation is similar in other dimensions.

In 2D, we consider a line impurity proportional to δ(x). Since the spatial translational invariance along y is not broken, k y is still a good quantum number. The correction to the number of available electronic states at position x and having momentum k y is then given by:

δρ(x, k y , E) = - 1 π Im tr δG(x, x; k y ; E), (3.35)
where δG corresponds to the perturbed part of the Green's function.

Averaging this over the direction perpendicular to the impurity, for which translational invariance was broken, we obtain:

δN(k y , E) ≡ dx δρ(x, k y , E) = - 1 π Im tr dx G 0 (x, k y , E)T (k y , E)G 0 (-x, k y , E) = - 1 π Im tr dx dk x L BZ,x dk x L BZ,x e ik x x e -ik x x G 0 (k x , k y , E)T (k y , E)G 0 (k x , k y , E) = - 1 π Im tr dk x L BZ,x G 0 (k x , k y , E)T (k y , E)G 0 (k x , k y , E) = - 1 π Im tr dk x L BZ,x δG(k x , k y ; k x , k y ; E) ≡ dk x L BZ,x δA(k x , k y , E), (3.36)
where δA(k x , k y , E) is the correction to the perturbed spectral function in momentum space.

We now have all the tools we need to study edge related phenomena in 1D, 2D and 3D systems. The following chapters are dedicated to the results obtained using the impurity method on a variety of examples. The impurity method allows us to simulate edges and surfaces from any bulk model. In this chapter we apply it to study several systems which exhibit edge states or surface states. We first investigate a 2D topological insulator: the Kane-Mele model. We will explain what a topological insulator is, obtain the edge states of the Kane-Mele model from a usual tight-binding calculation, and then show that we can reproduce these states with our method. After that, we move on to Weyl semimetals, 3D topological systems which display surface states known as Fermi arcs. As we did for the topological insulator, we will give a general introduction to Weyl semimetals, then study two given models and prove that our technique can recover their Fermi arcs.

| A topological insulator: the Kane-Mele model

What is a topological insulator?

Simply put, a topological insulator is an insulator which displays symmetry-protected conducting edge states (in the case of 2D systems) or surface states (in 3D systems). Similar edge states were first discovered through the integer quantum Hall effect (IQHE) (Klitzing et al., 1980) : in 2D systems subject to low temperatures and strong magnetic fields, the Hall conductance takes the quantized values νe 2 /h, where e is the elementary charge, h is the Planck constant and ν is an integer. Depending on the value of ν, a certain number of conducting edge states appear. It was understood later on [START_REF] Hatsugai | [END_REF] that ν is a topological invariant which, in particular, counts the number of conducting channels at the edges. Afterwards, models which reproduce the phenomena found in the IQHE, but without magnetic field, were introduced (Haldane, 1988;Kane and Mele, 2005;Bernevig et al., 2006), then discovered experimentally, both in 2D (König et al., 2007) and in 3D (Hsieh et al., 2008). They are now commonly referred to as "topological insulators".

We will give some more explicit details in the next section by taking the example of the Kane-Mele model (Kane and Mele, 2005). More general information can be found in books and reviews on the subject (for example see Andrei Bernevig and Hughes 2013). The Kane-Mele model describes a topological insulator created from graphene: we take the simple tight-binding model of graphene, described in 2.1, and add a next-nearest-neighbor (n.n.n.) term involving spin-orbit interactions. The Hamiltonian is written as:

The Kane-Mele model

H = i j ,α tc † i,α c j,α + i j ,α,β it 2 ν i j s z αβ c † i,α c j,β , (4.1)
where c † i,α creates an electron on the lattice site i, with spin α =↑, ↓ and s z αβ is the αβ component of the z Pauli matrix. In 2.1 we labeled the operators with the sublattice indices A, B. Here we can refrain from using them as the two sublattices are equivalent, and it is clear that i j relates to sites from different sublatttices while i j points to sites belonging to the same sublattice.

The parameter ν i j is a spin-dependent amplitude defined by the orientation of the n.n.n. hopping direction, as depicted in Fig. 4.1: a n.n.n. hopping consists in two successive n.n. hoppings, and depending on the orientation of the second n.n. hopping with respect to the first one, ν i j = +1 (if the electron hops to the left) or ν i j = -1 (if the electron hops to the right). To preserve complex conjugation, we have ν i j = -ν ji .

We plot the bulk energy spectrum in Fig. 4.2 for t = 1 and t 2 = 0.03. As we can see, the new term in the Hamiltonian opens a gap at the Dirac points between the valence band and the conduction band: the bulk is now insulating. We now consider an infinite ribbon, so as to have a system with edges. The energy spectrum for a ribbon with zigzag edges is shown in Fig. 4.3. As we can see, dispersing in-gap states appear in the same k intervals as the regular non-dispersing zigzag edge states of graphene (see Fig. 2.4). These states disperse between two inequivalent K points, and connect the valence and the conduction bands: they are conducting. In the intervals where they are defined, these states are localized on the zigzag edges, and start to delocalize into the bulk when approaching the Dirac points. In Sec. 2.1, we pointed out that the edge states of normal graphene are not topological because their appearance depends on the type of edge considered: they exist for zigzag and bearded edges, but not for armchair edges. In Fig. 4.4 we show the ribbon bands for the Kane-Mele model with the other types of edges. As we can see, the behavior observed with the zigzag edges is reproduced in the case of bearded edges: dispersing edge states appear in the same intervals as the usual bearded edge states. For armchair edges we now also see dispersing edge states. Unlike the zigzag and bearded ones, they do not disperse between Dirac points of different momentum, but they still connect the conduction and valence bands. This respects the previously mentioned condition on topological edge states: their presence is independent of the boundary considered.

For the Kane-Mele model, the topological invariant is the Z 2 invariant: it can only take the value 0 or 1. As we made explicit for the SSH model in Sec. 2.2, the invariant can only change if the bulk gap closes. This can be achieved for example with t 2 = 0, where we recover the normal graphene linear dispersion at low energy. However, this does not correspond to a topological phase transition: when the gap reopens for t 2 < 0, the dispersing edge states are again present. In order to change the topological invariant, one must consider other terms such as a Rashba term and a staggered sublattice potential (Kane and Mele, 2005). Nevertheless, since we are only interested in reproducing the edge states found in this topological phase, we will not take these additional terms into account.

-π -π √ 3 0 π √ 3 π k y a -3 -2 -1 0 1 2 3 E (a) Two bearded edges -π -π √ 3 0 π √ 3 π k y a -3 -2 -1 0 1 2 3 E (

Impurity method

x y A B u.c. Conversely, if the u.c. vector is parallel to the x axis, we get two effective bearded edges.

We will now show that by using our impurity scheme, we can recover the edge states without considering actual edges. The Kane-Mele model is defined on a 2D lattice, therefore, as explained in Sec. 3.3, we must introduce a line impurity in the infinite system. By setting the impurity potential to infinity, we should then reproduce the formation of the edge states.

The starting point is to rewrite the Hamiltonian in Eq. (4.1) in momentum space. In order to compare our results with the spectrum of Fig. 4.3, we must make sure that the impurity will recreate zigzag edge states. This can be done by defining an impurity which covers a whole unit cell. Depending on the unit cell considered however, we will create either two effective zigzag edges or two effective bearded edges, as shown in Fig. 4.5. Let us focus on the zigzag case for now.

We choose the basis (c A↑ , c A↓ , c B↑ , c B↓ ), where c iσ is an electron operator with spin σ = ↑/↓ on the sublattice i = A/B. The Kane-Mele Hamiltonian is then expressed as:

H(k) =                  h NNN 0 h NN 0 0 -h NNN 0 h NN h * NN 0 -h NNN 0 0 h * NN 0 h NNN                  , (4.2)
with

h NN = t 1 + e i √ 3k y a + e i √ 3 2 k y a e i 3 2 k x a , h NNN = 2t 2       sin √ 3k y a -2 cos 3 2 k x a sin       √ 3 2 k y a            
being the nearest-neighbor and the next-nearest-neighbor terms, with a the distance between nearest sites and * simply denoting the complex conjugation.

We then define the impurity matrix as V = UI, where I is the identity matrix, and U is the impurity strength. This represents an impurity localized on both sites of the unit cell, with no discrimination between different spins.

Following the steps described in Sec. 3.3, we now compute the average correction to the spectral function. The unperturbed Green's function is defined as:

G 0 (k, E) = [E + i0 + -H(k)] -1 , (4.3)
and the T-matrix for the line impurity is:

T (k 1y , k 2y , E) = δ(k 1y -k 2y ) I -V dk x L BZ,x G 0 (k x , k 1y , E) -1 V. (4.4)
For a square lattice of length a, we would simply integrate k x from -π/a to π/a, and would have L BZ,x = 2π/a. However in this case we have a honeycomb lattice, thus the limits of integration are not so simply defined. In App. 4.A we show that the integration must be performed from -2π/(3a) to 2π/(3a), with L BZ,x = 4π/(3a).

The average correction to the spectral function is then obtained through the equation:

δN(k y , E) = - 1 π dk x L BZ,x Im{tr[G 0 (k, E) T k y , k y , E G 0 (k, E)]}. (4.5)
In Fig. 4.6 we plot δN for U = 1, 10 and 100. To see the evolution of the impurity states more clearly, we plot in a reduced energy and momentum range in the first subfigures. As proven in Sec. 3.3.2, the average correction to the spectral function represents the correction to the number of available states at given E and k. We thus observe that the impurity changes some bulk states (δN < 0) into impurity states (δN > 0). For a weak impurity i.e. U = t = 1, the impurity states appear at distinct energies, still concentrated mostly out of the gap. When increasing U these states gradually evolve into edge states. Indeed for U = 100 i.e. U t, t 2 , we recover subgap states with the same linear dispersion as the ones derived via exact diagonalization (Fig. 4.3).

In Sec. 3.3 we also introduced the line spectral function as a tool to study states induced by a line impurity. The results obtained with this quantity are given in App. 4.B. They are similar to those obtained using the average correction to the spectral function, in that we recover the edge states dispersion for large U, but also show that the impurity states are localized around the impurity: they disappear when we move spatially away from the impurity.

In a similar way, we can reproduce the bearded edge states of Fig. 4.4a. To do this, we must define the unit cell as shown in Fig. 4.5b i.e. with the A and B sites at same y. In the Hamiltonian matrix of Eq. (4.2), only the nearest neighbor hopping term changes. We now have:

h NN = t       1 + 2 cos       √ 3 2 k y a       e -i 3 2 k x a       . (4.6)
The equations for G 0 , T and δN also remain unchanged and we obtain the results shown in Fig. 4.7a, which correspond to those obtained using exact diagonalization. The case of one bearded and one zigzag edges, as well as the case of two armchair edges, can also be reproduced as we show in Fig. 4.7b and 4.7c.

To obtain an effective bearded edge with an effective zigzag edge, we must define a line impurity localized on an odd number of sites along x (see Fig. 4.8a). There is however an added restriction due to the consideration of n.n.n. hopping. In order to effectively separate the system into two parts, and thus recreate edges, the impurity must restrict hoppings between the two subsystems. To restrict n.n.n. hopping, we must consider an impurity which extends to at least two sites along x: we can thus simulate both a bearded and a zigzag edge with an impurity localized over three sites. The easiest way to do this is to double the unit cell, thus considering 8×8 matrices instead of 4 × 4, and define V as a diagonal matrix with value U for three sites and value 0 for the last site. An alternative is to add a phase factor corresponding to the third site we must consider, which is out of the unit cell (see Eq. (4.24) in App. 4.C). However this phase factor depends on momentum, which means we would no longer be able to write V as independent of k, making the T-matrix equation more complicated. The armchair edges can also be reproduced using a doubled unit cell (see Fig. 4.8b). This time the line impurity breaks translational symmetry along y and the invariance remains only along x. The equations for G 0 , T and δN are obtained by integrating over k y instead of k x .

| Weyl semimetals and Fermi arcs

A gapless 3D topological system

We now wish to apply the impurity technique to a different topological system, this time in 3D: the Weyl semimetal. Similar to graphene, Weyl semimetals exhibit linear dispersion in the bulk energy spectrum at low energy: these systems are therefore gapless, in contrast to the topological systems we've introduced so far (the SSH chain in Sec. 2.2 and the Kane-Mele model in Sec. 4.1).

The points in momentum space where the bulk bands reach zero energy are called Weyl nodes or Weyl points, and always exist in pairs. The topological invariant for a Weyl semimetal is the Chern number, which takes integer values and is defined for each Weyl point. Again unlike the previous topological systems, the overall value of the topological invariant in Weyl semimetals is always zero: the Chern numbers for the two Weyl points of a pair are opposite, thus the net Chern number, defined as the sum of Chern numbers for each Weyl node, vanishes. What we are interested in is to reproduce the surface states appearing in Weyl semimetals: the so-called Fermi arcs. These states appear as open lines on the surface Brillouin Zone, which end at the projection of the two Weyl nodes of a pair on this surface Brillouin Zone (see Fig. 4.9). This fixes a condition on the surface to consider: the Weyl nodes must not project on the same point for this surface, otherwise no Fermi arc will be created. Experimentally, these Fermi arcs can be measured using angle-resolved photoemission spectroscopy (ARPES). An example of such a measurement is given in Fig. 4.10.

Application of the impurity method

We focus on two different models, described in Kourtis et al. 2016 andLau et al. 2017, which we denote by H 1 and H 2 , respectively. The Bloch Hamiltonians for these two systems are given by:

H 1,2 = k ψ † (k)H 1,2 (k)ψ(k), (4.7)
where ψ(k) = (c kA↑ , c kA↓ , c kB↑ , c kB↓ ) is a spinor with the index A/B denoting a generic unspecified orbital component, and the ↑ / ↓ the physical spin.

For the model in Kourtis et al. 2016 written in the basis above we have

H 1 (k) = g 1 (k)τ 1 σ 3 + g 2 (k)τ 2 σ 0 + g 3 (k)τ 3 σ 0 + g 0 (k)τ 0 σ 0 + βτ 2 σ 2 + α sin k y τ 1 σ 2 , (4.8)
where

g 0 (k) = 2d(2 -cos k x -cos k y ) g 1 (k) = a sin k x g 2 (k) = a sin k y g 3 (k) = m + t cos k z + 2b(2 -cos k x -cos k y ). (4.9)
and α, β are real parameters. The 2 × 2 identity matrices σ 0 /τ 0 and the Pauli matrices σ i /τ i , i = 1, 2, 3 act in the spin and the orbital spaces, correspondingly, and the multiplication of the σ and τ matrices indicates a tensor product.

We will consider two sets of parameters for this Hamiltonian: a) a = b = 1, t = -1, m = 0.5, d = 0.8, α = β = 0 and b) a = b = 1, t = -1.5, d = m = 0, β = 0.9, and α = 0.3. We can compute the position of the Weyl nodes by combing the BZ for zero energy points in the bulk spectrum. For both sets of parameters we find that two Weyl nodes of a pair are always at the same k x and k y . To see the Fermi arcs, we must thus consider a surface which is not parralel to the xy plane, otherwise the Weyl nodes' projections will superpose and no Fermi arc will be created on that surface.

We choose to take a surface perpendicular to the y axis. We can compute the spectrum using exact diagonalization by applying a partial inverse Fourier transform along y to the Hamiltonian (4.8). We thus end up with a system infinte along x and z and finite along y, with two infinite xz surfaces. In Fig. 4.11 we plot the minimum positive energy band obtained as a function of k x and k z , along with the Weyl points, for both sets of parameters. For the set a) we find two Weyl points, along with two Fermi arcs connecting them: one for k x < 0 and another for k x > 0. For the set b) we find two pairs of Weyl points, with one Fermi arc per pair. We can clearly see in Fig. 4.11b how the zero energy lines end at the Weyl nodes.

To recover the arcs with the impurity method, we add a plane-like impurity at y = 0 to the infinite system. For a strong impurity potential we thus create two effective xz surfaces: one at y = 1, the other at y = -1. Following the steps described in Sec. 3.3.1, we compute the surface Green's function and the corresponding surface spectral function at E = 0, for an impurity Using the surface Green's function, we can also easily describe the spin and orbital distribution of the Fermi arc states. Indeed, we have for the different spin and orbital components:

S x A (k x , k z ) = - 1 π Im G s 12 (k x , k z ) + G s 21 (k x , k z ) S x B (k x , k z ) = - 1 π Im G s 34 (k x , k z ) + G s 43 (k x , k z ) S y A (k x , k z ) = - 1 π Re G s 12 (k x , k z ) -G s 21 (k x , k z ) S y B (k x , k z ) = - 1 π Re G s 34 (k x , k z ) -G s 43 (k x , k z ) S z A (k x , k z ) = - 1 π Im G s 11 (k x , k z ) -G s 22 (k x , k z ) S z B (k x , k z ) = - 1 π Im G s 33 (k x , k z ) -G s 44 (k x , k z ) , (4.10)
where G s i j labels the i j component of the surface Green's function (we omit the energy dependence for the sake of brevity).

In We now perform the same analysis on the model from Lau et al. 2017: We now use the impurity method by considering a plane-like impurity at y = 0 on the infinite system. The surface spectral functions obtained at y = -1 are shown in Fig. 4.15. We successfully recover the Fermi arcs, as well as the Fermi pocket. We also compute the spin and orbital properties for this model (see Fig. 4.16). This time, for the set of parameters a), we find that the spins are not exactly opposite. This is most explicit when looking at the S z B component. We have shown that the impurity technique can be applied to topological systems such as topological insulators and Weyl semimetals. Using, respectively, a line-like and a surface-like impurity for 2D and 3D systems, we can derive the boundary states along with their properties, for example the spin and orbital distributions. Next, we will show that we can also use this technique to obtain states which arise due to the presence of boundaries, but are not localized at those boundaries, taking as example the Andreev bound states found in SNS junctions.

H 2 (k) = g 1 (k)τ 1 σ 3 + g 2 (k)τ 2 σ 0 + g 3 (k)τ 3 σ 0 + dτ 2 σ 3 + βτ 2 σ 2 + α sin k y τ 1 σ 2 + λ sin k z τ 0 σ 1 .
-π -π 2 0 π 2 π k x 0 π 2 π 3π 2 2π k z 1 2 3 -π -π 2 0 π 2 π k x -π -π 2
-π -π 2 0 π 2 π k x -π -π 2 0 π 2 π k z S z A -2 -1 0 1 2 -π -π 2 0 π 2 π k x -π -π 2 0 π 2 π k z S z B -2 -1 0

Appendix

4.A | Derivation of the T-matrix momentum limits of integration for a honeycomb lattice

While for a square or cubic lattice there is no subtlety concerning the integration limits, for the honeycomb lattice this is much more subtle. We will consider a line impurity localized on A atoms, but the conclusions do not depend on the type of impurity we consider. Let us start by writing down the form of the contribution of the impurity potential to the Hamiltonian, such as it has to be written in the continuum:

δH imp = dr V(r) (4.12)
where the integral is performed over the entire space, with

V(r) = Uρ(r), ρ(r) = j δ(r -R A j )a † j a j (4.13)
where the sum is performed over a line of lattice unit cells at x = 0, each denoted by the index j and R j = (0, j √ 3a). The a † j and a j operators describe the formation and annihilation of electrons at site j, they no longer live in the continuum but on the lattice and as such they are defined as:

a j = k∈BZ dk e -ik•R A j a k (4.14)
where k∈BZ ≡ BZ d 2 k S BZ with S BZ = 8π 2 3 √ 3a 2 . In order to use the momentum space T-matrix formalism we need to write δH imp in momentum space:

δH imp = U dr j δ(r -R A j ) k∈BZ k ∈BZ dkdk e i(k-k )•R A j a † k a k = U j k∈BZ k ∈BZ dkdk e i(k-k )•R A j a † k a k = U j k∈BZ dk x dk y k ∈BZ dk x dk y e i(k y -k y )( j √ 3a) a † k a k = U 2π √ 3a n k∈BZ dk x dk y k ∈BZ dk x dk y δ(k y -k y + n 2π √ 3a )a † k a k (4.15)
Since both k and k are in the first BZ the only possibilities for n are 0, 1 and -1. It appears that implementing this constraint is quite subtle, however things get much simpler if we deform the first BZ and instead we consider a rectangle with -2π 3 √ 3a < k y < 4π 3 √ 3a and -2π 3a < k x < 2π 3a . This is allowed since a k+Q µν = a k , where Q µν are all the reciprocal basis vectors (Bena 2009; in the tight-binding basis considered here the same relation is valid also for the B atoms (Bena and Montambaux, 2009)). Under this construction it is clear that the only possible solution for n is n = 0, and thus we have

δH imp = U 4π 3 √ 3a -2π 3 √ 3a dk y 2π/( √ 3a) 2π 3a -2π 3a dk x L k 2π 3a -2π 3a dk x L k a † k x ,k y a k x ,k y (4.16)
where L BZ,x = 4π/(3a).

4.B | Recovery of the Kane-Mele zigzag states with the line spectral function

The line spectral function introduced in Sec. 3.3.1 can also be used to visualize the transformation of the impurity states into edge states. The unperturbed Green's function G 0 and the T-matrix T are the same as the ones used to compute the average correction to the spectral function (see Eq. ( 4.3) and (4.4)). From there we compute the perturbed Green's function:

G(k 1 , k 2 , E) = G 0 (k 1 , E)δ(k 1 -k 2 ) + G 0 (k 1 , E)T (k 1 , k 2 , E)G 0 (k 2 , E).
(4.17)

The line Green's function at a given x is defined as:

G l (x, k y , E) = dk 1x dk 2x L 2 BZ,x G(k 1x , k y , k 2x , k y , E)e ik 1x x e -ik 2x x , (4.18)
and the corresponding line spectral function is:

A(x, k y , E) = - 1 π Im tr G l (x, k y , E) . (4.19)
For a system on a honeycomb lattice we define the tight-binding basis such that, in the exponentials of Eq. (4.18), x corresponds to the unit cell position (and not the site position, see App. 4.C for more information). The T-matrix considered here corresponds to a line impurity at x = 0, the effective edges (see Fig. 4.5) are thus situated on the two nearest unit cells i.e. at x = ±3a/2.

In Fig. 4.B.1 we plot the line spectral function at x = 3a/2 for U = 1 and U = 100, and at x = 3a and x = 12a for U = 100. At x = 3a/2 we find the same dispersion for the impurity states as in the results obtained via the average correction to the spectral function (Fig. 4.6). In this case the quantity we compute is always positive and is similar to a local density of states at given energy and position. As we can expect, we thus have at x = 3a/2 a high contribution from the impurity-induced states and a weak contribution from the bulk states. When we then move away from the impurity, the intensity of the impurity states decreases while that of the bulk states increases, until only the bulk contribution remains as we can see in the plot for x = 12a: the impurity states are indeed localized near the impurity.

4.C | Subtlety on the definition of the tight-binding basis for a honeycomb lattice

The honeycomb lattice has a unit cell of two sites, there are therefore two possible ways of defining the tight-binding basis (Bena and Montambaux, 2009). In one basis the positions of the two atoms in the unit cell are used as the origin for the Bloch wavefunctions. The Hamiltonian eigenstates are then written as linear combinations of two Bloch functions corresponding respectively to the A and B atoms, but with a different phase factor attached to each Bloch function:

Ψ k I (r) = c A I (k)Ψ A,k I (r) + c B I (k)Ψ B,k I (r) = 1 √ N j e ik•R A j c A I (k)φ(r -R A j ) + e ik•R B j c B I (k)φ(r -R B j ) , (4.20)
where N is the number of elementary cells, the functions φ(r) are the wavefunctions of the p z orbitals of the carbon atoms and R A/B j are the positions of the A and B atoms of the unit cell j. In the other basis, only one point per unit cell is used as the origin for the Bloch wavefunctions. It thus also consists of two p z orbital wavefunctions centered on the two sites of the unit cell, but with the same phase factor (determined by the position of the "origin" of the unit cell):

Ψ k II (r) = c A II (k)Ψ A,k II (r) + c B II (k)Ψ B,k II (r) = 1 √ N j e ik•R j c A II (k)φ(r -R A j ) + c B II (k)φ(r -R B j ) , (4.21)
where R j = na 1 + ma 2 specifies the position of the unit cell j = (n, m), with a 1,2 the lattice vectors. Depending on the basis considered, certain quantities are expressed differently. For example the component corresponding to nearest neighbor hopping in the Hamiltonian matrix (4.2) is written in the first basis as:

h I NN = t       2e -i 1 2 k x a cos       √ 3 2 k y a       + e ik x a       (4.22)
and in the second basis as (with the unit cell defined as the one drawn in the lowest half of Fig. 4.5a):

h II NN = t 1 + e i √ 3k y a + e i √ 3
2 k y a e i 3 2 k x a .

(4.23)

Of particular relevance to us is the expression for a delta-function impurity potential. For an impurity located on sublattice A, we have:

V A = U A a † j a j = U A k,k ∈BZ e i(k-k )•R A imp j a † I (k)a I (k ) = U A k,k ∈BZ e i(k-k )•R imp j a † II (k)a II (k ), (4.24)
and for an impurity on sublattice B:

V B = U B b † j b j = U B k,k ∈BZ e i(k-k )•R B imp j b † I (k)b I (k ) = U B k,k ∈BZ e i(k-k )•R imp j b † II (k)b II (k ), (4.25)
where the second quantization operators depend on the basis considered:

a I (k) = j e ik•R A j a j , b I (k) = j e ik•R B j b j , a II (k) = j e ik•R j a j , b II (k) = j e ik•R j b j .
Eq. (4.24) and (4.25) clearly show that if we must consider an impurity on both the A and B sublattices i.e. V = V A + V B , as is the case if we want to create two equivalent effective edges (see Fig. 4.5), the second basis is more convenient as there is no phase difference to consider: the exponential factor is the same for both sublattice terms. Moreover since we consider an impurity at x = 0, the exponential becomes 1 and we thus have no phase factor to consider.

For this reason, we always consider the second basis when applying the impurity method to systems on a honeycomb lattice, or on any lattice with more than one atom per unit cell, for which the same reasoning can be applied.

| From the BCS theory to the Andreev Bound States

In 1957, Bardeen, Cooper and Schrieffer published a paper presenting their theory of superconductivity (Bardeen et al., 1957), now known as the Bardeen-Cooper-Schrieffer (BCS) theory. According to this theory, the transition to the superconducting (SC) state is marked by the dominance of some attractive potential, often attributed to electron-phonon interactions, over the repulsive Coulomb interaction between electrons near the Fermi surface. It then becomes favorable for those electrons to combine into what we call Cooper pairs or BCS pairs: bound pairs of electrons with opposite spin and momentum.

The fermionic excitation of a superconductor is no longer a charge eigenstate, but a coherent superposition of an electron and a hole, named Bogoliubov quasi-particle (Bogoljubov, 1958). Here we take interest in a specific phenomenon originating from this particle-hole coherence: Andreev scattering, or Andreev reflection (Andreev, 1964). This novel scattering process occurs at interfaces between a superconductor and a normal conductor, which we refer to as normalsuperconductor (NS) junctions.

We show in Fig. 5.1.1 the difference between the usual reflection which can occur for example at a normal-insulator (NI) interface, and an Andreev reflection at an NS interface. In a normal scattering process, conduction electrons and holes are separately conserved. Therefore, as shown in the left figure, an electron with momentum p, velocity v p p and energy below the energy gap of the insulator is reflected by the insulating gap into an outgoing electron with velocity v p = v p -2 n( n • v p ), with n the normal to the interface. In a similar way at an NS interface, for energies below the SC gap, the electron is forbidden to propagate into the superconducting region. However, upon reflection, the outgoing excitation is a hole: the electron is converted into a hole. Furthermore, the reflected hole is not specularly reflected, i.e. it does not follow the change in velocity given previously, but retro-reflected: v p = -v p . Similarly, a hole will be converted into an electron upon Andreev scattering. Because of this, Andreev scattering is also called branch conversion scattering, where branch refers to the the electron and hole branches of excitation in a normal metal. As one can see, contrary to the case of normal scattering, charge is transported across the NS interface: a Cooper pair of charge -2e is produced in the superconducting region. Note that, since a Cooper pair binds two electrons of opposite spin, the retro-reflected hole has a spin opposite to the incoming electron. Let us now consider a superconductor-normal-superconductor (SNS) junction. Andreev scattering can take place on both sides of the normal section, and quantum interference between those scatterings can occur. This leads to the formation of Andreev bound states (ABS) in the normal region (Andreev, 1966). These states are a superposition of electron and hole excitations, with energies below the SC gap.

In what follows, we study the formation of ABS with the impurity technique. We found that analytical results can be obtained for ABS energies much smaller than the SC gap. This corresponds to the limit of long SNS junctions, thus the work presented here does not apply to the case of short SNS junctions. The impurity we consider is a little different from the ones presented in the previous chapter, as it must be able to model Andreev reflection. As we will explain, this can only be done by also considering imperfect SNS junctions, i.e. where both normal and Andreev scattering can happen at the interfaces. Interestingly, we found that this method is the first to give an analytical description of ABS under those two considerations (long and imperfect junctions), analytical tools being usually used to describe the two limits of perfect junctions (no normal scattering) and of junctions with very large normal scattering.

| Studying ABS with the impurity method

The "Andreev" impurity

To model an NS interface, we consider what we call an "Andreev" impurity. The goal being to model a full SNS junction, the system we will consider is similar to the ones introduced before (see Fig. 5.2.1): we will start with an infinite 1D normal system, and add two Andreev impurities to model the scattering at NS interfaces. This means that we do not model what happens inside the SC regions, but only what occurs inside the middle normal region. This is fine since we are only interested in what happens in that normal region, where the ABS will form. The Andreev impurity must scatter an electron into a hole with opposite spin (and vice versa). Our starting point is the Blonder-Tinkham-Klapwijk (BTK) theory (Blonder et al., 1982), which allows to compute the transmission and reflection coefficients at an NS interface (for both Andreev and normal scattering). The presence of normal scattering is introduced through a barrier strength Z: if Z = 0, the junction is perfect (only Andreev scattering). Here we focus on energies smaller than the SC gap, and an arbitrary barrier strength (neither zero nor too high). For these conditions, the BTK theory indicates that the normal reflection and Andreev reflection coefficients r and r A , resulting from injecting an electron from the normal side of the junction (see Fig. 5.2.2 for the schematics of the junction), are given by:

|r A | 2 = ∆ 2 E 2 + (∆ 2 -E 2 )(1 + 2Z 2 ) 2 ,
(5.1)

|r| 2 = 1 -|r A | 2 , (5.2)
where ∆ is the SC gap. Since we consider subgap energies, the transmission coefficients are equal to zero.

Figure 5.2.2 Schematics of an NS junction where we consider E to be inferior to the SC gap, i.e., the transmission coefficients t, t A = 0.

Here we focus only on energies much smaller than the gap, thus we have:

|r A | 2 = 1 (1 + 2Z 2 ) 2 ,
(5.3)

|r| 2 = 1 -|r A | 2 .
(5.4)

For Z = 0 we have a perfect Andreev junction, i.e., r A = 1, while for Z 1 we have a bad junction with |r| 2 |r A | 2 . In order to verify the equivalence between an Andreev impurity and an NS junction we will calculate the normal reflection and transmission coefficients, as well as the Andreev reflection and transmission for the Andreev impurity, and check in which regime they correspond to the traditional BTK values mentioned above. For a delta-function Andreev impurity (see Fig. 5.2.3) we can write down the Schrödinger equation:

- 2 2m d 2 ψ dx 2 (x) + V imp (x)ψ(x) = Eψ(x),
(5.5)

where ψ(x) = (ψ e (x), ψ h (x)) T is a two-component wave function with the upper and lower components standing for the electron and hole wave functions, respectively. We assume a quadratic dispersion with m denoting the quasiparticle mass. The impurity potential is a 2 × 2 matrix in the electron-hole space, and can be written as:

V imp (x) = Vδ(x) ≡ V N V A V * A -V N δ(x).
(5.6)

In the most general case a normal reflection component V N 0 should be present. However, as we will show, V A alone is already generating both an Andreev and a normal scattering component when taking the equivalence with the BTK formalism, therefore we will set V N = 0.

For x 0 the solution of Eq. (5.5) is just a linear combination of two-component vectors multiplied by right-moving and left-moving plane waves, e ikx and e -ikx , where the wavevector k ≡ √ 2mE/ . The most general solution in the presence of the impurity is then:

ψ L (x) = A r e ikx + A l e -ikx , x < 0, (5.7) ψ R (x) = B r e ikx + B l e -ikx , x > 0.
(5.8)

We focus on the simple case of an injected electron incoming on the barrier from the left, thus A r = (1, 0) T , and no injected electron or hole from the right, i.e., B l = (0, 0) T . Hence A l = (r, r A ) T , where r is the regular reflection coefficient, and r A is the Andreev reflection coefficient, and B r = (t, t A ) T with t being regular transmission and t A Andreev transmission.

By assuming that the wave function ψ(x) = ψ L (x)Θ(-x) + ψ R (x)Θ(x) is continuous at x = 0, we obtain 1 + r = t, r A = t A . We then write down the continuity equation for the derivative of the wave function. The delta-function potential gives rise to a discontinuity at the origin, proportional to the value of the impurity potential, in occurence here:

- 2 2m
[ψ R (0) -ψ L (0)] + Vψ(0) = 0, (5.9) which yields:

-

2 2m ik(-A r + A l + B r -B l ) + V(A r + A l ) = 0.
(5.10)

We note that the solution to this problem is energy-dependent, but since we are interested in energies very close to the Fermi level, we can take k to be constant k = k F , and we denote α ≡

2 k 2 F m .
The above equations yield:

r = - |V A | 2 α 2 + |V A | 2 ,
(5.11)

r A = -i αV * A α 2 + |V A | 2 .
(5.12)

In Fig. 5.2.4 we plot |r| 2 + |r A | 2 and |r A |/|r| as a function of V A , while setting α = 1. This shows that for large enough Andreev potentials the impurity models the NS junction asymptotically well, i.e., |r| 2 + |r A | 2 ≈ 1. On the other hand, we also see that with increasing V A the ratio between the Andreev and the regular reflection is decreasing, e.g., for V A = 3.5, for which |r| 2 +|r A | 2 ≈ 0.92, we barely have |r A |/|r| ≈ 0.3. This corresponds to a value of Z ≈ 1 for the NS junction (see Eq. (5.3)).

We can conclude that, while our impurity cannot model perfect Andreev reflection it can still recreate the Andreev bound states of an imperfect junction.

Recreating the ABS

Double impurity formalism

To model the SNS junction we must consider two Andreev impurities, which requires a slightly different T-matrix formalism from the one for single impurities. This time we will translate the formalism in real space, instead of writing everythin in momentum space like we did for the previous systems, as it makes the equations easier to obtain. Furthermore we are mostly interested in the evolution of the ABS energies as a function of the applied chemical potential and of the phase difference between the two SCs, which does not require us to be in momentum space.

For a single impurity in a 1D system, we had in momentum space (see Sec. 3.2):

G(k 1 , k 2 , E) = G 0 (k 1 , E)δ(k 1 -k 2 ) + G 0 (k 1 , E)T (k 1 , k 2 , E)G 0 (k 2 , E),
(5.13) with:

T (E) = I -V dk V BZ G 0 (k, E) -1
V.

(5.14)

The real space equivalents can be written as:

G(x, x , E) = G 0 (x -x , E) + G 0 (x, E)T (E)G 0 (-x , E) (5.15)
and (5.16) where G 0 (x = 0, E) is equivalent to integrating the Green's function G 0 (k, E) over all momenta (cf. Eq. (5.14)). Now we consider two delta-function impurities with amplitudes V i , localized at x = X i , i ∈ {1, 2}. The full impurity potential can be written as:

T (E) = [I -V • G 0 (x = 0, E)] -1 • V,
V imp (x) = V 1 δ (x -X 1 ) + V 2 δ (x -X 2 ) .
(5.17)

The perturbed Green's function is now expressed as: (5.18) where we now have four T-matrix terms: T 11 , T 12 , T 21 and T 22 , which can be found from the following equation (Choi, 2004(Choi, , 2005)):

G(x, x ) = G 0 (x-x )+ i j G 0 (x-X i )T i j G 0 (X j -x ),
T i j = V i δ i j +V i G 0 (X i -X 1 )T 1 j +V i G 0 (X i -X 2 )T 2 j .
(5.19)

For brevity we omit energy-dependence in all functions. Solving the system of four equations in Eq. ( 5.19), we get:

T 11 = I-V 1 G 0 (0)-V 1 G 0 (X 1 -X 2 )T (0) 2 G 0 (X 2 -X 1 ) -1 V 1 , T 12 = T (0) 1 G 0 (X 1 -X 2 )T 22 , T 21 = T (0) 2 G 0 (X 2 -X 1 )T 11 , T 22 = I-V 2 G 0 (0)-V 2 G 0 (X 2 -X 1 )T (0) 1 G 0 (X 1 -X 2 ) -1 V 2 ,
where we defined

T (0) i ≡ [I -V i G 0 (0)] -1 V i .
Using the Green's function of Eq. ( 5.18) we can compute the correction to the electronic local density of states due to the impurities:

∆ρ(x) = - 1 π Im         i j G 0 (x-X i )T i j G 0 (X j -x ) x =x         11 , (5.20)
where we only take the (1, 1) component of the perturbed part of the Green's function in order to see the electrons' behavior.

In what follows, we will evaluate this quantity as a function of different parameters in order to see the evolution of the ABS energies.

Results

The normal part of the junction is described using a simple lattice Hamiltonian H k = ξ k τ z , where τ z is the Pauli matrix acting in the particle-hole subspace, ξ k ≡ µ -2t cos ka, where a is the lattice constant, t denotes the hopping parameter and µ is the chemical potential. The unperturbed retarded Green's function in real space is computed as the Fourier transform of its momentumspace representation, G 0 (E, k) = E + iδ -H k -1 , and is given by:

G 0 (E, x) = G 11 (E, x) 0 0 G 22 (E, x) , (5.21)
where

G 11 (E, x) = -(E t -µ t + iδ) + sgn(E t -µ t ) (E t -µ t + iδ) 2 -4 |x| a at • sgn(E t -µ t ) • 2 |x| a • (E t -µ t + iδ) 2 -4 (5.22) G 22 (E, x) = + (E t + µ t + iδ) -sgn(E t + µ t ) (E t + µ t + iδ) 2 -4 |x| a at • sgn(E t + µ t ) • 2 |x| a • (E t + µ t + iδ) 2 -4
(5.23) if E 2 + µ 2 0, where we expressed the energy and the chemical potential in terms of the hopping amplitude, i.e., E t ≡ E/t and µ t = µ/t, and

G 11 (E, x) = - i 2at
e +i π|x| 2a , (5.24) (5.25) for E 2 t + µ 2 t = 0. The expressions in Eqs. (5.24) and (5.25) are obtained directly from the integral defining the Fourier transform of the momentum-space Green's function. In other words, we first set in the integrand E t = µ t = 0 and then we perform the integral over momenta to obtain the real-space formulae. Note that since the expressions above are obtained within the lattice model, the results are valid only for x = na, where n ∈ Z.

G 22 (E, x) = - i 2at e -i π|x| 2a ,
The results for the Green's function and the local density of states (LDOS) are evaluated using this full tight-binding model. However, in order to test the validity of our approximations, we also establish a correspondence between the continuum model (used in Sec. 5.2.1 to make the connection between an NS junction and an "Andreev impurity") and the lattice model. We thus expand ξ k in a quadratic form: .26) This allows us to extract:

ξ k = µ -2t cos ka ≈ ta 2 k 2 + µ -2t = 2 k 2 2m - 2 k 2 F 2m . ( 5
t = 2 2ma 2 , α ≡ 2 k 2 F m = 2(2t -µ). (5.27)
In what follows we set by default the hopping parameter t = 1, therefore making µ t and E t equivalent to µ and E respectively. We also use values of the chemical potential around µ = 1.5 so that α ≈ 1 (as set in Sec. 5.2.1). The broadening is set to δ = 0.001 and the lattice constant to a = 1.

-100 0 100 200 300 We now introduce the two Andreev impurities, described by the potential of Eq. (5.6), with V N = 0. In order to take into account the phase difference between the superconductors we replace V A → V A e ±iϕ for one of the impurities, choosing the signs differently for the 12 and 21 components, to preserve the Hermitian character of the Hamiltonian. Using the notations from Eq. ( 5.17), we can thus write:

V 1 = 0 V A V A 0 and V 2 = 0 V A e -iϕ V A e +iϕ 0 .
(5.28)

As mentioned in Sec. 5.2.1 we focus on some intermediate value of V A ≈ 3.5 corresponding to a ratio of Andreev and regular reflection of 0.3, and we study the formation of Andreev bound states using the formalism described in Sec. 5.2.2.1, the T -matrix for two delta-function impurities. This yields an exact formula for the perturbed Green's function of the system (see Eq. (5.18)), which allows to extract the perturbation to the local density of states (see Eq. (5.20)).

We first plot ∆ρ as a function of energy and position in Fig. 5.2.5, for fixed values of the chemical potential µ and the phase difference ϕ. We can see that bound states form at certain energies in the normal region, i.e. for x ∈ [0,200].

In order to compare the behavior of these bound states with that of previously studied ABS (Bena, 2012;Pillet et al., 2010), which we show examples in Fig. 5.2.6, we focus on a given position and we plot in Fig. 5.2.7 the dependence of the LDOS as a function of energy and chemical potential. In the plots we choose to average the LDOS over a few sites (i.e., 6 sites) to avoid any fluctuation effects due to the short wavelength oscillations visible in Fig. 5.2.5. As expected, we recover the oscillatory bell-shaped dependence of the ABS energy on the chemical potential. Furthermore, we plot in Fig. 5.2.8 examples that show the ABS dependence on the phase difference ϕ between the two SCs, and in Fig. 5.2.9 the dependence of the LDOS as a function of energy and ϕ, for a fixed chemical potential. 1 We recover a membrane-like form, with a certain amplitude. While not very large, this amplitude is still significant, marking the presence of nonzero Andreev reflection and of the corresponding phase coherence. This observed dependence of the ABS energies on the phase and chemical potential is connected to the fact that we are studying long SNS junctions in the imperfect long junction limit, and that we focus on the limit of small energies with respect to the gap, E

∆. In what follows we write down the solutions for the two extreme limits characterizing this junction, that have both been studied previously in the literature (Bena, 2012;Hurd and Wendin, 1994;Furusaki and Tsukada, 1991;Kulik, 1969).

First, for the long perfect SNS junctions (no normal backscattering at the leads), the corresponding results are described by the Andreev approximation. Following the results of Hurd and Wendin 1994;Kulik 1969, the energies of the ABS are given by: 2 arccos (E/∆) -(q eq h )L ± ϕ = 2nπ, (5.29)

where

q e = k F √ 1 + E/E F and q h = k F √ 1 -E/E F , with E F ≡ 2 k 2 F /2m = α/2
, are the momenta for the electron and hole of respective energy E and -E, and n = 0, ±1, ±2, ....

Here we consider E ∆, thus arccos (E/∆) ≈ π/2. Also E E F , thus q eq h ≈ k F E/E F and we end up with:

E ± n = 2 k F 2mL π (1 -2n) ± ϕ .
(5.30)

Using the equivalence between the tight-binding and continuum models (Eq. (5.27)), we get:

E ±A n = a L αt 2 π (1 -2n) ± ϕ , n ∈ Z.
(5.31)

In this limit the energies of the ABS are thus linear with the phase and independent of chemical potential.

To the other extreme, in the limit of zero Andreev scattering and infinite normal scattering at the leads, our setup is equivalent to the particle-in-a-box problem. For quadratically-dispersed quasiparticles, the momentum is quantized and takes the values:

k n = nπ L , (5.32)
with n = 1, 2, .... The energies are then: .33) which in the tight-binding model translates to:

E n = 2 π 2 n 2 2mL 2 - 2 k 2 F 2m , ( 5
E B n = ta 2 π 2 n 2 L 2 -2t + µ, n ∈ N * .
(5.34)

In this limit, the energy is thus independent of phase and linear with the chemical potential. In the left panel we fix ϕ = 0 and in the right panel µ = 1.5. The black dashed lines correspond to the energy spectrum in the particle-in-the-box limit (see Eq. ( 5.34)).

The plots in Fig. 5.2.7 and 5.2.9 represent a situation in between these two limits: both Andreev scattering and normal scattering coexist at the interfaces between the N and S regions. To compare with the limit cases, we plot the LDOS for V A = 0.5 in Fig. 5.2.10 and for V A = 10 in Fig. 5.2.11.

According to Fig. 5.2.4, V A = 0.5 corresponds to a non-physical situation, as |r| 2 +|r A | 2 1, but we get a high value for |r A |/|r|, so we could still use it to compare with the Andreev approximation spectrum. And indeed we see that the red lines in Fig. 5.2.10, which represent the energy spectrum in the Andreev approximation, fit quite well the LDOS results. In the left plot, the LDOS shows energy levels which are not quite constant as a function of µ. This is due to the approximation we made of α ≈ 1: α is actually function of µ (see Eq. 5.27) and so the energy levels are not exactly constant as a function of µ. However we can see that on this small range of chemical potential values, the levels still correspond pretty well to the Andreev limit energies. In the right plot we also recover the linear dispersion as a function of ϕ. The correspondence looks exact at zero energy, and a difference is seen when moving away from E = 0. This is not so surprising, as we used some low energy approximations in order to both compute the LDOS and make a comparison with the continuum model.

Going back to Fig. 5.2.4, we see that V A = 10 corresponds to a small |r A |/|r| ratio, getting close to the particle-in-the-box limit. In both plots of Fig. 5.2.11 we again see a pretty good correspondence between our results and the limit case. In the right plot, just like in the Andreev approximation case, we find that the coincidence is best at zero energy. In the left plot, the same thing is observed: the central line, which corresponds best to the limit spectrum, corresponds to the zero energy line of the right plot as we have set ϕ = 0 in the left plot and µ = 1.5 in the right plot.

We can see that by varying the parameter V A , we can effectively model the ABS behaviour in cases ranging from the Andreev approximation limit, where only Andreev scattering is considered, to the particle-in-the-box limit where only normal scattering exists. It is also important to note that these results are obtained very quickly (of the order of a minute for each figure presented here), due to the existence of an analytical solution for the Green's function (Eq. ( 5.22) -(5.25)).

We now change systems again to examine a different way of using the impurity scheme. By adding another layer to the technique, we can also model a phenomenon called quasiparticle interference. The impurity scattering mixes states of different momenta and same energy. Let us consider the scattering of a state at momentum k i to one at momentum k f of the same energy. The quantum interference between the initial and final states will create a standing wave pattern (see Fig. 6.1.1) with a characteristic momentum q = k f -k i (Zhang et al., 2009). By drawing a constant energy contour, for example the Fermi surface, we can determine the available k i and k f at a given energy. The difference q between these two momenta would then be represented by arrows between the two corresponding momentum points (see right panel of Fig. 6.1.2 for an example). The QPI pattern for the same energy will show features for momenta equal to the different q's (see left panel of Fig. 6.1.2).

Since the QPI pattern depends on the states' momenta, it is a way to obtain momentum space information from STM experiments, which can only give real space maps unlike angle-resolved photoemission spectroscopy (ARPES) experiments. QPI patterns can also be used to study scattering processes, which can be valuable in transport experiments (Chang et al., 2016), or to obtain knowledge on the system itself: the states' spin dependence (Huang et al., 2020), information about states inaccessible to ARPES (McElroy et al., 2003), the effective mass of excitations (Aynajian et al., 2012), the symmetries protecting topological states (Zhang et al., 2009)...

Common theoretical methods to compute QPI patterns

Joint Density of States

The most traditional method to obtain the QPI features theoretically is to approximate them by computing the joint density of states (JDOS) for the surface states observed: if we take a constant energy contour, for example the Fermi surface in Fig. 6.1.2, we can pick out the regions of high DOS and assume that all the states within these regions scatter with one another, in proportion to the associated DOS value (Roushan et al., 2009;Wang and Lee, 2003). Mathematically this consists in an autocorrelation of the surface spectral function A(k):

JDOS (q) = A(k)A(k + q)d k. (6.1)
For complicated systems, the spectral function is often computed using density functional theory (Kourtis et al., 2016;Gyenis et al., 2016). Alternatively, one can take the intensity obtained directly from an ARPES experiment instead of the spectral function (Roushan et al., 2009).

Spin-Dependent Scattering Probability

As we just saw, the JDOS does not discriminate between the different surface states and the computed QPI only depends on the spectral function intensity at the different momenta. This is a problem when studying spin-polarized states, as states with opposite spin orientations do not scatter with one another. A common solution is to compute the spin-dependent scattering probability (SSP) instead of the JDOS (Roushan et al., 2009). The SSP is similar to the JDOS, but also includes a spin-dependent scattering matrix element T (k, q)1 :

S S P(q) = A(k)T (k, q)A(k + q)d k. (6.2)
A common form for the matrix element is Roushan et al., 2009;Gyenis et al., 2016), where S(k) is the spin of the state at momentum k. In this form, when the states at k and k + q have aligned spins, T (k, q) = 1 and the term in the integral is the same as for the JDOS; when the spins are nonaligned, T (k, q) decreases, reflecting the weaker scattering between the two states, and eventually goes to zero for opposite spins, completely suppresing scattering.

T (k, q) = | S(k)|S(k + q) | 2 (

Using our method to compute the Fourier Transform of the Local Density of States

The problem, as underlined by Derry et al. 2015, is that a faithful description of interfering scattering processes typically requires relative phase information; and both the JDOS and the SSP lack information about overlap matrix elements between states in the impurity-coupled system. In consequence, the QPI may be small even when the JDOS/SSP is large, in particular it has been shown that these methods fail to capture basic scattering physics in the case of graphene (Simon et al., 2011).

In order to calculate the QPI features accurately, one would need to compute the Fourier transform of the LDOS (FTLDOS). This can be done using the Green's functions and T-matrix formalism. In this case, one needs to first obtain the surface Green's function of the system, then apply the T-matrix formalism to take into account the effect of impurities. As we have shown, our impurity technique allows us to compute this surface Green's function in an intuitive way. Other techniques exist, such as recursion methods, however they are harder to implement numerically (Velev and Butler, 2004).

To obtain the QPI with our technique, we must follow two steps: first compute the surface Green's function for the surface we want to study; then add one or multiple impurities, using the T-matrix formalism, and compute the resulting FTLDOS on the surface. More accurately, since we are interested in the LDOS variations caused by the impurities, we will only compute the correction to the FTLDOS i.e. the part which is modified when adding the impurities.

The first step, computing the surface Green's function, is explained in Sec. 3.3.1: if we consider a 3D system, we compute the bulk Green's function, then add an infinite plane impurity and compute the resulting surface Green's function for a surface parallel to the impurity plane. If the system we consider is 2D, then the surface Green's function is simply the Green's function when the system is taken to be an infinite surface.

In the second step, we add impurities to the surface. Here we will detail the math for a single localized impurity; we have also studied the case of two localized impurities, and give the details in App. 6.A. To simplify, we'll suppose that the BZ extends from -π to π, so that the integrals over k = (k x , k y ) are all normalized by (2π) 2 . We start with the surface Green's function computed in the first step, which we will call G 0 s , so that the T-matrix for a single impurity, described by a delta potential Vδ(r), is given by:

T (E) = I -V d k (2π) 2 G 0 s (E, k) -1 V. (6.3)
To obtain the correction δρ(k) of the FTLDOS, we will first express the correction δρ(r) of the LDOS in real space. Conveniently, the T-matrix in Eq. ( 6.3) only depends on energy, so that we can choose to compute it from the Green's function in momentum space G 0 s (k) even when using the formalism in real space. The perturbed surface Green's function G s in real space can be expressed using the T-matrix:

G s (E, r, r ) = G 0 s (E, r -r ) + G 0 s (E, r)T (E)G 0 s (E, -r ), (6.4)
which is the real space equivalent of Eq. (3.22). The correction to the LDOS due to the impurity can be written in terms of the Green's function:

δρ(E, r) ≡ - 1 π Im tr G s (E, r, r) -G 0 s (E, 0) = - 1 π Im tr G 0 s (E, r)T (E)G 0 s (E, -r) . (6.5)
By applying a Fourier transform (FT) to Eq. ( 6.5), we can calculate the correction to the FTLDOS δρ(E, k) = dr δρ(E, r)e -ik•r . If we rewrite the imaginary part using the identity Im z = 1 2i (zz * ), we get:

δρ(E, k) = - 1 2πi tr dr G 0 s (E, r)T (E)G 0 s (E, -r) -G 0 * s (E, r)T * (E)G 0 * s (E, -r) e -ik•r . (6.6)
Next, we express the unperturbed Green's functions in the integrand as

G 0 s (E, r) = dq G 0 s (E, q)e iq•r : δρ(E, k) = - 1 2πi tr dq (2π) 2 d q (2π) 2 G 0 s (E, q)T (E)G 0 s (E, q) -G 0 * s (E, q)T * (E)G 0 * s (E, q) dr e iq•r e -i q•r e -ik•r = - 1 2πi tr dq (2π) 2 d q (2π) 2 G 0 s (E, q)T (E)G 0 s (E, q) -G 0 * s (E, q)T * (E)G 0 * s (E, q) × (2π) 2 δ(q -q -k) = - 1 2πi tr dq (2π) 2 G 0 s (E, q)T (E)G 0 s (E, q -k) -G 0 * s (E, q -k)T * (E)G 0 * s (E, q) .
If we define:

g(E, q, k) ≡ G 0 s (E, q)T (E)G 0 s (E, q -k) -G 0 * s (E, q -k)T * (E)G 0 * s (E, q), (6.7)
then we obtain:

δρ(E, k) = i 2π dq (2π) 2 tr g(E, q, k) .

(6.8)

This quantity is complex. In all of our results, we show the absolute value of δρ, as is common in QPI analysis.

| Application to Weyl semimetals

In Sec. 4.2.1, we introduced Weyl semimetals (WSM) and their peculiar surface states called Fermi arcs. A large amount of work, both theoretical (to cite a few: Mitchell and Fritz 2016;Gyenis et al. 2016;Lambert et al. 2016;Chang et al. 2016;Kourtis et al. 2016;McCormick et al. 2017;Lau et al. 2017;Xu et al. 2018;Zheng and Hasan 2018 ) and experimental (Deng et al., 2016;Zheng et al., 2016a;Inoue et al., 2016;Batabyal et al., 2016;Zhang et al., 2017;Deng et al., 2017;Lin et al., 2017;Yuan et al., 2018;Li et al., 2018), has been devoted to the impurity scattering and the corresponding QPI features, allowing to extract information on the Weyl nodes, the arcs or the topology of Weyl semi-metals in general.

Here we use the impurity technique to compute the QPI for two different Weyl models. These two models, from Kourtis et al. 2016 andLau et al. 2017, were already introduced in Sec. 4.2.2: they describe Weyl systems with four Weyl nodes and two Fermi arcs. The results we obtain differ from those of Kourtis et al. 2016, where SSP calculations were used, and we argue that it is proof that the impurity technique can capture certain Fermi arc features better than the JDOS and SSP methods.

We give again, in Fig. 6.2.1, the surface spectral function obtained from the Kourtis model with the second set of parameters. If we add an impurity to the surface, the Fermi arcs we see here will scatter and produce a QPI pattern. The surface spectral function was obtained from the surface Green's function, through the formula

A(E, k x , k z ) = -1 π Im{tr[G s (E, k x , k z )]}.
We now use this surface Green's function G s again, in order to obtain the FTLDOS δρ (see Eq. (6.8)). We will set the impurity potential strength V = I. 6.2)), which can be obtained with good precision using a tight-binding calculation, the impurity method computes the full surface Green's function, which is more complete but requires a lot of points when evaluating the different integrals in order to be precise.

We can see two main differences between the two plots: the central feature has a different shape and the wing-like features around k x = ±π/2 observed in the SSP results are absent from our FTLDOS results. The central feature corresponds to intra-arc scattering, i.e. scattering between states belonging to the same Fermi arc. This can be seen by evaluating the extent of a Fermi arc: two states of to the same arc will have a maximum momentum difference of about π/4 in k x and about π/2 in k z , thus the intra-arc scattering pattern should extend to k x = ±π/4 and k z = ±π/2. The fact that the intra-arc feature we obtain extends beyond k x = ±π/4 is most likely due to the broadening issue mentioned above.

In contrast, the wing-like features correspond to inter-arc scattering, i.e. scattering between states belonging to different arcs. The wings' position matches the arcs' relative position: the distance between the Fermi arcs is about π/2 in k x , which creates inter-arc scattering around k x = ±π/2. So, contrary to SSP calculation results, we do not observe inter-arc scattering when obtaining the QPI pattern with the impurity method.

In Sec. 4.2.2, we studied the spin structure of the Fermi arcs, and saw that they had opposite spin (see Fig. 4.13). This means that scattering between the two arcs would require spin flips, which cannot occur when considering a scalar potential as we did here and as was done in Kourtis et al. 2016 (it was shown for instance that in superconductors, scattering between states with opposite spins does not yield any QPI features (Kaladzhyan et al., 2016)). It is mentioned in Kourtis et al. 2016 that the SSP calculations show a diminished intensity of the inter-arc features compared to the JDOS calculation results, which also indicates that taking spin into account tends to restrain scattering between the two arcs; however, their results still show non-zero inter-arc scattering. We interpret this as the SSP method yielding inaccurate results, contrary to the impurity method which takes into account all-order scattering processes through the T-matrix formalism. The tight-binding model introduced in Kourtis et al. 2016 aimed to describe, in a simple way, Weyl materials such as MoTe 2 . In Fig. 6.2.3, we give an experimentally obtained QPI pattern for MoTe 2 (Zheng et al., 2016b). As we can see, inter-arc scattering "wings" are actually observed, similar to those found in Kourtis et al. 2016. One might then question our criticism of the SSP method, since the SSP results do reproduce experimentally observed features. However, this experiment, and others which measure the QPI on MoTe 2 (Deng et al., 2016;Li et al., 2018), do not measure the spin-polarized LDOS (they use non-spin-polarized STM), which means that they cannot measure scattering between states of opposite spins. Indeed, both the presence of magnetic impurities (to allow spin-flips) and the use of spin-polarized STM are required in order to see such scattering features (Kaladzhyan et al., 2016). We conclude that the Fermi arcs in MoTe 2 should not have opposite spins, and thus the tight-binding model of Kourtis et al. 2016 does not model the spin structure of this material correctly. Consequently, the SSP results shown in Fig. 6.2.2 cannot be compared to experimental results, and we maintain that the SSP method itself yields inaccurate results with respect to the model it was applied to.

We can perform a similar analysis for the other Weyl model studied in Sec. 4.2.2, the model from Lau et al. 2017. In This time, as the distance between the two arcs is rather small, it is hard to differentiate intrafrom inter-arc scattering. Even though the pattern is different from what can be obtained with JDOS, we assume that both intra-and inter-arc scattering are present in our results. This is in agreement with the study of the spin structure done in Sec. 4.2.2 (Fig. 4.16), where we saw that the arcs do not have fully opposite spin, thus allowing inter-arc scattering to occur.

We have also studied the effects of changing the impurity strength, taking a magnetic impurity instead of a scalar one, and considering double-impurity scattering on the QPI pattern of the Kour-tis model, which we show in App. 6.B. In all cases, the shape of the intra-arc scattering pattern changes, but no inter-arc scattering pattern emerges.

In this chapter, we showed that we can use the impurity technique to obtain the QPI patterns of different Weyl semimetal models. For the Kourtis model, we found that no inter-arc interference pattern arise, which differs from original predictions based on SSP calculations. Since the impurity technique uses the T-matrix formalism to compute the perturbed Green's function, it can take into account relative phase information, whereas the JDOS and SSP techniques throw this information away when computing the QPI pattern from the spectral function only. For this reason, we believe that the impurity technique is better suited to compute QPI patterns than the commonly used JDOS and SSP calculations.

G s (E, r i , r j ) = G 0 s (E, r i -r j ) + α,β=0,1 G 0 s (E, r i -r α )T αβ (E)G 0 s (E, r β -r j ), (6.9)
where α and β run over the two impurities, and T αβ (E) is the real-space T-matrix:

T αβ (E) = I -VG 0 s (E, r α -r β ) -1 V = I -V d k (2π) 2 G 0 s (E, k)e ik•(r α -r β ) -1
V.

(6.10)

As we can see, this time the T-matrix depends on r, which makes the formula for δρ a little more complicated. However, we will show that it can keep the simple form of Eq. (6.8).

In real space, the correction to the FTLDOS is now:

δρ(E, r) = - 1 π Im tr         α,β=0,1 G 0 s (E, r -r α )T αβ (E)G 0 s (E, r β -r)         .
(6.11)

In momentum space, using again the identity Im z = 1 2i (zz * ), we obtain:

δρ(E, k) = - 1 2πi tr dr e -ik•r α,β=0,1 G 0 s (E, r -r α )T αβ (E)G 0 s (E, r β -r) -G 0 * s (E, r -r α )T * αβ (E)G 0 * s (E, r β -r) = - 1 2πi tr dq (2π) 2 d q (2π) 2
dr e iq•r e -i q•r e -ik•r α,β=0,1 G 0 s (E, q)T αβ (E)G 0 s (E, q)e -iq•r α e i q•r β -G 0 * s (E, q)T * αβ (E)G 0 * s (E, q)e i q•r α e -iq•r β .

If we now define the momentum space T-matrix as:

T (E, q, q) = α,β=0,1 e -iq•r α e i q•r β T αβ (E), (6.12)

then we can rewrite δρ as:

δρ(E, k) = - 1 2πi tr dq (2π) 2 d q (2π) 2 G 0 s (E, q)T (E, q, q)G 0 s (E, q) -G 0 * s (E, q)T * (E, q, q)G 0 * s (E, q) dr e iq•r e -i q•r e -ik•r = - 1 2πi tr dq (2π) 2 G 0 s (E, q)T (E, q, q -k)G 0 s (E, q -k) -G 0 * s (E, q -k)T * (E, q -k, q)G 0 * s (E, q) ,
where we recover Eq. (6.8), with the expression for g now modified to:

g(E, q, k) = G 0 s (E, q)T (E, q, qk)G 0 s (E, qk) -G 0 * s (E, qk)T * (E, qk, q)G 0 * s (E, q). (6.13)

6.B | Study of different scattering impurities for the Kourtis model

In this Appendix we show the effects of various scattering impurities on the QPI pattern for the Kourtis model.

6.B.1 Single scalar or magnetic impurity

We first study the QPI patterns for a single scalar impurity V = UI with different values for the impurity strength U. The results are shown in Fig. 6.B.1. We find that changing U modifies the shape of the intra-arc scattering features. However, no inter-arc scattering appears.

-ππ We now consider a magnetic impurity: V = U diag{1, -1, 1, -1}, where "diag" stands for a diagonal matrix. The QPI patterns obtained are shown in Fig. 6.B.2. For a strong enough impurity i.e. U = 10, the intra-arc pattern resembles the one obtained with a scalar impurity. For lower values of U, the pattern's shape is quite different, but we still do not find any inter-arc scattering.

6.B.2 Double impurity

QPI patterns for two impurities localized at r 0 = (x = 0, z = 0) and r 1 = (x = 0, z = 1) or r 1 = (x = 1, z = 0) are shown in Fig. 6.B.3. We find slightly different shapes for the intra-arc scattering features, nevertheless we observe no inter-arc scattering. Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of trigonal warping 71

| Different stackings in multi-layer graphene and graphite

Graphite, or multi-layer graphene, is a system composed of 2D graphene layers stacked on top of each other. The term "stacking" refers to the arrangement of the successive layers with regard to each other. Here we focus on the two most stable forms of stacking in graphite: Bernal stacking, also referred to as ABA stacking, and rhombohedral stacking, also known as ABC stacking. An illustration of these two stackings is given in Fig. 7.1.1. In both cases, as depicted in panel B, the first two layers are stacked on top of each other with a relative translation of a 0 x, where a 0 is the distance between two sites of a given layer, and x designates the horizontal axis. In the xy plane (the plane parallel to the layers) we refer to the position of the first layer (the top, yellow layer) as A, and that of the second layer (in blue) as B. The discrimination between the two stackings appears with the third layer: in Bernal stacking the third layer will be at position A, it is thus exactly below the first layer (panel C), while in rhombohedral stacking the third layer is at a different position from the other two layers, which we label C (panel D). From the fourth layer on, the stacking then repeats: for Bernal stacking we would get the sequence ABAB and for rhombohedral stacking ABCA, etc.

Regarding hoppings, we will consider three different terms, shown in the lower left picture of Fig. 7.1.1: an intra-layer hopping γ 0 (which is denoted as t in the graphene study of Sec. 2.1), and two inter-layer hoppings γ 1 and γ 3 . The γ 1 term connects sites which belong to two nearest layers and have the same (x, y) coordinates; due to the relative translation of these layers, it thus always connects a site from an A sublattice to one from a B sublattice. The γ 3 term also connects sites from opposite sublattices, but which have different (x, y) coordinates. As we will show , this last term is responsible for an effect called trigonal warping (Charlier et al., 1991). We will fix γ 0 = 3.3, γ 1 = 0.42, and will consider either γ 3 = 0 or γ 3 = 0.3. These are the values, in eV, previously used by the experimental teams when comparing their results to a theoretical model (Joucken et al., 2020).

Hamiltonians

For an easier understanding, let us first give the tight-binding Hamiltonian for bilayer graphene (BLG). We consider the momentum space basis (c k,l -,A , c k,l -,B , c k,l + ,A , c k,l + ,B ), where c k,l,α is the electron operator at momentum k = (k x , k y ), layer l and sublattice α, with l -and l + corresponding to the lower and upper layer of the BLG respectively. In this basis, we have:

H BLG (k) =                  0 h 0 (k) 0 h 3 (k) h * 0 (k) 0 γ 1 0 0 γ 1 0 h 0 (k) h * 3 (k) 0 h * 0 (k) 0                  , (7.1)
where we defined

h 0 (k) = γ 0       1 + 2e -i 3 2 a 0 k x cos       √ 3 2 a 0 k y             , (7.2) h 3 (k) = γ 3 ε(k), (7.3)
and

ε(k) ≡ 2e -i 3 2 a 0 k x cos       √ 3 2 a 0 k y       + e -i3a 0 k x . (7.4)
Here a 0 = 1.42 is the distance (in Å) between two nearest-neighbor atoms in the honeycomb lattice. As can be seen in the matrix of Eq. (7.1), both the γ 1 and the γ 3 terms connect sites from opposite sublattices. However, the γ 1 term connects a B site from the lower layer to an A site from the upper layer, while the γ 3 term connects an A site from the lower layer to a B site from the upper layer. Moreover, γ 1 is added to the matrix as is, while the γ 3 term (labeled h 3 ) includes ε, function of momentum, which translates the difference between the (x, y) coordinates of the two connected sites.

Let us now add a third layer to obtain trilayer graphene (TLG). The new basis we consider is (c k,l -,A , c k,l -,B , c k,l 0 ,A , c k,l 0 ,B , c k,l + ,A , c k,l + ,B ), where l -, l 0 and l + correspond to the lower, middle and upper layer. For rhombohedral stacking, the Hamiltonian in matrix form is then:

0 h 0 (k) 0 h 3 (k) 0 0 h * 0 (k) 0 γ 1 0 0 0 0 γ 1 0 h 0 (k) 0 h 3 (k) h * 3 (k) 0 h * 0 (k) 0 γ 1 0 0 0 0 γ 1 0 h 0 (k) 0 0       0 h 0 (k) 0 h 3 (k)χ(k) h * 0 (k) 0 γ 1 χ(k) 0 0 γ 1 χ * (k) 0 h 0 (k) h * 3 (k)χ * (k) 0 h * 0 (k) 0                  , (7.11)
with χ(k) = 1 + e -2id 0 k z . (7.12)

Here d 0 = 3.35 is the distance (in Å) between two neighboring layers. As we can see, these Hamiltonians are very close to the TLG Hamiltonian for rhombohedral stacking in Eq. (7.5) and the BLG Hamiltonian in Eq. (7.1), except for the terms with ξ(k) and χ(k). These two functions, which are the only ones with a dependence on k z , embody the periodicity of the stacking.

In what follows, we compare the behavior of graphite with that of few-layer graphene in both stackings. In coherence with the definition of the unit cells, we choose BLG for the Bernal case and TLG for the rhombohedral case, so that we compare the graphite behavior with that of a system composed of a single unit cell. The Hamiltonians used are thus those of Eq. (7.1, 7.5, 7.9, 7.11).

| Surface spectral functions

We can obtain the surface spectral function for graphite the same way we obtained it in the previous chapters. Here we create a surface at fixed z, thus the impurity potential is V = U I δ(z), where we will set U = 100 000. The small δ, which comes from the analytical continuation iω n → E + iδ used to obtain the retarded Green's function from the Matsubara Green's function, is set to 0.005.

The surface spectral function must be computed at a z which is one unit cell away from the impurity plane. This originates from the definition of the phase factors in the Hamiltonian (see App. 4.C): here the same phase factor is applied to the sites of a same unit cell, so that the effective surfaces created by the impurity are one unit cell away from it. Since the unit cells defined for the two stackings are different, the spectral function must be computed at different z's: at z = 2d 0 for Bernal stacking and z = 3d 0 for rhombohedral stacking1 .

To compare with the graphite results, we also compute the surface spectral function for BLG and TLG with rhombohedral stacking. In these two systems, the number of layers is fixed so that we do not need the impurity technique to create a surface artificially. Instead we directly compute the spectral function by first obtaining the Green's function from the Hamiltonian, then calculating the quantity -1 π Im tr G where the trace is taken on the sites of the upper layer only. We compute the surface spectral function at low energy, for two different values E = 0 and E = 0.05. As is characteristic of graphene and graphite, the states at these energies only exist around the Dirac point, therefore we restrict our calculation to momentum values around one of these points: we choose K = 2π 3a 0 , 2π 3 √ 3a 0

. We thus plot in Fig. 7.2.1 and 7.2.2 the surface spectral functions of the different graphene and graphite systems, with and without trigonal warping and at the two different energy values.

For Bernal stacking, in Fig. 7.2.1, we first observe that in the absence of the warping term γ 3 , the spectral function seems to be rotationally symmetric in the (k x , k y ) plane, both for graphite and bilayer graphene and for both energy values. Then, once the trigonal warping is turned on, the shape clearly acquires a three-fold symmetry. The resulting pattern seems to be much less sharp for graphite than for BLG: the equal energy lines are much wider and less well defined than for BLG. This is however also true, even if less obvious in particular at E = 0, when γ 3 = 0. This suggests that quasiparticles have shorter lifetimes on the surface of graphite than in BLG, and thus a decreased coherence length.

Looking now at the rhombohedral-stacked systems, in Fig. 7.2.2, we see some differences. Let us first compare the results to those of Fig. 7.2.1. The patterns for BLG and TLG look pretty much the same, except for γ 3 = 0.3, E = 0, however the momentum space range is different: the TLG states extend further than the BLG ones (the plot ranges for k x and k y are different). This is also the case when comparing ABA graphite with ABC graphite. This wider extension is due to the flat band, characteristic of rhombohedral stacking (see App. 7.A). This flat band also explains the saturation of the blobs for ABC graphite at E = 0 (panels b and d), in contrast with the finer features of ABA graphite at the same energy. As for the trigonal warping effect, we can see that the effect is somehow less pronounced for ABC graphite than it was for ABA graphite: the blobs transform into something mildly triangular, while the three-fold symmetry was much more obvious for ABA graphite. We will now calculate the QPI patterns for the systems described above, by considering a localized scalar impurity on the surface. For ABA and ABC graphite we start from the surface Green's functions and we compute the FTLDOS, as described in Sec. 6.1.3. In this case though, compared to the Weyl semimetal study, we make one small adjusment to the FTLDOS formula:

δρ(k, E) = i 2π dq S BZ tr sur f g(q, k, E) . (7.13)
Here, instead of taking the usual trace, we take the trace tr sur f only over the matrix components corresponding to the two top-layer atoms. This must be done in order to compare with experimental results, since the contribution to the density of states that is usually measured experimentally is the one from the top layer. We will thus compute G s at negative z and trace over the last two components, corresponding to the upper layer of the unit cell. In the case of Bernal graphite, this layer will then be a "B" layer, while it will be a "C" layer for rhombohedral graphite.

For bilayer graphene and trilayer graphene the same formalism applies, but with the surface Green's function G s replaced by the unperturbed Green's function G 0 (k, E) = [(E +iδ)I-H(k)] -1 , with H(k) the unperturbed BLG and TLG Hamiltonians.

Unlike the WSM study, we will also not compute the full QPI pattern. Here we focus on the intra-scattering pattern created by the states around the Dirac points, in accordance with the spectral function results we showed above.

In graphene and graphite only two inequivalent Dirac points exist; to compute the QPI, we can thus take into consideration the surface Green's function around two complementary Dirac points only. This is supported by the plots in Fig. 7.3.1, where we show the surface spectral function for TLG in the whole BZ: only two different patterns can be seen, each corresponding to one type of Dirac point.

However we can also notice, by comparing Fig. 7.3.1 (b) and (c), that the surface states around the two different Dirac points are related to each other: one pattern can be obtained from the other by applying a rotation of π2 around the corresponding K point, and a translation from one K point to the other. Let us now consider the integral in Eq. (7.13). Since surface states only exist in small areas around the Dirac points, we can separate this integral into two: one over an area close to K, the other over an area close to K . For each of these sub-integrals, the scattering between states at two different momenta q and q is computed "both ways": the scattering momentum k, over which the QPI is computed, takes both the values qq and qq. The translation from K to K , or from K to K, leaves k unchanged: [q -(K -K )] -[q -(K -K )] = qq . A rotation by π changes qq into qq, however since both terms are already included in the integral, this does not change the value of δρ(k). We thus see that even though the surface states around K and K are different, they will yield the same intra-scattering QPI pattern,3 which means we can compute the QPI from the surface states around a single K point. In order to see a more extended pattern, here we focus on E = 0.05, and we set γ 3 = 0.3 to consider trigonal warping. For the point-like impurity situated on the surface, different configurations are possible depending on where the impurity is situated in the unit cell: A or B sublattice, upper, middle or lower layer. We can thus consider a buried impurity if we localize it on the middle or lower layer. The associated impurity potential V p will only have one non-zero component, corresponding to the site on which the impurity is considered. For example, for an impurity on the sublattice A of the upper layer in Bernal graphite, the associated component is the 3rd component, and we have V p = diag(0, 0, U, 0). In what follows, we will fix U = 10. In Fig. 7.3.2 we plot the QPI pattern for the different systems, considering an impurity on the lowest layer and on the B sublattice (2nd matrix component). As the main conclusions remain the same, we collected the results for impurities localized on the other sites in App. 7.B.

One thing we notice is that, while in BLG and TLG the intensity is highest along an outer hexagon of the pattern, in graphite we instead see highest intensity in the center, particularly for rhombohedral stacking where, for this impurity, it is concentrated in a small circle-like area. Looking back at the spectral function for ABA graphite in Fig. 7.2.1 (h), where we saw that the states were more concentrated on the sides of a triangle, we infer that scattering mostly occurs inside those sides rather than between different sides. As for ABC graphite, which spectral function is shown in Fig. 7.2.2 (h), it seems that even though the states concentrate around a uniform triangular shape they don't necessarily scatter together.

The patterns for BLG and TLG also appear sharper and more detailed than the ones for graphite. This fits with the spectral function results, where we found more blurry patterns for graphite. Like the spectral function motif, the QPI is also wider for graphite, since the allowed states are more widely distributed.

Real space QPI

We can also study the QPI pattern in real space, in the vicinity of the impurity (located at x = 0, y = 0). To do this, we simply have to Fourier transform the momentum QPI δρ(k) from Eq. (7.13) In Fig. 7.3.3 we plot the results obtained for two types of impurities: one on the sublattice A of the lower layer, the other on the sublattice B of the upper layer. We refer to App. 7.C for the results obtained with the other types of impurities.

First we observe that the interference pattern is more extended for BLG and TLG than for graphite. This is consistent with the reduction of the quasiparticle lifetime observed in graphite, i.e. the blurring of the surface spectral function compared to BLG and TLG. This can be understood as graphite having one more dimension in which the quasiparticles can scatter (BLG and TLG can be considered close to 2D systems). This effect of shorter quasiparticle lifetime due to an added dimension has been observed previously in the case of Yu-Shiba-Rusinov states in two-versus three-dimensional superconductors (Yazdani et al., 1997;Ji et al., 2008;Ménard et al., 2015).

The other feature that we observe, is that we obtain a triangular pattern for buried impurities (situated on the lower layer of the unit cell), while we see a circular one for impurities directly on the surface (upper layer), even though trigonal warping is present in equal measure in both cases. We also observe a faster decay for the circular pattern. Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of trigonal warping 83

Dependence on the number of layers and energy

We now study the dependence of the real-space QPI patterns on the number of layers and energy.

For the dependence on the number of layers, we focus on a Bernal-stacked system with an impurity on the sublattice A of the second last layer. The parameters are the same as before: we take E = 0.05, γ 3 = 0.3 and U = 10. We plot in Fig. 7.3.4 the QPI for two, three, four and eight layers. We note that the pattern for three layers is oriented differently from the others: this orientation is linked to whether the number of layers is odd or even, since the unit cell of Bernal-stacked graphite consists of two layers. Overall, the results coincide with what we observed when comparing the QPI patterns of few-layer graphene versus graphite: when increasing the number of layers, the oscillations decay faster when moving away from the impurity. Already for eight layers we recover the same pattern as the one for graphite in Fig. 7.3.3 (b).

-300 -150 0 150 300

x/a 0 For the dependence on energy, we take BLG with again an impurity on the sublattice A of the lower layer and the same parameters as before, except that we vary E. We plot in Fig. 7.3.5 the QPI for E = 0, 0.02, 0.05 and 0.25. We observe that when increasing the energy, the oscillations have higher frequency: the space between two positive or two negative peaks of the QPI decreases.

Correspondence with experiments

We will now compare the results obtained through this study with that of experimental results obtained by the previously mentioned experimental teams. Information on how those results were obtained can be found in [START_REF] Joucken | Direct visualization of native defects in graphite and their effect on the electronic properties of bernal-stacked bilayer graphene[END_REF]Kaladzhyan et al. 2021a. All the maps presented here were acquired at low-tip sample bias (around 5mV), so that they are essentially proportional to the LDOS at the Fermi level. In Fig. 7.3.6 we show the real and momentum space QPI patterns obtained for BLG at different positive gate voltages. The change in gate voltage induces a variation of the Fermi level, and we indicate the distance between the Fermi level and the charge neutrality point (Dirac point) for each gate voltage value. We can see that, even if not so clear for V G = +20V, for V G = +40V and +60V the momentum space pattern shows stronger intensity on an external ring, which looks like a hexagon. This is in agreement with what we found when computing the FTLDOS for BLG. As for the real space images, we clearly see that when increasing the gate voltage, the frequency of the oscillations increases. This also corresponds to the energy-dependent study of the real-space QPI where we found that the space between two positive or two negative peaks decreases when moving away from E = 0, i.e. the charge neutrality point.

In Fig. 7.3.7 we show spatial maps of graphite and of a 10-layer graphene film along with the computed FFT image. For the momentum-space image, we see that the intensity is this time great-Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of trigonal warping 85 est in the center of the pattern. This confirms the tendency observed in the simulations where the FTLDOS showed overall stronger values in the center for graphite. As for the real-space experimental maps, they clearly show a triangular symmetry. From the results we obtained previously, this thus suggests that the impurities at the origin of those triangular patterns are embedded in the material and not directly at its surface. We also observe that, contrary to the results on BLG, the oscillations do not extend far from the impurities, consistent with the fast decay seen in the computed real-space QPI for graphite.

Conclusion

Using the impurity technique, we studied the surface spectral function as well as the impurityinduced quasiparticle interference patterns in real and momentum space of bernal-and rhombohedralstacked graphene and graphite. We found that the spectral function for graphite surfaces shows a decrease in coherence and quasiparticle lifetime compared to that of few-layer graphene. Then we observed that, independent of the type of stacking, surface impurity-induced oscillations decay much faster in graphite than in BLG and TLG. Furthermore, the shape of the oscillations seem to depend on the location of the impurity with respect to the layers: a buried impurity creates triangular shaped oscillations, while an impurity situated on the top layer will create a rather circular one. By comparing these results with experimental measurements, we found that they are in overall good agreement with the experimental data.

Appendix 7.A | Recovering the predicted surface states of Bernal-and rhombohedral-stacked graphite

In this appendix, we recover the surface states of graphite for both stackings using the impurity technique. In particular, we obtain the surface flat bands predicted to appear in multi-layer ABC graphene [START_REF] Zhang | [END_REF]Yelgel, 2016), and are able to determine their extent in a fully-infinite system. We will consider no trigonal warping here, i.e. we set γ 3 = 0.

For comparison, we start by calculating the energy spectrum of finite slabs of graphene using tight-binding calculations. To do this, we use the Hamiltonians in Eq. (7.7, 7.8). In Fig. 7.A.1 we plot the spectrum around the K point for a 100-layer ABA graphite slab. The spectrum shows multiple bands with parabolic dispersion. We then calculate the spectrum of ABC graphite, which we plot in Fig. 7.A.2 for different numbers of layers. As we can see, the surface states in this case form a flat band at zero energy, which extends further away from the K point when the number of layers is increased. It seems that when the number of layers tends to infinity, the extension of the flat band reaches a limit value in k space.

We now compute the surface spectral function for semi-infinite graphite, over the same range of energies and momenta, using the impurity technique. In Fig. 7.A.3 we consider both stackings and plot A s for the corresponding surfaces: z = 2d 0 for Bernal stacking and z = 3d 0 for rhombohedral stacking.

For ABA graphite we recover a parabolic-shaped band, which is consistent with the tightbinding spectrum presented in Fig. 7.A.1. Indeed the energy spectrum should contain both the surface and the bulk bands for the finite graphite slab, and we find only part of that spectrum in Fig. 7.A.3: the exterior parabolas, with higher intensity close to zero energy corresponding to a higher density of surface states.

For ABC graphite we observe once more a flat band, consistent with the results of the tightbinding calculations of the finite-size graphite slabs. The extension of the flat band observed here seems close to the one observed in the energy spectrum for 100 layers, confirming that it is converging to a finite value when the width of the slab is going to infinity.

Conclusion

In this thesis, we took interest in a new method to study edge and surface phenomena. First we showed that it can correctly recover known edge states, such as those characteristic of topological insulators and Weyl semimetals; then we saw that it can model more complicated boundary effects, recreating for instance Andreev bound states in certain SNS junctions; and finally we found that the technique actually makes a good foundation for the study of quasiparticle interference, both in momentum and real space, showing examples on Weyl semimetal models and on graphene and graphite systems.

Through the study of all these systems, we observed several advantages to the technique. For one, it could be considered semi-analytical, in the way that it is exact up to some integrals' calculation. The coding of the technique is also rather simple: starting from a small Hamiltonian matrix, one needs only follow the successive equations to compute the Green's functions and T-matrices, and the resulting spectral function or LDOS for example. It works in all physical dimensions (1D, 2D and 3D), and with both topological and non-topological systems. Its greatest strength most probably lies in its construction: it can model edges while bypassing the need to consider finite systems.

While this manuscript was being written, other applications of the impurity technique were pursued. Fully analytical results have been obtained for the QPI of bilayer graphene at lowenergy (Kaladzhyan et al., 2021b), granting in particular more insight into the effects of trigonal warping. The method has also been used to study Shiba chains and their topology [START_REF] Sedlmayr | New tools to determine the topological character of shiba chains[END_REF], where it can ignore the size-effect problems known to occur in tight-binding calculations. Among possible future applications, the quasiparticle interference seems like a promising choice, since it can take into account important contributions which are suppressed when using the usual approximations of JDOS and SSP. Abstract: We present here a new method to help study boundary phenomena, which we call "impurity method". The technique consists in simulating a boundary, such as an edge or a surface, by adding impurities along a line or a plane, with a high impurity potential (compared to the typical energy values of the system). Through a series of example systems, we then explore the possibilities which this method offers, along with its advantages and its limits. In particular, we show that it can be applied for any physical system dimension (1D, 2D and 3D), that it works for both topological and non-topological systems, and that it can bring several advantages compared to certain conventional techniques, such as tight-binding calculation or the approximations of JDOS and SSP used in the study of quasi-particle interference.

Figure 1

 1 Figure 1 Illustration de deux phases différentes : une phase ordonnée (cristal) où la position des atomes est fixe, et une phase désordonnée (liquide) dans laquelle les atomes n'ont qu'une certaine probabilité d'exister en un point donné de l'espace. Crédit à Trung Van Phan.

Figure 2

 2 Figure 2 Représentation de la technique des impuretés sur un système 2D : à un système infini (en gris) on ajoute une ligne d'impuretés (en rouge), coupant effectivement le système en deux. Chaque sous-système possède ainsi un bord effectif, où des états de bord peuvent se former (en bleu).

Figure 3

 3 Figure 3 Représentation de l'impureté que nous avons utilisée pour créer des bords dans les systèmes 1D, 2D et 3D. Dans chaque cas, le système infini est soumis à une impureté, qui le divise effectivement en deux sous-systèmes semi-infinis. Si le système possède des états de bord, ceux-ci sont recréés de chaque côté de l'impureté. Ces états de bord sont représentés en bleu en (a) et (b), ou existent sur la surface représentée en rouge en (c).

Figure 1 . 1

 11 Figure 1.1 Illustration of two different phases: an ordered phase (crystal) where the position of the atoms is fixed, and a disordered phase (liquid) for which the atoms only have a certain probability of existing at a given point in space. Credit to Trung Van Phan.

Figure 2 . 1

 21 Figure 2.1 The honeycomb lattice in graphene. The unit cell, in red, is composed of two sites labeled A and B. The black lines connect nearest-neighbor atoms.

Figure 2 . 2

 22 Figure 2.2 Left: upper energy band of bulk graphene with the BZ drawn on top. K points are indicated by circles and crosses, equivalent points are drawn with the same marker. Right: energy spectrum of bulk graphene along the left plot's black dashed line.

Figure 2 . 3

 23 Figure 2.3 The three possible edges of a graphene ribbon.

  bearded edge states exist for k y a I. In Fig.2.4c, both a zigzag and a bearded edge state are present: the flat band extends to all values of k y .

Figure 2 . 4

 24 Figure 2.4 Energy spectrum of a graphene ribbon depending on the edges considered. In each case, we consider a ribbon width of 50 cells in the open direction.

E

  Figure 2.7 SSH chain and its energy spectrum for t 1 = 0.

Figure 2 . 8

 28 Figure 2.8 Spectrum of an SSH chain of 100 sites for t 2 = 1 and different values of t 1 .

Figure 2 . 9

 29 Figure 2.9 Modulus square of the zero-energy eigenstates of an SSH chain of 100 sites as a function of the position x along the chain, for t 2 = 1 and different values of t 1 .

  4) where θ is the Heaviside function, A and B are two operators and [A, B] ζ = AB -ζBA with ζ = 1 if we consider bosons and ζ = -1 if we consider fermions.

( a )

 a Figure 3.1 Reprensentation of the impurity we used to create edges in 1D, 2D and 3D systems.In each case, the infinite system is subjected to an impurity, which effectively splits it into two semi-infinite subsystems. If the system possesses edge states, these are recreated on each side of the impurity. These edge states are represented in blue in (a) and (b), or exist on the surface represented in red in (c).

  Appendix 4.A Derivation of the T-matrix momentum limits of integration for a honeycomb lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Appendix 4.B Recovery of the Kane-Mele zigzag states with the line spectral function 41 Appendix 4.C Subtlety on the definition of the tight-binding basis for a honeycomb lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Figure 4 . 1

 41 Figure 4.1 The honeycomb lattice with the ν i j convention.

  Figure 4.2 Bulk energy spectrum for the Kane-Mele model along a line crossing the two inequivalent Dirac points (see Fig. 2.2), with t = 1 and t 2 = 0.03.

Figure 4 . 3

 43 Figure 4.3 Energy spectrum for the Kane-Mele model on a 50-cells-wide ribbon with zigzag edges. t = 1 and t 2 = 0.03.

Figure 4 . 4

 44 Figure 4.4 Energy spectrum for the Kane-Mele model on a 50-cells-wide ribbon with different edges. t = 1 and t 2 = 0.03. In the case of armchair edges, we have restricted the energy values to show the low energy behavior more clearly.

  Two effective bearded edges

Figure 4 . 5

 45 Figure4.5 Infinite graphene lattice where we add a line impurity, localized on a single unit cell along x (periodicity is conserved along y). Impurity sites, represented by larger dots, are effectively disconnected from the rest of the lattice, as indicated by the dotted lines. If the vector which goes between the two sites of the u.c. has a component along y (two possibilities exist for the u.c. definition), we then obtain effective zigzag edges on each side of the impurity (blue lines). Conversely, if the u.c. vector is parallel to the x axis, we get two effective bearded edges.

  Figure 4.6 Average correction to the spectral function due to the line impurity for two effective zigzag edges, for U = 1, 10 and 100. For better comparison with Fig. 4.3, we also plot the U = 100 case for a larger range of k y . We have set t = 1 and t 2 = 0.03.

  Figure 4.7 Average correction to the spectral function due to the line impurity for different types of effective edges, for U = 100. We have set t = 1 and t 2 = 0.03.

  Two effective armchair edges

Figure 4 . 8

 48 Figure 4.8 Infinite graphene lattice with a doubled unit cell where we add a line impurity such as to reproduce different edge configurations.

Figure 4 . 9

 49 Figure 4.9 Depiction of a Weyl semimetal, with two Weyl nodes in the bulk and a Fermi arc extending between the projections of the two nodes on the surface Brillouin Zone. (Modified from Lv et al. 2015)

Figure 4 .

 4 Figure 4.10 High-resolution ARPES Fermi surface map of crescent-shaped Fermi arcs on the surface of tantalum arsenide (TaAs), modified from Xu et al. 2015)

-Figure 4

 4 Figure 4.11 Minimum positive energy band, obtained via exact diagonalization, for a system defined by the Hamiltonian (4.8) with periodic boundary conditions in x and z and open boundary conditions in y (length of 100 sites). We take two different sets of parameters which generate a different number of Weyl nodes.

Figure 4 .

 4 Figure 4.13 Spin components on the A and B orbitals at E = 0 for the Hamiltonian (4.8) with the two sets of parameters: set a) in the first line of plots, and set b) in the second line.

-

  Figure 4.14 Minimum positive energy band, obtained via exact diagonalization, for a system defined by the Hamiltonian (4.11) with periodic boundary conditions in x and z and open boundary conditions in y (length of 100 sites). We take two different sets of parameters which generate two Fermi arcs, and one Fermi pocket in the second configuration.

Figure 4 .

 4 Figure 4.15 Surface spectral function at E = 0, y = -1, for the Hamiltonian (4.11) with a surface impurity at y = 0. The plots are obtained with the same set of parameters as in Fig. 4.14. We set U = 100.

1 2

 1 Figure 4.16 Spin components on the A and B orbitals at E = 0 for the Hamiltonian (4.11) with the two sets of parameters: set a) on the left and set b) on the right.

  U = 100, x = 12a

Figure 4 .

 4 Figure 4.B.1 Line spectral function for two effective zigzag edges, for U = 1 or U = 100 and at different positions away from the impurity. We have set t = 1 and t 2 = 0.03.

Figure 5

 5 Figure 5.1.1 Left: usual reflection at a normal-insulator boundary. Right: Andreev reflection at a normal-superconductor boundary (Sauls, 2018).

Figure 5

 5 Figure 5.2.1 Schematics of the equivalence between "Andreev" impurities and NS junctions.

Figure 5

 5 Figure 5.2.3 Schematics of an "Andreev" impurity with the incoming, reflected and transmitted plane waves.

  Figure 5.2.4 In red the ratio |r A |/|r| and in blue |r| 2 + |r A | 2 as a function of V A (in units of α).

  Figure 5.2.5 Correction to the electronic LDOS as a function of position and energy, with Andreev impurities at x = 0 and x = 200. The chemical potential and phase difference are fixed to µ = 1.5 and ϕ = 0, with V A = 3.5.

Figure 5

 5 Figure 5.2.6 Examples of the ABS profile as a function of chemical potential obtained theoretically (upper plot, from Bena 2012) and experimentally (lower plot, from Pillet et al. 2010).

  Figure 5.2.7 Correction to the electronic LDOS, averaged over 6 sites x ∈ [90, 95], as a function of chemical potential and energy. µ varies between 1.46 and 1.54, so that α ≈ 1 does not vary with µ. The phase difference is set to ϕ = 0, with V A = 3.5.

Figure 5

 5 Figure 5.2.8 Examples of the ABS profile as a function of the phase difference ϕ obtained theoretically (left plot, from Bena 2012) and experimentally (right plot, from Pillet et al. 2010).

  Figure 5.2.9 Correction to the electronic LDOS, averaged over 6 sites x ∈ [90, 95], as a function of the phase difference ϕ and energy. We set µ = 1.5, with V A = 3.5.

Figure 5 Figure 5

 55 Figure 5.2.10 Correction to the electronic LDOS, averaged over 6 sites x ∈ [90, 95], as a function of the chemical potential µ (left plot) or the phase difference ϕ (right plot) and energy, for V A = 0.5. In the left panel we fix ϕ = 0 and in the right panel µ = 1.5. The red dashed lines correspond to the energy spectrum in the Andreev limit (see Eq. (5.31)).

Figure 6

 6 Figure 6.1.1 STM image showing standing wave patterns around impurities (Crommie et al., 1993).

Figure 6

 6 Figure 6.1.2 Left: part of a QPI pattern obtained experimentally from a Weyl semi-metal surface at the Fermi level. Right: calculated Fermi surface for the same system. On the Fermi surface, different momenta q = k f -k i , are designated by arrows. The shape of the QPI pattern follows these arrows, showing interference between the different states observed in the Fermi surface image (Zheng et al., 2016a).

  Figure 6.2.1 Surface spectral function at E = 0, y = -1, for the Kourtis model with a surface impurity at y = 0 (same plot as in Fig. 4.12).

Figure 6

 6 Figure 6.2.2 Left: Quasiparticle interference pattern corresponding to the surface spectral function of Fig. 6.2.1, with V = I. Right: SSP obtained in Kourtis et al. 2016 for the same model and parameters.

Figure 6

 6 Figure 6.2.3 Experimental QPI pattern for the Weyl semi-metal MoTe 2 (Zheng et al., 2016b).

Figure 6

 6 Figure 6.2.4 Left: Surface spectral function at E = 0, y = -1, for the Lau model with a surface impurity at y = 0 (same plot as in Fig. 4.15).

Figure 6

 6 Figure 6.2.5 Left: Quasiparticle interference pattern corresponding to the surface spectral function of Fig. 6.2.4, with V = I. Right: JDOS obtained in Lau et al. 2017 for the same model and parameters.

Figure 6

 6 Figure 6.B.1 QPI patterns for the Kourtis model calculated for a scalar impurity with varying strength U. From left to right we set: U = 0.5, U = 1, U = 3, and U = 10.

Figure 6

 6 Figure 6.B.2 QPI patterns for the Kourtis model calculated for a magnetic impurity with varying strength U. From left to right we set: U = 0.5, U = 1, U = 3, and U = 10.

Figure 6

 6 Figure 6.B.3 QPI patterns for the Kourtis model calculated for two impurities, considering different values of the impurity strength U. In the left column we set U = 1, in the central column U = 3 and in the right column U = 10. For the first row we set r 1 = (x = 0, z = 1), and for the second one, r 1 = (x = 1, z = 0).

  Figure 7.1.1 Illustration of the ABA and ABC stackings (modified from Shan et al. 2018). A: single layer of graphene. B: bilayer graphene with Bernal stacking. C: Bernal-stacked graphite. D: rhombohedral-stacked graphite.

Chapter 7 .----

 7 Figure 7.2.1 Surface spectral function of BLG and ABA graphite around the K point (corresponding to k x = 0, k y = 0), with and without trigonal warping.

Chapter 7 .

 7 Figure 7.3.1 Surface spectral function of TLG for E = 0.05 with γ 3 = 0.3 for the entire BZ and around the two inequivalent Dirac points K and K .

Chapter 7 .

 7 Figure 7.3.2 Quasiparticle interference pattern computed for the different graphene and graphite stackings, for an impurity situated on the lowest layer and sublattice B.

Chapter 7 .

 7 Figure 7.3.3 Real space QPI for the different systems around two distinct impurities: one on the sublattice A of the lower layer, the other on the sublattice B of the upper layer.

  Figure 7.3.4 Real space QPI for a Bernal-stacked system with different numbers of layers. The impurity is set on the sublattice A of the second last layer.

Chapter 7 .

 7 Figure 7.3.5 Real space QPI for bilayer graphene at different energies. The impurity is set on the sublattice A of the lower layer. In some of the plots we have saturated the colors, as indicated by the arrow-like ends of the colorbars, so as to see the oscillations better. Note that the we have doubled the extent of the plot in space compared to the previous figures.

Figure 7

 7 Figure 7.3.6 Experimental measurements of the real-space and Fourier-space quasiparticle interference patterns for impurities in BLG for different gate voltages V G . Left panels: dI/dV S maps. Right panels: zoom-ins around the center of the corresponding FFT images. The induced difference between the Fermi level E F and the charge neutrality point E CNP is indicated for each gate voltage value (Joucken et al., 2021).

  Figure 7.A.1 Energy spectrum for Bernal-stacked graphite with 100 layers, along a line with fixed k x and crossing the K point (corresponding to k y = 0 on the horizontal axis).

Chapter 7 .

 7 Figure 7.A.2 Energy spectrum for ABC graphite with 10, 50 and 100 layers, along a line with fixed k x and crossing the K point (corresponding to k y = 0 on the horizontal axis).

Figure 7 .

 7 Figure 7.A.3 Surface spectral function for both stackings of graphite, along a line with fixed k x and crossing the K point (corresponding to k y = 0 on the horizontal axis). For ABA stacking we choose a surface at z = 2d 0 , and for ABC stacking one at z = 3d 0 .

Chapter 7 .

 7 Figure 7.B.1 Quasiparticle interference pattern computed for BLG and Bernal graphite, for an impurity situated on different layers and sublattices.

Chapter 7 .

 7 Figure 7.B.2 Quasiparticle interference pattern computed for TLG and rhombohedral graphite, for an impurity situated on different layers and sublattices.

Titre:

  Etude de phénomènes de bords via l'utilisation d'impuretés Mots clés: Bord, Impuretés, Matrice T, Fonction de Green, Topologie Résumé: Nous présentons ici une nouvelle méthode permettant d'étudier des phénomènes dus aux bords, nommée "méthode des impuretés". Cette technique consiste à simuler un bord ou une surface en plaçant des impuretés le long d'une ligne ou d'un plan, et en faisant en sorte que le potentiel qui leur est associé soit grand (par rapport aux énergies caractéristiques du système). A travers une suite d'exemples de systèmes, nous explorons alors les possibilités qu'offre cette nouvelle méthode, ses avantages et ses limites. Nous montrons en particulier qu'elle peut s'appliquer indépendamment de la dimension du système (1D, 2D, 3D), qu'elle fonctionne pour des systèmes topologiques comme nontopologiques, et qu'elle peut apporter des avantages par rapport à certaines techniques usuelles telles que le calcul tight-binding ou les approximations JDOS et SSP utilisées pour étudier les interférences de quasi-particules. Title: Study of edge phenomena through the use of impurities Keywords: Boundary, Impurities, T-matrix, Green's function, Topology

An introduction to edge states and other edge-related phenom- ena

	Chapter 2 |

Contents 2.1 Simple graphene: a non-topological system with edge states . . . . . . . . . . . 10 2.2 A brief introduction to topology: the Su-Schrieffer-Heeger chain . . . . . . . . 13

Green's functions and the T- matrix formalism

  11 Bulk spectrum of the SSH chain for t 2 = 1 and several values of t 1 across the phase transition.

	Chapter 3 | Contents
	3.1	Green's functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	18
		3.1.1 The retarded Green's function . . . . . . . . . . . . . . . . . . . . . . .	19
		3.1.2 The Matsubara Green's function . . . . . . . . . . . . . . . . . . . . .	20
		3.1.3 The spectral function and DOS . . . . . . . . . . . . . . . . . . . . . .	21
	3.2	The T-matrix formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	21
	3.3	Using the T-matrix formalism to study edge states . . . . . . . . . . . . . . . .	23
		3.3.1 Point, line and surface Green's functions . . . . . . . . . . . . . . . . .	24
		3.3.2 Average correction to the spectral function . . . . . . . . . . . . . . . .	26

Studying edge and surface states with the impurity method

  Chapter 4 |

Contents 4.1 A topological insulator: the Kane-Mele model . . . . . . . . . . . . . . . . . . 27 4.1.1 What is a topological insulator? . . . . . . . . . . . . . . . . . . . . . . 27 4.1.2 The Kane-Mele model . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Weyl semimetals and Fermi arcs . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2.1 A gapless 3D topological system . . . . . . . . . . . . . . . . . . . . . 34 4.2.2 Application of the impurity method . . . . . . . . . . . . . . . . . . . 35

  Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of

	trigonal warping			81
	back into real space, i.e.:
	δρ(r) =	d k S BZ	e ik•r δρ(k)	(7.14)

La rupture de symétrie peut en fait être utilisée pour décrire certaines transitions de phase topologiques, et ce n'est donc pas la raison pour laquelle les phases topologiques ne peuvent pas être expliquées par la théorie de Landau. La différence réelle réside dans le paramètre d'ordre : la théorie de Landau décrit les transitions de phase à l'aide de paramètres d'ordre locaux, dont la valeur moyenne est nulle dans une phase désordonnée et non nulle dans une phase ordonnée ; alors qu'il est impossible d'en définir un pour les transitions de phase topologiques, nous définissons donc à la place des invariants topologiques. Cependant, ces paramètres d'ordre local étant directement liés à la brisure de symétrie, la brisure de symétrie reste un concept très pertinent. Pour des raisons pédagogiques nous avons donc décidé de le mentionner dans le texte principal, à la place du paramètre d'ordre local.

L'existence des états liés de Majorana est à ce jour encore considérée comme théorique, leur preuve expérimentale restant très controversée.

Tous les résultats présentés dans cette thèse ont d'ailleurs été calculés de cette manière.

Certaines comparaisons seront néanmoins faites avec la méthode de diagonalisation exacte, d'une part parce qu'elle est couramment utilisée, d'autre part parce que nous l'utiliserons nous-mêmes pour vérifier que la technique de l'impureté peut retrouver les mêmes résultats, et ainsi confirmer sa validité.

Symmetry-breaking can be used to describe some topological phase transitions, and so this is not the reason why topological phases cannot be explained by the Landau theory. The actual difference lies in the order parameter: the Landau theory describes phase transitions using local order parameters, which average value is zero in a disordered phase and non-zero in an ordered one; while it is impossible to define one for topological phase transitions, instead we define topological invariants. However, these local order parameters being directly connected to the breaking of symmetry, symmetry-breaking is still a very relevant concept. For pedagogical purposes we thus decided to mention it in the main text, instead of the local order parameter.

The existence of Majorana bound states is to this date still considered theoretical, their experimental proof remaining very controversial.

All the results presented in this thesis have actually been computed this way.

Some comparisons will still be done with the exact diagonalization method, for one because it is commonly used, for another because we will use it ourselves to verify that the impurity technique can recover the same results, to confirm its validity.

The integral in Eq. (3.26), as well as the other ones running over momentum space presented in the rest of this chapter, are all approximated numerically by either Riemann sums or by using the trapezoidal rule.

Note that the ϕ range in the left plot of Fig.5.2.8 is shifted by a value of π compared to the plot of Fig.5.2.9.

This T must not be confused with the T-matrix that we have been working with so far, we are just reusing the common notation for this matrix element. Luckily we will not be using this T outside of this section, avoiding any possible confusion.

Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of trigonal warping

This technique also allowed us to recover the predicted surface states of graphite, including the flat-band states in rhombohedral stacking. We show this in App. 7.A

Since the patterns possess a C

symmetry, this is also valid for a rotation of π/3 or 5π/3. 3 This has been verified through computations.

| Introduction

Quasiparticle interference patterns and how they can be used

When impurities are present in a system, the quasiparticle states are scattered by them. This can be seen as a modulation in the LDOS. When moving to momentum space, this modulation generates what we call a quasiparticle interference (QPI) pattern. The QPI pattern is thus a picture in momentum space, generated by the presence of impurities (it is however to note that in some articles, the momentum space picture is referred to as the "Fourier transform of the QPI": they refer to the QPI as the real-space oscillations around the impurities).

Experimentally, the QPI pattern can be obtained using Fourier-transform scanning tunneling spectroscopy (FT-STS). Scanning tunneling spectroscopy (STS) is the measurement of the differential conductance dI/dV of a tunnel junction in a scanning tunneling microscopy (STM) experiment (Crommie et al., 1993). The dI/dV map for a given bias voltage V can be assimilated to a map of the LDOS at a specific energy E F + |e|V. Applying a Fourier transform to this map, for example by using a fast Fourier transform (FFT) algorithm, thus yields the QPI pattern at a given energy.

Appendix

6.A | QPI pattern derivation for double-impurity scattering

Here we give the derivation to compute the FTLDOS when we add two impurities to the surface. The T-matrix formalism can be modified to include scattering from both impurities (Derry et al., 2015;Mitchell et al., 2015). Assuming that the impurities are localized at r 0 and r 1 , the perturbed surface Green's function can be written as: We now use the impurity technique to study one final set of systems: few-layer graphene and graphite. The aim is to study the difference and similarities between graphite and few-layer graphene for two different types of stacking: Bernal (ABA) stacking and rhombohedral (ABC) stacking. In particular we examine the effects of trigonal warping, a symmetry distortion which originates from a particular interlayer hopping term. Part of this study originates from a joint work with experimental teams from the University of California in Santa Cruz, USA and the National Institute for Materials Science in Tsukuba, Japan [START_REF] Joucken | Direct visualization of native defects in graphite and their effect on the electronic properties of bernal-stacked bilayer graphene[END_REF], who investigated the role of defects on graphite samples commonly used to prepare graphene devices. As we will present in more detail in this chapter, we were able to use the impurity technique to provide a theoretical understanding for the scattering effects observed in the experiments.

We will first say a few words on stacking types in multi-layer graphene and graphite, introducing the two stacking which are studied here. Then, similar to the study of the Weyl semimetals, we will evaluate the surface spectral function for the different systems, then compute the QPI pattern, both in momentum and real space. Finally, we will compare our results with the experimental measurements mentioned above.

In this case, the type of sites connected between the two upper layers is opposite to the one between the two lower layers: the γ 1 term connects an A site from the middle layer to a B site from the upper layer, while the γ 3 term connects a B site from the middle layer to an A site from the upper layer. Also the term corresponding to the operators c † k,l 0 ,B c k,l + ,A (fourth line, fifth column of the matrix) is h * 3 (k) instead of h 3 (k) as the difference between the (x, y) coordinates of the upper layers is opposite to the one of the lower layers (as can be seen in the lower left panel of Fig. 7.1.1).

By repeating the layer-related patterns highlighted in the TLG Hamiltonians, we can write equations for graphite Hamiltonians in both stackings, for any number of layers. For rhombohedral stacking, the pattern repeats every layer (the γ 1 term always connects a B site from a layer to an A site from the upper layer, while the γ 3 term always connects an A site from a layer to a B site from the upper layer), thus we get:

In contrast, the pattern in Bernal stacking repeats only every two layers:

where we consider the first layer to be l = 0, thus even.

When increasing the number of layers, the characteristics of graphite converge to that of an infinite system. This includes for example the energy spectrum and the scattering behavior on the surface, which we will examine in this chapter. The Hamiltonian for infinite graphite does not depend on l, but instead on k z ; its matrix form for a given momentum point (k x , k y , k z ) is thus smaller than the one for a finite number of layers (of size 2l × 2l). Since the impurity technique is also more efficient when using small Hamiltonian matrices, we will here model graphite with an infinite number of layers (except for one section near the end of this chapter, where we will briefly study the dependence of the real space QPI on the number of layers).

Chapter 7. Study of few-layer graphene versus graphite for different stackings, and the effect of trigonal warping 74

For the infinite system, we now have k = (k x , k y , k z ). Since we have periodic boundary conditions in the z direction, we must redefine the unit cell to include that direction (up to now, the unit cell for a given layer consisted only of an A and a B site). This unit cell depends on the stacking. For rhombohedral stacking, the unit cell must consist of three layers: as we can see in the panel D of Fig. 7.1.1, the position of the layers in the (x, y) plane repeats every three layers (ABCABC...). For the same reason, the unit cell for Bernal stacking only contains two layers (ABABAB...). The size of the Hamiltonian matrices are thus respectively 6 × 6 and 4 × 4: