$W^{3}L$ as can be seen from the equations of the gluing sequence (one
$W\times W$ matrix inversion per slice added). It is therefore very efficient
for systems that are close to being one-dimensional $L\gg W$.

The original RGF algorithm has been generalized to different strategies for a
parallel implementation [78] and adding slices of different shapes in an
optimized way. This can be done using algorithms for analyzing the
connectivity graph of the Hamiltonian matrix [307, 189, 197, 167], using
minimum slices with just a single site [130] or circular slicing [277]. These
more recent algorithms have the advantage of being easy to deploy on multi-
terminal systems, arbitrary geometries and lattices but have a scaling similar
to RGF. Other approaches take advantage of a structure of the problem, such as
some parts of the device being defect free, to speed up the calculation [276].
[193] extends RGF to calculate other quantities than the transmission
probability such as the local density of states.

Figure 8: Schematic of one step of the original recursive Green’s Function
(RGF) algorithm. the gluing sequence is used iteratively to add new sites
slice by slice. The site coloring is the same as in Fig. 6

For systems whose different dimensions are of similar sizes ($L\approx W$), an
efficient way to build the system is to use the nested dissection algorithm
[147, 165]. In nested dissection one recursively double the system size by
using the gluing sequence according to the schematic of Fig. 9. For a $L\times
L$ square system, the computing time reduces to $L^{3}$ which is
parametrically faster than the $L^{4}$ of RGF. However, implementing the
nested dissection algorithm in a stable way is not easy as it implies gluing
blocks that are not connected to the electrodes (hence finite with badly
conditioned Green’s functions, see the discussion above). It can thus be
beneficial to use one of the existing direct sparse solvers as discussed in
Sec. VI.2.

There have been many algorithms related to RGF that all rely on some form of
the gluing sequence. Those include the modular Green’s function algorithm
[241, 240], the (parallel) patchwork algorithm [63] or matching methods [31]
that glue together various precalculated subparts of the system. Other
applications include multiscale modeling [52] where a central part of the
system is described accurately within DFT while the periphery is modeled by a
parametrized tight-binding model.

A very closely related class of algorithm uses “decimation” to integrate out
sites that are not of direct interest. The sites integrated out appear through
self-energy terms. Despite being apparently different, these algorithms are
very close to RGF kind of algorithms. Early algorithms include [103, 217].

Figure 9: Schematic of one step of the nested dissection algorithm. The gluing
sequence is used to recursively double the size of a $L\times L$ square system
to a $2L\times 2L$ system. The site coloring is the same as in Fig. 6

###  VI.5 Other approaches beyond this review, discussion.

In the early days of RGF, another competing algorithm was written in terms of
the transfer matrix that relates the wave-function in one layer to the next
[285]. The transfer matrix algorithm requires to be stabilized against
accumulation of errors as it does not naturally preserve unitarity [286]. It
has been mostly abandoned in favor of RGF.

The RGF class of algorithms have a computational effort that scales as
$L^{3d-2}$. These algorithms are very efficient in one dimension of for very
elongated systems. For systems where all dimensions have similar sizes,
approached based on nested dissection are parametrically faster and in
practice outperform RGF for $L\geq 100$ [104]. In two dimensions a calculation
of a squared system of $10^{6}$ sites takes less than 1 hour on a single core,
making 2D calculations easily tractable.

In three dimensions, the situation is more complex since it is difficult to
calculate systems large enough for the calculations not to be dominated by
finite size effects. An alternating idea to RGF kind of calculations is to use
a different starting point: instead of starting from vacuum and adding sites
layers by layers, one may start from a pristine material and use Dyson
equation to modify a finite number of Hamiltonian matrix elements [208, 256].
This approach allows one to start directly in the thermodynamic limit hence
address difficult 2D or 3D systems. In some simple cases, the pristine Green’s
function may be obtain analytically through complex contour integration [208,
253, 256]. The method has been generalized by [119] to arbitrary pristine
systems by a numerical calculation of the poles and residues of the problem.

When the current is forced to go through a small constriction with very few
propagating channel, one may take advantage of this fact to build an effective
quasi-1D problem that is much faster to solve [69, 70] even if the original
electrodes are fully 3D.

For quantities that are mostly local such as the electronic density, local
density of states or the conductivity (as opposed to the full conductance of a
device that is essentially global), linear scaling methods such as the Kernel
Polynomial Method (KPM) can be very efficient [304]. Linear scaling methods
for transport have been recently reviewed in [83].

This review has focused on discrete models. In some specific cases, it is also
possible to perform wave-matching directly in the continuum [274, 286, 64] or
by treating part of the system (e.g., the leads) directly in the continuum
[159].

Last, let us mention two alternative techniques that are slightly different
from the mainstream approach. The Complex Absorbing Potential (CAP) technique
can be used to bypass solving the lead problem by including a finite part of
the lead in the calculation in presence of a smoothly varying complex
potential [77]. In the Contact Block Reduction (CBR) approach, one
approximates the Green’s function of the isolated scattering region part by
selecting only a finite fraction of the eigenstates of the closed system
[183].

##  VII Physical observables in the Landauer-Büttiker approach

So far, we have been concerned solely with one-body quantum mechanics. We now
turn to expressing the out-of-equilibrium observables, such as the
conductance, in terms of the solution of the quantum mechanical problem, i.e.,
we now incorporate statistical physics.

Introducing the destruction $c_{n}$ (creation $c^{\dagger}_{n}$) operators on
site $n$, we are now considering Hamiltonians in second-quantized form,

${\cal H}=\sum_{nm}\hat{H}_{nm}c^{\dagger}_{n}c_{m}.$ (134)

where the matrices $\hat{H}$ correspond to the various infinite matrices of
the one-body problem.

###  VII.1 Many-body steady states

The solutions of the scattering problem $\hat{\psi}(E)$ computed in Sec. III.3
together with the bound states from Sec. III.5 form a complete basis of
single-particle wave functions of the Hamiltonian (46).

Because we focus on time-independent quantum transport, we limit our
consideration to steady many-body states, where each of the single particle
states has a specific occupation number $f_{\alpha E}$. Requiring additionally
that the fermion occupation numbers are uncorrelated, we obtain the general
form of the density matrix of such states

$\displaystyle\rho$ $\displaystyle=\prod_{\alpha
E}\left([1-f_{\alpha}(E)](1-\Psi_{\alpha E}^{\dagger}\Psi_{\alpha
E})+f_{\alpha}(E)\Psi_{\alpha E}^{\dagger}\Psi_{\alpha E}\right),$ (135)
$\displaystyle\Psi_{\alpha E}$ $\displaystyle=\sum_{n}\hat{\psi}_{\alpha
E,n}c_{n}.$ (136)

The single particle density matrix for such a state is

$\langle
c^{\dagger}_{n}c_{m}\rangle=\sum_{\alpha}\int\frac{dE}{2\pi}\hat{\psi}_{\alpha
E,n}^{*}\hat{\psi}_{\alpha E,m}f_{\alpha}(E)$ (137)

where we remember that the density of state of the one dimensional channels
has been already incorporated (up to the factor $2\pi$) in the normalization
of the wave function. Indeed, we normalize the states to carry unit current
and in one dimension the velocity $(1/\hbar)dE/dk$ is the inverse of the
density of state. When the sites $n,m$ belong to the scattering region, we
arrive at

$\langle
c^{\dagger}_{n}c_{m}\rangle=\sum_{\alpha}\int\frac{dE}{2\pi}(\Psi_{\text{sr}}^{*})_{n\alpha}(\Psi_{\text{sr}})_{m\alpha}f_{\alpha}(E)$
(138)

which is expressed in terms of quantities that we learned to calculate. We
narrow our scope further and consider one particularly important class of
steady states where each of the leads has fixed chemical potential voltage
$\mu_{a}$ and temperature $\Theta_{a}$. Then the modes $\hat{\psi}_{\alpha E}$
incoming from the lead $a$ have an occupation number given by the Fermi
distribution of that lead

$f_{a}(E)=\frac{1}{e^{\beta_{a}(E-\mu_{a})}+1},$ (139)

with $\beta=1/k_{b}\Theta$. It is worth noting that, although such a steady
state is not a thermal state (the different leads can have different
temperatures and/or chemical potentials), Eq. (135) is the thermal state of
the effective nonlocal Hamiltonian $H_{\mathrm{eff}}$,

$\rho=\frac{1}{Z}e^{-H_{\mathrm{eff}}},\quad H_{\mathrm{eff}}=\sum_{\alpha
E}\beta_{a}(E-\mu_{a})\Psi^{\dagger}_{\alpha E}\Psi_{\alpha E}.$ (140)

This form makes it explicit that the Wick theorem is applicable to this
density matrix. ($Z$ is the partition function that ensures normalization).

###  VII.2 Landauer formula and its generalizations

All the observable quantities follow from the solution of the scattering
problem, namely from the quantities $\Phi$, $\Lambda$, $S_{\mathrm{tp}}$,
$\Psi_{\text{sr}}$. The most important object is the scattering matrix
$S=S_{\mathrm{pp}}$ between propagating modes. We use the abbreviated notation
$S_{\alpha\beta}=[S_{\mathrm{pp}}]_{\alpha\beta}$ for the corresponding matrix
elements between channel $\alpha$ and $\beta$. Below we review the formulae
used for the most common observables calculated in quantum transport.

####  VII.2.1 Electric current and conductance

In analogy to the derivation of Eq. (13), the out-of-equilibrium current that
flows through lead $a$ reads,

$I_{a}=\frac{e}{h}\sum_{b}\int dE\sum_{\alpha\in a,\beta\in
b}\left|S_{\alpha,\beta}\right|^{2}\left[f_{b}(E)-f_{a}(E)\right]\,.$ (141)

Here, we use the convention that current $I_{a}>0$ means current flowing
towards the scattering region and the electron charge is $-e$ with $e>0$.
Using the unitarity of the scattering matrix $S$, this expression can be
rewritten as

$I_{a}=\frac{e}{h}\sum_{b}\int dE\,g_{ab}(E)\,f_{b}(E)\,,$ (142)

where $g_{ab}(E)$ is defined as

$g_{ab}(E)=\sum_{\alpha\in a,\beta\in
b}(|S_{\alpha,\beta}|^{2}-\delta_{\alpha,\beta}\delta_{a,b})\,.$ (143)

For $a\neq b$, $g_{ab}(E)$ is the total transmission from lead $b$ to lead
$a$. $g_{aa}(E)$ is negative, with $|g_{aa}(E)|$ equal to the total
transmission from lead $a$ to all other leads $b$.

While Eq. (142) is now very well established and has been validated in
particular by the observation of the quantization of conductance in quantum
point contacts, several other formula were initially proposed [51].

When all the leads have equal temperatures and the voltage applied $V_{a}$ are
small $\mu_{a}=E_{F}-eV_{a}$, we get in linear response,

$I_{a}=\sum_{b}\mathcal{G}_{ab}V_{b}\,,$ (144)

with the conductance matrix defined as,

$\mathcal{G}_{ab}=\frac{e^{2}}{h}\int\ dEg_{ab}(E)\ \frac{\partial f}{\partial
E}(E)\,.$ (145)

At zero temperature, this expression further reduces to what is known as the
multi-terminal Landauer formula [49].

$\mathcal{G}_{ab}=-\frac{e^{2}}{h}g_{ab}(E_{F})\,.$ (146)

The $N_{\mathrm{leads}}\times N_{\mathrm{leads}}$ conductance matrix
$\mathcal{G}_{ab}$ satisfies current conservation,
$\sum_{a}\mathcal{G}_{ab}=0$. Additionally gauge invariance requires that
adding a constant to all voltages does not change currents, therefore
$\sum_{b}\mathcal{G}_{ab}=0$. In other words the vector
$\mathbb{1}=(1,1,\ldots,1)^{T}$ is both a left and a right eigenvector of the
$\mathcal{G}_{ab}$ matrix with eigenvalue zero. Note that both relations are
due to the unitarity of the scattering matrix $S$.

A commonly occurring measurement technique of multiterminal samples is to send
a current between two electrodes, while measuring all the possible voltage
differences. In particular, the four terminal resistance is

$R_{ab,cd}=\frac{V_{a}-V_{b}}{I}\,,$ (147)

with $V_{a}$ and $V_{b}$ corresponding to $I_{c}=-I_{d}=I$, and all the other
currents equal to 0. The four terminal resistance can correspond to
longitudinal or Hall resistance depending on the measurement configuration.
For $c=a$ and $d=b$ it corresponds to a two-terminal resistance. Because the
conductance matrix has a zero eigenvalue, straightforward inversion of Eq.
(144) to calculate $R_{ab,cd}$ is impossible. A possible workaround is
eliminating a single row (the corresponding current can be calculated at a
later stage using current conservation) and a column (this is equivalent to
fixing the voltage on the corresponding terminal to 0) of the conductance
matrix and invert the truncated conductance matrix. An alternative solution is
to use our knowledge of the zero eigenvector of $\mathcal{G}$ to construct the
resistance matrix as a pseudoinverse. Introducing the projector
$\mathbb{P}=\mathbb{1}\mathbb{1}^{T}/N_{\mathrm{leads}}$ onto the kernel of
$\mathcal{G}$, one can simply shift the matrix $\mathcal{G}$ by an arbitrary
constant $a\neq 0$ times $\mathbb{P}$ to make it invertible. This allows one
to define the voltage $V$ as,

$V=(\mathcal{G}-a\mathbb{P})^{-1}I\,.$ (148)

Indeed, when the vector $I$ satisfies current conservation (i.e.,
$\mathbb{P}I=0$), Eq.(148) implies that the average of $V$ vanishes (i.e.,
$\mathbb{P}V=0$ since, i.e., $\mathbb{P}\mathcal{G}=0$) and it follows that
$V$ satisfies $I=\mathcal{G}V$.

Note that the conductance matrix from Eq. (144) only takes into account the
device described by the scattering region. In an experiment, the device will
be connected to electric lines that have a resistance that can be comparable
to the one of the device. In that situation, reconstructing the conductance
matrix needs additional care as one needs to embed it inside a classical
electric circuit. This can be further complicated if nonlinear effects start
to play a role.

####  VII.2.2 Conductance in normal/superconductor hybrids

If one grounded superconducting electrode is present in the system, the
conductance can still be described with the scattering formalism using the
Bogoliubov-De Gennes wave function (see also Sec. IX.5). In this case, charges
may enter the superconductor in Cooper pairs resulting in Andreev reflection a
process when an incoming electron $e$ converts into an outgoing hole $h$ in
the same lead via Andreev reflection, or in a different lead via crossed
Andreev reflection.

The current in a normal lead $a$ is then given as [152]

$I_{a}=\frac{e}{h}\sum_{b}\int
dEg_{ab}(E)\left[f_{b}(E)-f_{\text{S}}(E)\right]\,,$ (149)

where the summation is over all normal leads $b$, $f_{\text{S}}(E)$ the Fermi-
Dirac distribution function of the grounded superconductor, energy $E$ is
defined with respect to the Fermi level of the superconductor, and

$g_{ab}(E)=\sum_{\alpha\in a,\beta\in b}|S_{\alpha e,\beta e}|^{2}-|S_{\alpha
h,\beta e}|^{2}-\delta_{\alpha e,\beta e}\delta_{a,b}\,.$ (150)

This conductance matrix is only valid for the normal leads, and alone does not
exhibit current conservation. The reason for this is the presence of the
superconductor, that not only allows current in the quasi-particle continuum
but also in the Cooper pair condensate.

In typical experiments, bias voltages are on the order of the superconducting
gap. Since the conductance of a normal/superconductor hybrid changes
drastically over this energy scale [36], conductance is usually not considered
in linear response in this case. Instead, it is obtained by taking the
derivative of Eq. (149) with respect to the voltage $V_{b}$. In the limit of
temperatures much lower than the superconducting gap, the Fermi-Dirac
distribution is well-approximated by a step function and we find

$\frac{dI_{a}}{dV_{b}}=-\frac{e^{2}}{h}g_{ab}(-eV_{b})\,.$ (151)

Cases with more than one superconducting lead are beyond the scope of this
review. Then, currents can flow in equilibrium (DC Josephson effect) and
voltage differences can lead to time-dependence (AC Josephson effect, multiple
Andreev reflection [22]).

####  VII.2.3 Thermoelectric effects

Equation (142) can also be used to calculate Seebeck effects by introducing
temperature gradients [261]. In the limit of small gradients,
$k_{B}\Theta_{a}=k_{B}\Theta_{\mathrm{av}}+k_{B}\delta\Theta_{a}$, we find

$I_{a}=\sum_{b}\mathcal{G}^{S}_{ab}\,\delta\Theta_{b}\,,$ (152)

with the Seebeck conductance matrix defined as

$\mathcal{G}^{S}_{ab}=\frac{e}{h\Theta_{\mathrm{av}}}g_{ab}^{(1)}\,,$ (153)

where we have introduced

$g_{ab}^{(n)}=\int dE\ g_{ab}(E)\left(-\frac{\partial f}{\partial
E}\right)(E)(E-E_{F})^{n}\,.$ (154)

At low temperatures $\Theta_{\mathrm{av}}$, the expression simplifies after a
Bohr-Sommerfeld expansion to

