A numerical method to efficiently calculate the transport properties of large systems: an algorithm optimized for sparse linear solvers

TL;DR

This paper introduces an efficient numerical method for calculating electronic transport in large nanostructures by leveraging sparse linear algebra, significantly reducing computational complexity.

Contribution

The paper develops a wave-function matching method optimized for sparse solvers, enabling scalable transport calculations for large systems with improved efficiency.

Findings

Operation count scales with N_S x N_P for large systems

Method effectively handles large nanostructures

Sparse linear solver accelerates calculations

Abstract

We present a self-contained description of the wave-function matching (WFM) method to calculate electronic quantum transport properties of nanostructures using the Landauer-B\"uttiker approach. The method is based on a partition of the system between a central region ("conductor") containing $N_S$ sites and an asymptotic region ("leads") characterized by $N_P$ open channels. The two subsystems are linearly coupled and solved simultaneously using an efficient sparse linear solver. Invoking the sparsity of the Hamiltonian matrix representation of the central region, we show that the number of operations required by the WFM method in conductance calculations scales with $\sim N_S\times N_P$ for large $N_S$.

Figures

Equations

$$H=\sum _{j,j^{\prime }}{H}_{j,j^{\prime }}|j⟩⟨j^{\prime }|,$$

$$H|{\Psi }_m^{\pm }(E)⟩=E|{\Psi }_m^{\pm }(E)⟩,$$

$$Q=\sum _{\mu }|{\varphi }_{\mu }⟩⟨{\varphi }_{\mu }|$$

$$P=\sum _{m\in \alpha }\int dE|{\chi }_m(E)⟩⟨{\chi }_m(E)|,$$

$$H=H_{PP}+H_{QQ}+(H_{PQ}+H_{QP}),$$

$$P|{\Psi }_m^{\pm }(E)⟩$$

$$Q|{\Psi }_m^{\pm }(E)⟩$$

$$S_{m{m}^{\prime }}(E)={\delta }_{m{m}^{\prime }}-2i\pi {\rho }_m^{1/2}(E)\sum _{\mu {\mu }^{\prime }}[H_{PQ}{]}_{m\mu }{\left[\frac{1}{E-H_{QQ}-{\Sigma }^{+}(E)}\right]}_{\mu {\mu }^{\prime }}[H_{QP}{]}_{{\mu }^{\prime }{m}^{\prime }}{\rho }_{m^{\prime }}^{1/2}(E),$$

$${\mathcal{G}}_{\alpha \beta }=\frac{e^2}{h}{\int }_{-\infty }^{\infty }dE\left(-\frac{\partial f}{\partial E}\right)T_{\alpha \beta }(E),$$

$$T_{\alpha \beta }(E)=\sum _{\begin{array}{c}n\in \alpha \\ m\in \beta \end{array}}{\left|S_{nm}(E)\right|}^2,$$

$$H=\left(\begin{array}{ccc}{H}_L& V_{LS}& 0\\ V_{SL}& H_S& V_{SR}\\ 0& V_{RS}& H_R\end{array}\right).$$

$$H_L=\left(\begin{array}{cccc}\ddots & \ddots & & \\ \ddots & H_l& V_l& \\ & V_l^{†}& H_l& V_l\\ & & V_l^{†}& H_l\end{array}\right)H_R=\left(\begin{array}{cccc}{H}_r& V_r^{†}& & \\ V_r& H_r& V_r^{†}& \\ & V_r& H_r& \ddots \\ & & \ddots & \ddots \end{array}\right),$$

$${\small H=\left(\begin{array}[]{c:c c:c c c}H_{S}-E&V_{SL}^{\dagger}&V_{SR}^{\dagger}&&&\\ \hdashline V_{SL}&H_{l}-E&0&V_{l}^{\dagger}&0&\\ V_{SR}&0&H_{r}-E&0&V_{r}^{\dagger}&\\ \hdashline&V_{l}&0&H_{l}-E&0&\\ &0&V_{r}&0&H_{r}-E&\ddots\\ &&&&\ddots&\ddots\end{array}\right).}$$

$$\left(\begin{array}{ccccc}{H}_S-E& V_{TS}^{†}& & & \\ V_{TS}& H_T-E& V_T^{†}& & \\ & V_T& H_T-E& \ddots & \\ & & \ddots & \ddots & \end{array}\right)\left(\begin{array}{c}{\psi }_S\\ {\psi }_0\\ {\psi }_1\\ ⋮\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ ⋮\end{array}\right),$$

$$V_T{\varphi }_{n-1}+(H_T-E){\varphi }_n+V_T^{†}{\varphi }_{n+1}=0,$$

$${\varphi }_n=\chi {\lambda }^n,$$

$$V_T\chi +(H_T-E)\chi \lambda +V_T^{†}\chi {\lambda }^2=0.$$

$${\chi }^{\prime }\equiv {\lambda }^{-1}{V}_T\chi$$

$$\left(\begin{array}{cc}{H}_T-E& \hat{1}\\ V_T& \hat{0}\end{array}\right)\left(\begin{array}{c}\chi \\ {\chi }^{\prime }\end{array}\right)=\lambda \left(\begin{array}{cc}-V_T^{†}& \hat{0}\\ \hat{0}& \hat{1}\end{array}\right)\left(\begin{array}{c}\chi \\ {\chi }^{\prime }\end{array}\right).$$

