Pushing the limit of quantum transport simulations

TL;DR

This paper introduces a new class of algorithms for quantum transport simulations that efficiently access the thermodynamic limit in systems with translational invariance, enabling studies of larger and more complex quantum materials.

Contribution

The authors develop a novel Green's function-based algorithm that analytically integrates over momentum, allowing direct simulation of infinite or semi-infinite quantum systems with reduced computational cost.

Findings

Efficient simulation of infinite graphene electrodes.

Analysis of Friedel oscillations in 2D systems.

Study of topological surface states and their disorder resilience.

Abstract

Simulations of quantum transport in coherent conductors have evolved into mature techniques that are used in fields of physics ranging from electrical engineering to quantum nanoelectronics and material science. The most efficient general-purpose algorithms have a computational cost that scales as $L^{6 \dots 7}$ in 3D, which on the one hand is a substantial improvement over older algorithms, but on the other hand still severely restricts the size of the simulation domain, limiting the usefulness of simulations through strong finite-size effects. Here, we present a novel class of algorithms that, for certain systems, allows to directly access the thermodynamic limit. Our approach, based on the Green's function formalism for discrete models, targets systems which are mostly invariant by translation, i.e. invariant by translation up to a finite number of orbitals and/or quasi-1D electrodes and/or the presence of edges or surfaces. Our approach is based on an automatic calculation of the poles and residues of series expansions of the Green's function in momentum space. This expansion allows to integrate analytically in one momentum variable. We illustrate our algorithms with several applications: devices with graphene electrodes that consist of half an infinite sheet; Friedel oscillation calculations of infinite 2D systems in presence of an impurity; quantum spin Hall physics in presence of an edge; the surface of a Weyl semi-metal in presence of impurities and electrodes connected to the surface. In this last example, we study the conduction through the Fermi arcs of the topological material and its resilience to the presence of disorder. Our approach provides a practical route for simulating 3D bulk systems or surfaces as well as other setups that have so far remained elusive.

Figures

Equations

$${\hat{H}}_{\mathrm{tot}}=\sum _l{\hat{H}}_l^{\infty }+\sum _a{\hat{W}}_a,$$

$$\hat{H}^{\infty}=\sum_{\mu\nu}\sum_{x,y,z=-\infty}^{\infty}{H}_{0}^{\mu\nu}|x,y,z,\mu\rangle\langle x,y,z,\nu|\\ \hskip 56.9055pt+\big{(}V_{x}^{\mu\nu}|x-1,y,z,\mu\rangle\langle x,y,z,\nu|\\ \hskip 56.9055pt+V_{y}^{\mu\nu}|x,y-1,z,\mu\rangle\langle x,y,z,\nu|\\ \hskip 28.45274pt+V_{z}^{\mu\nu}|x,y,z-1,\mu\rangle\langle x,y,z,\nu|+\text{h.c.}\big{)},$$

$${\hat{W}}_0=\sum _{ij}{W}_{ij}|i⟩⟨j|,$$

$${\hat{W}}_1=\sum _{z=-\infty }^{\infty }\sum _{i,j}{W}_{ij}|z,i⟩⟨z,j|$$

$$+W_{ij}^z|z+1,i⟩⟨z,j|+h.c.,$$

$${\hat{W}}_2$$

$$\hat{G}(E)\equiv \frac{1}{E-{\hat{H}}_{\mathrm{tot}}+i\eta },$$

$${\hat{g}}_l(E)=\frac{1}{E-{\hat{H}}_l^{\infty }+i\eta }.$$

$$⟨\text{r’}|{\hat{g}}_l(E)|\text{r}⟩={\int }_{-\pi }^{+\pi }\frac{d{k}_{x}d{k}_{y}d{k}_z}{(2\pi {)}^3}{e}^{i\text{k.(r - r’)}}⟨\text{k}|{\hat{g}}_l(E)|\text{k}⟩.$$

$$\begin{array}{rl}{\hat{H}}_{2D}=& t\sum _{x,y}\left(|x+1,y⟩⟨x,y|+|x,y+1⟩⟨x,y|\right)\\ & +\epsilon |0,0⟩⟨0,0|,\end{array}$$

$$H_0$$

$$V_x$$

$$V_y$$

$$\rho (\text{r},E)=-\frac{1}{\pi }\text{Im}(⟨\text{r}|G(E)|\text{r}⟩).$$

$${\mathcal{H}}_{\mathrm{BHZ}}(\text{k})=\left(\begin{array}{cc}h(\text{k})& 0\\ 0& h^{\ast }(-\text{k})\end{array}\right),$$

$$h(\text{k})=ϵ(\text{k})+\vec{d}(\text{k})\cdot \vec{\sigma },$$

$$ϵ(\text{k})$$

$$\vec{d}(\text{k})$$

$$H_0$$

$$V_x$$

$$V_y$$

$${\hat{H}}_{\mathrm{G}}=\sum _{⟨i,j⟩}|i⟩⟨j|,$$

$$T_{l{l}^{\prime }}=\mathrm{Tr}(\hat{G}{\Gamma }_l{\hat{G}}^{†}{\Gamma }_{l^{\prime }}),$$

$${\Gamma }_l=i({\Sigma }_l-{\Sigma }_l^{†})$$

$${\Sigma }_l=W_l{G}_l{W}_l^{†}.$$

$$\begin{array}{r}{\mathcal{H}}_{\mathrm{WEYL}}(\text{k})={\tau }_z\left[t({\sigma }_x\sin (k_x)+{\sigma }_y\sin (k_y))+t_z{\sigma }_z\sin (k_z)\right]\\ +\mu (\text{k}){\tau }_x{\sigma }_0+\frac{1}{2}{b}_0{\tau }_z{\sigma }_0+\frac{1}{2}\beta {\tau }_0{\sigma }_z,\end{array}$$

$$\begin{gathered} \mu (\text{k})={\mu }_0+t(2-\cos (k_x)-\cos (k_y)) \\ +t_z(1-\cos (k_z))). \end{gathered}$$

$$H_0$$

$$V_x$$

$$V_y$$

$$V_z$$

$$\left(\begin{array}{c}{I}_0\\ I_1\\ I_2\end{array}\right)=\frac{e^2}{h}\left(\begin{array}{ccc}{T}_{02}+T_{01}& -T_{01}& -T_{02}\\ -T_{10}& N_1-R_{11}& -T_{12}\\ -T_{20}& -T_{21}& N_2-R_{22}\end{array}\right)\left(\begin{array}{c}{V}_0\\ V_1\\ V_2\end{array}\right),$$

$$\begin{split}R_{aa}={}&N_{\rm a}+\text{Tr}(\Gamma_{a}\hat{G}\Gamma_{a}\hat{G}^{\dagger})\\ &+i\big{(}\text{Tr}(\Gamma_{a}\hat{G})-\text{Tr}(\Gamma_{a}\hat{G}^{\dagger})\big{)},\end{split}$$

$$\rho (x,k_y,k_{z}E)=-\frac{1}{\pi }\text{Tr}\text{Im}⟨x,k_y,k_z|G(E)|x,k_y,k_z⟩,$$

$${\hat{g}}_l(E)=\frac{1}{E-{\hat{H}}_l^{\infty }+i\eta }.$$

Full text

I Introduction

Quantum nanoelectronics is changing from a domain of fundamental research into one of the main platforms for development of quantum technologies. The fabrication techniques that are used in this field are related to those employed at industrial scale in microelectronics and thus hold the promise of scalable superconducting or semi-conducting quantum bits. Yet, unlike in microelectronics, where a full simulation stack for the devices is available (from the device layout to its functionalities), the modeling of quantum systems is still a work in progress with various problems remaining unsolved.