$\mathcal{G}^{S}_{ab}=\frac{e\pi^{2}k_{B}^{2}\Theta_{\mathrm{av}}}{3h}\frac{\partial
g_{ab}}{\partial E}(E_{F})\,.$ (155)

In analogy with the charge current, one can further introduce the heat current
as the energy flowing in the lead, counted from the chemical potential,

$I_{a}^{H}=-\frac{1}{h}\sum_{b}\int
dE\,g_{ab}(E)(E-\mu_{a})\left[f_{b}(E)-f_{a}(E)\right]\,.$ (156)

Here we use sign conventions such that $I_{a}^{H}>0$ means heat flowing
towards the scattering region. In the limit of small bias voltages, one gets
the Peltier effect

$I_{a}^{H}=\sum_{b}\mathcal{G}^{P}_{ab}\,\,V_{b}$ (157)

with the Peltier conductance matrix defined as

$\mathcal{G}^{P}_{ab}=\frac{e}{h}g_{ab}^{(1)}=\Theta_{\mathrm{av}}\mathcal{G}^{S}_{ab}\,.$
(158)

The direct correspondence with the Seebeck effect is a consequence of the
thermodynamic Onsager relations. Last, in the limit of small temperature
gradients, we get

$I_{a}^{H}=\sum_{b}\mathcal{G}^{H}_{ab}\,\delta\Theta_{b}\,,$ (159)

with the Heat conductance matrix defined as,

$\mathcal{G}^{H}_{ab}=\frac{1}{h\Theta_{\mathrm{av}}}g_{ab}^{(2)}\,.$ (160)

In the limit of small temperature, the Heat conductance simplifies into

$\mathcal{G}^{H}_{ab}=-\frac{\pi^{2}k_{B}^{2}\Theta_{\mathrm{av}}}{3h}g_{ab}(E_{F})\,.$
(161)

We note that in this limit, we recover

$\frac{\mathcal{G}^{H}_{ab}}{\Theta_{\mathrm{av}}\mathcal{G}_{ab}}=\frac{\pi^{2}k_{B}^{2}}{3e^{2}}\,,$
(162)

the (mesoscopic version of the) Wiedemann-Franz law.

####  VII.2.4 Spin currents

A Landauer formula can also be derived for spin currents. It can be used to
calculate, e.g., the spin torque exerted on a magnetic layer by an injected
current [294, 37]. The definition of the spin current, which provides an
explicit link between spin current and spin torque, is obtained from the
continuity equation of the current [293]. Introducing the vector of Pauli
matrices $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})^{T}$, the spin
current $\vec{I}$ flowing through a surface $\Sigma$ takes the form
$\propto\int_{\vec{r}\in\Sigma}d\vec{r}\ \mathrm{Im}\
\sum_{\eta\eta^{\prime}}\Psi_{\eta}(\vec{r})\vec{\sigma}_{\eta\eta^{\prime}}\vec{\nabla}\Psi_{\eta^{\prime}}(\vec{r})$
where $\Psi_{\eta}(\vec{r})$ is the wave function for the spin component
$\eta$. In other words, the spin current is formally equivalent to the charge
current with an addition of Pauli matrices in its definition. In the limit of
zero temperature and small bias voltages, the spin current takes the form

$\vec{I}_{a}=\sum_{b}\vec{\mathcal{G}}_{ab}\,V_{b}\,,$ (163)

where the spin conductance matrix $\vec{\mathcal{G}}_{ab}$ between non-
magnetic leads $a$ and $b$ takes the form

$\vec{\mathcal{G}}_{ab}=\frac{e}{4\pi}\sum_{\eta\eta^{\prime}\eta^{\prime\prime}}\sum_{\alpha\in
a,\beta\in
b}\vec{\sigma}_{\eta\eta^{\prime}}S_{\alpha\eta^{\prime},\beta\eta^{\prime\prime}}(E_{F})S_{\alpha\eta,\beta\eta^{\prime\prime}}^{*}(E_{F})\,.$
(164)

Here, the scattering matrix is defined for states with a well-defined spin
component $\eta$, i.e., this expression is only applicable to leads where spin
is a conserved quantity. Spin currents, being transfer of spin quanta $\hbar$
per unit time, have the dimensions of energy.

Note that in contrast to charge current, there is no conservation of spin
current in a magnetic system. Hence, it may be needed to calculate the spin
current inside the system as well as in the lead. This can be done using Eq.
(138). There is also no “Gauge invariance” so that spin current can flow in
the system at equilibrium. In the absence of spin-orbit coupling, such a spin
current can be identified with magnetic exchange interaction [43, 293].

####  VII.2.5 Injectivities and Emissivities

Two secondary concepts of the Landauer-Büttiker formalism are the emissivity
and the injectivity [47, 50]. The injectivity $\partial\langle
c_{n}^{\dagger}c_{n}\rangle/\partial\mu_{a}$ expresses how much the electronic
density $\langle c_{n}^{\dagger}c_{n}\rangle$ varies when one raises the
chemical potential $\mu_{a}$ of lead $a$. The emissivity $\partial
I_{a}/\partial(H_{\text{sr}})_{nn}$ is the variation of current $I_{a}$ upon a
change of electric potential $(H_{\text{sr}})_{nn}$. These two concepts arise
naturally when one extends quantum transport to AC regimes or non-linear
regimes to introduce a minimum treatment of the effect of electron-electron
interactions, see the discussion around Eq. (230).

The injectivity arises naturally in the formalism and reads at zero
temperature,

$\frac{\partial\langle
c_{n}^{\dagger}c_{n}\rangle}{\partial\mu_{a}}=\frac{1}{2\pi}\sum_{\alpha\in
a}|\Psi_{\text{sr}}(\mu_{\alpha})|_{n\alpha}^{2}\,,$ (165)

where the scattering wave function $\Psi_{\text{sr}}$ is calculated at
$E=\mu_{\alpha}$. Other generalizations can be obtained from Eq. (138) in a
straightforward manner.

To obtain the emissivity, one differentiates Eq. (50) and arrives at the
following linear system,

$\begin{pmatrix}H_{\text{sr}}-E&P_{\text{sr}}^{T}V^{\dagger}\Phi_{\mathrm{t+}}\Lambda_{\mathrm{t+}}\\\
VP_{\text{sr}}&-V\Phi_{\mathrm{t+}}\end{pmatrix}\begin{pmatrix}\frac{\partial\Psi_{\text{sr}}}{\partial(H_{\text{sr}})_{nm}}\\\
\frac{\partial
S_{\mathrm{tp}}}{\partial(H_{\text{sr}})_{nm}}\end{pmatrix}=\begin{pmatrix}-\Xi\Psi_{\text{sr}}\\\
0\end{pmatrix}$ (166)

that must be solved numerically. Here, the matrix $\Xi$ is defined as
$\Xi_{ij}=\delta_{in}\delta_{jm}$. Equation (166) can be used to implement
automatic differentiation schemes in quantum transport calculations. Eq. (166)
has the same structure as Eq. (93) from which it follows that the emissivity
can also be obtained from the knowledge of the Green function,

$\displaystyle\frac{\partial\Psi_{\text{sr}}}{\partial(H_{\text{sr}})_{nm}}=G_{\text{sr}}\Xi\Psi_{\text{sr}}$
(167) $\displaystyle\frac{\partial
S_{\mathrm{tp}}}{\partial(H_{\text{sr}})_{nm}}=G_{\mathrm{t+}}\Xi\Psi_{\text{sr}}$
(168)

###  VII.3 Quantum noise

The observables defined in the preceding sections are all one-body
observables, i.e., mean observables. Another class of important observables
are the quantum fluctuations around the mean and in particular the quantum
fluctuation of the current [35]. The general expression of the noise
correlation $\bar{S}_{ab}(\omega)$ between lead $a$ and $b$ (noise power when
$a=b$) at finite frequency $\omega$ reads,

$\begin{split}\bar{S}_{ab}&(\omega)=\frac{e^{2}}{h}\sum_{c,d}\sum_{\gamma\in
c,\delta\in d}\int dE\\\
&A^{cd}_{\gamma\delta}(a,E,E+\hbar\omega)A^{dc}_{\delta\gamma}(b,E+\hbar\omega,E)\\\
&\\{f_{c}(E)[1-f_{d}(E+\hbar\omega)]+f_{d}(E+\hbar\omega)[1-f_{c}(E)]\\}\end{split}$
(169)

where $A^{cd}_{\gamma\delta}(a,E,E^{\prime})$ is defined in terms of the
scattering matrix as,

$A^{cd}_{\gamma\delta}(a,E,E^{\prime})=\delta_{ac}\delta_{ad}\delta_{\gamma\delta}-\sum_{\alpha\in
a}S_{\alpha\gamma}^{*}(E)S_{\alpha\delta}(E^{\prime})$ (170)

Eq. (169) has a number of simplified forms, notably at zero frequency, small
bias, zero temperatures, etc. We refer to [35] for a comprehensive review. See
also the next section for a simple route for deriving such formula.

##  VIII Physical observables: the Non-equilibrium Green’s function approach

In this section, we take an alternative, yet equivalent, route compared to the
one in the preceding section to calculate out-of-equilibrium physical
observables such as the electric conductance from the solution of the quantum
mechanical problem. While, in the preceding section, the latter was expressed
in terms of the scattering matrix and scattering states, here quantum
propagation is encoded in the retarded Green’s function.

The connection from the retarded Green’s function to the physical observable
is done through the so-called Keldysh formalism [133]. As we do not consider
electron-electron interactions, but simpler quadratic Hamiltonians, the
Keldysh formalism considerably simplifies into what is often known as NEGF
(Non Equilibrium Green Function’s) formalism [55, 191]. NEGF is mathematically
equivalent to the Landauer-Büttiker approach. It is perhaps less intuitive but
has the advantage that once a few basic equations have been established,
everything follows from straightforward algebra.

###  VIII.1 Keldysh formalism in a nutshell

In this subsection, we summarize the main definitions and results of the
Keldysh formalism that will be needed. We refer the reader to textbooks such
as [233, 266] for the derivation of the results that we state here without
proofs.

The Keldysh formalism introduces two independent Green’s functions, the
“lesser” ($<$) and “greater” ($>$) Green’s functions

$\displaystyle\hat{G}^{<}_{nm}(\tau)$ $\displaystyle=$ $\displaystyle i\langle
c^{\dagger}_{m}e^{i{\cal H}\tau}c_{n}e^{-i{\cal H}\tau}\rangle,$ (171)
$\displaystyle\hat{G}^{>}_{nm}(\tau)$ $\displaystyle=$ $\displaystyle-i\langle
e^{i{\cal H}\tau}c_{n}e^{-i{\cal H}\tau}c^{\dagger}_{m}\rangle.$ (172)

These Green’s functions are very convenient since at $\tau=0$ they correspond
directly to the observables (e.g., density or current) of interest. From these
two Green’s functions, one may construct the “time-ordered” ($T$) , “anti-time
ordered” ($\bar{T}$), “retarded” ($R$) and “advanced” ($A$) Green’s functions
as,

$\displaystyle\hat{G}^{T}_{nm}(\tau)$ $\displaystyle=$
$\displaystyle\theta(\tau)\hat{G}^{>}_{nm}(\tau)+\theta(-\tau)\hat{G}^{<}_{nm}(\tau),$
(173) $\displaystyle\hat{G}^{\bar{T}}_{nm}(\tau)$ $\displaystyle=$
$\displaystyle\theta(\tau)\hat{G}^{<}_{nm}(\tau)+\theta(-\tau)\hat{G}^{>}_{nm}(\tau),$
(174) $\displaystyle\hat{G}^{R}_{nm}(\tau)$ $\displaystyle=$
$\displaystyle\theta(\tau)\hat{G}^{>}_{nm}(\tau)-\theta(\tau)\hat{G}^{<}_{nm}(\tau),$
(175) $\displaystyle\hat{G}^{A}_{nm}(\tau)$ $\displaystyle=$
$\displaystyle\theta(-\tau)\hat{G}^{<}_{nm}(\tau)-\theta(-\tau)\hat{G}^{>}_{nm}(\tau).$
(176)

Note that, in the preceding sections, we encountered only one type of Green’s
function, the retarded Green’s function, and the subscript $R$ has been
omitted. These Green’s functions are linearly dependent. Below, we will focus
on the Green’s functions $\hat{G}^{<}$, $\hat{G}^{R}$ and $\hat{G}^{A}$ using
the fact that
$\hat{G}^{<}=\hat{G}^{T}-\hat{G}^{R}=\hat{G}^{\bar{T}}+\hat{G}^{A}$ as can be
proven trivially from the above definitions.

In this review, we only consider many-body Hamiltonians ${\cal H}$ restricted
to be quadratic in the destruction $c_{n}$ and creation $c^{\dagger}_{n}$
operators on site $n$,

${\cal H}=\sum_{nm}\hat{H}_{nm}c^{\dagger}_{n}c_{m}.$ (177)

The equivalence between the definition of the retarded Green’s function used
in the preceding sections and the one in Eq. (175) is not entirely trivial and
indeed valid only for quadratic Hamiltonians. In the following, we make use of
this simplification.

In the Keldysh formalism, these Green’s functions naturally appear in the form
of a $2\times 2$ matrix associated with the so-called Keldysh contour,

$\displaystyle\mathbf{\hat{G}}_{nm}(\tau)=\begin{pmatrix}\hat{G}^{T}_{nm}(\tau)&\hat{G}^{<}_{nm}(\tau)\\\
\hat{G}^{>}_{nm}(\tau)&\hat{G}^{\bar{T}}_{nm}(\tau)\end{pmatrix}.$ (178)

Since this review is concerned with d.c. quantum transport (i.e., time-
independent Hamiltonians), we’ll work chiefly in energy representation,

$\mathbf{\hat{G}}_{nm}(E)=\int d\tau\mathbf{\hat{G}}_{nm}(\tau)e^{iE\tau}.$
(179)

Now that we have defined the mathematical objects of interest, we can state
the very few results from non-equilibrium many-body theory that will be needed
to establish the basic equations of quantum transport. The equations of motion
read

$\displaystyle\sum_{p}\begin{pmatrix}E\delta_{mp}-\hat{H}_{mp}&0\\\
0&E\delta_{mp}-\hat{H}_{mp}\end{pmatrix}\mathbf{\hat{G}}_{pn}(E)=\begin{pmatrix}\delta_{mn}&0\\\
0&\delta_{mn}\end{pmatrix},$ (180)

or in compact notations,

$\displaystyle\begin{pmatrix}E-\hat{H}&0\\\
0&E-\hat{H}\end{pmatrix}\mathbf{\hat{G}}(E)=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}.$ (181)

Splitting the Hamiltonian matrix $\hat{H}$ into an “unperturbed” matrix
$\hat{H}_{0}$ and a “perturbation” $\hat{W}$,

$\hat{H}=\hat{H}_{0}+\hat{W},$ (182)

and denoting $\mathbf{\hat{G}}_{0}(E)$ the Green’s function associated with
$\hat{H}_{0}$, one can write a form of Dyson equation [233, 266] as

$\displaystyle\mathbf{\hat{G}}(E)=\mathbf{\hat{G}}_{0}(E)+\mathbf{\hat{G}}_{0}(E)\begin{pmatrix}\hat{W}&0\\\
0&-\hat{W}\end{pmatrix}\mathbf{\hat{G}}(E),$ (183)

and also

$\displaystyle\mathbf{\hat{G}}(E)=\mathbf{\hat{G}}_{0}(E)+\mathbf{\hat{G}}(E)\begin{pmatrix}\hat{W}&0\\\
0&-\hat{W}\end{pmatrix}\mathbf{\hat{G}}_{0}(E).$ (184)

Note that the partition of $\hat{H}$ into $\hat{H}_{0}$ and $\hat{W}$ is
arbitrary, although in what follows the $\hat{W}$ matrix will contain the
matrix elements that connect the leads to the scattering region and
$\hat{H}_{0}$ the rest of the Hamiltonian. Identifying the different
components of the $2\times 2$ Dyson equation Eq. (183), one may write

$\displaystyle\hat{G}^{R}$ $\displaystyle=$
$\displaystyle\hat{G}^{R}_{0}+\hat{G}^{R}_{0}\hat{W}\hat{G}^{R},$ (185)
$\displaystyle\hat{G}^{A}$ $\displaystyle=$
$\displaystyle\hat{G}^{A}_{0}+\hat{G}^{A}_{0}\hat{W}\hat{G}^{A},$ (186)
$\displaystyle\hat{G}^{<}$ $\displaystyle=$
$\displaystyle\hat{G}^{<}_{0}+\hat{G}^{R}_{0}\hat{W}\hat{G}^{<}+\hat{G}^{<}_{0}\hat{W}\hat{G}^{A}.$
(187)

The first equation (or equivalently the second one) encodes the quantum
mechanics and forms a close set of equations for the retarded Green’s function
alone. The third equation encodes the statistical physics. Similarly, Eq.
(184) leads to