$${\chi }_p=\left(\cdots ,0,{\chi }_p^{\alpha },0,\cdots \right),$$

$$j_p$$

$${\psi }_n=\sum _{q=1}^{N_P}{A}_q^{-}{\chi }_q^{-}({\lambda }_q^{-}{)}^n+\sum _{p=1}^{N_P}{A}_p^{+}{\chi }_p^{+}({\lambda }_p^{+}{)}^n,$$

$$\left(\begin{array}{c}{A}_1^{+}\\ A_2^{+}\\ ⋮\\ A_{N_P}^{+}\end{array}\right)=\tilde{S}\left(\begin{array}{c}{A}_1^{-}\\ A_2^{-}\\ ⋮\\ A_{N_P}^{-}\end{array}\right).$$

$${\psi }_{nq}$$

$$(H_S-E){\psi }_{Sq}+V_{TS}^{†}{\psi }_{0q}=0,$$

$$V_{TS}{\psi }_{Sq}+(H_T-E){\psi }_{0q}+V_T^{†}{\psi }_{1q}=0,$$

$$V_{TS}{\psi }_{Sq}=V_T{\psi }_{-1q},$$

$${{\chi }^{\prime }}_q^{\pm }={\left({\lambda }_q^{\pm }\right)}^{-1}{V}_T{\chi }_q^{\pm },$$

$$V_{TS}{\psi }_{Sq}={{\chi }^{\prime }}_q^{-}+\sum _{p=1}^{N_P}{{\chi }^{\prime }}_p^{+}{\tilde{S}}_{pq}={{\chi }^{\prime }}_q^{-}+{{\chi }^{\prime }}^{+}{\tilde{S}}_q,$$

$${\psi }_{0q}={\chi }_q^{-}+{\chi }^{+}{\tilde{S}}_q.$$

$$\left(\begin{array}{cc}{H}_S-E& V_{TS}^{†}{\chi }^{+}\\ V_{TS}& -{{\chi }^{\prime }}^{+}\end{array}\right)\left(\begin{array}{c}{\psi }_{Sq}\\ {\tilde{S}}_q\end{array}\right)=\left(\begin{array}{c}-V_{TS}^{†}{\chi }_q^{-}\\ {{\chi }^{\prime }}_q^{-}\end{array}\right).$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Advances in the fabrication of high-quality samples at the micro and nanoscale paved the way for the discovery of unusual electronic transport properties. As a consequence, the demand for numerical methods to realistically describe such systems on an atomistic/microscopic basis has dramatically increased. At low temperatures, experiments have reported electronic coherence lengths as long as tens of microns, typically larger than the characteristic sample size [1, 2, 3, 4]. In such mesoscopic systems, the electronic transport is dominated by quantum interference. In this paper, we critically analyze the Wave Function Matching (WFM) method [5, 6, 7], whose numerical implementations allow to efficiently compute the quantum transport properties of electrons in nanostructures, modeling realistic sample sizes and non-trivial geometries.

Quantum electronic transport in mesoscopic systems is usually described by the Landauer-Büttiker approach [8], that gives a simple relation between the conductance and the quantum transmission coefficients of a single-particle scattering problem. In other words, the problem is reduced to solving a Schrödinger equation for an open quantum system. We show in this paper that the WFM method is one of the most efficient numerical tools for this task. The latter introduces a partition between a central or scattering region (“conductor”) and the asymptotic one (“leads” or terminals) and, by matching the corresponding wave function at the partition boundaries, gives the system scattering matrix $S$ [9].

Alternatively, transport properties in mesoscopic systems can be calculated using non-equilibrium Green’s function (NEGF) techniques [8]. This formalism is widely used due to its successful combination with Density Function Theory [10, 11, 12, 13, 14, 15, 16]. At the single-particle level, NEGF is equivalent to the Landauer-Büttiker approach (see, for instance, Ref. [17]). The standard method to compute transport properties in large systems using NEGF is the Recursive Green’s Function (RGF) method [18, 19, 20]. The latter takes advantage of the sparsity of the system Hamiltonian to partition the scattering region into conveniently chosen small domains [21, 22]. The corresponding Green’s functions are recursively combined using the Dyson equation to obtain matrix elements of the full system Green’s function that are relevant for transport calculations. The RGF method is robust, accurate, has a simple implementation, and has been widely used [23, 24, 25, 26, 27, 28, 29].

A recent open source implementation of the WFM method, the Kwant package [5], has significantly increased its usage. Kwant is developed under a user-friendly platform coded in Python and handles general-shaped scattering regions, multiple orbitals, and multi-probes [30]. Furthermore, extensions to the Kwant package can be easily joined [31]. Kwant also explores the sparsity of the system Hamiltonian by using the MUMPS libraries, a forefront package for sparse linear algebra [32]. Ref. [5] shows that Kwant significantly outperforms the RGF method in a wide range of applications.

In this paper, we show that the number of operations required by the WFM method to compute the conductance of a given system is much smaller than previously claimed [5]. To explain this finding, we first give a self-contained presentation of the method – whose documentation is scattered and scarce – critically analyzing its main features. Next, we numerically study a number of systems to corroborate our analytical findings.