The simulation of quantum transport is one of the areas that has a long history, almost as old as quantum nanoelectronics itself PhysRevLett.47.882 ; Thouless_1981 ; MacKinnon1985 . Early techniques were mostly based on the recursive Green function algorithm MacKinnon1985 that was later generalized to various geometries, materials and multi-terminal systems Bruno_RGF_2005 ; Baranger_1991 . Other techniques employed wave-functions and the scattering matrix approach PhysRevB.44.8017 . Such simulations are now rather mature and open source simulation packages are becoming available kwant_article . While the computational cost of these calculations is easily affordable in one ($\sim L$) and two ($\sim L^{3\dots 4}$) dimensions, it quickly reaches prohibitive levels in three ($\sim L^{6\dots 7}$) dimensionsagrawal1993cutting , a case of high practical interest ($L$: typical length of the device). There are many situations where one must simulate large three dimensional systems: it is the case in particular for any realistic geometry; or in the presence of different characteristic length scales (such as Fermi wave length and mean level spacing); or in the case of topological materials such as 3D topological insulators or Weyl semi-metalsXu613 where well-separated surfaces are important to avoid surface states mixing.

The need to simulate 3D quantum systems has provoked the development of new methods that scale linearly with the system size. These techniques include the Kernel polynomial expansionkpm_method ; fan2018linear and variants of the Lanczoslanczos1950iteration method like Lanczos recursion methodlow_T_Lanczos ; finite_T_lanczos . See Ref. fan2018linear, for a recent review. The accuracy of these approaches increases with the number of iterations and very often they focus on local properties such as (semi-classical) conductivities or the local density of state (with exceptions such as the noticeable Ref. Mayou_transport, ). The conductance of a coherent quantum conductor is, however, intrinsically a global quantity since it stems from the interference of spatially separated paths and new approaches must be developed to tackle this problem.

Here, we present a set of algorithms for quantum transport that work directly in the thermodynamic limit. Our approach addresses mostly translationally invariant systems (MTIS), i.e. systems that are translationally invariant up to a finite number of modifications. MTIS are infinite in size, and hence do not suffer from the finite size effects that often affect traditional approaches. They allow to describe many setups of practical interest such as bulk systems with impurities; bulk systems with line or surface defects; multilayer quantum wells; surfaces (half of the 3D space is filled with the material) with impurities and/or surface (topological) states and or electrodes attached to the surface;ai etc. Our approach takes advantage of the structure of MTIS: instead of starting from vacuum and adding new parts to the systems (the usual “bottom-up” approach), we first calculate the properties of the true translationaly invariant system and then, in a second step, include the modifications (“top-down” approach). A very appealing aspect of this technique is that it naturally handles systems with very different characteristic scales as the computing cost increases with the number of modifications made with respect to a pristine material but is independent of distances.

The article is organized as follows. In section II, we provide a mathematical formulation of MTIS and explain our strategy for computing their properties. Section III describes applications to a few concrete systems of interest together with some benchmarks. We discuss (i) Friedel oscillations in a Two-dimensional electron gas; (ii) the edge state properties of a 2D quantum spin Hall model; (iii) The conductance of an infinite sheet of graphene with a constriction and (iv) the multi-terminal differential conductance of the surface of a Weyl semi-metal to which has been attached quasi-1D electrodes (and impurities). Readers not interested in the technical aspects can stop at the end of Sec. III. The rest of the article explains our algorithm which can be broken down into a collection of separate subproblems. The main technical challenge lies in performing a momentum integral of matrices that contain Dirac and principal part distributions. This is addressed in section IV using an novel algorithm for calculating poles and residues of Green’s function matrices. The other subproblems involve delicate quadrature rules for performing Fourier transform in presence of integrable singularities (section V), using the Dyson equation (section VI) and calculating the contributions due to bound states (section VII).

II Mostly translationally invariant systems

In this section we introduce the concept of mostly translationally invariant systems (MTIS) and its mathematical formulation. MTIS allow to describe a great variety of quantum systems of high interest that are out of scope of traditional simulation techniques. Examples of MTIS are given in Fig. 1 with additional MTIS shown in the application section III. The most advanced system that is studied in this article is the disordered surface of a Weyl semi-metal (that spans an entire 3D half-space) to which two semi-infinite quasi-one-dimensional electrodes are attached from above, as shown in Fig. 7.

II.1 Problem formulation

The method presented in this paper allows to simulate infinite systems whose Hamiltonians can generically be written as the sum of two terms,

[TABLE]

where the first one contains translationally invariant Hamiltonians, and the second one contains perturbations that break translational invariance along one or several directions and/or connect several infinite systems together. The translationally invariant systems (along one, two or three dimensions) will be called TIS hereafter.

Although not fully general, Eq. (1) covers a very large subgroup of tight-binding systems, including multilayer systems or the surface of a bulk material (3D half-space filled with a material). A few examples of MTIS in two dimensions are given in Fig. 1d-f: Fig. 1d shows an infinite 2D sheet where a finite number of sites have defects (vacancies, missing bounds and/or modified hoppings or on-site energies); Fig. 1e shows a semi-infinite 2D sheet (a system of this kind will be used to form graphene electrodes in the application section); Fig. 1f shows a semi-infinite sheet with a reconstruction of the edge attached to a molecule that in turn is attached to a quasi-one dimensional electrode. A three dimensional example that generalizes Fig. 1f is shown in Fig. 2h. Other examples are presented in the application sections.

The Hamiltonian matrix of TIS takes the generic form (hereafter, we drop the suffix $l$ unless needed explicitly)

[TABLE]

where a ket $|x,y,z,\mu\rangle$ is labeled by the spatial positions $x,y,z$ on the lattice as well as an orbital degree of freedom $\mu$ that accounts for spin, particle-hole, different atoms in the unit cell and/or different orbitals. The Hamiltonian matrix Eq. (2) contains the on-site matrix $H_{0}$ and the nearest-neighbour hopping matrices $V_{x}$, $V_{y}$ and $V_{z}$ along $x,y$ and $z$ respectively. In general, it should also include nine different diagonal hoppings such as $\hat{V}_{xy}=V_{x^{-}y^{+}}^{\mu\nu}|x-1,y+1,z,\mu\rangle\langle x,y,z,\nu|$. We have omitted them for clarity but the algorithms presented in this article do account for these terms. The case of second-nearest-neighbor hoppings can be represented in the above model by merging two unit cells into a larger unit cell with first-neighbor hoppings only, etc. The matrices $H_{0}$, $V_{x}$, $V_{y}$ and $V_{z}$ account for $N_{\rm o}$ orbitals per unit cell. Summation over the orbital degrees of freedom will often be implicit in the following. The structure of Eq. (2) is shown schematically in Fig. 1a,b (1D) and Fig. 1c (2D).

The second type of terms, $\hat{W}_{a}$ in Eq. (1), breaks translational invariance along one, two or all directions. Terms that break translational invariance in all directions take the form

[TABLE]

with $i=(x,y,z,\mu)$. These terms include e.g. impurities or hoppings between two infinite systems. Terms that break translational invariance along all directions but one take the form

[TABLE]

with $i=(x,y,\mu)$. These terms describe e.g. an edge reconstruction as in Fig. 1f. They are also used to cut an infinite 2D sheet into two separate parts to create systems like the one shown in Fig. 1e (by adding the negation of the hopping $V_{x}$ to the Hamiltonian at the bound to be severed). Terms that break translational invariance along all directions but two are defined similarly,