$\displaystyle\hat{G}^{R}$ $\displaystyle=$
$\displaystyle\hat{G}^{R}_{0}+\hat{G}^{R}\hat{W}\hat{G}^{R}_{0},$ (188)
$\displaystyle\hat{G}^{A}$ $\displaystyle=$
$\displaystyle\hat{G}^{A}_{0}+\hat{G}^{A}\hat{W}\hat{G}^{A}_{0},$ (189)
$\displaystyle\hat{G}^{<}$ $\displaystyle=$
$\displaystyle\hat{G}^{<}_{0}+\hat{G}^{R}\hat{W}\hat{G}^{<}_{0}+\hat{G}^{<}\hat{W}\hat{G}^{A}_{0}.$
(190)

For a system at thermal equilibrium, one has

$\hat{G}^{<}(E)=-f(E)\left[\hat{G}^{R}(E)-\hat{G}^{A}(E)\right],$ (191)

which completes the basic formalism that we shall need in this section.

###  VIII.2 Main NEGF equation

Within the NEGF formalism, computing out-of-equilibrium amounts to performing
two tasks. First one computes the retarded Green’s function (i.e., solve the
quantum mechanical problem) using, e.g., any algorithm from Sec. VI . Second,
one must compute the lesser Green’s function (i.e., the observables). In this
subsection, we derive the main result of NEGF that relates the lesser Green’s
function to the retarded one in the context of quantum transport. The
derivation is a generalization of the original work by Caroli et al. [55].
These results have later been generalized and popularized by Wingreen and Meir
[191].

We start by splitting and applying the Dyson equation to
$\hat{H}_{\text{sys}}$. The perturbation $\hat{W}$ contains the matrix
elements that connect the leads to the scattering region. $\hat{H}_{0}$
contains the Hamiltonian of the disconnected elements of the system: the semi-
infinite leads and the isolated scattering region. Let us write Dyson equation
in block form where the block index indicates the lead ($a=1,2,\dots$) or the
scattering region ($0$). By construction, the disconnected Green’s functions
$(\hat{G}_{0})_{kl}=\delta_{kl}(\hat{G}_{0})_{l}$ are diagonal with respect to
the block index. Eq. (187) for the scattering region block
${G}^{<}_{\text{sr}}=\hat{G}^{<}_{00}$ reads

$\hat{G}^{<}_{\text{sr}}=(\hat{G}^{<}_{0})_{0}+(\hat{G}^{R}_{0})_{0}\hat{W}_{0a}\hat{G}^{<}_{a0}+(\hat{G}^{<}_{0})_{0}\hat{W}_{0a}\hat{G}^{A}_{a0},$
(192)

with an implicit summation over leads $a>0$. Note that in the rest of this
section we have absorbed the projection matrices $P_{\text{sr}}$ into
$\hat{W}$, such that, e.g., $\hat{W}_{a0}=V_{a}P_{\text{sr}}$, where $V_{a}$
is the hopping matrix of lead $a$. Multiplying Eq. (192) on the left by
$E+i\eta-H_{\text{sr}}$ and using $[E-H_{\text{sr}}](\hat{G}^{<}_{0})_{0}=0$
and $[E-H_{\text{sr}}](\hat{G}^{R}_{0})_{0}=1$ we obtain

$\left[E-H_{\text{sr}}\right]{G}^{<}_{\text{sr}}=\hat{W}_{0a}\hat{G}^{<}_{a0}.$
(193)

To close the above equation, we write Eq. (187) and (190) for the $a0$ and
$0a$ blocks respectively,

$\hat{G}^{<}_{a0}=(\hat{G}^{R}_{0})_{a}\hat{W}_{a0}{G}^{<}_{\text{sr}}+(\hat{G}^{<}_{0})_{a}\hat{W}_{a0}{G}^{A}_{\text{sr}},$
(194)
$\hat{G}^{<}_{0a}={G}^{R}_{\text{sr}}\hat{W}_{0a}(\hat{G}^{<}_{0})_{a}+{G}^{<}_{\text{sr}}\hat{W}_{0a}(\hat{G}^{A}_{0})_{a}.$
(195)

Inserting Eq. (194) into (193), we arrive at

$\left[E-H_{\text{sr}}-\hat{W}_{0a}(\hat{G}^{R}_{0})^{a}\hat{W}_{a0}\right]{G}^{<}_{\text{sr}}=\hat{W}_{0a}(\hat{G}^{<}_{0})_{a}\hat{W}_{a0}{G}^{A}_{\text{sr}}.$
(196)

We recognize Eq. (96) on the left hand side, so that,

${G}^{<}_{\text{sr}}={G}^{R}_{\text{sr}}\hat{W}_{0a}(\hat{G}^{<}_{0})_{a}\hat{W}_{a0}{G}^{A}_{\text{sr}}.$
(197)

To conclude, we use Eq. (191) supposing each lead to remain at equilibrium
with its Fermi function $f_{a}(E)$. We arrive at the main result of NEGF
theory,

${G}^{<}_{\text{sr}}(E)={G}^{R}_{\text{sr}}(E)\Sigma_{\text{sr}}^{<}(E){G}^{A}_{\text{sr}}(E),$
(198)

with the lesser self-energy $\Sigma_{\text{sr}}^{<}$ defined as

$\Sigma_{\text{sr}}^{<}(E)=-\sum_{a}f_{a}(E)\left[\Sigma_{a}^{R}(E)-\Sigma_{a}^{A}(E)\right],$
(199)

where

$\Sigma^{R}_{a}=\hat{W}_{0a}(\hat{G}^{R}_{0})_{a}\hat{W}_{a0}$ (200)

is the partial contribution to the retarded self-energy
$\Sigma^{R}(E)=\sum_{a}\Sigma^{R}_{a}$. Since
${G}^{A}_{\text{sr}}=({G}^{R}_{\text{sr}})^{\dagger}$ and
${\Sigma}^{A}_{a}=({\Sigma}^{R}_{a})^{\dagger}$, we can equivalently write Eq.
(199) as

$\Sigma_{\text{sr}}^{<}(E)=i\sum_{a}f_{a}(E)\Gamma_{a}(E)\,,$ (201)

using definition of the linewidth (106) for lead $a$,
$\Gamma_{a}=i(\Sigma^{R}_{a}-\Sigma^{A}_{a})$.

Eq. (198) provides a direct expression of the Lesser Green’s function in terms
of the Retarded one and the occupation function in the leads. From there, one
can directly calculate physical observables as $\langle
c^{\dagger}_{m}c_{n}\rangle=-iG^{<}_{nm}(\tau=0)$ which leads to the one-body
density matrix,

$\langle
c^{\dagger}_{m}c_{n}\rangle=-i\int\frac{dE}{2\pi}(G^{<}_{\text{sr}})_{nm}(E),$
(202)

which generalizes the concept of local density of states to off-diagonal
observables and out-of-equilibrium situations.

Inserting Eq. (111) into Eqs. (201) and (198) as well as using Eq. (114) one
can relate the lesser Green’s function to the scattering wave-functions; one
arrives at [306, 97]

$G^{<}_{\rm sr}(E)_{nm}=i\sum_{a}f_{a}(E)\sum_{\alpha\in
a}\Psi_{\text{sr}}(E)_{n\alpha}\Psi_{\text{sr}}^{*}(E)_{m\alpha}.$ (203)

At equilibrium, this equation simplifies into $G^{<}_{\rm
sr}(E)=i\Psi_{\text{sr}}(E)\Psi_{\text{sr}}^{\dagger}(E)f(E)$ which can be
derived straightforwardly using the definition of $G^{<}$ [266].

It is worth noting that there is a subtle point in the derivation above: the
limit $\eta\rightarrow 0$ has to be taken at the end of the calculation, not
at the beginning, whenever the energy $E$ matches a bound states of the
system, owing to the emergence of (diverging) Dirac functions at these
energies. It follows that, strictly speaking, one should add the bound states
contributions to Eq. (198) or (203) which only contain contributions from the
scattering states. For example, Eq. (203) then reads

$\begin{split}G^{<}_{\rm sr}(E)_{nm}&=i\sum_{a}f_{a}(E)\sum_{\alpha\in
a}\Psi_{\text{sr}}(E)_{n\alpha}\Psi_{\text{sr}}^{*}(E)_{m\alpha}\\\
&+i\sum_{j}f_{j}(E)\delta(E-E_{j})\left(\Psi_{\text{sr}}\right)_{nj}\left(\Psi_{\text{sr}}^{*}\right)_{mj},\end{split}$
(204)

where the index $j$ now labels bound states at energy $E_{j}$ with
corresponding wave function $(\Psi_{\text{sr}})_{nj}$. These contributions are
usually omitted because (i) they do not contribute to transport properties and
(ii) the filling factor $f_{i}(E)$ depends on the history of the system and
must be decided based on physics not described within our framework [266].
These bound states may play a role nevertheless in certain situations [229].

At equilibrium, the expression for ${G}^{<}_{\text{sr}}$ can be simplified. To
do so, we use relation for the Retarded and Advanced Green’s function obtained
by (i) multiplying the equation
${G}^{R}_{\text{sr}}(E-H_{\text{sr}}-\Sigma^{R})=1$ by ${G}^{A}_{\text{sr}}$
on the right, (ii) multiplying its Hermitian conjugate
$(E-H_{\text{sr}}-\Sigma^{A}){G}^{A}_{\text{sr}}=1$ by ${G}^{R}_{\text{sr}}$
on the left and (iii) subtracting the second equation from the first. One
arrives at

${G}^{R}_{\text{sr}}[\Sigma^{R}-\Sigma^{A}]{G}^{A}_{\text{sr}}={G}^{R}_{\text{sr}}-{G}^{A}_{\text{sr}}.$
(205)

Inserting Eq. (205) into Eq. (198), we get

${G}^{<}_{\text{sr}}=-f(E)[{G}^{R}_{\text{sr}}-{G}^{A}_{\text{sr}}],$ (206)

recovering Eq. (191). Calculating the number of electrons on one site at
equilibrium using Eq. (202), we get back the usual expression in terms of the
local density of state:

$\langle
c^{\dagger}_{n}c_{n}\rangle=-\int\frac{dE}{\pi}f(E)\mathrm{Im}(G^{R}_{\text{sr}})_{nn}(E).$
(207)

###  VIII.3 Landauer formula within NEGF

Let us now define the current flowing towards the scattering region and
establish the Green’s function version of the Landauer formula. Introducing
the total number of electrons in the scattering region
$\mathcal{Q}_{\text{sr}}(t)=\langle{\cal Q}\rangle$ with ${\cal
Q}=\sum_{n\in\text{sr}}c^{\dagger}_{n}c_{n}$, the evolution of ${\cal Q}$ (in
Heisenberg representation) $\partial{\cal Q}/\partial t=i[{\cal H},{\cal Q}]$
provides the continuity equation

$\frac{\partial\mathcal{Q}_{\text{sr}}}{\partial t}=\sum_{a}I_{a}(t)\,.$ (208)

In steady state, i.e., time-independent current flow $I_{a}(t)=I_{a}$, we find
the current $I_{a}$ from lead $a$ towards the scattering region as

$I_{a}=\frac{e}{h}\int
dE\operatorname{Tr}[G^{<}_{a0}\hat{W}_{0a}-\hat{W}_{a0}G^{<}_{0a}].$ (209)

Using Eq. (194), Eq.(195) and the definitions of the self energies we get

$I_{a}=\frac{e}{h}\int
dE\operatorname{Tr}\left[G^{<}_{\text{sr}}(\Sigma^{R}_{a}-\Sigma^{A}_{a})-(G^{R}_{\text{sr}}-G^{A}_{\text{sr}})\Sigma^{<}_{a}\right].$
(210)

Further using Eqs. (205) and Eq. (201), we finally arrive at

$I_{a}=\frac{e}{h}\sum_{b}\int
dE\operatorname{Tr}\left[G^{R}_{\text{sr}}\Gamma_{b}G^{A}_{\text{sr}}\Gamma_{a}\right]\left[f_{b}(E)-f_{a}(E)\right].$
(211)

This equation is often (somewhat abusively) referred to as the Landauer
formula. It automatically satisfies two important requirements: the absence of
incoming currents at equilibrium (note that there may be persistent currents
inside the scattering region), and the global conservation of current
$\sum_{a}I_{a}=0$.

To conclude this section, we end with the proof that Eq. (211) is indeed
equivalent to the (original) scattering matrix version of the Landauer
formula. Inserting Eq. (109) and Eq. (111) in Eq. (211) and using the
invariance of the trace under cyclic permutation, we get

$\begin{split}&\operatorname{Tr}\left[G^{R}_{\text{sr}}\Gamma_{b}G^{A}_{\text{sr}}\Gamma_{a}\right]\\\
&=\operatorname{Tr}\left[G^{R}_{\text{sr}}\Gamma_{b}\Phi_{\mathrm{p}-}^{b}\Phi_{\mathrm{p}-}^{\dagger
b}\Gamma_{b}G^{A}_{\text{sr}}\Gamma_{a}\Phi_{\mathrm{p}+}^{a}\Phi_{\mathrm{p}+}^{\dagger
a}\Gamma_{a}\right]\\\
&=\operatorname{Tr}\left[\left(i\Phi_{\mathrm{p}+}^{\dagger
a}\Gamma_{a}G^{R}_{\text{sr}}\Gamma_{b}\Phi_{\mathrm{p}-}^{b}\right)\left(-i\Phi_{\mathrm{p}-}^{\dagger
b}\Gamma_{b}G^{A}_{\text{sr}}\Gamma_{a}\Phi_{\mathrm{p}+}^{a}\right)\right]\\\
&=\operatorname{Tr}\left[S^{ab}(S^{ab})^{\dagger}\right],\end{split}$ (212)

where the $\Phi_{\mathrm{p}+}^{a(b)}$ matrix is restricted to the modes of
lead $a$ ($b$). we have introduced a new object $S^{ab}$,

$S^{ab}=i\Phi_{\mathrm{p}+}^{\dagger
a}\Gamma_{a}G^{R}_{\text{sr}}\Gamma_{b}\Phi_{\mathrm{p}-}^{b}\,.$ (213)

From the Fisher-Lee relation (116), we can indeed identify

$S^{ab}=S^{ab}_{\text{p}\text{p}},$ (214)

i.e., the scattering matrix between lead $a$ and lead $b$, restricted to the
propagating channels, concluding our of equivalence between the scattering
approach and the NEGF approach. Note that in the sum of Eq. (211) only $b\neq
a$ contribute, and hence Eq. (212) is valid only for $a\neq b$, i.e.,
describes transmission from lead $b$ to $a$.

###  VIII.4 Higher-order observables

The observables defined in the preceding sections are all one-body observables
(except for the quantum noise), i.e., quadratic in the creation/destruction
operators. To calculate higher-order observables, one can make use of the fact
that, even though the system is out of equilibrium, the density matrix is the
exponential of a (highly non-local) quadratic Hamiltonian, see Eq.(140). It
follows that Wick theorem applies, allowing expression of, e.g., quartic terms
in terms of quadratic ones,

$\displaystyle\langle c^{\dagger}_{n}(t)c_{m}(t)$ $\displaystyle
c^{\dagger}_{p}(t^{\prime})c_{q}(t^{\prime})\rangle=$ $\displaystyle\langle
c^{\dagger}_{n}(t)c_{m}(t)\rangle$ $\displaystyle\langle
c^{\dagger}_{p}(t^{\prime})c_{q}(t^{\prime})\rangle+\langle
c^{\dagger}_{n}(t)c_{q}(t^{\prime})\rangle\langle
c_{m}(t)c^{\dagger}_{p}(t^{\prime})\rangle,$

i.e., in terms of sum of products of Lesser and Greater Green’s functions.
This approach can be used directly to calculate quantities such as shot noise
[96].

The same approach can be used to derive expressions relating higher-order
observables to the scattering wavefunction. For this, one uses the following
dictionary [97] :

$[G^{<}(\tau)]_{nm}=i\sum_{\alpha}\int\frac{dE}{2\pi}(\Psi_{\text{sr}})_{n\alpha}(\Psi_{\text{sr}}^{*})_{m\alpha}e^{-iE\tau}f_{\alpha}(E),$
(216)

and

$[G^{>}(\tau)]_{nm}=-i\sum_{\alpha}\int\frac{dE}{2\pi}(\Psi_{\text{sr}})_{n\alpha}(\Psi_{\text{sr}}^{*})_{m\alpha}e^{-iE\tau}[1-f_{\alpha}(E)].$
(217)

Note that likewise, the retarded Green’s function can be expressed as

$[G^{R}(\tau)]_{nm}=-i\theta(\tau)\sum_{\alpha}\int\frac{dE}{2\pi}(\Psi_{\text{sr}})_{n\alpha}(\Psi_{\text{sr}}^{*})_{m\alpha}e^{-iE\tau}.$
(218)

##  IX Selected applications to model systems

We have now completed the technical part of this review where we presented the
mathematical formalisms and the associated algorithms used to make actual
numerical calculations. We are left to discuss how these techniques are used
for real world applications. Numerical quantum transport calculations are by
now very common in a wide variety of fields ranging from molecular electronics
through spintronics to nanoelectronics. They are applied in the study of
devices made of semiconductors, graphene, topological materials,
superconductors, and various combinations thereof. We shall not attempt to
review all these applications but will confine ourselves to selecting a few
examples that we find either important or illustrative. Such a choice is
necessarily subjective.