The paper is organized as follows: In Sec. 2 we provide a short review of the relation between quantum transport and the scattering theory. Next, we adapt the theory for the tight-binding approximation and cast the scattering problem as the solution of a linear system. In Sec. 3 we describe the WFM method and discuss its computational cost. In Sec. 4 we benchmark the WFM method comparing its CPU time, memory usage and precision with a standard RGF implementation. Next, we present an application of the WFM method for a realistic-sized multi-probe graphene Hall bar system. In Sec. 5 we summarize our conclusions.

2 Theoretical background

The WFM method is suited to calculate the scattering properties of a system with arbitrary geometry and dimension $d=1,2,$ or 3. It is aimed to describe the non-interacting on-shell scattering processes in mesoscopic samples or crystalline structures coupled to multiple terminals.

Figure 1 illustrates a generic multi-terminal two-dimensional $(d=2)$ system. A central scattering region is coupled to electrodes represented by semi-infinite leads labeled by $\alpha=1,\cdots,\Lambda$, where incoming and outgoing electrons propagate coherently. Due to the transverse confinement, the leads states are quantized in open modes (scattering channels) labeled by $n=1,\dots,N_{\alpha}$. The index $n$ labels both the transverse modes and the electron spin projection. The mesoscopic sample corresponds to the central or scattering region, while the leads are associated to the asymptotic domain.

Let us describe the system single-particle Hamiltonian by a tight-binding model. This approximation is suited to model both an atomistic system represented by a linear combination of atomic orbitals and a continuous system in a finite element representation [8]. The system Hamiltonian is written as

[TABLE]

where the index $j=({\bm{r}}_{i},\sigma)$ labels both the position in the lattice and the internal degrees of freedom $\sigma$ such as spin, atomic orbital, etc., of the state $\left|{j}\right\rangle$.

2.1 Quantum transport and scattering theory

The Schrödinger equation of the scattering system reads

[TABLE]

where $|\Psi_{m}^{+}(E)\rangle$ ($|\Psi_{m}^{-}(E)\rangle$) stands for the outgoing (incoming) scattering state at channel $m$. Here $m$ labels both $\alpha$ and $n$. The S-matrix is defined by the scattering amplitudes $\langle\Psi_{m}^{-}(E)|\Psi_{m^{\prime}}^{+}(E)\rangle=S_{mm^{\prime}}(E)\delta(E-E^{\prime})$.

The scattering matrix $S$ can be formally written in terms of projection operators that decompose the Hilbert space in the partition described by Fig. 1 [33]. Let us assume, for instance, normal boundary conditions at the interface ${\cal B}$ between the scattering and asymptotic regions. One defines the projection operator

[TABLE]

in terms of the complete set of discrete orthonormal states $\left\langle{\phi_{\mu}}\middle|{\phi_{\mu^{\prime}}}\right\rangle=\delta_{\mu\mu^{\prime}}$ defined in the scattering (or central) region and obeying the boundary conditions at ${\cal B}$. In turn, at the asymptotic region, one defines

[TABLE]

where $\left|{\chi_{m}(E)}\right\rangle$ form a complete set of continuous orthogonal states, $\left\langle{\chi_{m}(E)}\middle|{\chi_{m^{\prime}}(E^{\prime})}\right\rangle=\delta_{mm^{\prime}}\delta(E-E^{\prime})$, defined in the asymptotic (or leads) region. Since the asymptotic region is not compact, the projection operator $P$ is continuous. By construction, $P$ and $Q$ span the system Hilbert space and, hence, $P+Q=1$.

The system Hamiltonian is conveniently decomposed into three pieces

[TABLE]

where we introduced the notation $AHB=H_{AB}$ 111It is convenient to use as the channel basis in the asymptotic region the eigenfunctions of $H_{PP}$, namely, $H_{PP}\left|{\chi_{m}(E)}\right\rangle=E\left|{\chi_{m}(E)}\right\rangle$ with normal boundary conditions at ${\cal B}$..

The projection operators allow one to write Eq. (2) as a Lippmann-Schwinger equation, namely

[TABLE]

where $E^{\pm}=E\pm i\eta$, with $\eta$ an infinitesimal positive number. After some algebra, one writes the $S$-matrix as [33, 9]

[TABLE]

where $\Sigma^{\pm}(E)=H_{QP}(E^{\pm}-H_{PP})^{-1}H_{PQ}$ is the embedding self-energy, which accounts for coupling the to the continuum and describes the resonance processes, while $\rho_{m}(E)$ stands for the electronic density of states at the channel $m$. Here we explicitly neglect direct tunneling processes between different electrodes [33, 9]. This approximation is accurate provided the central region is sufficiently large to prevent direct tunneling processes across the system. This condition is true for most mesoscopic systems, except for small molecular junctions (for more details, see, for instance, Ref. [34]).

The Landauer-Büttiker theory [8] relates the linear conductance of an electronic sample to the transmission probability as

[TABLE]