[TABLE]

with $i=(x,\mu)$. The restriction that we impose on the $\hat{W}$ matrices is that only a finite numbers of matrix elements $W_{ij}$, $W_{ij}^{z}$ and $W_{ij}^{y}$ may be non-zero. In practice, up to a few tens of thousand sites can be involved by each of these terms. An important aspect is that these sites need not be placed close to each other spatially. The computational cost for handling a system with e.g. 1000 impurities is independent of the distance between the impurities.

II.2 Principle of the technique

We now turn to the description of the method that we have developped to address the MTIS. The approach takes full advantage of the decomposition shown in Eq. (2) into an infinite translationally invariant and a finite arbitrary part. While the global algorithm may appear somewhat complex, it decomposes into well-defined subproblems. These subproblems have been (at least partially) resolved in the past except for one that we call the “residue problem”. The chief result of this article is an algorithm that solves the residue problem thereby unlocking the development of the present algorithmic suite. Below, we describe the principle of our method. The mathematical details are given in later sections.

The main mathematical object studied in this article is the retarded Green’s function of the system,

[TABLE]

where $\eta$ is an infinitely small positive number. $\hat{G}$ captures the single-particle propagation of the problem. In the absence of electron-electron interactions (or in a mean field treatement), the knowledge of $\hat{G}$ is sufficient to calculate all the observables including out-of-equilibrium PRL_Wingreen ; elke2017molecular ; datta_1995 . Note that $\hat{G}$ is an infinite matrix. However, only a few of its elements need to be computed, typically the matrix elements between the electrodes at the Fermi energy (conductance) or its diagonal elements at sites of interest (local density of states)knit .

The starting point of our calculation is the Green’s function of the TIS parts,

[TABLE]

Since the TIS are invariant by translational, $\hat{g}_{l}$ can be obtained easily in momentum $\mathbf{k}$-space. To calculate $\hat{g}_{l}$ in real $\mathbf{r}$-space, one must perform a Fourier transform which formally reads

[TABLE]

Performing this Fourier transform cannot be done using standard FFT or other quadrature approaches since the integrand contains Dirac and principal value distributions. Even using a small finite value of $\eta$ to regularize the integrand, direct numerical approaches are bound to fail. We solve this problem by using integration in the complex plane and specially design tools to calculate the poles and residues of $\hat{g}(k_{x},k_{y},k_{z},E)$ where $k_{x}$ is considered as a complex variable while $E$, $k_{y}$ and $k_{z}$ are real parameters. We call the corresponding problem the residue problem. The associated “residue solver” provide the function $\bra{x^{\prime},k_{y},k_{z},E}\hat{g}_{l}\ket{x,k_{y},k_{z},E}$ which is now a well behaved regular function. Subsequent integration over $k_{y}$ and eventually $k_{z}$ can be performed using quadrature rules. However, we shall see that the presence of cusps and kinks make these numerical integrals somewhat delicate and they must be handled with care. Once $\bra{\textbf{r'}}\hat{g}_{l}\ket{\textbf{r}}$ has been obtained, one can use standard Green’s function techniques to combine these elements with the $\hat{W}_{0}$ matrix and calculate the desired elements $\bra{\textbf{r'}}\hat{G}\ket{\textbf{r}}$. Following Ref. knit, , we call this step the glueing sequence. The matrices $\hat{W}_{1}$ and $\hat{W}_{2}$ are dealt with in a similar way before the integration over $k_{z}$ ($\hat{W}_{1}$) or $k_{y}$ ($\hat{W}_{2}$). During the calculation, for instance when a bond is cut between two parts of the system, new (bound) states may appear. These bound states include in particular edge or surface states present at the boundary of topological systems. We use the bound state algorithm developed in Ref. Bound_state_algo, to address this problem. Altogether, our MTIS algorithm consists of a non-trivial combination of the residue solver, the numerical Fourier transform solver, the glueing sequence solver and the bound state solver.