###  IX.1 A hierarchy of discrete models

c|c[cell-space-limits=8pt] 3-1Optionally

self-consistent [fill=[HTML]98db65]1-1Discretized $k\cdot p$ Hamiltonians

[fill=[HTML]c7db65]1-1Tight-binding models (TBM)

[fill=[HTML]dbc365]1-1Semi-empirical tight-binding models

(SE-TBM)

2-1Self-consistent [fill=[HTML]dba065]1-1First principles tight-binding models

(FP-TBM)

[fill=[HTML]db8165]1-1Density functional theory (DFT)

Table 2: Schematic of the different levels of modeling used in numerical
quantum transport, see text.

The first step in this regard is to select the appropriate level of modeling
for the Hamiltonian of the scattering region as well as of the (semi-infinite)
electrodes and of the system-electrode interface. There exists a hierarchy of
discrete Hamiltonians that can be used to describe a system, and physicists
must select an appropriate level of modeling, requiring a tradeoff between
accuracy and numerical tractability. Fig. IX.1 shows a schematic of this
hierarchy. At the bottom are calculations performed entirely within the
framework of density functional theory. In the absence of strong electronic
correlations, DFT provides an accurate description of the system from first
principles and has _de facto_ become the reference approach to the calculation
of realistic materials. However, strictly speaking, the Kohn-Sham equations of
DFT are valid for the ground state of the system and often need to be
corrected to properly account for dynamical properties such as the band gap.
Next come the First Principle Tight-Binding Models (FP-TBM) where a DFT
calculation is followed by the projection of the corresponding (self-
consistent) Hamiltonian onto a set of localized states such as maximally
localized Wannier functions. Semi-Empirical Tight-Binding Models (SE-TBM) are
also models that are derived from a DFT calculation. However, in contrast to
FP-TBM where the model is constructed in a systematic way, the SE-TBM are
parametrized and then fitted to properly reproduce the DFT calculations
(typically the band structure) and/or experimentally known properties of the
material (such as the band gap). The orbitals and terms in the Hamiltonian
included in the SE-TBM are typically dictated by symmetry and physics
considerations. When the Tight-Binding Model (TBM) is not tightly bound to a
first principle calculation, we simply call it a general TBM. These TBM’s can
be studied from a purely theoretical perspective or in situations where a
simple model is sufficient to account for the physics. These model systems can
also come from the discretization of a continuum model such as $\mathbf{k.p}$
model or a simple effective mass model. All these models describe essentially
the quantum mechanics of the conducting electrons of the system.

Another important ingredient is the electromagnetic and in particular
electrostatic environment in which these electrons propagate. In fact, the
electrostatic energy is usually the largest energy of the problem. In first
principle calculations, the electrostatic problem is solved self-consistently
with the quantum mechanical one. In model systems (general TBM and
$\mathbf{k\cdot p}$ models) one often assumes that the electric potential is
known and has a simple form.

In this section and the next, we will give a few selected examples to
illustrate the use of computational quantum transport. In the remaining part
of this section, we focus on “model” systems (the two upper types of models in
Fig. IX.1) while the next section will focus on more “realistic” models (the
three lower levels in Fig. IX.1).

###  IX.2 Applications to the two-dimensional electron gas

####  IX.2.1 The early days: Disordered systems and Anderson localization

An important application of numerics to quantum transport is the study of the
effect of static disorder on conductance as well as of the associated Anderson
localization phenomena. This is perhaps the historically first application of
the techniques described in this review and also one of the most successful
ones.

The model is the discretized effective mass Schrödinger equation Eq. (22) in
the presence of static disorder. In two dimensions the model reads

$\displaystyle\Psi_{n+1,m}$ $\displaystyle+$
$\displaystyle\Psi_{n-1,m}+\Psi_{n,m+1}+\Psi_{n,m-1}$ (219) $\displaystyle+$
$\displaystyle V_{n,m}\Psi_{n,m}=E\Psi_{n,m},$

where the potential $V_{n,m}$ takes different random values on different
sites. Typically, one chooses $V_{n,m}$ randomly with a flat distribution
$V_{n,m}\in[-W/2,W/2]$ where the parameter $W$ is the strength of the
disorder. In a seminal paper, the “gang of four” [1] had predicted that the
conductance $g(W,L)$ (where $L$ is the size of the sample) followed a simple
scaling law from which one could predict the insulating or metallic behavior
of the sample. The scaling hypothesis, which was supported by theoretical
arguments both in the strong and weak disorder regimes, was that the function
$\beta\equiv\partial\log g/\partial\log L$ was a function of $g$ only (i.e.,
$\exists f,\ \beta(L,W)=f[g(W,L)]$ ). One of the most striking conclusions of
the scaling theory was that in three dimensions there was a critical value of
disorder associated to a metal-insulator transition, but in one and two
dimensions the system was always an insulator in the sense that the
conductance would become vanishingly small in the $L\rightarrow\infty$ limit.

While the scaling hypothesis was plausible, it was merely a hypothesis. It was
put on firm grounds in a series of works combining numerical calculations with
finite-size scaling analysis. These calculations were pioneered by [223, 224]
who, however, wrongly concluded the existence of a metal-insulator transition
in two dimensions. Subsequent more precise numerics [182] could indeed confirm
the finite-size scaling hypothesis, as shown in Fig. 10. A large activity
followed in the following decade with the study of the effect of various
parameters such as spin-orbit coupling, magnetic field, type of disorder, as
reviewed in [142].

Figure 10: Scaling function $\beta(g)$ computed numerically for a 2D
(triangles) and 3D (pluses) tight-binding model, confirming the scaling
hypothesis. From [182].

Besides the Anderson localization regime, numerics have also been used to
study the quantum corrections to the classical limit in the diffusive regime.
An early work in this direction is the study of the universal conductance
fluctuations [270]. More recent applications include simulations of diffusive
metal-superconductor junctions [186] or disproving of the conjecture of the
existence of a two-parameter scaling regime in a system with partially broken
time reversal symmetry [255].

####  IX.2.2 Ballistic systems and Quantum billiards

In the early nineties, the opposite, ballistic regime (no disorder) attracted
wide interest in the context of “quantum billiards”. Quantum billiards consist
of a high mobility two-dimensional electron gas (typically GaAs/GaAlAs
heterostructures) where a billiard is patterned using electrostatic gates
deposited on top of the heterostructure. A study of such billiards (also known
as quantum dots in the regimes where they are almost closed and charging
energy plays an important role) was performed in [122] using the recursive
Green’s function algorithm. 4-terminal (X-shaped) and 3-terminal (T-shape)
ballistic junctions were studied in [24, 262]. One of the successes of these
ballistic calculations was the finding that the statistics of the height of
the resonances (in the almost closed regime) matched the Porter-Thomas
distribution prediction of random matrix theory [123]. Numerics also clarified
_weak localization_ , i.e., the interference contribution to the magneto-
conductance [25]. Numerics could in particular clearly distinguish the
difference in behavior between the magneto-conductance of _integrable_ and
_chaotic_ billiards: a chaotic billiard gives rise to a Lorentzian lineshape
for the magnetoconductance, whereas an integrable billiard features a
triangular lineshape, as shown in the left panels of Fig. 11. This prediction
was soon confirmed experimentally [57], shown in the right panels of Fig. 11.
In billiards where there exists a short straight path between the entrance and
the exit, the observation of Fano resonances [105] was also well reproduced in
numerics [62]. A review of quantum billiards can be found in [7]. Some recent
ballistic devices focus on building electronic interferometers analogous to
those known in optics [48]. Those systems can also be simulated both with
[130] (quantum Hall regime) and without [29] magnetic field.

Figure 11: Left: Numerically computed weak localization correction to the
conductance for a chaotic stadium billiard (top) and an integrable circular
billiard (bottom). Solid lines correspond to specular scattering at the
boundary, other lines to a situation with disorder at the boundary. Right:
Corresponding experimental data obtained by averaging over 48 billiards
defined in a high-quality GaAs electron gas. The lineshape in the chaotic case
is Lorentzian and for the integrable billiard triangular, in agreements with
simulation. Adapted from [57].

####  IX.2.3 Scanning gate microscopy

In a scanning gate microscope, one measures how the conductance through a
constriction (“a quantum point contact”) is affected by a conducting tip that
scans the surface of the heterostructure. Numerical simulations are very handy
for these systems as the signal is highly dependent on any residual disorder
as well as on the sample geometry. Numerics were shown to correctly capture
the experimental observation, at least on a qualitative level [141, 268].
Combined with machine learning techniques, numerical calculations can be used
to unveil the underlying disorder potential seen by the electrons [220].

###  IX.3 Quantum transport in graphene

####  IX.3.1 Single $p_{z}$-orbital tight-binding model

The discovery of graphene triggered a very large number of numerical studies
that is far beyond the scope of the present review [221]. Many of these works
use a simple next-nearest neighbor model

$t\sum_{j=N(i)}\Psi_{j}=E\Psi_{i},$ (220)

where the site $i$ sits on the hexagonal lattice of graphene and the sum over
$j=N(i)$ stands for the three neighbors $j$ of site $i$. The value $t=2.7$ eV
is the hopping integral between two neighboring $p_{z}$ orbitals. Recursive
Green’s Function algorithms were adapted from the square lattice to the
hexagonal case of graphene [162]. They were used very early to check the
prediction that the conductivity of pristine graphene at the Dirac point was
$4e^{2}/(\pi h)$ while its Fano factor had the same value $1/3$ as in a
diffusive wire [283] due to transport through evanescent states in the limit
of wide samples. These predictions were later observed experimentally [67].
Recursive Green’s function calculations were also used extensively to study
the effect of disorder which makes the conductivity dependent on the disorder
strength and correlation length [161] or to assess the regime of validity of,
e.g., semi-classical Boltzmann approaches [148].

####  IX.3.2 Massless Dirac equation

At low energy, pristine graphene is described by 2 independent relativistic
Weyl (massless Dirac) equations

$H_{\mathrm{D}}=v_{F}\begin{pmatrix}0&p_{x}-ip_{y}\\\
p_{x}+ip_{y}&0\end{pmatrix}$ (221)

corresponding to two different valleys K and K’, where the $2\times 2$
structure arises from graphene having two atoms per unit cell. Relativistic
equations of this type play an important role in modern condensed matter
physics as they arise in a number of topological materials, either in the bulk
or for the description of their surface states [303].

Attempts to simulate only one of the two valleys (for instance in order to
distinguish between intervalley and intravalley scattering) with any
discretization scheme faces the so-called “fermion doubling” problem: in any
tight-binding model, Dirac points always come in pairs [273, 200].

This problem can be circumvented by borrowing solutions from the lattice gauge
theory [263], such as by evaluating finite differences on a lattice that is
displaced symmetrically from the original lattice [282] or by adding a small
term which breaks the time-reversal symmetry [116, 107] or $k$-space
discretization [27]. See [213] for a recent discussion of various ways of
addressing Fermion doubling. Similar approaches can be used for other Dirac-
like systems such as 3D Weyl semi-metals [252]

###  IX.4 Quantum Spin Hall

Topological materials are attracting a lot of attention and have prompted many
numerical studies of their transport properties. Reviews of these materials
can be found in [110, 230, 23]. Computational methods have been useful in
particular for 2D models where large systems can be studied.

####  IX.4.1 Kane-Mele model

One of the most important models of topological insulators is the Kane-Mele
model that describes a two-dimensional quantum spin Hall insulator [125, 126].
This model corresponds to a single orbital model for graphene Eq. (220) to
which one adds a spin-orbit coupling term of strength $\lambda$. The tight-
binding equation reads

$\displaystyle
t\sum_{j=N(i)}\Psi_{j}+i\lambda\sum_{j=N^{\prime\prime}(i)}\nu_{ij}\sigma_{z}\Psi_{j}=E\Psi_{i},$
(222)

where $\Psi_{i}=(\Psi_{i\uparrow},\Psi_{i\downarrow})$ is a two component
spinor, $\sigma_{z}$ a Pauli matrix, $N(i)$ [resp. $N^{\prime\prime}(i)$]
stands for the [resp. second] nearest neighbors of $i$, and $\nu_{ij}=-1$
(resp. $\nu_{ij}=1$) for clockwise (resp. counter-clockwise) hoppings. The
spin-orbit term opens an energy gap in the spectrum. Inside this gap are
topologically protected chiral states that appear at the edge the system.

In actual pristine graphene, the value of $\lambda$ is expected to be very
small in the $10\ \mu$eV range [139]. Various strategies have been designed to
engineer higher values of $\lambda$ including depositing heavy adatoms on top
of graphene [301, 138]. Simulations of devices based on quantum spin Hall
effect have shown interesting behavior or thermoelectric effects [258, 58] or
electric switches [259].

####  IX.4.2 Bernevig-Hughes-Zhang model

A second important model that leads to a quantum spin Hall phase is the
Bernevig-Hughes-Zhang (BHZ) model [32]. The BHZ Hamiltonian is a four band
model that reads

$\displaystyle
H=C+M\sigma_{0}\tau_{0}-B(p_{x}^{2}+p_{y}^{2})\sigma_{0}\tau_{z}$
$\displaystyle-D(p_{x}^{2}+p_{y}^{2})\sigma_{0}\tau_{0}+A(p_{x}\sigma_{z}\tau_{x}-p_{y}\sigma_{z}\tau_{y}),$
(223)

where $\sigma_{a}$ and $\tau_{a}$, $a\in\\{0,x,y,z\\}$ are Pauli matrices (or
the identity matrix for $a=0$) that act, respectively, on the spin degree of
freedom or an extra orbital degree of freedom. The interest for the BHZ model
stems from the fact that it arises naturally as an effective 2D model for a
quantum well made in inverted semi-conductor heterostructures such as
HgTe/HgCdTe. The parameters $A,B,C,D$ and $M$ are material-specific and depend
on heterostructure geometry parameters such as well thickness. Discretized
versions of the BHZ model are a natural playground for numerical transport
calculations in devices [5]. In particular, the presence of disorder has been
shown to stabilize a “topological Anderson insulator” phase [163].

###  IX.5 Mesoscopic superconductivity

It is natural to apply computational quantum transport methods to Mesoscopic
superconductivity problems and this has been done in this field since the
early days. A natural framework to do so is the mean field Bogoliubov De
Gennes (BdG) method [98] that extends an initial normal (non-superconducting)
Hamiltonian $H$ to a larger Hilbert space with a $2\times 2$ Nambu block
structure with an “electron” and a “hole” sector:

$H_{\mathrm{BdG}}=\begin{pmatrix}H&\Delta\\\
\Delta^{\dagger}&-H^{*}\end{pmatrix}.$ (224)

The underlying normal Hamiltonian $H$ can be either in the continuum or
already discretized into some tight-binding model Hamiltonian [315] so that
the technique can be applied to a wide class of devices where some part
becomes superconducting. The anti-symmetric matrix $\Delta$ accounts for the
superconducting pairing. For conventional $s$-wave superconductors, $\Delta$
is a diagonal matrix $\Delta_{nm}=\Delta_{n}\delta_{nm}$ where $|\Delta_{n}|$
is the local superconducting gap on site $n$. Using various forms of off-
diagonal $\Delta$, it is also possible to account for more exotic pairings
such as $p$-wave or $d$-wave superconductivity [21]. In principle $\Delta$
must be calculated self-consistently from the mean-field treatment of the
effective (attractive) electron-electron interaction. However, in some
situations (such as a bulk superconductor in contact with a small mesoscopic
region) it is sufficient to consider $\Delta_{n}$ to be equal to its bulk
value inside the superconductor and zero in the normal part. This
approximation is often used in practice.

####  IX.5.1 Andreev reflection

In the presence of a single superconductor (or several at the same
electrochemical potential) and normal electrodes, the superconductor acts as
an “Andreev mirror” that reflects electrons into holes. The conductance of
such a system can be obtained through a simple extension of the Landauer
formula [152]. In the simple case where the superconducting gap is much
smaller than the Fermi energy and one considers voltage biases smaller than
the gap, it is sufficient to study solely the scattering matrix of the normal
part of the system and infer the properties of the full normal-superconducting
system from it [187, 30]. For example, [30] showed that Andreev reflection
gives rise to an enhancement of Andreev back-scattering even in disordered
systems, as shown in Fig. 12. In more general situations, one simulates the
full BdG Hamiltonian [111, 91] with a fixed gap or even the full self-
consistent problem with models inherited from first principle calculations
[296]. Recent applications include the search for Majorana bound states in
semiconducting nanowires proximitized with a superconductor [4, 149].

Figure 12: Numerical simulation of a $300\times 300$ tight-binding model for a
disordered normal metal, in series with a superconductor (inset). The
histograms give the modal distribution for reflection of an electron, with
mode number 1 indicating normal incidence. The top two panels give the
distribution of reflected holes (in the absence and presence of a magnetic
field $B$), the bottom panel of reflected electrons (for $B=0$). The arrow
indicates the ensemble-averaged height of the giant backscattering peak for
Andreev reflection, predicted from theory. From [30].

.

####  IX.5.2 Supercurrents

Another important question that can be considered is the current-phase
relation for the supercurrent that can flow in-between two superconductors
[101]. This is not strictly speaking quantum transport since the supercurrent
actually flows at equilibrium.