where $f(E)=\left[1+e^{(E-\mu)/kT}\right]^{-1}$ is the Fermi-Dirac distribution with $\mu$ and $T$ giving the equilibrium chemical potential and temperature of the reservoirs 222In general, the reservoirs have different chemical potentials and temperatures, thus, $f_{\alpha}(E)=[1+e^{(E-\mu_{\alpha})/kT_{\alpha}}]^{-1}$. For simplicity we take all temperatures equal to $T$ and since we restrict ourselves to linear response, the small differences between $\mu_{\alpha}$ and the equilibrium $\mu$ lead to Eq. (9). . The transmission $T_{\alpha\beta}(E)$ is given by

[TABLE]

where $S_{nm}(E)$ is given by Eq. (8). The WFM method also gives local properties such as local currents and the local density of states (LDOS), as discussed in Sec. 4.

2.2 The scattering problem in tight-binding approximation

Let us now write the system Hamiltonian in a suitable form to implement the WFM method. For the sake of simplicity, we discuss in detail the two-terminal case and, at the end, we generalize the results to the multi-terminal case.

Let us consider a mesoscopic system attached to semi-infinite leads, $\alpha=R,L$, as illustrated by Fig. 2a. Following the partition operators presented in the previous section, we introduce the standard matrix representation: (i) $H_{QQ}\leftrightarrow H_{S}$ for the scattering region Hamiltonian; (ii) $H_{PP}\leftrightarrow H_{L}+H_{R}$, for the leads Hamiltonian; (iii) $H_{QP}\leftrightarrow V_{SL}+V_{SR}$, for the coupling term connecting the mesoscopic system to the leads.

The full Hamiltonian is written in a block matrix form as

[TABLE]

$H_{L}$ and $H_{R}$ can be written in the block-diagonal structure

[TABLE]

where $H_{l(r)}$ stands for suitable $L(R)$-lead unit cell Hamiltonian of dimension $M_{L(R)}$ (represented by boxes in Fig. 2). $V_{l(r)}$ are the hopping matrices between nearest-neighboring unit-cells and the unwritten matrix elements are identically zero.

It is advantageous to use the structure of the leads matrices $H_{L}$ and $H_{R}$ to group them into an effective single-lead with disjoint sections. The rearranged layout is depicted in Fig. 2b. The modified $H$ reads

[TABLE]

The effective lead, hereafter denoted by $T$, compacts the eigenvalue problem to a single semi-infinite partition, namely

[TABLE]

where $\psi_{S}$ corresponds to the scattering wave function at the central region and $\psi_{n}$, to the lead wave function at the $n$-th slice, with $n=0,0,2,\dots$ (see Fig. 2). The generalization to a multi-terminal setup is straightforward. In this case, $H_{T}$ accounts for all $H_{\alpha}$, with $\alpha=1,\dots,\Lambda$ and has dimension $M_{T}=\sum_{\alpha=1}^{\Lambda}M_{\alpha}$.

3 The wave function matching method

Let us now solve Eq. (14). For that purpose we introduce the eigenmode basis $\phi_{n}$:

[TABLE]

which corresponds to rows of Eq. (14) far from the scattering region. Due to translational symmetry, one can use Bloch’s theorem to conveniently write $\phi_{n}$ as

[TABLE]

where $\chi$ is the lead unit cell eigenfunction (independent of $n$) and $\lambda$ is a complex constant. Hence,

[TABLE]

The standard procedure to solve this quadratic eigenvalue problem (QEP) in $\lambda$ is to introduce an auxiliary vector

[TABLE]

and to linearize Eq. (17)[35] as

[TABLE]

The advantage of casting Eq. (17) as a Generalized eigenvalue problem (GEP) is that one can calculate the eigenvalue $\lambda$, which is associated to the crystal momentum $k$ (using $\lambda=e^{ika_{\alpha}}$, where $a_{\alpha}$ is the $\alpha$-lead lattice constant), and the eigenvector $\chi$ as a function of the electronic energy $E$. The QEP is translated into a linear problem at the expense of doubling the equation dimension. Hence, the number of eigenvalues is twice the rank $M_{T}$ of the matrices $V_{T}$ and $H_{T}$.

One can solve the GEP in Eq. (19) by means of well-known numerical algorithms [35, 36, 37, 38]. Given an electronic energy $E$, we calculate the eigenvectors $(\chi_{p},{\chi^{\prime}}_{p})$ and the corresponding eigenvalues $\lambda_{p}$, where $p=1,\cdots,2M_{T}$.

We can infer from the scattering problem that the $2M_{T}$ solutions correspond to $M_{T}$ incoming modes and $M_{T}$ outgoing modes, as depicted in Fig. 3. Since the terminals are uncoupled, the eigenstate $\chi_{p}$ has a block structure

[TABLE]

where each block $\chi_{p}^{\alpha}$ describes the eigenstate of the $\alpha$-lead with eigenvalue $\lambda_{p}$ for $p=1,\cdots,M_{\alpha}$ and $\alpha=1,\cdots,\Lambda$.

The modes can be propagating $|\lambda_{p}|=1$ or evanescent $|\lambda_{p}|<1$ ($|\lambda_{p}|>1$ gives a non-physical behavior). The probability current for the $p$-th propagating mode reads [39]

[TABLE]