To make the above discussion more concrete, Fig. 2 shows the workflow of a typical 3D calculation. the starting point in Fig. 2a is the Green’s function of a 3D TIS given in momentum space. We use our residue solver to integrate over $k_{x}$ and obtain (b) the Green’s function as a function of space (along $x$) and momentum ($k_{y},k_{z}$). The glueing sequence is used to cut this infinite system into two and obtain a semi-infinite system (c). In this step, we also use our bound state solver to incorporate possible surface states. The semi-infinite system is modified on its surface (green) as well as within its bulk (blue) to account for the various layers that form the material. We obtain (d). We use our numerical integrator to perform the Fourrier transform along $k_{y}$ and obtain (e) which now depends on two spatial variables ($x,y$) and one momentum ($k_{z}$). We use the glueing sequence to add a quasi-one dimensional wire deposited on the surface (red points), we obtain (f). After integration over the last momentum variable $k_{z}$, we finally obtain a 3D system in real space (g). This describes a multilayer system on the surface of which a one-dimensional wire has been deposited. This system is infinite in all the directions where shaded dashed lines are shown. We can use the glueing sequence a last time to modify this system and include impurities (purple circles) or attach the system to another MTIS (here the orange quasi-one dimensional electrode. In practice it is often useful to store the result of (g) in memory or on disk so that different kinds of steps (g)-(h) (e.g. average over disorder) can be performed at very low computational cost. A computation that applies a sequence of steps similar to the one shown in the above schematic to the surface of a Weyl semi-metal is presented in the application section III.

In the next section, we proceed directly to actual numerical calculations for concrete systems and postpone detailed mathematical derivations of the method to Sec. IV and onwards.

III Applications

In this section, we present four applications that demonstrate the usefulness of the MTIS approach:

(A) Friedel oscillations in a two-dimensional electron gas. In this pedagogical example that could be performed fully analytically we study the effect of an impurity embedded in an infinite sheet of 2D material (described by a simple effective-mass Hamiltonian). In a corresponding experiment, the resulting Friedel oscillations could be observed with a scanning tunneling microscopecrommie1993imaging ; harrison1980solid .

(B) Quantum spin Hall effect. This example involves a 2D topological insulatorBernevig1757 that possesses edge states on its boundaries. We calculate the local density of states of a half-plane filled with this material. In our calculation, the system has a unique boundary in contrast to more traditional calculations of finite-width ribbons that would feature two boundaries (with possible overlap between the edge states due to the finite width).

(C) Graphene nanoribbon. The conductance of a graphene nanoribbon is computed. The electrodes are truly two-dimensional, in contrast to standard calculations where quasi-one-dimensional electrodes are used.

(D) Weyl semimetal three-terminal device. The last example is a full-fledged three dimensional computation of the disordered surface of a topological Weyl semi-metal to which quasi-one dimensional electrodes are attached. Our conductance calculation provides direct evidence for the role in transport of Fermi-arcs which are present at the surface of Weyl semi-metals, a topic that has attracted a lot of attention recentlyQuantum_transport_in_Weyl ; Weyl_t_matrix_impurity ; Das_sarma_disorder_weyl .

The flowchart of Fig. 3 illustrates the course of calculations for examples A, C, and D.

III.1 Friedel oscillations in a two-dimensional electron gas

Let us consider an impurity surrounded by a 2D electron gas modeled on an infinite square lattice by the tight-binding Hamiltonian

[TABLE]

where $t$ is the hopping integral between nearest neighbors and $\varepsilon$ the on-site energy shift of the impurity at site $\ket{0,0}$. The TIS part of this model is defined by one on-site and two hopping matrices of size $1\times 1$:

[TABLE]

We calculate the local density of states (LDOS) $\rho(\textbf{r},E)$ in a region surrounding the impurity. This quantity can be measured directly by a scanning tunneling microscope in the tunneling regimecrommie1993imaging . It is related to the Green’s function by

[TABLE]

Fig. 4a shows the LDOS for a finite portion of space as calculated by our method. Note that the underlying physical system is infinite – the result does not suffer from any finite size effects. As expected, the LDOS exhibits Friedel oscillations $\rho(\textbf{r})\sim\cos(2k_{F}r)/r^{2}$ where $k_{F}$ is the Fermi momentum.

To demonstrate the merits of the MTIS approach we contrast it with a more traditional approach: Fig. 4b shows an approximation of the same LDOS function, calculated with Kwantkwant_article for the “cross” geometry shown in panel (d) that consists of a square of size $L\times L$ to which four quasi-one dimensional electodes of width $L$ are attached. Boundary effects that distort the Friedel oscillations are clearly visible near the edges in the colormap (b). In addition, this calculation is much more computationally intensive, with the computational effort scaling as $L^{3}$. As a benchmark, Fig. 4e compares both calculations at the two line cuts shown in Fig. 4a and b. In order to obtain a quantitative match between the two approaches, a very large value of $L$ must be used, as shown in Fig. 4f.

III.2 Quantum spin Hall effect

The second application demonstrates the quantum spin Hall effect within the BHZ model Bernevig1757 for a two-dimensional topological insulator. The continuous BHZ model is described by the Hamiltonian

[TABLE]

with

[TABLE]

and

[TABLE]

where $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of Pauli matrices, and $A$, $B$, $C$, $D$ and $M$ are scalar parameters of the model.

The discretized tight-binding model that corresponds to $h(\textbf{k})$ in the small momentum limit is given by the following onsite and hopping matrices,

[TABLE]

Since the interesting physics occurs at the boundary (topological insulators are gapped in the bulk but have conducting edge states), we consider a semi-infinite 2D sheet.

We consider a sample made of a pristine semi-infinite sheet of this topological material, as shown in the inset of Fig. 5a . Fig. 5a shows the local density of states (LDOS) as a function of energy $E$, $\rho(x,E)=-\text{Tr}\ \ \text{Im}\bra{\textbf{r}}G(E)\ket{\textbf{r}}/\pi$ at three different distances from the edge. We also plot the bulk density of states (DOS) of the infinite sheet (i.e. the value of the LDOS infinitely far from the edge). The latter shows a vanishing DOS inside the gap $[-0.7,0.7]$ of the system. The non-zero values of the LDOS inside the gap for the semi-infinite sheet is due to the presence of the propagating edge states. Note that all quantities here are independant of $y$ as all systems are invariant by translation along that direction. As the distance from the boundary is increased, the edge states amplitude decrease exponentially which leads to a corresponding decrease of the LDOS inside the gap.

Fig. 5b shows a calculation of the LDOS versus $x$ for the three energies marked by crosses in Fig. 5a. We observe a slow (algebraic) convergence to the (independently calculated) bulk limit for the two values of energy inside the spectrum, characteristic of the Friedel oscillations created by the presence of the edge. Inside the gap, we observe a much faster exponential convergence towards zero due to the exponential localization of the edge states. We emphasize again that the computational time to obtain this data is totally independent of the value of $x$ used in the calculation, in sharp contrast with standard approaches.

III.3 Graphene nanoribbon

We now turn to a quantum transport problem: the calculation of the conductance of a graphene nanoribbon connected to two semi-infinite graphene sheets. The sample geometry corresponds to the one shown in Fig. 6d with the width $W$ approaching infinity. For simplicity, we use the standard nearest-neighbor Hamiltonian on a honeycomb lattice,

[TABLE]

but the method applies to any other tight-binding model as well.

Obtaining the conductance of this setup with the help of a traditional quantum transport code (that is limited to quasi-1D leads of finite width) involves a series of calculculations where the semi-infinite sheets are approximated by leads of increasing finite width. The conductance in the thermodynamic limit $W=\infty$ is then extrapolated from these finite-size results. In contrast, our new method works directly in the thermodynamic limit.

We note that this particular system could also be addressed with the technique of Ref. Fabry_perot_graphene, that uses an iterative Lanczos-like approach that takes advantage of the existence of the nanoribbon constriction to restrict the Hamiltonian to the states that are connected to it. This approach is orthogonal (and possibly complementary) to our method which exploits translational invariance.

In order to directly compute observables of the infinitely wide sample, we first calculate the Green’s functions of the two semi-infinite 2D graphene sheets and then proceed to connect them to the central nanoribbon using the Dyson equation. Following the standard formalism datta_1995 , the differential conductance $g_{ll^{\prime}}=dI_{l}/dV_{l^{\prime}}$ between electrode $l^{\prime}$ and $l$ is then given in terms of the total transmission probability $g_{ll^{\prime}}=(e^{2}/h)T_{ll^{\prime}}$ with

[TABLE]

where $\hat{G}$ designates the Green’s function of the total system, and

[TABLE]

the rate matrix expressed in terms of the lead self-energy

[TABLE]

In the above definition, $G_{l}$ is the surface Green’s function of one graphene electrode (here a semi-infinite 2D graphene sheet). $W_{l}$ corresponds to the matrix elements of the Hamiltonian that connect the central ribbon to electrode $l$.

Fig. 6 shows the results of a calculation for a particular geometry of the central graphene ribbon (inset d). The main panel displays the transmission probability $T_{ll^{\prime}}$ between the two leads as a function of energy. We observe that the results of finite-width computations (colored dots) are scattered around the infinite-width values (solid line) that have been obtained with the MTIS method.

For perspective, the dashed line represents the result of a calculation for $W=0$, i.e. for a simple infinite ribbon of fixed width. The conductance is quantized and the total transmission simply counts the number of conducting sub-bands at the Fermi level. In contrast, we note the presence of fluctuations in the other curves. They are caused by the abrupt widening on both ends of the ribbon that induces reflexions which in turn create a Fabry-Perot cavity.Fabry_perot_graphene The associated interference pattern is at the origin of the fluctuations.

At lower energies, we observe a clear difference between the finite-width data and the infinite-width case. This is to be expected since, for small energies, the associated wave length $\lambda\sim 1/E$ is comparable to the ribbon width and hence finite size effects occur. The apparent noise in the finite-width data is due to rapid oscillations in energy. Fig. 6b shows the convergence of the result as a function of the width for a fixed energy. We find the convergence towards infinite-width limit (blue horizontal curve) slow and oscillatory. Consequently, accessing the large-width physics with standard finite-width techniques turns out to be extremely difficult at low energies.

For large values of the energy (typically $E>1$), we find that the finite-width data is almost indistinguishable from the $W=\infty$ curve. A plot of the convergence at $E=1.3$ is shown in Fig. 6c. In this regime, the MTIS approach has nevertheless a clear computational advantage because the computational cost of the finite-width quasi-1D electrodes scales as $W^{3}$. With our current implementation, we find that computing directly the $W=\infty$ limit is faster than the finite-width calculation as soon as $W\geq 100$ for $E=0.5$.

III.4 Weyl semimetal three-terminal device

We conclude this section with an application that showcases the full power of the MTIS approach. We consider a 3D topological Weyl semi-metal with impurities on its surface and calculate the differential conductance of a three-terminal geometry. This computaton combines the difficulties of addressing an infinite 3D system (half of the space is filled with the 3D material) with disorder, in a multiterminal geometry and in presence of surface states (the so-called Fermi arcs). The geometry of the device is shown in Fig. 7a with a top view shown in 7b.

We use a 3D 4-band model of a Weyl semimetal that is defined as Vazifeh2013Electromagnetic

[TABLE]

with

[TABLE]

The Pauli matrices $\tau_{i}$ and $\sigma_{i}$ ($i=x,y,z$) act, respectively, on the orbital and spin degrees of freedom. We use the parameters $t=2$, $t_{z}=1$, $\mu_{0}=-0.1$, $b_{0}=0$ and $\beta=1$. After discretization, the on-site matrix $H_{0}$ and hopping matrices $V_{x}$, $V_{y}$ and $V_{z}$ are given by

[TABLE]

The geometry of Fig. 7a consists of the Weyl semimetal occupying the $x\geqslant 0$ half of the 3D space to which two quasi-1D electrodes are attached. We also deposit some impurities on the surface. The complex sequence of submethods used to calculate the differential conductance of this geometry is shown in Fig. 3. The quasi-1D electrodes are considered to be made from a normal metal described, respectively, by on-site and hopping parameters $t=1$ and $\mu=0$. The electrodes have a square cross section of $5\times 5$ sites. In the following, the quasi-1D electrodes are labeled as 1 and 2, and the infinite Weyl semi-metal is labeled as 0. The conductance matrix that relates the current $I_{a}$ to the voltage $V_{b}$ applied to lead $b$ takes the form

[TABLE]

where $N_{a}=\sum_{j}T_{lj}=\sum_{j}T_{jl}$ is the number of channels in the quasi-1D electrode $a$. The conductance matrix conserves current and is “gauge invariant” in the Büttiker sense, i.e. raising simultaneously all voltages by the same amount does not change the currentButtiker1995 . Notice the special treatment of electrode [math]: since $R_{00}$ and $N_{0}$ are both infinite, only their difference is a well-defined quantity. The total reflection terms $R_{aa}$ are given byWimmer_thesis

[TABLE]

where $G_{a}$ is the Green’s function of the full system at the surface of a quasi-1D electrode.

We now present the results of the calculation for a clean surface. We assume the geometry of Fig. 7b where one electrode (serving as Ohmic contact) is fixed on the surface while the second one scans the surface modeling the tip of a STM (scanning tunneling microscope). The letters ABCDEFG indicate the trajectory of the STM tip on the surface. Fig. 7c shows the transmissions between the different electrodes as a function of the position of the STM tip on its trajectory for a value of energy close to the Weyl points. At this energy, the bulk density of states is very small and transport is entirely dominated by the surface states, the so-called Fermi arcs. These arcs are extremely anisotropic with a dispersion relation $E=vk_{y}$ so that they only propagate along the $y$ axis in the positive direction. It follows that transport between electrodes $1$ and $2$ is only possible when they are aligned along the $y$-direction; and only from $1$ to $2$, not from $2$ to $1$. This is illustrated by the curves between A and B where $T_{21}\approx 0.55$ and $T_{12}\approx 0$. As soon as this alignment is lost, both $T_{21}$ and $T_{12}$ vanish and all the current disappears into the substrate. Note that the above results are closely linked to the band structure of the model considered here which only contains one pair of Weyl points. In real materials, there might be several pairs of such pointsxu2015discovery and surface reconstruction can also lead to a bending of the Fermi arcsdeng2016experimental which would lead to a decrease of the extreme anisotropy observed in our numerics.

Intermediate steps in the calculation of the above transport curve involve the computation of the Fermi arcs wavefunction, first in momentum space and then in real space. This information is necessary to analyze the final calculation of observables. The color plot of Fig. 8 shows the local density of states of the pristine surface,

[TABLE]

while the blue line indicates the position of the Fermi arc. From such plots, one can extract for example the dispersion relation of the Fermi arcs.

Fig. 9 shows the different transmissions as a function of energy $E$ when the STM tip is positioned at point A’. The transmission towards the bulk increases quadratically which is consistent with a linear density of states in the bulk (dashed lines). More interestingly, the transmissions $T_{21}$ and $T_{12}$ are very weakly affected by the presence of the bulk states and remain very anisotropic ($T_{12}\approx 0$). Hence, transport on the surface is very weakly affected by the presence of the bulk states and remains dominated by the Fermi arcsPhysRevB.93.235127 .

Fig. 7d shows the differential conductance measured across the two top contacts in the following setup: one injects a current $I$ into electrode $1$ and retrieves the same current from electrode $2$ (i.e. vanishing current leaves electrode 0). The strong anisotropy observed in the two non-local resistances $(V_{1}-V_{0})/I$ and $(V_{2}-V_{0})/I$ bears the signature of the Fermi arc anisotropy.

We now turn to the study of the resilience of the above picture in the presence of disorder. While disorder in bulk Weyl semimetals has been the focus of several studiesPhysRevLett.113.026602 ; annurev-conmatphys-033117-054037 ; PhysRevLett.115.076601 ; PhysRevLett.115.246603 , the effect on the surface has received much less attentionDas_sarma_disorder_weyl . We consider discrete impurities randomly scattered at the surface in the region around the leads (see Fig. 10a).

We model the impurities by modifying $5\%$ of the on-site Hamiltonian matrices located around the leads (20 to 30 sites) to which we add a shift in energy $w\cdot\mathbb{1}_{4}$ (Fig. 10b and c or a fully random on-site matrix $w\sum_{i,j}h_{i,j}\sigma_{i}\tau_{j}$ where $h_{i,j}$ are random numbers from the $[-1,1]$ interval (Fig. 10d). Different traces correspond to different samples.

In the low energy regime of Fig. 10b, at energy $0.02$, the effect of the impurities is only significant for a large value of $w$ comparable to the bandwidth of the model (equal to approximately $5$), which indicates that the anisotropic transport on the surface is resistant to the presence of disorder. This is not unexpected since at this energy the density of bulk states to scatter to is very low. However, the same situation is shown in Fig. 10b for an energy ten times larger, $E=0.2$. The resilience to disorder persists when we use the random matrix disorder of Fig. 10d that breaks all possible symmetries. We conclude that the topological protection of the Fermi arcs can be observed directly in the differential conductance which is not perturbed by the presence of the bulk states. Our observations are compatible with the statement that disorder only renormalizes the dispersion relationDas_sarma_disorder_weyl , which remains non-zero even close to Anderson localization.

IV The Residue problem

We turn to the detailed derivation of the formalism and algorithms. In a first step we calculate a finite set of matrix elements for the Green’s function of a TIS, i.e. $\bra{\textbf{r'},\mu}\hat{g}_{l}(E)\ket{\textbf{r},\nu}$ for a set of values $(\textbf{r'},\mu,\textbf{r},\nu)$. The Green’s function of the TIS is formally defined as the inverse of the Hamiltonian (shifted in energy):

[TABLE]

Since the TIS is invariant by translation, $\hat{H}_{l}^{\infty}$ can be diagonalized in momentum space. The momentum states $|\alpha\textbf{k}\rangle$ satisfy $\hat{H}_{l}^{\infty}|\alpha\textbf{k}\rangle=E|\alpha\textbf{k}\rangle$ and take the form

[TABLE]

where the eigenvectors $\Psi_{\alpha\textbf{k}}$ are eigenstates

[TABLE]

of the momentum Hamiltonian

[TABLE]

We arrive at

[TABLE]

Eq. (26) forms the starting point of our calculation. The momentum Hamiltonian $\mathcal{H}(k)$ is a small $N_{o}\times N_{o}$ matrix that can easily be diagonalized numerically. Hence, at first glance, the Fourrier transform in Eq. (26) could be performed numerically using stanrd quadrature rules such as Simpson. Such an approach is too naive however, since the fraction $1/[E-E(\textbf{k})+i\eta]$ diverges when $E(\textbf{k})$ crosses $E$. As $\lim_{\eta\rightarrow 0}1/(X+i\eta)$ equals $P(1/X)-i\pi\delta(X)$, i.e. the principal part and a Dirac distribution, the integral over at least one of the momentum variables (here $k_{x}$) must be calculated analytically. This is the residue problem which we discuss in this section.

For a fixed value of $k_{y}$ and $k_{z}$, we are left with a 1D TIS problem described by an on-site matrix $H$ and a hopping matrix $V$:

[TABLE]

IV.1 Formulation of the residue problem

The 1D translationally invariant Hamiltonian $\hat{H}_{\rm\mathrm{1D}}$ of the system,

[TABLE]

consists of repeated blocks as in Fig. 1a,b. Each block is described by $N_{o}\times N_{o}$ matrices $H$ (on-site) and $V$ (hopping to the next cell). We seek a finite number of elements of the 1D Green’s function,

[TABLE]

Introducing the $N_{o}\times N_{o}$ matrix $g^{\mathrm{1D}}_{x,x^{\prime}}(E)$ whose elements are given by

[TABLE]

the residue problem consists of calculating the integral

[TABLE]

or, in a more compact form,

[TABLE]

with

[TABLE]

IV.2 Solution to the residue problem

Here we provide without justification the solution of the residue problem. The proof of this result is given in the following section.

The first step of the algorithm is to solve the generalized eigenvalue problem

[TABLE]

where $E$ is an input and we seek the values of $k$ for which the above equation has a solution $\phi$. We emphasize that this is a very different problem from diagonalizing $\mathcal{H}(k)$ for a given value of $k$. Introducing explicitly the form of $\mathcal{H}(k)$ and designating $e^{ik}$ as $\lambda$, Eq. (34) takes the form

[TABLE]

which is a generalized eigenproblem. Solving this problem numerically can be handled by the Kwant numerical packagekwant_article for example.

The second step is to sort the resulting set of eigenvectors/eigenvalues $(\phi_{a},\lambda_{a})$ into two categories: left-goers and right-goers. Eigenstates with $\lambda_{a}<1$ ($\lambda_{a}>1$) are attributed to the right- (left-) goers category. Eigenstates with $|\lambda_{a}|=1$ are classified according to the velocity

[TABLE]

States with positive (negative) velocity are right- (left-) goers.

The third step is to construct the projectors onto the different subspaces spanned by the eigenvectors $\phi_{a}$:

[TABLE]

The special case of projector $P_{0}$ corresponds to the projector onto the kernel of $V$ ($VP_{0}=0$). Introducing the matrix

[TABLE]

the result of the calculation reads

[TABLE]

where the sum is extended to the right-goers for $x\geq 0$ and to the left-goers for $x<0$. Eq. (39) is the chief analytical result of this article. It provides a fast and stable solution to the residue problem. The computation of a set of different values of $x$ involves only little numerics with small ($N_{o}\times N_{o}$) matrices. The matrices in front of the $\lambda^{x}$ terms can be cached so that the computational effort to obtain $N_{x}$ different values of $x$ (for a fixed energy) scales as $N_{o}^{4}+N_{o}^{2}N_{x}$.

IV.3 Proof of Eq. (39)

In a first step, we extend the integral of Eq. (32) onto the complex $k$-plane. For $x\geqslant 0$, we use the red contour shown in Fig. 11 which consists of four branches: the integral over $\Gamma_{1}$ is the orginal integral, the integrals over $\Gamma_{2}$ and $\Gamma_{4}$ compensate each other by symmetry of the integrand and the integral over $\Gamma_{3}$ is evaluated in the limit where the horizontal segment $\Gamma_{3}$ is shifted towards $i\infty$. For $x<0$, we use the mirror-symmetric green contour. For $x\neq 0$, the factor $\lambda^{x}$ in the integrand is exponentially small which makes the integral over $\Gamma_{3}$ to vanish in the limit $\kappa\rightarrow\infty$, where $\kappa$ is the imaginary part of the momentum. Applying Residue theorem to this contour integral, we arrive at

[TABLE]

where the first terms accounts for the residues of the integrand evaluated at its poles and the second to the integral over the $\Gamma_{3}$ segment.

Let us focus on the first part of this expression, the residues. The first step consist in finding the poles of the integrand of Eq. (32). These poles corresponds to the values $k_{a}$ solution of Eq. (34). These poles are either in the upper half plane $\Im m\ k_{a}>0$ (right goers) or lower half plane $\Im m\ k_{a}<0$ (left goers). The fate of the non-evanescent solutions with $\Im m\ k=0$ is found by remembering the existence of the infinitely small positive $\eta$. To obtain the residues, we need to calculate the expansion of the integrand close to the pole. We start by expanding ${\cal H}(k)$ around $k_{a}$,

[TABLE]

where we note that high derivatives of ${\cal H}(k_{a})$ are simply related to the lower ones: $\partial_{k}^{n}{\cal H}(k_{a})=-\partial^{n-2}_{k}{\cal H}(k_{a})$ for $n\geq 3$. We now seek the expansion of $1/[E-{\cal H}(k_{a}+q)]$ around $q=0$. The key point to notice here is that since $[E-{\cal H}(k_{a})]$ is not invertible, this expansion has a term proportional to $1/q$ that we need to calculate to apply the residue theorem. The algebra to obtain in a systematic way the coefficients of the developpement of $1/[E-{\cal H}(k_{a}+q)]$ is given in Appendix A. We obtain

[TABLE]

with

[TABLE]

Applying residue theorem with the above expression for the matrix residue $D_{-1}$ provides the first term of Eq. (40). Altough we did not encounter the case in practice, it is possible that the matrix $E-{\cal H}(k_{a})-\partial_{k}{\cal H}(k_{a})P_{\lambda_{a}}$ above is not invertible. In that case the sum in Eq. (42) starts at $-2$. Appendix A provides a systematic iterative construction of the residue in that case. Appendix A also provides a general algorithm to calculate the other coefficients $D_{n}$ if ever needed.

Let us now calculate the second term of Eq. (40), the integral over the $\Gamma_{3}$ segent of the contour which is only present when $x=0$. When $\kappa$ is large, one gets from the definition of $\mathcal{H}(k)$,

[TABLE]

We need to integrate $1/[E-\mathcal{H}(k+i\kappa)]$ over $k$. If the matrix $V$ is invertible, then the second term of the development is not necessary: $1/[E-\mathcal{H}(k+i\kappa)]\approx-e^{-\kappa+ik}V^{-1}$ and the integrand (hence the integral) vanishes in the limit $\kappa\rightarrow+\infty$. However, if $V$ is not invertible, the expansion of $1/[-V+(E-H_{0})e^{-\kappa+ik}]$ in power of $e^{-\kappa+ik}$ has a pole $\propto e^{\kappa-ik}$ which provides a term of order $\mathcal{O}(1)$ in the integrand. Using the result of Appendix A, we obtain

[TABLE]

It follows that the integrand of the second term of Eq. (40) is $k$-independent and can be calculated explicitely. This concludes the proof of Eq. (39).

V Numerical Fourier transform

Once the integration over $k_{x}$ has been performed using the complex integration technique of the previous section, the TIS Green’s function of Eq. (26) takes the form

[TABLE]

where $g^{\rm 1D}_{x}(k_{y},k_{z})$ is the result of the residue solver. In this section, we focus on the case where the system covers the full 3D space, $-\infty\leq x\leq+\infty$. To study, e.g., a surface with $-\infty\leq x\leq 0$, one introduces intermediate steps before the integration over $k_{y}$ and $k_{z}$. These steps are discussed in the sections further below.

To integrate over $k_{y}$ and $k_{z}$ we use regular numerical integration (quadrature) where the interval $[-\pi,+\pi]$ is broken up into small subintervals and the integrand is approximated by a polynomial in each subinterval. A typical example of the integrand is shown in Fig. 12. We observe that it possesses a number of cusps and kinks, so that integration over momentum is somewhat delicate. To obtain reliable results in a robust way, we have found that the doubly adaptative scheme of Ref. quad_pedro, is particularly efficient. Of prime importance is the algorithm used to evaluate the residual error of the integral, that allows one to compute $g^{\mathrm{3D}}_{x,y,z}(E)$ to a given precision. We plan to release an efficient implementation of the resulting algorithm as open source softwareC.Groth-in-preparation .

VI The glueing sequence

The Dyson equation can be used to introduce modifications to the TIS Green’s function. This approach, that we call pictorially “the glueing sequence”, has been discussed in different contexts Bruno_RGF_2005 ; knit ; molecular_electronics ; fetter2012quantum . We follow the presentation of Ref. knit, . An alternative, equivalent, formalism uses the so-called T-matrix approachPhysRevB.69.134517 ; PhysRevB.86.115433 .

Let us suppose that the Hamiltonian of the system splits into two contributions $\hat{H}$ and $\hat{W}$. Typically $\hat{H}$ is the Hamiltonian of the TIS system and $\hat{W}$ a perturbation that involves only a finite number of sites. (However, we shall see below that we can be slightly less restrictive.) We suppose that we have already calculated the Green’s function of the TIS,

[TABLE]

and we seek the Green’s function of the full Hamiltonian,

[TABLE]

The starting point is a simple equation that is valid for any two matrices $A$ and $B$:

[TABLE]

In the present context, it translates into

[TABLE]

which is known as the Dyson equation. Eq. (50) involves infinite matrices. Let us introduce the set of points $\Omega$ that is the support of the matrix $\hat{W}$: $\forall i,j\notin\Omega,\ W_{ij}=0$. The crucial property of the Dyson equation is that its projection on any finite set of points $\Theta$ that contains $\Omega$ forms a close set of equations,

[TABLE]

that can be solved using standard linear algebra routines. More formally, we introduce the projectors $P_{\Omega}$ and $P_{\Theta}$ defined as

[TABLE]

with a similar definition for $P_{\Omega}$. The finite matrix $\hat{G}_{\Omega\Theta}=P_{\Omega}\hat{G}P_{\Theta}$ is then given by

[TABLE]

In a second step, one calculates

[TABLE]

The above sequence allows to introduce any finite number of modifications to a system whose corresponding Green’s function elements are known. This is a straightforward generalization of the basic sequence at the core of the recursive Green’s function algorithm. Denoting be $N_{\Omega}$ and $N_{\Theta}$ the number of elements of the respective sets, the computational cost of the glueing sequence scales as $N_{\Omega}N_{\Theta}^{2}$ (for non-local properties such as transport) or $N_{\Omega}^{2}N_{\Theta}$ (for local properties such as local density of states). This cost is independent of the system size (which is infinite) and only depends on the number of modifications and points of interest.

In this article, we use the glueing sequence in the following situations.

*Point-like impurities. * The simplest application of the glueing sequence is to modify selected on-site energies or hopping elements. In this case, one strictly follows the sequence given above. An example of such a calculation is the inclusion of on-site disorder in the application to Weyl semi-metals discussed in Sec. III.4. If one is interested in calculating averages over impurity configurations, then only the last stage of the calculation, the glueing sequence, needs to be recalculated for each sample. Hence, these averages can be relatively cheap computationally. A similar approach has been used previously in the context of disordered grapheneOstrovsky_2010 ; Schelter_2011 . In Ref. Ostrovsky_2010, ; Schelter_2011, the TIS Green’s function has been obtained analytically through a route that is similar to the route that we have followed numerically.

*Attaching two systems together * The glueing sequence can also be used to connect two systems together, e.g. the semi-infinite 2D graphene sheet and the graphene nanoribbon shown in Fig. 6. In this case the Dyson equation has an additional structure, as the unperturbed Greeen’s function is block-diagonal (the two subsystems are initially unconnected). This structure may be used to simplify the algebra of the glueing sequence, see Ref. knit, .

Creating a multilayer system. The matrix $\hat{W}_{2}$ introduced in Sec. II has an infinite number of terms, but it conserves momentum along the two directions $y$ and $z$. It can be cast in the form

[TABLE]

Hence, for a given value of $k_{y},k_{z}$, i.e. before we have performed the numerical integration on these variables, $\hat{W}_{2}$ takes the form of a finite perturbation and we may apply the glueing sequence. Such a step allows to create, e.g., multilayer systems. The above equation is actually not fully general. A straightforward extension is to consider matrices $W_{ij}(k_{y},k_{z})$ that include an explicit dependance on momentum. Terms of the class $\hat{W}_{1}$ are treated in a similar fashion.

Slicing a system in two. A particular example of perturbation of the $\hat{W}_{2}$ class is when the perturbation is exactly opposite to one bond $V_{x}$.

[TABLE]

In that case the perturbation slices the system and creates two disconnected systems. This is what we use to create the semi-infinite 2D graphene sheet as well as the surface of the Weyl semi-metal. It is important to notice that upon slicing, one may create bound states (1D), edge states (2D) or surface states (3D). In particular in topological materials, these states (that we refer collectively as bound states) are always present. Bound states may also appear in the multilayer example above. Properly dealing with bound states is the topic of the next section. In the applications shown in this article, we focus on surfaces in the [100] direction. The same technique can be extended to other surfaces, such as [110], by simply enlarging the unit cell so that [110] becomes effectively [100] for the large unit cell.

VII The bound state problem

The last problem that remains to be addressed is the appearance of bound states in the system. These states appear in a wide variety of situations that include slicing the system (creation of a surface), multilayers (quantum wells), Josephson junctions or around impurities. In the present article, the edge states of the quantum spin Hall effect and the surface states (Fermi arcs) of the Weyl semi-metals are (generalized) bound states. Bound states do not hybridize with the continuum and must therefore be handled as separate contributions. In MTIS, bound states may be invariant by translation upon zero (true bound states), one (edge states) or two directions (surface states). To simplify the discussion, we focus below on surface states. The results can be straightforwardly extended to the other situations.

Suppose that we have used the residue solver to construct the TIS Green’s function $g_{x}(E,k_{y},k_{z})$ for fixed values of the transverse momentum $(k_{y},k_{z})$. In a second step, we have used the glueing sequence to slice the system and obtain $G_{x,x^{\prime}}(E,k_{y},k_{z})$, where the system terminates at $x=0$ (Fig. 2c). The contribution of a bound state $\psi_{\alpha}(k_{y},k_{z})$ with energy $E(k_{y},k_{z})$ to $G_{x,x^{\prime}}(E,k_{y},k_{z})$ takes the form

[TABLE]

Since $\lim_{\eta\to 0}1/(X+i\eta)$ equals $P(1/X)-i\pi\delta(X)$, the presence of a bound state leads to two complications: (i) The principal value cannot be integrated with numerical routines as the integral is formally divergent. (ii) The Dirac function associated with the divergence at $E=E_{\alpha}(k_{y},k_{z})$ is not captured by $G_{x,x^{\prime}}(E,k_{y},k_{z})$. Numerically, one only observes a numerical instability of $G_{x,x^{\prime}}(E,k_{y},k_{z})$ for values of $(k_{y},k_{z})$ such that $E=E_{\alpha}(k_{y},k_{z})$. An example of the integrand that results after the slicing sequence is shown in Fig. 13 where we indeed observe the behaviour discussed above.

To proceed, one needs first to calculate the bound states independently. We have devoted an earlier article to this problem to which we refer for all technical detailsBound_state_algo . We suppose that we have computed the energy $E_{\alpha}(k_{y},k_{z})$ and the associated state $\psi_{\alpha}(x,k_{y},k_{z})$ of the semi-infinite problem. Once this is done, we numerically inverse the function $E_{\alpha}(k_{y},k_{z})$ to obtain $k_{\alpha}=\bar{k}_{y}(E,k_{z})$. Our goal is to perform the integration over $k_{y}$ and calculate

[TABLE]

(i) Regularizing the principal value. The problem of the presence of the principal value is easily handled once one has evaluated the function $k_{\alpha}=\bar{k}_{y}(E,k_{z})$. We follow the standard procedure for handling principal values and select a small interval $[k_{\alpha}-\epsilon,k_{\alpha}-\epsilon]$. In this interval, we use the symmetrized integrand

[TABLE]

where $F(k_{y})=e^{ik_{y}(y-y^{\prime})}G_{x,x^{\prime}}(E,k_{y},k_{z})$. The symmetrized integrand that gives the principal value of the real part plotted in Fig. 13 is displayed in the inset of the same figure. We find that although the original integral was divergent, this procedure correctly regularizes the integral. Substracting $\psi_{\alpha}\psi_{\alpha}^{\dagger}/(k_{\alpha}-k_{y})$ from $G_{x,x^{\prime}}(E,k_{y},k_{z})$ could also regularize the integral. Indeed, the dashed line of Fig. 13 perfectly fits the divergence of the Green’s function at the bound state.

(ii) Contribution from the Dirac function. We now evaluate the contribution of the bound state due to the Dirac function,

[TABLE]

where we have droped the dependance on $k_{z}$ for compactness. Performing this integral, we arrive at

[TABLE]

The last step is to calculate $\partial E_{\alpha}/\partial k_{y}$. In the case of a simple slicing where the system is invariant by translation away from the boundary, it is given by

[TABLE]

where $\lambda_{\alpha}$ is the evanescent momentum of the bound state. For more complicated cases, e.g. a coated surface, there is a contribution from the part that is not invariant by translation, that plays the role of the scattering region in Ref. Bound_state_algo, .

To summarize this section, the result of the integral over $k_{y}$ is given by

[TABLE]

Each of these integrals can now be performed using the numerical techniques discussed in Sec. V.

VIII Conclusion

Despite being ubiquitous, quantum transport simulations face severe limitations in a number of situations where large systems are to be simulated. This is in particular the case for most 3D systems. The Mostly Translational Invariant Systems (MTIS) that we have discussed in this article cover a significant fraction of systems of interest that were unaccessible to simulations up to now. We have presented a general method that can handle arbitrary MTIS, and demonstrated the power of the approach in a number of situations. The advantage of the MTIS approach stems from the fact that one works directly in the thermodynamic limit. Not only does the computational cost not depend on the size of the (infinite) system, but it is also independent of distances such as the distance between two impurities. It is therefore a natural tool for systems that possess several different characteristic length scales. The accuracy of the method is determined by that of the numerical integration and can be tuned at will.

We have studied the transport properties of the surface of a Weyl semi-metal. In our calculation, there is a single surface present as the opposite second surface that would be present in a finite size calculation is sent infinitely far away. Removing the second surface greatly simplifies the calculations and their interpretations. Our approach could be used for other surface problems such as the formation of Majorana bound states around magnetic impurities on top of a superconductorCren_2015 or the study of the Dirac cones at the surface of a 3D topological insulatorshasan_2010 or more generally for simulating scanning tunneling microscope (STM) experiments.

For quantum transport, MTIS provides the possibility to use 2D or 3D electrodes such as the semi-infinite graphene sheet that we have presented. These electrodes are often more realistic than the quasi-1D electrodes used in almost all approaches. They are also less computationally intensive in most situations. We expect these electrodes to quickly become a standard feature of quantum transport toolkits.

MTIS also covers many applications that we have barely discussed so far. For pure bulk systems, they could be used for the study of defects, (RKKY) interaction between magnetic impurities or to calculate the collision integrals that enter the semi-classical Boltzman equation. The MTIS approach can be trivially combined with the recursive Green’s function approach and its generalizations to build complex geometries.

Acknowledgments

Warm thank to Arthur Waintal for making Fig.7a. We acknowledge interesting discussions with Anton Akhmerov and Michael Wimmer.

This project has been funded by the US Office of Naval Research (ONR) and the graphene flagship through the ANR GRANSPORT.

Appendix A Laurent expansion of the inverse of a matrix

Let $A(x)$ be a matrix that depends on a continuous parameter $x$. $A(x)$ is defined in terms of its expansion in power of $x$,

[TABLE]

Let $G$ be the inverse of $A$,

[TABLE]

The aim of this appendix is to calculate the terms $G_{k}$ of the Laurent expansion of $G(x)$,

[TABLE]

in terms of the expansion of $A(x)$. If $A_{0}$ is inversible, this problem is trivial: one can simply insert the expansions of $A(x)$ and $G(x)$ into Eq. (65) and identify the terms one by one. One obtains

[TABLE]

In this appendix, we focus on the case where $A_{0}$ is not invertible which leads to the appearance of a term $G_{-1}$. Let $P$ be the projector on the Kernel of $A_{0}$, i.e. $A_{0}P=0$, and $Q=\mathds{1}-P$. To proceed, we write Eq. (65) in block form in the $P$ and $Q$ subblocks:

[TABLE]

This expands to

[TABLE]

Introducing the new variable

[TABLE]

and the new series

[TABLE]

with

[TABLE]

allows to map the problem of calculating $\overline{G}$ to a problem that has the same structure as Eq. (65):

[TABLE]

If $B_{0}=A_{1}P+A_{0}Q=A_{1}P+A_{0}$ is invertible, then a term-by-term identification of Eq. (73) leads to

[TABLE]

Eq. (70) can be inverted to

[TABLE]

and we finally arrive at

[TABLE]

Eq. (78) is the central result used for the construction of the residue solver. In all the calculations that we have performed, we have always found $B_{0}$ to be inversible. If $B_{0}$ is not inversible, the same construction that has been applied to Eq. (65) can be performed on Eq. (73). In that case, one would introduce the projector $P^{\prime}$ onto the kernel of $B_{0}$ and proceed to expand $\overline{\overline{G}}=xP^{\prime}\overline{G}+(1-P^{\prime})\overline{G}$. In that case, the developpement of $G(x)$ would start at $G_{-2}$. The same procedure can be extended iteratively to the case where the development starts at $G_{-n}$.