Numerical calculations of supercurrents are particularly useful for the self-
consistent problem which is usually difficult to address analytically. There
are many examples of such calculations including [160, 202, 34]

An interesting aspect of superconducting devices is that when the different
superconducting electrodes do not share the same electrochemical potential,
the problem becomes effectively time-dependent even though the initial
formulation is not. For instance, a DC voltage $V_{b}$ across a
superconducting-normal-superconducting junction gives rise to an AC current at
frequency $2eV_{b}/h$ – this is the AC Josephson effect. To simulate such
problems, one need to leave the DC framework of this article and uses time-
dependent techniques [66, 222, 267, 305].

####  IX.5.3 Search for Majorana bound states in nanowires

Majorana bound states (MBSs) are the simplest example of non-Abelian anyons
[198]. They have been predicted to be realized for instance in one-dimensional
hybrid superconducting/semiconducting nanowire. The simplest model describing
such a system is the Bogoliubov-de Gennes Hamiltonian

$H=\left(\frac{p_{x}^{2}}{2m}-\mu+V(x)+\frac{\alpha_{R}}{\hbar}p_{x}\sigma_{y}\right)\tau_{z}+E_{Z}\sigma_{x}\tau_{0}+\Delta\tau_{x},$
(225)

where $\sigma_{i}$ ($\tau_{i}$) are Pauli matrices in spin (electron-hole)
space, $i=0,x,y,z$. $p_{x}$ is the momentum along the nanowire, $m$ the
effective mass of the semiconductor, $\mu$ the global chemical potential, and
$V(x)$ a potential induced for example by external gates, allowing to change
the electron density in different parts of the wire and defining a tunnel
barrier. $\alpha_{R}$ is the strength of Rashba spin-orbit coupling in the
semiconductor, $E_{Z}$ the Zeeman splitting induced by an external magnetic
field along the nanowire, and $\Delta$ the induced superconducting gap. This
model exhibits MBSs when $E_{Z}^{2}>(\mu-V_{\text{wire}})^{2}+\Delta^{2}$, in
a region with potential $V_{\text{wire}}$ [178, 207]. Arising through
topological phase transition, the appearance of a MBS is always accompanied by
a gap closing at the transition point.

Many transport simulations rely on this simple model, but it has also been
extended to higher spatial dimensions, including additional physical effects
such as the orbital effects of a magnetic field [201], a more realistic
modeling of the proximity effect [17, 235], and including all these effects
together with electrostatic simulations [308].

Experimental and theoretical work in this field is numerous, and previous
experimental results are nowadays discussed controversially. Here we only
highlight selected work employing numerical transport simulations that shaped
the understanding in this field: MBSs are states at zero-energy, and give rise
to a zero-bias peak in transport. However, the particular appearance may
depend on device details. Moreover, zero-bias peaks are not specific to MBSs.
Numerical transport simulations have been crucial in shaping the understanding
of these aspects.

Early work focused on studying the appearance of the zero-bias peak in
different device geometries [228, 231]. It was also shown that the gap closing
accompanying the appearance of MBS may not be visible in transport
measurements [264, 225].

Numerical transport simulations have been essential in shaping the
understanding that zero-bias peaks are not a smoking-gun signature of MBSs.
Early work showed that disorder can lead to trivial low-lying Andreev bound
states giving rise to zero-bias peaks in the conductance [172, 226, 231]. Fig.
13 shows one example of such zero-bias peaks sticking to zero energy even in a
trivial system. Later, it was shown that even in the absence of disorder
particular Andreev bound states can appear very similar to MBS [171],
mimicking for example quantized conductance signatures [194], that previously
were considered a smoking gun. Even more, it was shown that any local
measurement signal associated with MBSs can be mimicked by trivial states
[292]. For this reason, new measurement paradigms have been proposed, in
particular non-local measurements to detect the topological phase transition
[239, 68].

Figure 13: Numerical simulation of a two-dimensional version of Eq.(225) for a
nanowire including a tunnel barrier and disorder. The differential conductance
is shown for a single disorder realization as a function of bias voltage $V$
and Zeeman splitting $E_{Z}$ due to a parallel magnetic field. The system is
in the trivial regime, but still exhibits a conductance peak pinned to zero
voltage (green circles) that could be mistaken for a Majorana bound state. The
two panels correspond to two parameter regimes. From [226].

###  IX.6 Spintronics

Quantum transport numerics is also used widely in spintronics. The techniques
can be used for mainstream metallic multilayers [37] but have been
particularly popular in semiconducting spintronics [81].

A typical Hamiltonian considered in model calculations includes a Rashba spin-
orbit coupling term of strength $\alpha$ and ferromagnetic regions with an
exchange coupling $J$ and a magnetization oriented around $\vec{m}$:

$\displaystyle
H=\frac{1}{2M}(p_{x}^{2}+p_{y}^{2})+\alpha(p_{y}\sigma_{x}-p_{x}\sigma_{y})$
$\displaystyle+J(\vec{r})\vec{m}\cdot\vec{\sigma}.$ (226)

This type of models can be used to study a variety of phenomena such as the
spin Hall effect for electrons [204] or holes [109].

Spintronics studies are not limited to 2D but also include various 3D
calculations such as, e.g., 3D topological insulators [59].

###  IX.7 Self-consistent quantum-electrostatic simulations of nanoelectronics

The model calculations discussed previously assume a known Hamiltonian, i.e.,
that the electric potential seen by the conducting electrons is an input of
the model. This level of modeling is sufficient in many situations but suffers
from a number of limitations: the chosen potential may be unrealistic, the
relations between the electric potential and, e.g., the voltage applied to
electrostatic gate may be unknown, or they may ignore important effects such
as edge state reconstruction or non-linear corrections to the conductance. The
next level in modeling thus consists in calculating the electric potential
self-consistently by solving the quantum problem (the focus of this review)
together with the Poisson equation. Below, we review the main aspects of the
self-consistent problem.

####  IX.7.1 Formulation of the self-consistent problem

The self-consistent problem is formulated as a set of three equations that
account, respectively, for quantum mechanics, statistical physics, and
electrostatics. We denote the electric potential as $U$ and for simplicity
suppose that the tight-binding model is formulated in real space (otherwise
$U$ has to be replaced by the corresponding matrix elements in the quantum
problem). We write $H_{0}$ for the tight-binding Hamiltonian without electric
potential and $\Psi_{\alpha E}$ for the scattering state in mode $\alpha$ and
energy $E$ (ignoring for simplicity possible bound states present in the
system). The problem reads

$\displaystyle\left[H_{0}-eU\right]\Psi_{\alpha E}=E\Psi_{\alpha E},$ (227)
$\displaystyle n=\int\frac{dE}{2\pi}\sum_{\alpha}|\Psi_{\alpha
E}|^{2}f_{\alpha}(E),$ (228)
$\displaystyle\vec{\nabla}\cdot(\epsilon\vec{\nabla}U)=en+en_{0},$ (229)

where Eq. (227) is the scattering problem formulated at the beginning of this
review. Eq. (228) relates the electronic density $n$ to the local density of
states $|\Psi_{\alpha E}|^{2}/(2\pi)$ in a non-equilibrium situation where
different leads may have different chemical potentials or temperatures. Eq.
(229) is the Poisson equation in a (possibly spatially dependent) dielectric
constant $\epsilon$, $e$ is the electron charge, and $n_{0}$ an additional
electronic density due to, e.g., dopants.

The self-consistent quantum-electrostatic problem in the presence of
electrodes is a non-linear integro-differential problem. Its solution is
usually obtained by some sort of iterative scheme. One starts with an initial
electronic density, calculates the potential from the Poisson equation, then
solves the quantum problem to recalculate the density. This sequence is
iterated until convergence. Many techniques have been developed to obtain
better convergence properties for such schemes. We refer to [19] for an
entrance point to the corresponding literature. Besides the convergence of the
self-consistent loop, the calculation of the density itself, Eq. (228) can be
computationally intensive. To accelerate the computation, it is useful to
separate the density into the sum an equilibrium contribution and an out-of-
equilibrium one. The calculation of the equilibrium density can be accelerated
using integration in the complex plane [211, 128] while the integration of the
out-of-equilibrium contribution spans a small range of energies close to the
Fermi energy [251].

The importance of taking into account electrostatics can be illustrated by the
non-linear part of the current-voltage characteristics $I(V)$ of a device
[61]. Let us write $g_{2}$ for the second-order contribution such that
$I(V)=gV+g_{2}V^{2}+\mathcal{O}(V^{3})$. The non-interacting Landauer theory
then gives $g_{2}=\frac{1}{2}\frac{\partial g}{\partial E_{F}}$ with
$g_{2}(B)=g_{2}(-B)$ in presence of magnetic field. When taking electrostatics
into account, one needs to calculate the _emissivity_ $\frac{\partial
n_{i}}{\partial V}$ (change of density due to the bias voltage, also known as
Landauer dipole), then compute the change of potential $\frac{\partial
U_{i}}{\partial n_{j}}$ due to the change of density $n_{j}$ (through Poisson
equation), then finally compute the _injectivity_ , i.e., the change of
conductance $\frac{\partial g}{\partial U_{i}}$ due to the change of
potential. One arrives at

$g_{2}=\frac{1}{2}\frac{\partial g}{\partial E_{F}}+\sum_{ij}\frac{\partial
g}{\partial U_{i}}\frac{\partial U_{i}}{\partial n_{j}}\frac{\partial
n_{j}}{\partial V},$ (230)

where the second terms is entirely due to electron-electron interactions [16,
227, 113]. In particular the antisymmetric part $g_{2}(B)-g_{2}(-B)$ vanishes
in the absence of this term.

####  IX.7.2 Examples of applications

Quantum transport within the self-consistent quantum electrostatic model has
been studied in many situations. (Self-consistent calculations within the
density functional theory framework will be discussed in the next section.)
They are particularly important for the quantum Hall effect as the
electrostatic interaction leads to a reconstruction of the edge states into
compressible and incompressible stripes [247, 19, 20], as shown in Fig. 14.
The technique has also been applied to address silicon devices including
quantum dots [94] and transistors [134, 205, 150]. Other calculations have
focused on semiconductor nanowires used in the search for Majorana bound
states [291, 17, 310].

Figure 14: Numerical calculation of the electron density (top) and band
structure (bottom) for a top gate-defined wire in a two-dimensional electron
gas in the quantum Hall regime. The device is invariant by translation along
$y$, the electrons are confined by two gates along $x$ with the magnetic field
applied along $z$. Blue lines correspond to a self-consistent Thomas-Fermi
calculation, orange lines to a full solution of the Poisson-Schrödinger
problem. Gray regions indicate the compressible stripes while the white
regions are incompressible. Adapted from [20].

##  X Selected applications to realistic systems

In this section, we turn to a few applications of computational quantum
transport for “realistic” models, i.e., the DFT, FP-TBM and SE-TBM discussed
in Fig. IX.1. Once again, we do not aim at being exhaustive, the corresponding
literature being very large, but just point to a few illustrative cases.

###  X.1 DFT models

First-principles density functional theory (DFT) can provide predictive
calculations in the absence of any prior experimental data. The combined
approach is often referred to as DFT + NEGF. There exist many reviews and
pedagogical materials that explain DFT in great detail. We refer the
interested reader to the articles that accompany DFT software packages, a list
of which can be found towards the end of Sec. I.2.

####  X.1.1 Non-linear conductance in short junctions

We start with a use case where density functional theory (DFT) is used to
compute the non-linear current voltage characteristics of nanodevices. As we
have seen in the preceding section, it is crucial to treat the electron-
electron interaction at least at the mean field level to calculate non-linear
corrections to the conductance. DFT calculations do that by design and have
been used in many situations [251].

Early calculations and later benchmarks focused on short wires made of a few
atoms, e.g., C, Al or Au [154, 42, 214, 238]. The calculations soon extended
to molecular nanoelectronics [76, 131, 82, 86, 246]. Despite the sensitivity
of the conductance to the nature of the molecule-electrode contact, it was
found in [75] that the qualitative shape of the experimental $I-V$
characteristics could be recovered. Another topic that attracted a strong
interest is carbon nanotubes [275] as well as graphene-based devices with
either a transistor-like geometry [212, 18, 216] or multi-terminal cross bar
geometry [38, 106, 243]. Possible applications of these devices include DNA
sensing using nanopores in graphene [245]. Recent works have extended the
studies to other 2D materials such as MoS2 [95]

Despite their widespread use, it should be noted that these DFT techniques
have serious limitations for the study of small junctions. Perhaps the most
important one is that these techniques do not account for the ubiquitous
phenomena of Coulomb blockade. Another limitation is that most studies ignore
phonons and more generally atomic displacements that may be triggered by the
out-of-equilibrium situation. In the presence of bound states, an additional
prescription must be given to determine their occupations [164]. In some
situations, it is important to properly account for the charge redistribution
rather deeply inside the electrodes [192].

####  X.1.2 Applications to spintronics

A very successful example of the DFT + NEGF approach is the theoretical
prediction [46, 190] of the existence of a large tunneling magneto-resistance
of a Fe/MgO/Fe magnetic tunnel junction. There, the coherent spin-polarized
tunneling and the symmetry of electron wave functions plays an essential role
that can only be captured by quantum transport calculations. This prediction
has motivated large experimental efforts [313] to fabricate such devices.
These magnetic tunnel junctions now play a central role in both basic [298]
and applied [173] spintronics research.

Figure 15: (a) Schematic view of a magnetic tunnel junction where 6 monolayers
of MgO(001) insulating barrier are sandwiched between two semi-infinite
Fe(001) leads. The junction is infinite in the transverse direction, so that
the depicted supercell is periodically repeated in the direction transverse to
current flow from the left to the right Fe lead. (b) Electronic band structure
of a periodic Fe/MgO/Fe/MgO…Fe/MgO lattice obtained by DFT calculations using
LCAO basis (solid line) versus those using full potential linearized augmented
plane-wave basis (circles). Panel (b) is adapted from [295].

Fig. 15a shows a schematic of a Fe-MgO-Fe junction simulated in [295]. The
system is supposed to be invariant by translation along the direction
perpendicular to the junction (no disorder) so that the transmission is
calculated as a function of the transverse momentum which must be integrated
upon at the end of the calculation. The Hamiltonian of the scattering region
was written in a linear combination of atomic orbitals (LCAO). An accurate
calculations requires as many as 10 000 orbitals for this calculation while
the supercell only contains a total of 24 atoms (Fe Mg and O). To control this
accuracy, one compares the LCAO calculation with a more accurate one done in a
plane wave basis, as illustrated in Fig. 15b. In addition to checking the
convergence of the calculation with the number of LCAO orbitals, one must also
obtain convergence of the momentum integration which requires $\sim 10^{4}$
momentum points at zero bias (only $\sim 10^{2}$ for the band structure, up to
$\sim 10^{5}$ at finite bias voltage, up to $\sim 10^{7}$ for spin-torque
[299]).

Later works tried to optimize the magneto-resistance of the tunnel junction.
It has been predicted, for instance, that a few layers of graphene with Ni
electrodes act as a perfect spin filter [127] with regimes showing a negative
differential resistance [244] at finite bias. Another line of work used
DFT+NEGF to address metallic magnetic systems [254, 269]. An important
evolution of the field has been the discovery of the phenomena of spin-
transfer torque in non-collinear systems. Spin torque has enabled fully
electric control of a magnetic configuration [232]. DFT studies of spin
transfer torque include [269, 108, 299, 79]. Recent studies also capture a
related effect, the spin orbit torque [203].

###  X.2 First-principles tight-binding models (FP-TBM)

The direct DFT approach to quantum transport is very computationally
intensive. These calculations are usually restricted to one-dimensional
systems or quasi-one-dimensional systems that are invariant by translation
along the perpendicular direction. For more complex geometries or in the
presence of disorder [132, 312] one needs to develop effective models that
reduce the complexity. The FP-TBM approach projects the DFT Hamiltonian onto a
smaller basis set, typically atomic or pseudo-atomic orbitals. FP-TBM are
designed to describe the physics in a finite energy window centered around the
Fermi level. We distinguish FP-TBM (where there exists a systematic procedure
to construct the basis set) from the semi-empirical SE-TBM discussed in the
next subsection, where the construction is less systematic. The systematic
approach of FP-TBM methods is particularly important when several bands get
hybridized at the Fermi level.

Perhaps the most well-known FP-TBM basis are the _maximally localized Wannier
functions_ which can be constructed out of Bloch wave functions [188]. These
orbitals are rather compact, retain the original symmetries of the bands and
reproduce well the original DFT results. A popular approach uses a plane wave
code such as VASP [143] followed by a projection on Wannier orbitals using,
e.g., the code Wannier90 [196]. An illustrative example is the case of
graphene. Using a minimal basis set of a single $p_{z}$-like orbital per site,
[84] find that including hoppings up to the 8th neighbor (i.e., 400
parameters) is required to reproduce the bands predicted by DFT. This can be
contrasted with the single parameter model used in the preceding section.
Examples of quantum transport calculations performed with Wannier functions
include the calculation of the conductance of disordered nanowires [53, 257]
and of single molecule junctions [279, 272].