The incoming modes correspond to $j_{p}>0$ and the outgoing ones to $j_{p}<0$. We label those two sets of solutions as $\lambda_{p}^{\pm},\chi_{p}^{\pm}$ for $p=1,\dots,N_{P}$, where $\pm$ indicates the corresponding current direction and $N_{P}\leq M_{T}$ is the number of incoming/outgoing propagating channels at the electronic energy $E$. See Fig. 3. Since we are interested in transmission coefficients, we restrict ourselves to the analysis of the propagating modes. The evanescent modes ($j_{p}=0$) can be treated straightforwardly as a generalization of this method.

Using the sets $\chi_{p}^{\pm}$ as basis, we write the wave functions $\psi_{n}$ as

[TABLE]

where $n=0,1,\cdots$, and $A_{p}^{\pm}$ are unknown amplitudes.

One defines the scattering matrix $\tilde{S}$ that relates incoming with outgoing amplitudes as

[TABLE]

Since the eigenchannel basis used by the WFM method is not normalized as the one introduced in Sec. 2, the matrix $\tilde{S}$ does not preserve the flux. As we discuss below, $S$ is obtained from $\tilde{S}$ by a simple relation.

To calculate the S-matrix, we consider the scattering process of a single incoming mode $q$, namely

[TABLE]

The corresponding S-matrix can be obtained by solving the first two lines of Eq. (14)

[TABLE]

where $\psi_{Sq}$ is the scattering region wave function upon injection from mode $q$. Substituting Eq. (24) into Eq. (26) and recalling that the basis functions $\chi^{\pm}_{p}$ satisfy Eq. (17), we find

[TABLE]

where $\psi_{-1q}$ is also given by Eq. (24). Note, however, that $\psi_{-1q}$ has no physical meaning, since in Eq. (14) there is no slice defined for $n=-1$. Here, $\psi_{-1q}$ is an auxiliary mathematical quantity designed to represent the contributions of the terms including $\psi_{0q}$ and $\psi_{1q}$ in Eq. (26).

Applying the definition of $\chi^{\prime}$, Eq. (18), to each propagating mode as

[TABLE]

Eq. (27) becomes

[TABLE]

where ${\chi^{\prime}}^{\pm}\equiv({\chi^{\prime}}^{\pm}_{1},{\chi^{\prime}}^{\pm}_{2},\cdots,{\chi^{\prime}}^{\pm}_{N_{P}})$ with dimension $M_{T}\times N_{P}$ and $\tilde{S}_{q}$ is the column $q$ of the S-matrix with dimension $N_{P}\times 1$. Analogously, using Eq. (24) we write $\psi_{0q}$ as

[TABLE]

The linear system composed by Eqs. (25) and (29) reads

[TABLE]

Let us now generalize Eq. (31) to account for different $q$-modes

[TABLE]

We cast this result into the compact form

[TABLE]

where $\tilde{S}$ is the full S-matrix and $\psi_{S}=\begin{pmatrix}\psi_{S1}&\psi_{S2}&\cdots&\psi_{SN_{P}}\end{pmatrix}$ is the wave function of the scattering region. The S-matrix has dimension $N_{P}\times N_{P}$ while $\psi_{S}$ has dimension $N_{S}\times N_{P}$, since it is defined for all the $N_{S}$ sites in the central region upon injection from all the $N_{P}$ channels.

Hence, the solution of Eq. (33) has a computational cost that depends on the number of propagating channels $N_{P}$ at the electronic energy $E$. Due to the sparsity of $H_{S}$, we infer that CPU time required to compute a given system conductance scales as $N_{S}\times N_{P}$. In Sec. 4 we numerically verify that the WMF method indeed follows this prediction.

Note that Eq. (33) involves representations in different spaces, while the scattering wave function is given in the tight-binding basis, the S-matrix is expressed in eigenmode basis. The matrices $\chi^{\pm}$ give a connection between theses two basis [40]. For a sufficiently large system, $H_{S}$ and $V_{TS}$ are sparse matrices making the problem appropriate to the sparse solvers.

3.1 Connection to Green’s functions

The coupling with leads gives a finite line-width to the resonances in the scattering region via a so-called self-energy. In the non-equilibrium Green’s functions formalism (NEGF) (see, for instance, Refs. [8, 34]) the embedding self-energy modifies the scattering region Hamiltonian as $H_{S}\rightarrow H_{S}+\Sigma$. In what follows we demonstrate that $\Sigma$ can be calculated from the presented equations.

Let us define the dual space states $\tilde{\chi}_{p}^{\pm}$, where

[TABLE]

and identify the first and the second terms on the RHS of Eq. (24) with

[TABLE]

Introducing the translation operator $F_{\pm}$ [7]

[TABLE]

where

[TABLE]

one can write $\psi_{0q}$ and $\psi_{1q}$, respectively, as

[TABLE]

Substituting Eq. (39) into Eq. (26) and solving Eq. (25) for $\psi_{Sq}$ we find

[TABLE]

where

[TABLE]

is a source term dependent of which channel $q$ is injecting,

[TABLE]

is the embedding self-energy, and

[TABLE]

is the surface Green’s function of the semi-infinite leads. Since Eq. (43) involves outgoing states $F_{+}$, $G_{T}$ and $\Sigma$ correspond to retarded Green’s function and self-energy, respectively [8].

We stress that both $\Sigma$ and $G_{T}$ are independent of $q$ and contain information about all the propagating modes at the energy $E$. In A we show that Eq. (43) reproduces the exact analytical surface Green’s function of a 1D chain.

Notice that we can solve Eq. (40) for $\psi_{Sq}$ as

[TABLE]

where $G_{S}$ is the scattering region Green’s function given by

[TABLE]

Thus, knowing the full Green’s function matrix $G_{S}$, we can calculate $\psi_{Sq}$ for any $q$ using Eq. (44).

With the help of the dual vector $\tilde{\chi}_{p}^{+}$ defined in Eq. (34) and the definition of $\psi_{0q+}$ given by Eq. (35), we calculate the amplitudes $\tilde{S}_{pq}$ as

[TABLE]

The outgoing wave function $\psi_{0q+}$ is a superposition of states $\chi_{p}^{+}$ with amplitudes $\tilde{S}_{pq}$. Those states carry a probability current

[TABLE]

Here $j_{pq}$ depends on the injecting mode $q$ and $j_{p}$ is given by Eq. (21).

The transport coefficients $P_{pq}$ defined as the ratio between the incoming probability current $j_{q}$ and the outgoing probability current $j_{pq}$ at mode $p$ reads

[TABLE]

where we defined the scattering amplitudes $S_{pq}$ as [8]

[TABLE]

where $S$ is unitary and conserves the current probability [8].

3.2 Generalized Fisher-Lee expression for the transmission amplitudes

Let us use the WFM method elements introduced above to derive the relation between the transmission amplitudes as functions of the scattering region Green’s functions.

First we write $\psi_{0q+}$ in the RHS of Eq. (38) as a function of the scattering region wave function $\psi_{Sq}$ using Eqs. (39) and Eq. (26), namely

[TABLE]

Hence, the scattering amplitude $S_{pq}=\sqrt{{j_{p}}/{j_{q}}}\tilde{S}_{pq}$ reads

[TABLE]

Here we used Eqs. (41) and (44) to substitute the dependence on $\psi_{Sq}$ by a dependence on the scattering region Green’s function $G_{S}$.

We assume that the modes $q$ and $p$ belong to different leads $\alpha$ and $\beta$, respectively. Due to the block structure of Eq. (20) and to the absence of coupling between the leads, the matrices $G_{T}$, $V_{T}^{\dagger}$ and $\left(F_{-}-F_{+}\right)$ are block diagonal in the leads subspace. The two-contacts Hamiltonian in Eq. (13) illustrates the diagonal block structure of $V_{T}^{\dagger}$, for instance. Thus, the second term in Eq. (51) identically vanishes.

In this case, the scattering amplitude in Eq. (51) becomes

[TABLE]

where $t_{pq}^{\beta\alpha}$ is the current-normalized transmission amplitude for the scattering from mode $q$ in the lead $\alpha$ to the mode $p$ in the lead $\beta$.

Although one can calculate the transmission coefficients by means of $G_{S}$ from Eq. (52), only few Green’s functions matrix elements, such as the elements connecting sites belonging to the interface with the leads, are required to compute the transmission (see, for instance Ref. [20]). Therefore, a simplification of Eq. (52) is desirable. For that purpose we use a sub-block division of the scattering region similar to the one used in Ref. [20].

We divide the scattering region into $\Lambda+1$ blocks, where $C$ is the a central block, which has no connection with the leads, and $\alpha$ represents the $\alpha$-interface, which is connected to $C$ and only to the lead $\alpha$, where $\alpha=1,\cdots,\Lambda$. In this picture, $G_{S}$ and $V_{TS}$ read

[TABLE]

Since $G_{T}V_{T}^{\dagger}\left(F_{-}-F_{+}\right)$ is diagonal, where $[G_{T}]_{\alpha\beta}=\delta_{\alpha\beta}G_{\alpha}V_{\alpha}^{\dagger}\left(F_{-}^{\alpha}-F_{+}^{\alpha}\right)$ and $F_{\pm}^{\alpha}\equiv\sum_{p\in\alpha}^{N_{P}}\lambda_{p}^{\pm}\chi_{p}^{\pm\alpha}\left(\tilde{\chi}_{p}^{\pm\alpha}\right)^{\dagger}$, and the states $\chi_{q}^{-}$ and $\tilde{\chi}_{p}^{+}$ have different non-vanishing blocks given by Eq. (20), we find

[TABLE]

which is the generalized Fisher-Lee expression [8].

4 Benchmark and Application

Let us now demonstrate the efficiency of the sparse solvers associated with the WFM method implemented in the Kwant package. To this end, we compare the processing time and memory usage of the WFM method with the standard RGF approach for a two-dimensional model system as a function of its size and aspect ratio. We conclude this section by discussing an application of the WFM method, namely, the calculation of longitudinal and transverse resistance of a realistic-sized graphene Hall bar.

As mentioned in the introduction, nowadays the RGF method is one of the most standard technique to compute the conductance of nanoscale systems. This method is designed to compute only the system full Green’s function matrix elements related to transport properties [8]. For that purpose, the system is divided into partitions. The computational time necessary to calculate the transmission scales with the number of partitions times the cube of the typical number of sites within the partitions.