An alternative to Wannier function uses the linear-muffin-tin-orbitals (LMTO)
and its generalizations [13]. Other FP-TBM include the quasiminimal basis
orbitals [297] or predetermined basis of fixed localized wave functions [2,
3]. All these techniques require particular care in the treatment of
boundaries or interfaces if the model is to be transferred to other geometries
or structures.

###  X.3 Semi-empirical tight-binding models (SE-TBM)

Continuing with the hierarchy of Hamiltonians, we now consider semi-empirical
tight-binding models (SE-TBM). In this approach a set of relevant orbitals and
couplings is chosen a priori. The model is then fitted to ab-initio data
(e.g., often to reproduce the DFT band structure) and/or to reproduce
experimental data. SE-TBM models can have up to hundreds of parameters. The
simplest SE-TBM models, with very few parameters, essentially match the model
Hamiltonians discussed in the preceding section. The SE-TBM models are also
closely connected to the $\mathbf{k}\cdot\mathbf{p}$ models discussed in the
next subsection. The main difficulty in preparing SE-TBM for quantum transport
calculations is the same as the one discussed in the preceding section for FP-
TBM models: a proper treatment of the boundaries of the device and the
interface between different materials. For instance, the description of the
interface between graphene and metallic contacts can require as many as
hundred parameters for a faithful description [215, 28]. The models are
usually fitted with classical least mean square optimization methods, but
recently machine learning approaches have started to be used [112].

SE-TBM models are extremely popular and we cannot do justice to the literature
in this review. They have been applied to MoS2 devices [237], adatoms on
graphene [301], graphene-hBN devices [185], quantum dots [177], 3D topological
insulators Bi2Se3, Bi2Te3, Sb2Te3 [170], and many more. Let us also mention
the possibility to perform multiscale calculations where critical parts of the
system are described at the DFT level while other parts use SE-TBM models
[52].

####  X.3.1 Semiconductor devices

SE-TBM models for semiconductors go back to the eighties with minimum next
nearest-neighbor models accounting for the s and p orbitals [289]. More
accurate models also include the d bands: [54, 39]. Assuming that the hopping
matrix elements only depend on the distance between atoms [215], these models
can be parametrized as a function of atomic distance and then used for an
arbitrary atomic configuration. This is very useful for semiconductor devices
where stress can have a significant effect on transport properties [205, 80].
Another appealing aspect of these models is that they can correct typical
deficiencies of DFT approaches, e.g., the inability to correctly predict the
semi-conductor gap, by fitting the model to experimental data.

####  X.3.2 Carbon based nanodevices

SE-TBM models are also extremely popular for modeling graphene nanoribons or
carbon nanotube devices. It was recognized early on that the simple nearest-
neighbor TBM for a single $p_{z}$ orbital Eq. (220) is often inaccurate: it
predicts perfectly flat bands in zigzag nanoribons as well as a metallic
(gapless) band structure for armchair nanoribbons that are $3N+2$ carbon atoms
wide. A more accurate model uses second- and third-nearest-neighbor hoppings
[236]. Equivalently, one may use nearest-neighbor hoppings with more bands
[65, 41]. These additional bands are also important when describing the spin-
orbit interaction [139]. Again, the transport properties across these carbon-
based devices are very sensitive to the accuracy of model and in particular
the interfaces. For instance, SE-TBM modeling [180] that did not properly
address the interface between carbon nanotubes and the metallic Ohmic contact
could not predict the experimental findings quantitatively [89]. Another
example is the off current of graphene nanoribbon transistors which varies by
orders of magnitude upon moving from the simple single-orbital $p_{z}$-SETBH
model to a more elaborate multi-orbital $p/d$-SETBH model [41]. Similar
conclusions were drawn for the validity of simulations of non-local
resistances in graphene Hall bars [185].

###  X.4 Discretized $\mathbf{k}\cdot\mathbf{p}$ Hamiltonians

We end our hierarchy of realistic models with $\mathbf{k}\cdot\mathbf{p}$
Hamiltonians. The $\mathbf{k}\cdot\mathbf{p}$ Hamiltonians are low-energy
effective theories valid in the vicinity of some high-symmetry $k$-point. They
describe the physics at large scale, not at the atomic scale, which can be
advantageous for the simulation of large systems. There exists a large
literature on the $\mathbf{k}\cdot\mathbf{p}$ method [290]. [140] describes
the particular case of two-dimensional transition metal dichalcogenide
semiconductors while [184] focuses on graphene. $\mathbf{k}\cdot\mathbf{p}$
Hamiltonians can account for magnetic field as well as spin-orbit effects
[309].

In the simplest cases, $\mathbf{k}\cdot\mathbf{p}$ models often reduce to the
model Hamiltonians discussed in the preceding section. For instance, the
effective-mass approximation can be thought as arising from developing the
dispersion relation around the $\Gamma$ point of the conductance band of a
III-V semiconductor. In the case of graphene, the simple tight-binding model
Eq. (220) can be expanded around the two points $K$ and $K^{\prime}$ of the
Brillouin zone which leads to two copies of the 2D massless Dirac equation Eq.
(221). More complex examples arise in the description of semi-conductors and
in particular to account for the $p$-character of their valence band. For
instance, the $4\times 4$ Luttinger Hamiltonian [179],

$H=\frac{1}{2m}\left[\left(\gamma_{1}+\frac{5}{2}\gamma_{2}\right)\hat{\mathbf{p}}^{2}-2\gamma_{2}\left(\hat{\mathbf{p}}\cdot\hat{\mathbf{S}}\right)^{2}\right],$
(231)

describes the heavy and light hole bands of semi-conductors. Here $m$ is the
electron mass and $\hat{\mathbf{S}}=(S_{x},S_{y},S_{z})$ is a vector of
spin-3/2 operators. The dimensionless parameters $\gamma_{1}$ and $\gamma_{2}$
are material specific. Prior to quantum transport calculations, these models
must be discretized using, e.g., finite differences. When discretizing these
models, particular care needs to be taken to avoid spurious solutions, e.g.,
in $8\times 8$ $\mathbf{k}\cdot\mathbf{p}$ semiconductor models [88].
Additionally, Dirac materials are not amenable to regular discretization due
to the occurrence of the fermion doubling problem as explained in Sec. IX.3.2.
$\mathbf{k}\cdot\mathbf{p}$ models are often used for quantum transport
calculations such as [311, 260]. As for other models, special care must be
given to the treatment of boundaries and interfaces [45] so that they are
error prone when used for confined structures such as quantum dots [90] or
thin films [199].

##  XI Final remarks

###  XI.1 Simulation workflow

In practice, numerical simulations of quantum transport devices are often
carried out by approaching the complete problem step by step. One starts by
studying the bulk Hamiltonian(s) for the material(s) under consideration.
Written in $k$-space, these Hamiltonians are readily diagonalized and their
spectrums provides useful information such as the positions of gaps or Dirac
points.

In a second step, one constructs electrodes of finite width and analyzes them
in isolation. The associated $k$-space Hamiltonians provide information such
as the number of propagating modes in each electrode. Packages such as
KwantSpectrum [137] provide tools to facilitate this analysis.

It is only in a third step that one proceeds to treating scattering problems
that consist of a scattering region with attached electrodes. It it often
useful to start out with subparts of the final device. For instance, in an
electronic interferometer that contains several quantum point contacts, one
could study the transport properties of a system that contains a single
quantum point contact in order to calibrate the corresponding model. Another
example is the interface between two materials that could be studied
separately before being integrated into a larger device. When the final device
is very large, it can actually be computationally advantageous to stop at this
level and concatenate the scattering matrices of the different pieces
numerically (following, e.g., [241]). In most cases, however, one proceeds, at
last, to simulate the entire device directly.

While in former times research groups relied on in-house codes, nowadays in
most cases one will want to take advantage of already existing solutions, at
least as a starting point or foundation. Out of the codes that are listed at
the end of Sec. I.2, the following three target quantum transport specifically
and are under active maintenance: nextnano [33], Kwant [104], and Quantica.jl
[249].

###  XI.2 Outlook

Numerical calculations of quantum transport have become an ubiquitous tool in
many fields, from electrical engineering to materials science and quantum
nanoelectronics. In this review, we focused on providing a comprehensive
presentation of the two equivalent formalisms used in this field, the
scattering matrix approach and the Green’s function approach, with an emphasis
on explicitly providing the links between the two. We also provided an in-
depth review of the different algorithms that have been developed to perform
actual calculations.

Many topics fall beyond the scope of this review. One is the simulation of
time-dependent quantum transport [265, 145, 137], see [97] for a review, and
its special case of finite-frequency transport or Floquet scattering [195].
Another major aspect that we did not touch is the effect of electron-electron
interactions or the interaction of electrons with other degrees of freedom
such as phonons or magnons. It is indeed a very frustrating fact that even the
simplest and most common of these effects, for instance Coulomb blockade,
cannot be easily incorporated into the type of approach presented in this
review. The effects of decoherence are also beyond the scope of this review.

In terms of applications, the range is so vast that we could only point to a
few examples. Our choice of topics was necessarily subjective and biased by
our own interests.

With regard to future work, much remains to be done besides addressing the
above points. One important aspect is improving the accuracy of predictions
compared to experiments [60]. This requires better models, efficient tools to
address electrostatics, solvers capable to handle large systems (in particular
in three dimensions), and close collaboration with experimental teams with the
goal of producing systematic data sets dedicated to the calibration of models.

## Author contributions

AA, MW, and XW developed the formalism. AA and MW developed the stable
algorithm for solving the scattering problem. MI described the bound state
problem. TOR and DV implemented and described symmetry analysis of the modes.
BKN wrote the section on applications to realistic systems. AA, MW, CG, and XW
wrote the rest of the manuscript. XW organized the writing. All authors
reviewed the text.

###### Acknowledgements.

Many colleagues have given us important help in putting this review together.
We are particularly grateful to Piet Brouwer, Rodolofo Jalabert, Colin
Lambert, Caio Lewenkopf, Horacio Patawski, Klaus Richter, Stefano Sanvito,
Henning Schomerus, and Dietmar Weinmann.

##  Appendix A Derivation of Feynman-Hellman equation for the velocity

The Feynman-Hellman equation relates the velocity of a mode to the average of
the velocity operator. Considering normalized modes,

$\phi^{\dagger}\phi=1$ (232)

we have

$\frac{\partial\phi^{\dagger}}{\partial
k}\phi+\phi^{\dagger}\frac{\partial\phi}{\partial k}=0.$ (233)

From $E(k)=\phi^{\dagger}H(k)\phi$ we get,

$\displaystyle v$ $\displaystyle\equiv$ $\displaystyle\frac{\partial
E(k)}{\partial k}$ (234) $\displaystyle=$
$\displaystyle\frac{\partial\phi^{\dagger}}{\partial
k}H(k)\phi+\phi^{\dagger}H(k)\frac{\partial\phi}{\partial
k}+\phi^{\dagger}\frac{\partial H(k)}{\partial k}\phi$ $\displaystyle=$
$\displaystyle E(k)\left[\frac{\partial\phi^{\dagger}}{\partial
k}\phi+\phi^{\dagger}\frac{\partial\phi}{\partial
k}\right]+\phi^{\dagger}\frac{\partial H(k)}{\partial k}\phi$ $\displaystyle=$
$\displaystyle\phi^{\dagger}\frac{\partial H(k)}{\partial k}\phi$
$\displaystyle=$ $\displaystyle i\phi^{\dagger}\left[\lambda
V^{\dagger}-\lambda^{-1}V\right]\phi$

which is the Feynman-Hellman equation. Note that inside the matrices $\Phi$,
the modes $\phi$ are renormalized to carry a unit velocity $v=\pm 1$.

##  Appendix B Properties of the leads modes.

In this appendix, we establish a few technical results on the lead propagating
and evanescent modes that are needed in various proofs of the main text. We
begin by multiplying Eq. (37) by $\Phi^{\dagger}$:

$\Phi^{\dagger}V\Phi+\Phi^{\dagger}(H-E)\Phi\Lambda+\Phi^{\dagger}V^{\dagger}\Phi(\Lambda)^{2}=0.$
(235)

The complex conjugate of the above equation reads

$(\Lambda^{*})^{2}\Phi^{\dagger}V\Phi+\Lambda^{*}\Phi^{\dagger}(H-E)\Phi+\Phi^{\dagger}V^{\dagger}\Phi=0.$
(236)

Now, multiplying Eq. (235) by $\Lambda^{*}$ on the left and Eq. (236) by
$\Lambda$ on the right, we arrive after subtracting one equation from the
other at

$[\lambda_{\alpha}^{*}-(\lambda_{\alpha}^{*})^{2}\lambda_{\beta}]\left.\Phi^{\dagger}V\Phi\right|_{\alpha\beta}=[-\lambda_{\alpha}^{*}(\lambda_{\beta})^{2}+\lambda_{\beta}]\left.\Phi^{\dagger}V^{\dagger}\Phi\right|_{\alpha\beta}.$
(237)

which is valid for all types of modes, evanescent and/or propagating.

###  B.1 Application to evanescent modes

If one of the two modes $\alpha$ or $\beta$ is evanescent, we can simplify Eq.
(237) by $1-\lambda_{\alpha}^{*}\lambda_{\beta}\neq 0$ and arrive at

$\lambda_{\alpha}^{*}\left.\Phi^{\dagger}V\Phi\right|_{\alpha\beta}=\lambda_{\beta}\left.\Phi^{\dagger}V^{\dagger}\Phi\right|_{\alpha\beta}$
(238)

which can be written in matrix form as,

$\Lambda_{\mathrm{e}+}^{*}\Phi_{\mathrm{e}+}^{\dagger}V\Phi_{\mathrm{e}+}=\Phi_{\mathrm{e}+}^{\dagger}V^{\dagger}\Phi_{\mathrm{e}+}\Lambda_{\mathrm{e}+}$
(239)

Similarly,

$\Lambda_{\mathrm{p}+}^{*}\Phi_{\mathrm{p}+}^{\dagger}V\Phi_{\mathrm{e}+}=\Phi_{\mathrm{p}+}^{\dagger}V^{\dagger}\Phi_{\mathrm{e}+}\Lambda_{\mathrm{e}+}$
(240)

and

$\Lambda_{\mathrm{e}+}^{*}\Phi_{\mathrm{e}+}^{\dagger}V\Phi_{\mathrm{p}+}=\Phi_{\mathrm{e}+}^{\dagger}V^{\dagger}\Phi_{\mathrm{p}+}\Lambda_{\mathrm{p}+}$
(241)

###  B.2 Application to propagating modes

For propagating modes $1-\lambda_{\alpha}^{*}\lambda_{\beta}$ vanishes
whenever $\lambda_{\alpha}=\lambda_{\beta}=\lambda$, i.e., the two modes
$\alpha$ and $\beta$ have the same momentum. In that situation, we can form
the velocity matrix

$i\Phi_{\mathrm{p}+}^{\dagger}\left[\lambda
V^{\dagger}-\lambda^{*}V\right]\Phi_{\mathrm{p}+}|_{\alpha\beta}$ (242)

where $\alpha$ and $\beta$ lie inside the degenerate space. This matrix is
Hermitian hence can be diagonalized. The resulting linear combination of
propagating modes (for which we abusively use the same notation), normalized
to carry unit velocities remain solutions of the mode equation and further
satisfy,

$i\Phi_{\mathrm{p}+}^{\dagger}\left[\lambda
V^{\dagger}-\lambda^{*}V\right]\Phi_{\mathrm{p}+}|_{\alpha\beta}=\delta_{\alpha\beta}$
(243)

from which we finally get,

$\Lambda_{\mathrm{p}+}^{*}\Phi_{\mathrm{p}+}^{\dagger}V\Phi_{\mathrm{p}+}=\Phi_{\mathrm{p}+}^{\dagger}V^{\dagger}\Phi_{\mathrm{p}+}\Lambda_{\mathrm{p}+}+i$
(244)

A corresponding expression can be derived for incoming modes. In that case,
the sign of the current is opposite and one arrives at,

$\Lambda_{\mathrm{p}-}^{*}\Phi_{\mathrm{p}-}^{\dagger}V\Phi_{\mathrm{p}-}=\Phi_{\mathrm{p}-}^{\dagger}V^{\dagger}\Phi_{\mathrm{p}-}\Lambda_{\mathrm{p}-}-i$
(245)

##  Appendix C Proof of invertibility of mode matrices

In this section, we discuss the relation between the rank of $\Phi_{\bar{t}+}$
and the absence of bound state in the corresponding semi-infinite lead.

Let us suppose that $\Phi_{\bar{t}+}$ is not full rank. It follows that
$V\Phi_{\bar{t}+}$ is not full rank neither, hence that there is a vector
$c\neq 0$ such that $V\Phi_{\bar{t}+}c=0$. On the other hand, we know that the
wave function

$\psi(j)=\Phi_{\bar{t}+}\Lambda_{\bar{t}+}^{j}c$ (246)

is an eigenstate of the _infinite_ lead, i.e., it satisfies the infinite lead
equation

$V\psi(j-1)+(H-E)\psi(j)+V^{\dagger}\psi(j+1)=0.$ (247)

In addition, we now have $V\phi(0)=0$, which means that $\psi(j)$ also
fulfills the boundary condition,

$(H-E)\psi(1)+V^{\dagger}\psi(2)=0.$ (248)

In other words, $\psi(j)$ is an eigenstate of the _semi-infinite_ lead.

Let us now decomposes the vector $c$ into its propagating and evanescent
sector $c=(c_{\mathrm{p+}},c_{\mathrm{\bar{e}+}})$. Using Eq.(239) and
Eq.(244), we get,

$c^{\dagger}\Lambda_{\bar{t}+}^{*}\Phi_{\bar{t}+}^{\dagger}V\Phi_{\bar{t}+}c=c^{\dagger}\Phi_{\bar{t}+}^{\dagger}V^{\dagger}\Phi_{\bar{t}+}\Lambda_{\bar{t}+}c+ic_{\mathrm{p+}}^{\dagger}c_{\mathrm{p+}}$
(249)

and since
$V\Phi_{\bar{t}+}c=c^{\dagger}\Phi_{\bar{t}+}^{\dagger}V^{\dagger}=0$, it
follows that $c_{\mathrm{p+}}^{\dagger}c_{\mathrm{p+}}=0$ hence,

$c_{\mathrm{p+}}=0$ (250)

Hence, $V\Phi_{\mathrm{t}+}c=0$ reduces to $V\Phi_{\mathrm{e}+}c_{\bar{e}+}=0$
which in turns means that $\psi(j)$ is a _bound state_ solution of the semi-
infinite lead. This last paragraph is a technical way to reflect the fact that
the probability current of an eigenstate of the semi-infinite lead must be
zero everywhere - this is current conservation - hence the propagating
contribution must vanish since we include only modes that propagate along one
direction.

Conversely, if there are no bound states at energy $E$, then
$V\Phi_{\bar{t}+}$ is full rank. Since $V=AB^{\dagger}$, $V\Phi_{\bar{t}+}$
being full rank implies that $B^{\dagger}\Phi_{\bar{t}+}$ is full rank and
since it is a square matrix, it must therefore be invertible.

##  Appendix D Scattering problem with non-orthogonal basis set

This review has focused on tight-binding models with an orthogonal basis set.
However, non-orthogonal basis sets appear naturally for basis sets used in
Slater-Koster tight-binding methods or density functional theory, as, e.g., in
[124, 209, 210, 215]. In this case, different basis states
$\left|n_{i}\right>$ have finite overlap, $\left\langle
n_{i}|n_{j}\right\rangle=W_{ij}$. Note that usually the overlap matrix is
denoted by the letter $S$ which is already taken by a central object of this
review, the scattering matrix. Hence we resort to $W$ to denote the overlap
matrix.

In this appendix we will argue how the main formulas presented in this review
can be generalized to tight-binding models with a non-orthogonal basis set
with a simple replacement rule. More explicit and extensive derivations can be
found in the literature, such as for the Green’s function formalism in [174]
and its non-equilibrium version in [151].

The Schrödinger equation for the scattering system with a non-orthogonal basis
reads

$\hat{H}_{\text{sys}}\hat{\psi}=E\hat{W}_{\text{sys}}\hat{\psi}$ (251)

with the overlap matrix

$W_{\text{sys}}=\begin{pmatrix}W_{\text{sr}}&P_{\text{sr}}^{T}W_{V}^{\dagger}&\\\
W_{V}P_{\text{sr}}&W_{H}&W_{V}^{\dagger}\\\ &W_{V}&W_{H}&W_{V}^{\dagger}\\\
&&W_{V}&W_{H}&W_{V}^{\dagger}\\\ &&&\ddots&\ddots&\ddots\end{pmatrix},$ (252)

where $W_{\text{sr}}^{\dagger}=W_{\text{sr}}$ and $W_{H}^{\dagger}=W_{H}$ by
definition. Here we have also assumed that the overlap matrix $W_{\text{sys}}$
has the same sparsity structure as the Hamiltonian $H_{\text{sys}}$, which can
always be achieved by choosing the unit cell of the leads appropriately.

Most of the results in this review are obtained by direct manipulation of the
Schrödinger equation $\hat{H}_{\text{sys}}\hat{\psi}=E\hat{\psi}$. To
generalize to the case of a non-orthogonal basis, i.e., (251), it is thus
sufficient to replace

$\begin{split}H_{\text{sr}}-E&\rightarrow H_{\text{sr}}-E\,W_{\text{sr}},\\\
H-E&\rightarrow H-E\,W_{H}\text{, and}\\\ V&\rightarrow V-E\,W_{V}\end{split}$
(253)

in our main results such as Eq.(50).

The Bloch equation for the lead modes now reads

$H(k)\phi=E(k)W(k)\phi$ (254)

with $W(k)=W_{V}e^{-ik}+W_{H}+W_{V}^{\dagger}e^{ik}$. In analogy to Appendix
A, the expression of the velocity of a lead mode is obtained by projecting the
equation on $\phi^{\dagger}$ and taking a derivative with respect to $k$. This
yields

$v=\frac{1}{\hbar}\frac{dE}{dk}=\frac{i}{w}\phi^{\dagger}\left[\lambda(V^{\dagger}-EW_{V}^{\dagger})-\lambda^{-1}(V-EW_{V})\right]\phi,$
(255)

where $w=\phi^{\dagger}(\lambda W_{V}^{\dagger}+W_{H}+\lambda^{-1}W_{V})\phi$,
and again $\lambda=e^{ik}$. In our derivations we used the convention that
propagating modes are normalized such that $v=1$. Hence, for the case of a
non-orthogonal basis set, we need to absorb both $v$ and $w$ in the
normalization of $\phi$. With this convention, the properties of the lead
modes derived in Appendix B hold with the replacement rules (253), as they are
again obtained by direct manipulation of the Bloch equation in the leads.

Hence, with the replacement rules (253) the results presented in this review
can immediately be applied to systems with a non-orthogonal basis set.

##  Appendix E Notation