We recall that Ref. [5] draws conclusions by comparing the performance of the RGF and WFM methods for a square lattice system with $L\times L$ sites as a function of $L$. The authors [5] find that the CPU time required to compute the conductance using the RGF method scales with $L^{4}$, while the WFM implementation in Kwant scales with $L^{3}$. Here we explore more diverse situations to numerically verify that the WFM method is more efficient than $L^{3}$, as discussed in Sec. 3.

Let us begin considering a nearest neighbor (nn) tight-binding Hamiltonian in a two-dimension square-lattice of length $L$ and width $W$ in number of sites. We take $W^{\prime}$ as the width of the leads (see inset of Fig. 5a). We set $E=0$. In this case, we recall that for semi-infinite square lattice leads the number of open channels at the left and right leads $N_{L}=N_{R}=W^{\prime}$. This model stems for instance from a finite-difference discretization of the Schrödinger equation of a mesoscopic two-dimensional electron gas (2DEG) [40, 8]. For this model, the optimal partition of the RGF consists of $L$ partitions (slices) with $W$ sites each.

Figure 4 gives the CPU time (in arbitrary units) necessary to compute the conductance of the system, Eq. (9), as a function of $L$, $W$, and $W^{\prime}$. It should be emphasized that in both implementations, the linear algebra calculations are coded in lower level programming languages, making this comparison possible.

As discussed in Sec. 3 one has to solve $N_{P}$ times the sparse linear system of dimension $N_{S}+M_{T}$, Eq.(32). Since the number of operations to solve a sparse system scales as $\mathcal{O}(N)$ [35] and here $M_{T}=2W^{\prime}$, the WFM is expected scale as $(LW+2W^{\prime})W^{\prime}$. Figure 4a to 4c verify that this conjecture is indeed correct. As a consequence, the performance of the WFM method is much better than previously believed [5] for a realistic model of a nanostructure, $M_{T}=2W^{\prime}\ll W$.

Let us now examine a situation where $W=W^{\prime}$. Figure 5 clearly shows that the CPU time of the RGF (nn) method scales with $L^{4}$, as expected by the matrix multiplication and diagonalization operations involved. In distinction, the WFM shows a much better CPU performance scaling as $L^{3}$ (here $L=W=W^{\prime}$). However, the overall pre-factor is typically large, making the method clearly advantageous only for $W\gtrsim 10^{2}$ sites.

We use this setting to investigate the efficiency of the WFM method when dealing with tight-binding Hamiltonians that consider hopping matrix elements beyond nearest-neighbor sites. This is the case in tight-binding models based on Wannier wave functions [41], that are very practical and accurate tools to model large scale disorder systems. Let us consider a square lattice tight-binding Hamiltonian with up to the $3$-rd nearest-neighbor (3nn) hopping terms. Since for the RGF method, only neighboring partitions should be connected in this model, one has to double the size of each slice, $W\to 2W$, reducing the total number of slices by half $L\to L/2$. Hence the CPU time grows by a factor of $4$ (solid line of Fig. 5). In Fig. 5 we show that Kwant is practically insensitive to the coordination number of the lattice model, which represents a huge advantage over RGF.

Let us now analyze the memory usage of both methods. As already pointed out in Ref. [5], the memory usage in Kwant can be ten times larger than an RGF implementation which is a problem for computation of transport properties in large systems. In what follows we study this issue in more detail, examining the intermediate processes, such as the leads eigenmodes calculation, the linear system construction and factorization, and the solving stage, regions (ii)-(iv) of Fig. 6, respectively. This stage-by-stage information of the memory usage gives a clear view of the method advantages and bottlenecks.

Figure 6 shows the memory usage in a conductance calculation for both the WFM and the RGF implementations. A huge difference can be noted between the maximum memory used for each method. In Kwant, a preliminary time is spent in reading the input parameters, stage (i), which is negligible in the RGF Fortran 90 implementation and it is not displayed. The next stage in both methods, indicated by (ii) in Fig. 6, is related to the computation of the lead contribution, namely, the lead surface Green’s function in the RGF [8] and the eigenmode diagonalization in WFM. In both methods, this is done twice for our two-probe model and $\Lambda-$times in general systems. Kwant spends an extra time in the factorization of the linear system, Eq. (33).

At the solving stage, indicated by (iv) in Fig. 6, we observe that Kwant requires one order of magnitude more memory than the RGF method. This is the only feature where the RGF outperforms the WFM methods. We note however that WFM approach allows for the computation of local operators (such as local currents and LDOS) with no significant additional cost, which is not the case for the RGF method.

Both methods are very robust and accurate. In our extensive tests, the computed conductances agree within the numerical precision. Even in the cases where the Green’s function regularization factor $\eta$ is known to require a special choice in RGF, like transmission by evanescent modes in graphene [42], the WFM method gives reliable results without any particular adjustment.

4.1 Application: Graphene Hall Bar

Let us now show the results of the WFM method for the transport properties in a realistic size graphene Hall bar in the quantum Hall regime [43]. Despite the importance of such class of systems, few numerical studies have addressed the longitudinal and transverse resistances in Hall bar geometries due to the lack of an efficient multi-terminal electronic transport code.

We consider a graphene sample with $\sim 10^{6}$ atoms in a Hall bar geometry (inset of Fig. 8). The graphene tight-binding Hamiltonian [43] is $H=-\sum_{\langle i,j\rangle}(t_{ij}\left|{i}\right\rangle\left\langle{j}\right|+{\rm H.c.})+\sum_{i}\epsilon_{i}\left|{i}\right\rangle\left\langle{i}\right|$, where the sums run over the sites of a honeycomb lattice and $\langle\cdots\rangle$ restricts the pairs of sits to nearest-neighbors. The model includes a local (Anderson) scalar disorder by randomly choosing $\epsilon_{i}$ from an uniform distribution $[-\delta W,\delta W]$, where $\delta W=0.08t$.

The magnetic field ${\bf B}=B{\bf e}_{z}$ perpendicular to the graphene sheet is accounted for by Peierls substitution, namely, by taking $\phi_{ij}=\frac{e}{\hbar}\int_{{\bf r}_{i}}^{{\bf r}_{j}}{\bf A}({\bf r})\cdot d{\bf r}$ and ${\bf A}=B\left(-\beta y,(1-\beta)x,0\right)$. The gauge $0\leq\beta\leq 1$ is conveniently chosen according to the orientation of the leads. Since in our Hall bar we consider leads along both the $x$ and $y$ directions, we avoid discontinuities in the magnetic field by smoothly varying the vector potential according to Ref. [44, 29].

We calculate the system longitudinal and transverse resistances using the Landauer-Büttiker formula, Eq. (9), considering the terminals $\alpha=2,3,5,6$ as voltage probes, that is, $I_{2}=I_{3}=I_{5}=I_{6}=0$ (see inset of Fig. 8b). In this setting, $I_{1}=-I_{4}=I$. Hence, $R_{xx}=(V_{2}-V_{3})/I_{1}$ and $R_{xy}=(V_{3}-V_{5})/I_{1}$. Our results for a single disorder realization correspond to typical quantum Hall resistance curves for graphene samples [43]. Figure 8 shows quantized Hall plateaus at $R_{xy}=\frac{h}{2e^{2}}\frac{1}{2n+1}$ for integer values of $n$ and zero longitudinal resistance $R_{xx}$ at the $R_{xy}$ plateaus.

We also calculate the local density of states [5]. Figures 7a and 7b show localized edge states along the sample in the plateau region and a delocalized state in the transition region between two plateaus (solid and dashed vertical lines of Fig. 8, respectively). As expected, the local current [5] show quantum Hall edge states, see Fig. 7c.

The CPU time required for the resistance calculations is of the order of $80$ seconds for a single disorder realization and a single energy value in one core of an Intel® Xeon® X5650 processor. As discussed, such fast computation time in WFM relies on the $W^{\prime}\ll W$ condition.

5 Conclusion

We have reviewed the underlying theory of the (WFM) method applied to a tight-binding (finite element) Hamiltonian used to model the transport properties of mesoscopic systems. Our analysis revealed that the WFM method is computationally far superior than previously expected [5].

We numerically verify ou predictions in a number of settings, benchmarking the CPU time, memory usage and precision of the WFM versus the RGF method.

To illustrate the power of the method we calculate the longitudinal and transverse resistance of a realistic-sized disordered graphene sheet in the quantum Hall regime. We consider a sample patterned in a Hall bar geometry, corresponding to a multi-terminal setting difficult to treat with other numerical approaches.

We conclude mentioning that the WFM method allows for a straightforward generalization for multi-terminal systems with nontrivial sample geometries, while the RGF approach resorts on ingenious schemes to deal with such situations [25, 20]. In addition, the Kwant package also offers a set of implementation tools to facilitate the study of a wide range of settings, such as multi-orbital atomic states, general lattice connectivity and geometry, to name a few.

Acknowledgments

We thank Xavier Waintal and Bruno A. D. Marques for useful discussions. We acknowledge the financial support of the Brazilian funding agencies CNPq, CAPES, and FAPERJ.

Appendix A Surface Green’s function for a 1D chain

In this Section we explicitly derive the Green’s function of the first site in a semi-infinite linear chain of atoms using Eq. (43). We consider a system with one orbital per atom, where the hopping matrices are numbers given by $H_{T}=\epsilon_{0}$ and $V_{T}=-t$. We apply the definition of the leads eigenstates in Eq. (16) into Eq. (15) to find

[TABLE]

whose solutions are

[TABLE]

Notice that if we substitute $\lambda=e^{ika}$ into Eq. (64) we recover the well known dispersion relation for the 1D chain.

The eigenstates corresponding to $\lambda_{1}^{\pm}$ are $\chi_{1}^{\pm}=1$. The probability current $j_{1}$, given by Eq. (21), reads

[TABLE]

Thus, there is one state propagating forwards ($j_{1}^{+}>0$) and one propagating backwards ($j_{1}^{-}<0$), resulting in $N_{p}=1$.

Since the dual vectors are simply $\tilde{\chi}_{1}^{\pm}=1$, Eq. (36) leads to $F_{-}=\lambda_{1}^{+}$. Therefore, the surface Green’s function $G_{T}$ in Eq. (43) reads

[TABLE]