$\hat{H}_{\text{sys}}$ / $\hat{G}_{\text{sys}}$ / $\hat{\psi}$ |  Hamiltonian / Green’s function / Wave function of full (infinite) systems. The hat indicate an infinite system while subscript specifies the system (infinite lead alone $H_{\text{lead}}$ or semi-infinite lead + scattering region $H_{\text{sys}}$.
---|---
$\hat{H}_{\text{lead}}$ / $\hat{\phi}$ |  Hamiltonian / Wave function of the full infinite lead.
$H_{\text{sr}}$ / $G_{\text{sr}}$ |  Hamiltonian / Green’s function of scattering region
$H$, $V$ |  Hamiltonian and hopping matrix of lead
$\Sigma$ |  Lead self-energy
$H(k)$ |  Bloch Hamiltonian of the lead
$V_{j}$ |  hopping in lead between j-th nearest neighbor cell ($V=V_{1}+V_{2}+\dots$)
$P_{\text{sr}}$ |  projection from lead to scattering region
$\Phi$ / $\phi$ |  matrix of translation eigenvectors (transverse modes) / vector of a single translation eigenvector. The columns of $\Phi$ are made of the $\phi$ vectors.
$\Lambda$ / $\lambda$ |  matrix of translation eigenvalues / single eigenvalue
$\tilde{\Lambda}$ |  truncated $\Lambda$ matrix where the vanishing diagonal elements have been disregarded.
p,± / e,± |  subscript indicating a propagating ($\mathrm{p}$) / evanescent ($\mathrm{e}$) mode, with $\pm$ indicating whether it’s incoming ($-$, going toward the scattering region) or outgoing ($+$ going away from the scattering region).
t |  subscript ($\mathrm{t}$) indicating both propagating ($\mathrm{p}$) and evanescent ($\mathrm{e}$) modes.
${}_{\mathrm{\bar{t}}}$ |  subscript ($\mathrm{\bar{t}})$ similar to $(\mathrm{t})$ but the modes $\lambda_{\alpha}=0$ are excluded. See Table 1.
${}_{\mathrm{\bar{e}}}$ |  subscript ($\mathrm{\bar{e}})$ similar to $(\mathrm{e})$ but the modes $\lambda_{\alpha}=0$ are excluded. See Table 1.
$v$ |  velocity
$J$ |  current
$N$ |  number of modes or sites. $N_{\mathrm{t}}=N_{\mathrm{p}}+N_{\mathrm{e}}$ is the number of sites in each unit cell of the lead. $N_{\text{sr}}$ is the number of sites of the scattering region.
$V=AB$ |  decomposition of the lead hopping matrix
$S$ |  scattering matrix
$S_{\text{tp}}$ |  Generalized scattering matrix containing both outgoing propagating and evanescent modes.
$\psi_{\text{sr}}$ / $\Psi_{\text{sr}}$ |  wave function in scattering region / wave function matrix: the columns of $\Psi_{\text{sr}}$ are made of the vectors $\psi_{\text{sr}}$
$E$, $\varepsilon$ |  energy
$M$ |  transfer matrix
$T$ |  transmission
$g$ |  conductance
$\Theta$ |  temperature
$\mu$ |  chemical potential
$E_{F}$ |  Fermi energy
$\lambda_{F}$ |  Fermi wave length
$V_{b}$ |  Bias voltage
$f(E)$ |  Fermi function
$I_{a}$ |  Electric current in lead $a$ (counted positive when flowing towards the scattering region)
$g_{ab}(E)$ |  Dimensionless conductance matrix between lead $a$ and $b$
$\mathcal{G}_{ab}$ |  Conductance matrix between lead $a$ and $b$
$r,t$ |  reflection and transmission matrices
$\Sigma$ |  self-energy
$P$ |  projector onto conservation law blocks
$\mathcal{C}$, $\mathcal{P}$, $\mathcal{T}$ |  chiral, particle-hole, and time-reversal symmetry
$k$ |  $k$-vector
$q$ |  vector of amplitudes of the different (propagating/evanescent) modes.
$\alpha,\beta,\gamma,\delta$ |  Greek letters label mode index.
$n,m$ |  latin letters label site index.
$i,j$ |  labels the unit cells within a single lead
$a,b,c,d$ |  labels different leads
$c_{n}$ and $c^{\dagger}_{n}$ |  are the fermionic annihilation and creation operators on site $n$.
${\cal H}$ |  is the second quantized Hamiltonian
$\Theta(\tau)$ |  is the Heaviside function
$\tau$ |  is time.
$\eta$ |  is a small imaginary energy
$BS$ |  subscript, denote solutions of Eq. (56)
$\mathds{G}$, $\mathds{H}$ |  are generic Green’s function and Hamiltonian matrices

## References

  * Abrahams  _et al._ [1979] Abrahams, E, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan (1979), “Scaling theory of localization: Absence of quantum diffusion in two dimensions,” Phys. Rev. Lett. 42, 673–676.
  * Agapito  _et al._ [2016a] Agapito, Luis A, Marco Fornari, Davide Ceresoli, Andrea Ferretti, Stefano Curtarolo, and Marco Buongiorno Nardelli (2016a), “Accurate tight-binding Hamiltonians for two-dimensional and layered materials,” Phys. Rev. B 93, 125137.
  * Agapito  _et al._ [2016b] Agapito, Luis A, Sohrab Ismail-Beigi, Stefano Curtarolo, Marco Fornari, and Marco Buongiorno Nardelli (2016b), “Accurate tight-binding Hamiltonian matrices from ab initio calculations: Minimal basis sets,” Phys. Rev. B 93, 035104.
  * Akhmerov  _et al._ [2011] Akhmerov, A R, J. P. Dahlhaus, F. Hassler, M. Wimmer, and C. W. J. Beenakker (2011), “Quantized conductance at the Majorana phase transition in a disordered superconducting wire,” Phys. Rev. Lett. 106, 057001.
  * Akhmerov  _et al._ [2009] Akhmerov, A R, C. W. Groth, J. Tworzydło, and C. W. J. Beenakker (2009), “Switching of electrical current by spin precession in the first Landau level of an inverted-gap semiconductor,” Phys. Rev. B 80, 195320.
  * Akkermans and Montambaux [2007] Akkermans, E, and G. Montambaux (2007),  _Mesoscopic Physics of Electrons and Photons_ (Cambridge University Press).
  * Alhassid [2000] Alhassid, Y (2000), “The statistical theory of quantum dots,” Rev. Mod. Phys. 72, 895–968.
  * Allen [1979a] Allen, Roland E (1979a), “Green’s function and generalized phase shift for surface and interface problems,” Phys. Rev. B 19, 917–924.
  * Allen [1979b] Allen, Roland E (1979b), “Green’s functions for surface physics,” Phys. Rev. B 20, 1454–1472.
  * Amestoy  _et al._ [1996] Amestoy, Patrick R, Timothy A. Davis, and Iain S. Duff (1996), “An approximate minimum degree ordering algorithm,” SIAM Journal on Matrix Analysis and Applications 17 (4), 886–905, https://doi.org/10.1137/S0895479894278952 .
  * Amestoy  _et al._ [2019] Amestoy, PR, A. Buttari, J.-Y. L’Excellent, and T. Mary (2019), “Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures,” ACM Transactions on Mathematical Software 45, 2:1–2:26.
  * Amestoy  _et al._ [2001] Amestoy, PR, I. S. Duff, J. Koster, and J.-Y. L’Excellent (2001), “A fully asynchronous multifrontal solver using distributed dynamic scheduling,” SIAM Journal on Matrix Analysis and Applications 23 (1), 15–41.
  * Andersen and Saha-Dasgupta [2000] Andersen, O K, and T. Saha-Dasgupta (2000), “Muffin-tin orbitals of arbitrary order,” Phys. Rev. B 62, R16219–R16222.
  * Anderson  _et al._ [1999] Anderson, E, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen (1999),  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA).
  * Ando [1991] Ando, T (1991), “Quantum point contacts in magnetic fields,” Phys. Rev. B 44 (15), 8017–8027.
  * Angers  _et al._ [2007] Angers, L, E. Zakka-Bajjani, R. Deblock, S. Guéron, H. Bouchiat, A. Cavanna, U. Gennser, and M. Polianski (2007), “Magnetic-field asymmetry of mesoscopic $\mathit{dc}$ rectification in Aharonov-Bohm rings,” Phys. Rev. B 75, 115309.
  * Antipov  _et al._ [2018] Antipov, Andrey E, Arno Bargerbos, Georg W. Winkler, Bela Bauer, Enrico Rossi, and Roman M. Lutchyn (2018), “Effects of gate-induced electric fields on semiconductor Majorana nanowires,” Phys. Rev. X 8, 031041.
  * Areshkin and Nikolić [2010] Areshkin, Denis A, and Branislav K. Nikolić (2010), “Electron density and transport in top-gated graphene nanoribbon devices: First-principles Green function algorithms for systems containing a large number of atoms,” Phys. Rev. B 81, 155450.
  * Armagnat  _et al._ [2019] Armagnat, P, A. Lacerda-Santos, B. Rossignol, C. Groth, and X. Waintal (2019), “The self-consistent quantum-electrostatic problem in strongly non-linear regime,” SciPost Phys. 7, 31.
  * Armagnat and Waintal [2020] Armagnat, Pacome, and Xavier Waintal (2020), “Reconciling edge states with compressible stripes in a ballistic mesoscopic conductor,” Journal of Physics: Materials 3 (2), 02LT01.
  * Asano  _et al._ [2006] Asano, Yasuhiro, Yukio Tanaka, Takehito Yokoyama, and Satoshi Kashiwaya (2006), “Josephson current through superconductor/diffusive-normal-metal/superconductor junctions: Interference effects governed by pairing symmetry,” Phys. Rev. B 74, 064507.
  * Averin and Bardas [1995] Averin, D, and A. Bardas (1995), “ac Josephson effect in a single quantum channel,” Phys. Rev. Lett. 75, 1831–1834.
  * Bansil  _et al._ [2016] Bansil, A, Hsin Lin, and Tanmoy Das (2016), “Colloquium : Topological band theory,” Rev. Mod. Phys. 88, 021004.
  * Baranger  _et al._ [1991] Baranger, Harold U, David P. DiVincenzo, Rodolfo A. Jalabert, and A. Douglas Stone (1991), “Classical and quantum ballistic-transport anomalies in microjunctions,” Phys. Rev. B 44, 10637–10675.
  * Baranger  _et al._ [1993] Baranger, Harold U, Rodolfo A. Jalabert, and A. Douglas Stone (1993), “Weak localization and integrability in ballistic cavities,” Phys. Rev. Lett. 70, 3876–3879.
  * Baranger and Stone [1989] Baranger, Harold U, and A. Douglas Stone (1989), “Electrical linear-response theory in an arbitrary magnetic field: A new fermi-surface formation,” Phys. Rev. B 40, 8169–8193.
  * Bardarson  _et al._ [2007] Bardarson, J H, J. Tworzydło, P. W. Brouwer, and C. W. J. Beenakker (2007), “One-parameter scaling at the Dirac point in graphene,” Phys. Rev. Lett. 99, 106801.
  * Barraza-Lopez [2013] Barraza-Lopez, Salvador (2013), “Coherent electron transport through freestanding graphene junctions with metal contacts: a materials approach,” Journal of Computational Electronics 12, 145.
  * Bautze  _et al._ [2014] Bautze, Tobias, Christoph Süssmeier, Shintaro Takada, Christoph Groth, Tristan Meunier, Michihisa Yamamoto, Seigo Tarucha, Xavier Waintal, and Christopher Bäuerle (2014), “Theoretical, numerical, and experimental study of a flying qubit electronic interferometer,” Phys. Rev. B 89, 125432.
  * Beenakker  _et al._ [1995] Beenakker, C W J, J. A. Melsen, and P. W. Brouwer (1995), “Giant backscattering peak in angle-resolved Andreev reflection,” Phys. Rev. B 51, 13883–13886.
  * Bergeron  _et al._ [2005] Bergeron, D, N. Shtinkov, R. A. Masut, and P. Desjardins (2005), “Green’s function matching method for one- and zero-dimensional heterostructures,” Phys. Rev. B 72, 245308.
  * Bernevig  _et al._ [2006] Bernevig, B Andrei, Taylor L. Hughes, and Shou-Cheng Zhang (2006), “Quantum spin Hall effect and topological phase transition in HgTe quantum wells,” Science 314, 1757, arXiv:cond-mat/0611399 .
  * Birner  _et al._ [2007] Birner, Stefan, Tobias Zibold, Till Andlauer, Tillmann Kubis, Matthias Sabathil, Alex Trellakis, and Peter Vogl (2007), “nextnano: General purpose 3-d simulations,” IEEE Transactions on Electron Devices 54 (9), 2137–2142.
  * Black-Schaffer and Doniach [2008] Black-Schaffer, Annica M, and Sebastian Doniach (2008), “Self-consistent solution for proximity effect and Josephson current in ballistic graphene SNS Josephson junctions,” Phys. Rev. B 78, 024504.
  * Blanter and Buttiker [2000] Blanter, Ya M, and M. Buttiker (2000), “Shot noise in mesoscopic conductors,” Phys. Rep. 336, 1.
  * Blonder  _et al._ [1982] Blonder, G E, M. Tinkham, and T. M. Klapwijk (1982), “Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion,” Phys. Rev. B 25, 4515–4532.
  * Borlenghi  _et al._ [2011] Borlenghi, Simone, Valentin Rychkov, Cyril Petitjean, and Xavier Waintal (2011), “Multiscale approach to spin transport in magnetic multilayers,” Phys. Rev. B 84, 035412.
  * Botello-Méndez  _et al._ [2011] Botello-Méndez, A R, Eduardo Cruz-Silva, J. M. Romo-Herrera, Florentino López-Urias, Mauricio Terrones, Bobby G. Sumpter, Humberto Terrones, Jean-Christophe Charlier, and Vincent Meunier (2011), “Quantum transport in graphene nanonetworks,” Nano Letters 11 (8), 3058, http://pubs.acs.org/doi/pdf/10.1021/nl2002268 .
  * Boykin [2009] Boykin, Timothy B (2009), “Recent developments in tight-binding approaches for nanowires,” J. Comput. Electron. 8, 142.
  * Boykin  _et al._ [2008] Boykin, Timothy B, Mathieu Luisier, and Gerhard Klimeck (2008), “Multiband transmission calculations for nanowires using an optimized renormalization method,” Phys. Rev. B 77, 165318.
  * Boykin  _et al._ [2011] Boykin, Timothy B, Mathieu Luisier, Gerhard Klimeck, Xueping Jiang, Neerav Kharche, Yu Zhou, and Saroj K. Nayak (2011), “Accurate six-band nearest-neighbor tight-binding model for the pi -bands of bulk graphene and graphene nanoribbons,” Journal of Applied Physics 109 (10), 104304.
  * Brandbyge  _et al._ [2002] Brandbyge, Mads, José-Luis Mozos, Pablo Ordejón, Jeremy Taylor, and Kurt Stokbro (2002), “Density-functional method for nonequilibrium electron transport,” Phys. Rev. B 65, 165401.
  * Bruno [1995] Bruno, P (1995), “Theory of interlayer magnetic coupling,” Phys. Rev. B 52, 411–439.
  * Bruzzone  _et al._ [2014] Bruzzone, S, G. Iannaccone, N. Marzari, and G. Fiori (2014), “An open-source multiscale framework for the simulation of nanoscale devices,” IEEE Transactions on Electron Devices 61 (1), 48–53.
  * Burt [1992] Burt, M G (1992), “The justification for applying the effective-mass approximation to microstructures,” Journal of Physics: Condensed Matter 4 (32), 6651.
  * Butler  _et al._ [2001] Butler, W H, X.-G. Zhang, T. C. Schulthess, and J. M. MacLaren (2001), “Spin-dependent tunneling conductance of $\mathrm{Fe}|\mathrm{MgO}|\mathrm{Fe}$ sandwiches,” Phys. Rev. B 63, 054416.
  * Buttiker [1993] Buttiker, M (1993), “Capacitance, admittance, and rectification properties of small conductors,” Journal of Physics: Condensed Matter 5 (50), 9361–9378.
  * Bäuerle  _et al._ [2018] Bäuerle, Christopher, D. Christian Glattli, Tristan Meunier, Fabien Portier, Patrice Roche, Preden Roulleau, Shintaro Takada, and Xavier Waintal (2018), “Coherent control of single electrons: A review of current progress,” Reports on Progress in Physics 81 (5), 10.1088/1361-6633/aaa98a.
  * Büttiker [1986] Büttiker, M (1986), “Four-terminal phase-coherent conductance,” Phys. Rev. Lett. 57, 1761–1764.
  * Büttiker [1995] Büttiker, M (1995), “Time-dependent current partition in mesoscopic conductors.” Il Novo Cimento B 110, 509.
  * Büttiker  _et al._ [1985] Büttiker, M, Y. Imry, R. Landauer, and S. Pinhas (1985), “Generalized many-channel conductance formula with application to small rings,” Phys. Rev. B 31, 6207–6215.
  * Calogero  _et al._ [2019] Calogero, Gaetano, Nick Papior, Mohammad Koleini, Matthew Helmi Leth Larsen, and Mads Brandbyge (2019), “Multi-scale approach to first-principles electron transport beyond 100 nm,” Nanoscale 11 (13), 6153–6164.
  * Calzolari  _et al._ [2004] Calzolari, Arrigo, Nicola Marzari, Ivo Souza, and Marco Buongiorno Nardelli (2004), “ _Ab initio_ transport properties of nanostructures from maximally localized Wannier functions,” Phys. Rev. B 69, 035108.
  * Carlo [2003] Carlo, Aldo Di (2003), “Microscopic theory of nanostructured semiconductor devices: beyond the envelope-function approximation,” Semiconductor Science and Technology 18 (1), R1.
  * Caroli  _et al._ [1971] Caroli, C, R Combescot, P Nozieres, and D Saint-James (1971), “Direct calculation of the tunneling current,” J. Phys. C 4 (8), 916–929.
  * Chalker and Coddington [1988] Chalker, J T, and P D Coddington (1988), “Percolation, quantum tunnelling and the integer Hall effect,” Journal of Physics C: Solid State Physics 21 (14), 2665–2679.
  * Chang  _et al._ [1994] Chang, A M, H. U. Baranger, L. N. Pfeiffer, and K. W. West (1994), “Weak localization in chaotic versus nonchaotic cavities: A striking difference in the line shape,” Phys. Rev. Lett. 73, 2111–2114.
  * Chang  _et al._ [2014a] Chang, Po-Hao, Mohammad Saeed Bahramy, Naoto Nagaosa, and Branislav K. Nikolić (2014a), “Giant thermoelectric effect in graphene-based topological insulators with heavy adatoms and nanopores,” Nano Letters 14 (7), 3779–3784, pMID: 24932511, http://dx.doi.org/10.1021/nl500755m .
  * Chang  _et al._ [2014b] Chang, Po-Hao, Farzad Mahfouzi, Naoto Nagaosa, and Branislav K. Nikolić (2014b), “Spin-Seebeck effect on the surface of a topological insulator due to nonequilibrium spin-polarization parallel to the direction of thermally driven electronic transport,” Phys. Rev. B 89, 195418.
  * Chatzikyriakou  _et al._ [2022] Chatzikyriakou, Eleni, Junliang Wang, Lucas Mazzella, Antonio Lacerda-Santos, Maria Cecilia da Silva Figueira, Alex Trellakis, Stefan Birner, Thomas Grange, Christopher Bäuerle, and Xavier Waintal (2022), “Unveiling the charge distribution of a GaAs-based nanoelectronic device: A large experimental dataset approach,” Phys. Rev. Res. 4, 043163.
  * Christen and Büttiker [1996] Christen, T, and M. Büttiker (1996), “Gauge-invariant nonlinear electric transport in mesoscopic conductors,” Europhys. Lett. 35 (7), 523–528.
  * Clerk  _et al._ [2001] Clerk, A A, X. Waintal, and P. W. Brouwer (2001), “Fano resonances as a probe of phase coherence in quantum dots,” Phys. Rev. Lett. 86, 4636–4639.
  * Costa Girão and Meunier [2013] Costa Girão, Eduardo, and Vincent Meunier (2013), “Patchwork algorithm for the parallel computation of the Green’s function in open systems,” Journal of Computational Electronics 12 (2), 123–133.
  * Crawford and Brouwer [2002] Crawford, M G A, and P. W. Brouwer (2002), “Density of proper delay times in chaotic and integrable quantum billiards,” Phys. Rev. E 65, 026221.
  * Cresti  _et al._ [2008] Cresti, A, N. Nemec, B. Biel, G. Niebler, F. Triozon, G. Cuniberti, and S. Roche (2008), “Charge transport in disordered graphene-based low dimensional materials,” Nano Research 1, 361.
  * Cuevas  _et al._ [1996] Cuevas, J C, A. Martín-Rodero, and A. Levy Yeyati (1996), “Hamiltonian approach to the transport properties of superconducting quantum point contacts,” Phys. Rev. B 54, 7366–7379.
  * Danneau  _et al._ [2008] Danneau, R, F. Wu, M. F. Craciun, S. Russo, M. Y. Tomi, J. Salmilehto, A. F. Morpurgo, and P. J. Hakonen (2008), “Shot noise in ballistic graphene,” Phys. Rev. Lett. 100, 196802.
  * Danon  _et al._ [2020] Danon, Jeroen, Anna Birk Hellenes, Esben Bork Hansen, Lucas Casparis, Andrew P. Higginbotham, and Karsten Flensberg (2020), “Nonlocal conductance spectroscopy of Andreev bound states: Symmetry relations and BCS charges,” Phys. Rev. Lett. 124, 036801.
  * Darancet  _et al._ [2009] Darancet, Pierre, Valerio Olevano, and Didier Mayou (2009), “Coherent electronic transport through graphene constrictions: Subwavelength regime and optical analogy,” Phys. Rev. Lett. 102, 136803.
  * Darancet  _et al._ [2010] Darancet, Pierre, Valerio Olevano, and Didier Mayou (2010), “Quantum transport through resistive nanocontacts: Effective one-dimensional theory and conductance formulas for nonballistic leads,” Phys. Rev. B 81, 155422.
  * Datta [1995] Datta, Supriyo (1995),  _Electronic Transport in Mesoscopic Systems_ (Cambridge University Press).
  * Davis [2006] Davis, TA (2006),  _Direct Methods for Sparse Linear Systems_, Fundamentals of Algorithms (Society for Industrial and Applied Mathematics).
  * Davis [2004] Davis, Timothy A (2004), “Algorithm 832: UMFPACK v4.3—an unsymmetric-pattern multifrontal method,” ACM Trans. Math. Softw. 30 (2), 196–199.
  * Davis  _et al._ [2016] Davis, Timothy A, Sivasankaran Rajamanickam, and Wissam M. Sid-Lakhdar (2016), “A survey of direct methods for sparse linear systems,” Acta Numerica 25, 383–566.
  * Di Ventra  _et al._ [2000] Di Ventra, M, S. T. Pantelides, and N. D. Lang (2000), “First-principles calculation of transport properties of a molecular device,” Physical Review Letters 84, 979–982.
  * Di Ventra and Lang [2001] Di Ventra, Massimiliano, and Norton D. Lang (2001), “Transport in nanoscale conductors from first principles,” Phys. Rev. B 65, 045402.
  * Driscoll and Varga [2008] Driscoll, J A, and K. Varga (2008), “Calculation of self-energy matrices using complex absorbing potentials in electron transport calculations,” Phys. Rev. B 78, 245118.
  * Drouvelis  _et al._ [2006] Drouvelis, PS, P. Schmelcher, and P. Bastian (2006), “Parallel implementation of the recursive Green’s function method,” Journal of Computational Physics 215 (2), 741 – 756.
  * Ellis  _et al._ [2017] Ellis, Matthew O A, Maria Stamenova, and Stefano Sanvito (2017), “Multiscale modeling of current-induced switching in magnetic tunnel junctions using ab initio spin-transfer torques,” Phys. Rev. B 96, 224410.
  * Esseni  _et al._ [2017] Esseni, David, Marco Pala, Pierpaolo Palestri, Cem Alper, and Tommaso Rollo (2017), “A review of selected topics in physics based modeling for tunnel field-effect transistors,” Semiconductor Science and Technology 32 (8), 083005.
  * Fabian  _et al._ [2007] Fabian, Jaroslav, Alex Matos-Abiaguea, Christian Ertler, Peter Stano, and Igor Žutić (2007), “Semiconductor spintronics,” Acta Physica Slovaca 57, 565.
  * Faleev  _et al._ [2005] Faleev, Sergey V, François Léonard, Derek A. Stewart, and Mark van Schilfgaarde (2005), “Ab initio tight-binding LMTO method for nonequilibrium electron transport in nanosystems,” Phys. Rev. B 71 (19), 195422.
  * Fan  _et al._ [2021] Fan, Zheyong, José H. Garcia, Aron W. Cummings, Jose Eduardo Barrios-Vargas, Michel Panhans, Ari Harju, Frank Ortmann, and Stephan Roche (2021), “Linear scaling quantum transport methodologies,” Physics Reports 903, 1–69, linear scaling quantum transport methodologies.
  * Fang and Kaxiras [2016] Fang, Shiang, and Efthimios Kaxiras (2016), “Electronic structure theory of weakly interacting bilayers,” Phys. Rev. B 93, 235153.
  * Ferrer  _et al._ [2014] Ferrer, J, C J Lambert, V M García-Suárez, D Zs Manrique, D Visontai, L Oroszlany, R Rodríguez-Ferradás, I Grace, S W D Bailey, K Gillemot, Hatef Sadeghi, and L A Algharagholy (2014), “GOLLUM: a next-generation simulation tool for electron, thermal and spin transport,” New Journal of Physics 16 (9), 093029.
  * Finch  _et al._ [2009] Finch, C M, V. M. García-Suárez, and C. J. Lambert (2009), “Giant thermopower and figure of merit in single-molecule devices,” Phys. Rev. B 79, 033405.
  * Fisher and Lee [1981] Fisher, Daniel S, and Patrick A Lee (1981), “Relation between conductivity and transmission matrix,” Physical Review B 23 (12), 6851.
  * Foreman [1997] Foreman, Bradley A (1997), “Elimination of spurious solutions from eight-band $\mathbf{k}\cdot\mathbf{p}$ theory,” Phys. Rev. B 56, R12748–R12751.
  * Franklin  _et al._ [2012] Franklin, Aaron D, Mathieu Luisier, Shu-Jen Han, George Tulevski, Chris M. Breslin, Lynne Gignac, Mark S. Lundstrom, and Wilfried Haensch (2012), “Sub-10 nm carbon nanotube transistor,” Nano Lett. 12, 758.
  * Fu  _et al._ [1998] Fu, Huaxiang, Lin-Wang Wang, and Alex Zunger (1998), “Applicability of the $\mathbf{k}\cdot\mathbf{p}$ method to the electronic structure of quantum dots,” Phys. Rev. B 57, 9971–9987.
  * Fuchs  _et al._ [2021] Fuchs, Jacob, Michael Barth, Cosimo Gorini, İnanç Adagideli, and Klaus Richter (2021), “Crossed Andreev reflection in topological insulator nanowire T-junctions,” Phys. Rev. B 104, 085415.
  * Fujimoto and Hirose [2003] Fujimoto, Yoshitaka, and Kikuji Hirose (2003), “First-principles treatments of electron transport properties for nanoscale junctions,” Phys. Rev. B 67, 195315.
  * Fulga  _et al._ [2012] Fulga, I C, F. Hassler, and A. R. Akhmerov (2012), “Scattering theory of topological insulators and superconductors,” Phys. Rev. B 85, 165409.
  * Gao  _et al._ [2014] Gao, X, D. Mamaluy, E. Nielsen, R. W. Young, A. Shirkhorshidian, M. P. Lilly, N. C. Bishop, M. S. Carroll, and R. P. Muller (2014), “Efficient self-consistent quantum transport simulator for quantum devices,” Journal of Applied Physics 115 (13), 133707, https://doi.org/10.1063/1.4870288 .
  * Garcia-Lekue  _et al._ [2015] Garcia-Lekue, A, M.G. Vergniory, X.W. Jiang, and L.W. Wang (2015), “Ab initio quantum transport calculations using plane waves,” Progress in Surface Science 90 (3), 292 – 318.