\section{Introduction}
\label{sec:intro}

Edge states and the bulk-boundary correspondence are key features of systems with topological properties \cite{Hasan2010,Qi2011,Bernevig2013}. Particular experimental and theoretical interest has been devoted
to electronic two-dimensional (2D) topological systems, which are characterized by edge states that are robust against disorder and by a quantized conductance.
Currently, it is well established that  the current flow in integer quantum Hall (IQH)
and quantum spin Hall (QSH) systems  is associated with gapless chiral edge states \cite{Halperin1982} and helical edge
states \cite{Kane2005}, respectively.

Advances in device fabrication and detection techniques allow one to measure edge currents in quantum Hall \cite{vanHaren1995,McCormick1999,Yacoby1999,Suddards2012,Pasher2014} and quantum spin Hall systems \cite{Roth2009,Nowack2013,Spanton2014}, using a variety of methods. As a result, it is nowadays possible to experimentally assess local current maps of 2D samples \cite{Nowack2013,Tetienne2017}, whose interpretation calls for a state-of-the-art microscopic electronic transport theory.

Of particular interest is the investigation of dissipationless currents proposed in both in quantum Hall \cite{Geller1995} and quantum spin Hall systems \cite{Buttiker2009, Sonin2011, Ando2013, Maekawa2017}.
While recent experiments in graphene \cite{Uri2020}  has rekindled the interest in equilibrium currents \cite{Lee2004}
in IQH systems, to the best of our knowledge experimental evidence of equilibrium currents in QSH systems is still missing.

Equilibrium electronic currents in (non-topological) mesoscopic conductors with broken time-reversal symmetry have been theoretically addressed a long time ago \cite{Baranger1989}, but a microscopic study in the {IQH} regime in a multi-probe set-up is still lacking.
The situation is less clear for QSH systems.
Equilibrium currents in topological insulators have been studied for massless \cite{Mishchenko2014,Chen2020} as well as for massive Dirac electrons \cite{Silvestrov2019}. Interestingly, these works overlook the literature on spin currents in the spin Hall regime.
The existence of bulk spin currents in thermodynamic equilibrium in conductors lacking inversion symmetry have been first discussed more than 15 years ago by Rashba \cite{Rashba2003}. Application proposals of equilibrium spin currents \cite{Pareek2004} have been subject of controversy, until it was shown that both the Landauer-B\"uttiker approach \cite{Kiselev2005} and NEGF theory \cite{Souma2005,Nikolic2006} are incompatible with a net equilibrium spin transport.

We address the problem of equilibrium currents using the non-equilibrium Green's functions theory (NEGF) \cite{Jauho2008} combined with the Landauer-B\"uttiker approach \cite{Landauer1970, Buttiker1985} that allows one to obtain the local transport properties of a given multi-probe system with arbitrary geometry. We employ the formalism put forward in Refs.~\cite{Cresti2003,Nikolic2006} to compute local electronic and spin currents using the recursive Green's
functions (RGF) method  \cite{MacKinnon1985,Lewenkopf2013,Lima2016,Lima2018}.
We use this approach  to study the charge and the spin flow in topological systems.
We show that the equilibrium currents cannot be measured by the standard methods and discuss strategies to assess them. In addition, our analysis also shows the necessity of separating the equilibrium and the nonequilibrium contributions to the local currents for a correct theoretical interpretation of recent experimental results.

This paper is structured as follows. In Sec.~\ref{sec:multiprobe}, we define the problem and obtain close expressions for the local charge and spin currents in a multi-probe set-up using non-equilibrium Green's functions. We identify the equilibrium and non-equilibrium current components and discuss their properties in terms of the Landauer-B\"uttiker approach. In Sec.~\ref{sec:toymodel} we introduce an analytical solvable model to show the necessity of an external magnetic field to generate equilibrium charge currents. In Sec.~\ref{sec:IQH} we study the nonequilibrium and the equilibrium currents in the integer quantum Hall regime for a realistic multi-terminal setup. In Sec.~\ref{sec:QSH} we conduct a similar analysis for quantum spin Hall systems. We summarize our results and present our conclusions in Sec.~\ref{sec:conclusions}.

\section{Electronic current in a multiprobe mesoscopic system}
\label{sec:multiprobe}

We begin this section with a brief overview of the Landauer-B\"uttiker approach \cite{Landauer1970, Buttiker1985}
for the description of the transport properties of multiprobe quantum coherent mesoscopic systems.
These results, originally obtained using the scattering matrix approach \cite{Buttiker1986,Blanter2000}, can be derived using non-equilibrium Green's functions (NEGF) \cite{Datta1995,Hernandez2013}. In this paper we adopt the latter.
Next, we discuss the theory that describes local equilibrium and non-equilibrium electronic currents \cite{Todorov2002, Cresti2003, Nikolic2006, Zarbo2007}. We cast the corresponding expressions in a convenient way to implement an efficient recursive Green's function (RGF) method \cite{Lewenkopf2013,Lima2016,Lima2018} that is used to compute the results reported in Secs.~\ref{sec:IQH} and \ref{sec:QSH}.

We consider the standard model Hamiltonian
\begin{align}
	H=H_C+H_L+V,
	\label{hamiltonian}
\end{align}
where $H_C$ describes the mesoscopic system $C$, $H_L$ the leads that connect $C$ to electronic reservoirs in thermal and chemical equilibrium, and $V$ the coupling between the mesoscopic system $C$ and the leads.

The multiterminal Landauer-B\"uttiker formula for the electronic current $I_\alpha$ at the terminal
$\alpha$ reads \cite{Buttiker1986,Datta1995}
\begin{align}
	I_\alpha
	&= \frac{e}{h} \sum_\beta \int_{-\infty}^{+\infty} dE\ \mathcal T_{\alpha\beta}(E) \left[f_\alpha(E)-f_\beta(E)\right] ,
\label{totalcurrent0}
\end{align}
where the Greek letters label the terminals $\alpha=1,\cdots,\Lambda$.
The terminals, usually modeled by semi-infinite leads, are in contact with electronic reservoirs in thermal equilibrium and inject into the system electrons following a Fermi-Dirac distribution $f_\alpha(E)=[1+e^{(E-\mu_\alpha)/k_BT}]^{-1}$, where $\mu_\alpha=\mu_0+eV_\alpha$, $\mu_0$ is the equilibrium chemical potential and $V_\alpha$ is the voltage applied to the $\alpha$-terminal. For the sake of simplicity, we consider that all electronic reservoirs have the same temperature $T$.

The transmission $\mathcal T_{\alpha\beta} (E)$ is given by \cite{MeirWingreen1992}
\begin{align}
\label{eq:transmission}
\mathcal T_{\alpha\beta} (E)= {\rm tr} \left[ \mathbf \Gamma_\alpha(E) \mathbf G^r(E) \mathbf \Gamma_\beta(E) \mathbf G^a(E) \right]
\end{align}
where $\mathbf G^r=\left(\mathbf G^a\right)^\dagger$ is the retarded Green's function of the full system,
while $\mathbf \Gamma_\beta$ represents the decay width of the lead corresponding to the $\beta$-terminal.
For convenience we express $\mathbf G^r$ and $\mathbf \Gamma_\beta$ in a local basis
representation that renders simple expressions for the local currents. Here we consider a Wannier basis with the states as labelled by $\nu=(i,\ell, \sigma)$, where the corresponding wave function is centered at the lattice atom $i$ and has orbital and spin quantum numbers given by $(\ell, \sigma)$.
In this representation, $\mathbf G^r$ has the dimension of the number of Wannier states in the central region, the dimension of $\mathbf \Gamma_\beta$ is the number of states corresponding to sites at the $\beta$-terminal-central region interface.
The decay width reads
\begin{align}
\mathbf	\Gamma_\alpha = -2\,\text{Im}\left(\mathbf \Sigma_\alpha^r\right),
\end{align}
where $\mathbf \Sigma_\alpha^r$ is the retarded embedding self-energy, namely
\begin{align}
\mathbf{\Sigma}_\alpha^r = \mathbf V_{C \alpha} \mathbf G^r_\alpha \mathbf V^{}_{\alpha C},
\end{align}
$\mathbf V_{C \alpha}=\mathbf V_{\alpha C}^\dagger$
gives the coupling matrix elements between the terminal $\alpha$ and the central region $C$,
and $\mathbf G^r_\alpha$ is a contact Green's function that casts the electron dynamics in the leads.
The latter can be calculated in a number of ways \cite{LopezSancho1985, MacKinnon1985, Rocha2006, Wimmer2009thesis}.

In linear approximation, where $f_\alpha(E)=f_0(E)+(-\partial f_0/\partial E)eV_\alpha$ and $f_0(E)$ is the equilibrium Fermi-Dirac distribution, one obtains the familiar result \cite{Buttiker1986}
\begin{align}
I_\alpha =-\sum_{\beta=1}^{\Lambda} \mathcal G_{\alpha\beta} V_\beta = \sum_{\beta=1}^{\Lambda} \mathcal G_{\alpha\beta} \left(V_\alpha-V_\beta \right),
\label{current}
\end{align}
where $\mathcal G_{\alpha \beta}$ is the Landauer conductance given by
\begin{align}
\mathcal G_{\alpha \beta} = \frac{e^2}{h}\int_{-\infty}^{\infty} dE \left(-\frac{\partial f_0}{\partial E} \right) \mathcal T_{\alpha\beta}(E).
\label{conductance0}
\end{align}
Here $\mathcal T_{\alpha\beta}(E)$ accounts for the spin degrees of freedom and, hence, there is no spin degeneracy factor in Eq.~\eqref{conductance0}.
As standard, we chose the sign of $I_\alpha$, Eqs.~(\ref{totalcurrent0}) and (\ref{current}), to ensure that a positive current flows from terminal $\beta$ to terminal $\alpha$ when $V_\beta>V_\alpha$ or $\mu_\beta>\mu_\alpha$.
See Appendix~\ref{appendix:lbequations} for details.

Let us now calculate the local currents using NEGF.
We begin by considering the number operator $N_{\nu}(t) = d_{\nu}^\dagger(t) d^{}_{\nu}(t)$, where
$d_{\nu}^\dagger$ ($d^{}_{\nu}(t)$) creates (annihilates) an electron at the state $\nu$. As discussed
above, here $\nu$ labels a local basis of states. The local electron flow
is associated to the expectation value of the temporal rate of change of the number operator, namely,
$\langle \dot N_{\nu}(t)\rangle$. The average is taken over the grand canonical ensemble, namely, $\langle \cdots \rangle = \text{Tr}[e^{-\beta H}\cdots]/Z$, where $Z=\text{Tr}[e^{-\beta H}]$ is the partition function and $\beta = 1/k_BT$.
This can be done by different methods \cite{Stefanucci2013}.

We calculate $\langle \dot N_{\nu}(t)\rangle$ using the equations-of-motion method, in a similar way as it is done for the total electronic current \cite{Jauho2008,Lima2013}, and obtain
\begin{align}
	\langle \dot{N}_{\nu}(t)\rangle
	=  \int_{-\infty}^{+\infty} &\frac{dE}{h}  \bigg[
		\mathbf G^{<} \mathbf H	- \mathbf H	\mathbf G^{<} \nonumber\\
		& \!\!\!\!+\sum_\alpha
		\left(\mathbf G^{<}_{C\alpha} \mathbf V_{\alpha C} - \mathbf V_{C\alpha} \mathbf G^{<}_{\alpha C}\right) \bigg]_{\nu,\nu},
\label{current3_local}
\end{align}
where $\mathbf G^<$ is the lesser Green's function of sites in the central region and $\mathbf G^<_{\alpha C}$ is hybrid lesser Green's function containing propagation information between the $\alpha$-lead and the central region.

Recalling the identities $\mathbf H	\mathbf G^{<} = -\left(\mathbf G^{<} \mathbf H\right)^\dagger$ and $\mathbf V^{}_{C\alpha} \mathbf G^{<}_{\alpha C} = -\left(\mathbf G^{<}_{C\alpha} \mathbf V^{}_{\alpha C}\right)^\dagger$ we simplify the diagonal matrix element in Eq.~(\ref{current3_local}) to
\begin{align}
	\langle \dot N_{\nu}(t)\rangle
	&=  \frac{2}{h}\int_{-\infty}^{+\infty} \!\!{dE} \,{\rm Re}\bigg[
		\mathbf G^{<} \mathbf H	+ \sum_\alpha \mathbf G^{<}_{C\alpha} \mathbf V^{}_{\alpha C} \bigg]_{\nu,\nu}.
\label{current4_local}
\end{align}

It is straightforward to demonstrate that
\footnote
{
Using the standard NEGF relations
$\mathbf G^<_{C\alpha} =
\mathbf G^r \mathbf V_{C\alpha} \mathbf G^<_{\alpha} +
\mathbf G^< \mathbf V_{C\alpha} \mathbf G^a_{\alpha} $
and
$\mathbf G^< = \mathbf G^r \mathbf \Sigma^< \mathbf G^a $,
one writes
$
\sum_\alpha \mathbf G^<_{C\alpha} \mathbf V_{\alpha C}
= \mathbf G^r \mathbf \Sigma^< +\mathbf G^< \mathbf \Sigma^a
= \mathbf G^< \left[ \left(\mathbf G^a\right)^{-1}  + \mathbf \Sigma^a \right]
= \mathbf G^< \left(\mathbf E -\mathbf H \right).
$
Analogously,
$
\sum_\alpha \mathbf V_{C\alpha} \mathbf G^<_{\alpha C}
= \left(\mathbf E -\mathbf H \right)\mathbf G^<
$
}
\begin{align}
	\mathbf G^{<} \mathbf H	- \mathbf H	\mathbf G^{<}
	+\sum_\alpha
		\left(\mathbf G^{<}_{C\alpha} \mathbf V_{\alpha C} - \mathbf V_{C\alpha} \mathbf G^{<}_{\alpha C}\right) = 0,
\end{align}
which explicitly shows that in steady state the average rate of change of the number of particles in any site is zero.

In order to calculate the current flowing between any pair of states connected by the model Hamiltonian we write Eq.~(\ref{current4_local}) as
\begin{align}
	\left\langle \dot q_{\nu}(t)\right\rangle \equiv e\left\langle \dot N_{\nu}(t)\right\rangle = - \sum_{\nu'} \widetilde{I}_{\nu'\nu},
\label{localcurrent0}
\end{align}
where $\left\langle \dot q_{\nu}(t)\right\rangle$ is the rate of change of the charge in the state $\nu$ and ${\widetilde{I}}_{\nu'\nu}$ is the charge current flowing from the state $\nu$ to $\nu'$.

By comparing Eqs.~\eqref{current4_local} and \eqref{localcurrent0} the local electronic current $\widetilde{I}_{\nu'\nu}$ can be defined as
\begin{align}
	\widetilde{I}_{\nu'\nu}
	&\equiv  -\frac{2e}{h}\!\int_{-\infty}^{+\infty} \!\!dE\ \text{Re}\bigg[
		 G^{<}_{\nu\nu'}  H^{}_{\nu'\nu}	+ \sum_\alpha  G^{<}_{C\alpha,\nu\nu'}  V^{}_{\alpha C,\nu'\nu} \bigg] \nonumber\\
						&= \widetilde{I}_{CC,\nu'\nu} + \sum_\alpha \widetilde{I}_{\alpha C,\nu'\nu}.
\label{current5_local}
\end{align}
Here, the subscripts $CC$ and $\alpha C$ denote the partition to which the sites corresponding to $\nu'$ and $\nu$ belong (in that order).
Thus, $H_{\nu\nu'}\neq 0$ (or $V_{\nu\nu'}\neq 0$) is a necessary condition for a non vanishing local or bond current between the states $\nu$ and $\nu'$.
The sum of all the bond currents across the interface between the system and the $\alpha$-lead yields
\begin{align}
	\sum_{\nu'\nu} \widetilde{I}_{\alpha C,\nu'\nu}
	&=
	-\frac{e}{h}\int_{-\infty}^{+\infty} \!\!dE\
	{\rm Tr}\left[ \mathbf G^{<}_{C\alpha}  \mathbf V^{}_{\alpha C} - \mathbf V^{}_{C\alpha} \mathbf G^{<}_{\alpha C}\right],
\end{align}
which, due to charge conservation, is identical to $I_\alpha$, the total electronic current at terminal $\alpha$, given by Eq.~(\ref{totalcurrent0}) \cite{Lima2013}.

Let us assume that the bonds of interest are located sufficiently far from any lead interface so that all $V_{\alpha C,\nu'\nu}$ matrix elements are identically zero. Hence
\begin{align}
	\widetilde{I}_{\nu'\nu} = \widetilde{I}_{CC,\nu'\nu}
	&=  -\frac{e}{h}\int_{-\infty}^{+\infty} dE\
	2{\rm Re}\!\left( G^{<}_{\nu\nu'}  H^{}_{\nu'\nu} \right).
				\label{localcc}
\end{align}
Using the Keldysh equation \cite{Jauho2008} written as
\begin{align}
	\mathbf G^< = \mathbf G^r \mathbf \Sigma^< \mathbf G^a = \sum_\alpha if_\alpha\mathbf G^r \mathbf \Gamma_\alpha \mathbf G^a,
	\label{keldysh1}
\end{align}
Eq.~(\ref{localcc}) becomes
\begin{align}
	\widetilde{I}_{\nu'\nu}
	&=  \frac{e}{h}\sum_\alpha \int_{-\infty}^{+\infty} dE\ f_\alpha(E) \widetilde{\mathcal T}_{\nu'\nu}^{\alpha}(E),
\label{localcc2}
\end{align}
where
\begin{align}
	\widetilde{\mathcal T}_{\nu'\nu}^{\alpha}(E) \equiv
	2\,\text{Im}\left[ \left(\mathbf G^r \mathbf \Gamma_\alpha \mathbf G^a\right)_{\nu\nu'}  H^{}_{\nu'\nu} \right].
				\label{localtcc}
\end{align}
Recalling that $\mathbf G^r \mathbf \Gamma_\alpha \mathbf G^a$ is Hermitian,
it is straightforward to show that
$\widetilde{\mathcal T}_{\nu\nu'}^{\alpha}=-\widetilde{\mathcal T}_{\nu'\nu}^{\alpha}$.

Expanding $f_\alpha(E)$ up to linear order in the terminal applied voltages $\{V_\beta\}$, we write
the local current as
\begin{align}
		\widetilde{I}_{\nu'\nu}^{}		&= \widetilde{I}_{\nu'\nu}^{\rm neq}  + \widetilde{I}_{\nu'\nu}^{\rm eq},
		\label{localcurrent}
\end{align}
where the non-equilibrium component reads
\begin{equation}
		\widetilde{I}_{\nu'\nu}^{\rm neq}\equiv  \sum_\alpha \widetilde{\mathcal G}_{\nu'\nu}^{\alpha} V_\alpha
		\label{iklneq}
\end{equation}
with a local conductance
\begin{equation}
		\widetilde{\mathcal G}_{\nu'\nu}^{\alpha} \equiv \frac{e^2}{h}\int_{-\infty}^{+\infty} dE\ \left(-\frac{\partial f_0}{\partial E}\right) \widetilde{\mathcal T}_{\nu'\nu}^{\alpha}(E).
		\label{conductancelocal}
\end{equation}
In turn, the equilibrium local charge current is given by
\begin{equation}
		\widetilde{I}_{\nu'\nu}^{\rm eq}\equiv  \frac{e}{h}\int_{-\infty}^{+\infty} dE\ f_0(E) \left[\sum_\alpha \widetilde{\mathcal T}_{\nu'\nu}^{\alpha}(E)\right].
\label{ikleq}
\end{equation}

Let us project $\widetilde{I}_{\nu'\nu}^{}$ to real space. Recalling that  $\nu = (i, \ell, \sigma)$ and the atomic position can be represented by the index $i$, one can define the ``bond'' charge current as
\begin{equation}
\widetilde I_{i'i} \equiv \sum_{\ell,\ell',\sigma, \sigma'} \widetilde{I}_{\nu'\nu}^{}.
\end{equation}
Similarly, one can define a bond spin current \cite{Nikolic2006}
\begin{align}
\widetilde I_{i'i}^s \equiv \frac{\hbar}{2e} \sum_{\ell,\ell'} \left[\widetilde{I}_{(i',\ell',\uparrow),(i,\ell,\uparrow)} -
\widetilde{I}_{(i',\ell',\downarrow),(i,\ell,\downarrow)} \right],
\label{eq:I_spin_local}
\end{align}
where $\hbar/2e$ converts the units of a charge current into a spin one.
By suitable projection schemes, one can also obtain local valley currents \cite{Settnes2016} and orbital currents
\cite{Bernevig2005,Cysne2021} from $\widetilde{I}_{\nu'\nu}$.

Note that, similarly to the total currents $I_\alpha$, the local nonequilibrium current $\widetilde{{\bf I}}^{\rm neq}$ involves an energy integration over a typically small energy window of width $ \sim k_BT$ around $\mu_0$,
In distinction, the evaluation of the local equilibrium currents $\widetilde{{\bf I}}^{\rm eq}$ requires an energy integration of the local transmissions weighted by the equilibrium Fermi distribution over the entire band.

In the zero bias limit, $\{V_\alpha\} = 0$, Eqs.~(\ref{current}) and (\ref{localcurrent}) yield
\begin{align}
{I}_\alpha = 0 \qquad \mbox{and} \qquad
\widetilde{I}_{\nu'\nu} = 	\widetilde{I}_{\nu'\nu}^{\rm eq}.
\end{align}
These relations imply that if $\sum_\alpha \widetilde{\mathcal T}_{\nu'\nu}^{\alpha}(E) \neq 0$ it is possible to have a finite local electronic current in the central region, $\widetilde{I}_{\nu'\nu}^{\rm eq} \neq 0$, even in the absence of a bias voltage, but there is no net current flowing through the contacts ${I}_\alpha = 0$. The latter implies that there is no entropy production \cite{Bruch2018}, as required by equilibrium processes.
Equation (\ref{localcc2}) indicates that each lead $\alpha$ contributes by injecting electrons in the system for all energies up to its chemical potential $\mu_\alpha$.
Therefore, it is only natural that in the zero bias case ($\{V_\alpha\} = 0$) all the leads inject electrons for all energies up to the equilibrium chemical potential $\mu_0$.
In what follows we explore situations where the equilibrium steady state currents injected by all terminals do not cancel each other and there is a non-vanishing local current in the absence of bias voltage.

\section{Toy model}
\label{sec:toymodel}

Let us illustrate some of the main features of the charge equilibrium currents using an analytically solvable (non-topological) toy model.
We consider a central region consisting of 3 sites coupled to semi-infinite linear chains, see Fig.~\ref{fig:toymodel}.
We assume that the electrons are described by a nearest neighbor single-orbital tight-binding Hamiltonian with hopping matrix elements $-t_{kk'}$ with no spin dependence.
The central region sites lie in the $xy$-plane at $\mathbf r_1 = (0,0)$, $\mathbf r_2 = (a/2,a/2)$, $\mathbf r_3 = (a,0)$.
We account for a constant magnetic field $\mathbf B = B\hat z$ using the Peierls substitution, namely,
$t_{kk'} \rightarrow t_{kk'}e^{i\varphi_{kk'}}$, where
$\varphi_{kk'} = \frac{eB}{2\hbar} \left(x_{k'}-x_k\right) \left(y_{k'}+y_k\right)$ \cite{Lewenkopf2013}.
Thus,
\begin{align}
	\varphi_{12} = \varphi_{23} = \phi/2
	\quad \text{and} \quad
	\varphi_{13} = 0,
\end{align}
where $\phi \equiv 2\pi \Phi/\Phi_0$ is proportional to the ratio between the magnetic flux $\Phi = Ba^2/4$ enclosed by the triangular ``loop" in the central region (see Fig.~\ref{fig:toymodel}) and the magnetic flux quantum $\Phi_0 = h/e$.

\begin{figure}
	\centering
		\includegraphics[width=1.00\columnwidth]{fig01.pdf}
	\caption{Sketch of the toy model system: Three sites in the central region attached to two semi-infinite leads at left ($L$) and right ($R$).
	The dashed line represents the cross section used to evaluate the net local equilibrium current in the main text.
	}
	\label{fig:toymodel}
\end{figure}

Assuming that the sites $1$ and $3$ are attached to leads, the corresponding self-energies are $\Sigma_1^r=-i\Gamma_1/2$ and $\Sigma_3^r=-i\Gamma_3/2$, we write $\mathbf G^r(E)=(E\mathbf 1 - \mathbf H - \mathbf \Sigma^r)^{-1}$ as
\begin{align}
  \mathbf G^r(E) &=
	\left(\begin{array}{ccc}
		E+i\Gamma_1/2				& t_{12}e^{i\phi/2}		& t_{13}            \\
		t_{12}e^{-i\phi/2}	& E										& t_{23}e^{i\phi/2} \\
		t_{13}							& t_{23}e^{-i\phi/2}	& E+i\Gamma_3/2
	\end{array} \right)^{-1}.
	\label{ehstoymodel}
\end{align}

{
It is instructive to verify charge conservation by calculating the net electronic transmission through the central region sites.
Considering the injection from the left lead $L$, the net electronic transmissions read
\begin{align}
	\widetilde{\mathcal T}^{L}_{21}+\widetilde{\mathcal T}^{L}_{31} &=\frac{\Gamma_1\Gamma_3}{|D|^2} ( \tau_4-E\tau_3\cos\phi+ E^2\tau_2 ),\label{localt1}\\
	\widetilde{\mathcal T}^{L}_{21}+\widetilde{\mathcal T}^{L}_{23} &=0, \label{localt2}\\
	\widetilde{\mathcal T}^{L}_{31}+\widetilde{\mathcal T}^{L}_{32} &=\frac{\Gamma_1\Gamma_3}{|D|^2} ( \tau_4-E\tau_3\cos\phi+ E^2\tau_2 ),\label{localt3}
\end{align}
where $D = {\rm det}[ ({\bf G}^r)^{-1}], \tau_2\equiv t_{13}^{2}$, $\tau_3\equiv 2t_{12}t_{13}t_{23}$, and $\tau_4\equiv t_{23}^{2}t_{12}^{2}$. The derivation of the expressions for the transmissions are presented in  App.~\ref{app:gfsandtransmissions} and the sign convention is discussed in App. \ref{app:toymodel-sign-convention}.
The net electronic transmission at site $2$, given by Eq.~(\ref{localt2}), vanishes due to charge conservation, that is, the local transmission from site $1$ to site $2$ equals the local transmission from site $2$ to site $3$.
In contrast, the electronic transmissions given by Eqs.~(\ref{localt1}) and (\ref{localt3}), respectively, do not vanish since we are not taking into account the transmissions through the bonds with the leads, {\sl i.e.}, electrons entering site $1$ through the left lead $L$ and leaving site $3$ through the right lead $R$.
Moreover, assuming injection from the left lead $L$, the total transmission leaving site $1$, Eq.~(\ref{localt1}), must equal the transmission entering site $3$, Eq.~(\ref{localt3}).
}

{
The total transmission between the left and the right leads reads
\begin{align}
\mathcal T_{RL}
= \frac{\Gamma_1\Gamma_3}{|D|^2} ( \tau_4-E\tau_3\cos\phi+ E^2\tau_2 )
\label{trltoymodel}
\end{align}
This ensures charge conservation
for site $1$ since $\widetilde{\mathcal T}^{L}_{1L}+\widetilde{\mathcal T}^{L}_{12}+\widetilde{\mathcal T}^{L}_{13}=0$ and
for site $3$ since $\widetilde{\mathcal T}^{L}_{3R}+\widetilde{\mathcal T}^{L}_{31}+\widetilde{\mathcal T}^{L}_{32}=0$.
Here we used that $\widetilde{\mathcal T}^{L}_{1L}=-{\mathcal T}_{LR}$ and $\widetilde{\mathcal T}^{L}_{3R}={\mathcal T}_{LR}$.
The above arguments are also valid assuming injection from the right lead.
An analogous calculation can show that
$\widetilde{\mathcal T}^{R}_{3R}+\widetilde{\mathcal T}^{R}_{31}+\widetilde{\mathcal T}^{R}_{32}=0$ and
$\widetilde{\mathcal T}^{R}_{1L}+\widetilde{\mathcal T}^{R}_{12}+\widetilde{\mathcal T}^{R}_{13}=0$.
}

We are now ready to analyze the equilibrium current. Let us consider, for instance, the bond current flowing from site 2 to site 1, namely
\begin{equation}
\widetilde{I}_{12}^{\rm eq} = \frac{e}{h} \int_{-\infty}^\infty dE f_0(E) \left[ \widetilde{\mathcal T}^{L}_{12}(E) +\widetilde{\mathcal T}^{R}_{12}(E)\right]
\label{i12eq}
\end{equation}
with the sum of local transmissions given by
\begin{align}
\widetilde{\mathcal T}^{L}_{12}&+\widetilde{\mathcal T}^{R}_{12} =
-\frac{2}{|D|^2}t_{12}t_{23}t_{13} \times \nonumber\\
& \times  \left[ (\Gamma_1+\Gamma_3) E^2 - (\Gamma_1 t_{23}^2 + t_{12}^2 \Gamma_3)  \right] \sin\phi.
\label{3x3teq}
\end{align}

Interestingly, $\widetilde{\mathcal T}^{L}_{12}+\widetilde{\mathcal T}^{R}_{12}$ (and thus $\widetilde{I}_{12}^{\rm eq}$)
vanishes in three situations:
(i) if one breaks the loop by turning off any of the hopping matrix elements $t_{12}$, $t_{23}$ or $t_{13}$,
(ii) if there is no magnetic field ($\phi=0$) and
(iii) if the system is detached from the leads ($\Gamma_1=\Gamma_3=0$).

The presence of a magnetic field breaks time reversal symmetry inducing a preferential electronic flow through one of the two system branches, $1$-$3$ or $1$-$2$-$3$, depending on the injection direction.
This causes a current imbalance between the branches leading to a non-vanishing equilibrium electronic current given by Eq.~(\ref{i12eq}). Conversely, once the loop is broken by disconnecting one of the bonds, the electronic current flows through a single branch irrespective of the injection direction and there is no local current in the equilibrium.

The simplicity of the model allows us to analytically show that
\begin{align}
 \widetilde{\mathcal T}^{L}_{12}+\widetilde{\mathcal T}^{R}_{12}
=\widetilde{\mathcal T}^{L}_{31}+\widetilde{\mathcal T}^{R}_{31}
=\widetilde{\mathcal T}^{L}_{23}+\widetilde{\mathcal T}^{R}_{23}.
\label{eq:Tsums_toymodel}
\end{align}
The above expressions highlight two important properties of equilibrium currents:

{\sl (i) There is not net equilibrium current leaving or entering the system:}
The total equilibrium transmission from $2$ to $1$, given by the sum in Eq.~(\ref{i12eq}), namely $\widetilde{\mathcal T}^{L}_{12}+\widetilde{\mathcal T}^{R}_{12}$, equals the total equilibrium transmission from $1$ to $3$ given by $\widetilde{\mathcal T}^{L}_{31}+\widetilde{\mathcal T}^{R}_{31}$.
By invoking charge conservation and by accounting for the injection of both left and right terminals, we find that there is no net equilibrium transmission between $1$ and $L$.
The same happens for the net equilibrium transmission between the site $3$ and the right reservoir $R$.
In summary, both terminals inject the same equilibrium current into the system leading to a zero net {equilibrium} electronic current, that is, $\widetilde{I}^{\rm eq}_{1L} = \widetilde{I}^{\rm eq}_{3R} =0$. A non-vanishing net electronic current at the terminals requires a voltage bias. In this case, Eq.~(\ref{iklneq}) does not vanish and a non-equilibrium electronic current sets in.

{\sl ($ii$) The equilibrium current integrated over any system cross section is identically zero.}
This is nicely seen by analyzing the equilibrium current $\widetilde{I}^{\rm eq}_{\rm sec}$ flowing through the section
defined by the dashed line in Fig.~\ref{fig:toymodel}, namely,
\begin{equation}
		\widetilde{I}_{\rm sec}^{\rm eq}\equiv  \frac{e}{h}\int_{-\infty}^{+\infty} dE\ f_0(E)
		 \left[\widetilde{\mathcal T}^{L}_{12}+\widetilde{\mathcal T}^{R}_{12}
                  +\widetilde{\mathcal T}^{L}_{13}+\widetilde{\mathcal T}^{R}_{13} \right].
\label{iseceq}
\end{equation}
Recalling that $\widetilde{\mathcal T}^\alpha_{ij}=-\widetilde{\mathcal T}^\alpha_{ji}$ and Eq.~\eqref{eq:Tsums_toymodel},
it immediately follows that $\widetilde{I}^{\rm eq}_{\rm sec} =0$.

\section{Local currents in the quantum Hall regime}
\label{sec:IQH}

In this section we study the equilibrium and non-equilibrium local currents in the IQH regime \cite{Prange1990} for a multiprobe setup.
First, we discuss general qualitative aspects of the transport properties using the Landauer-B\"uttiker approach.
Next, we calculate the local and total currents using the formalism presented in Sec.~\ref{sec:multiprobe}
considering a graphene system as a case in point.

\begin{figure}[tbp]
	\centering
		\includegraphics[width=0.75\columnwidth]{fig02.pdf
\caption{Sketch of the $6$-terminal Hall bar used to discuss the electronic transport in the IQH and QSH regimes.
We assume that an electronic current $I$ is injected by the terminal $\alpha=1$ and collected at the terminal $\beta=6$.
The electronic current $I$ is driven by a bias $V_{\rm bias}=V_1-V_6$.}
	\label{fig:hall_bars_setup_A}
\end{figure}

Let us consider $6$-terminal setup, as sketched in Fig.~\ref{fig:hall_bars_setup_A},
where the applied bias voltage $V_{\rm bias}$ between terminals $\alpha=1$ and $\beta=6$ drives the electronic current $I$ from terminal $1$ ($I_1=I$) to terminal $6$ ($I_6=-I$).
For simplicity, we set $V_1=V_{\rm bias}$ and assume the terminal $6$ as grounded, $V_6=0$
We consider terminals $2$ through $5$ as voltage probes, namely, $I_2=I_3=I_4=I_5=0$.
The current $I$ and the voltages $V_2$, $V_3$, $V_4$ and $V_5$ are determined by the applied bias $V_{\rm bias}$.
We employ the same setup in Sec.~\ref{sec:QSH} to discuss the quantum spin Hall regime considering spin resolved currents.

\subsection{General discussion}
\label{sec:QH_general}

Here we address the charge transport properties of a Hall bar system in the IQH regime \cite{Prange1990,Datta1995}  paying particular attention to equilibrium currents.

We consider that the Hall bar, see Fig.~\ref{fig:hall_bars_setup_A}, is subjected to a perpendicular magnetic field $\mathbf B = B \mathbf{\hat z}$ that is sufficiently strong to give rise to quantized edge states.
The electrons injected at the terminal $\alpha=1$ along the $x$-direction are deflected towards their left-hand side flowing to the nearest terminal $\beta=2$ with edge states propagating along the $y$-direction.

Figure~\ref{fig:hall_bar_QH_local_t} shows the local transmissions for all possible single terminal injection possibilities.
In this section we neglect spin-dependent processes and, thus, the edge states are spin degenerate.

\begin{figure}[tbp]
	\centering
		\includegraphics[width=0.49\columnwidth]{fig03a.pdf}
		\includegraphics[width=0.49\columnwidth]{fig03b.pdf}
		\includegraphics[width=0.49\columnwidth]{fig03c.pdf}
		\includegraphics[width=0.49\columnwidth]{fig03d.pdf}
		\includegraphics[width=0.49\columnwidth]{fig03e.pdf}
		\includegraphics[width=0.49\columnwidth]{fig03f.pdf}
	\caption{
	Sketches of local transmissions/conductances for a system in the IQH regime.
	Panels (a) to (f) correspond to an electronic injection from a single terminal $\alpha=1$ to $6$, respectively.
	Spin dependent process are neglected and, thus, spin up (red) and spin down (blue) states are degenerate.
	}
		\label{fig:hall_bar_QH_local_t}
\end{figure}

For simplicity, we consider the zero-temperature limit and, further, assume that $E_F$ lies between the $N$th and $(N+1)$th Landau level energies.
In this situation, there are $N$ propagating channels per spin at any of the system edges.
The conductance matrix, given by Eq.~(\ref{conductance0}), reads
\begin{align}
	\boldsymbol{\mathcal{G}} =
	2N\frac{e^2}{h}
	\left(\begin{array}{cccccc}
-1 & 0 & 0 & 1 & 0 & 0\\
 1 &-1 & 0 & 0 & 0 & 0\\
 0 & 1 &-1 & 0 & 0 & 0\\
 0 & 0 & 0 &-1 & 1 & 0\\
 0 & 0 & 0 & 0 &-1 & 1\\
 0 & 0 & 1 & 0 & 0 &-1
\end{array}\right).
\label{conductanceqh}
\end{align}
The negative diagonal conductance values manifestly enforces current conservation.
Current conservation and gauge invariance imply that
$\sum_\alpha \mathcal G_{\alpha\beta} =\sum_\beta \mathcal G_{\alpha\beta} =0$ \cite{Buttiker1986,Ihn2010}.

Using Eqs.~(\ref{current}) and (\ref{conductanceqh}) we calculate the unknown voltages $\{V_\alpha\}$ for a bias voltage $V_{\rm bias}$ that drives an electronic current $I$, see Fig.~\ref{fig:hall_bars_setup_A}. We obtain
\begin{align}
\hspace{-0.04cm}
\left(
 V_1,  V_2,  V_3,  V_4,  V_5,  V_6 \right)
=
\left(
 1, 1, 1,  0, 0, 0\right)\frac{h}{2Ne^2}I,
\label{voltagesqh}
\end{align}
where $I=(2Ne^2/h)V_{\rm bias}$.
As expected, the ``top'' terminals $2$ and $3$ are in equilibrium with the source terminal $1$, the ``bottom'' terminals $4$ and $5$ are in equilibrium with the drain terminal $6$, and there is no voltage drop through the propagation direction ($V_2=V_3, V_4=V_5$) \cite{Datta1995}.

The longitudinal resistance \cite{Buttiker1986,Datta1995,Ihn2010} $R_{xx} \equiv R_{16,45} \equiv \left|V_2-V_3 \right|/I  = R_{16,23} = 0$ vanishes while the transverse resistance is $R_{xy}  \equiv R_{16,35} \equiv \left|V_3-V_5 \right|/I = R_{16,24} = h/2Ne^2$.
This is the standard B\"uttiker picture to describe IQH transport measurements \cite{Buttiker1986}.

Additionally, we calculate the contact resistance $R_c$ by subtracting the longitudinal resistance $R_{xx}$ from the resistance $R_{16,16}$, namely $R_{c} \equiv R_{16,16} - R_{xx} = h/2Ne^2$.
The contact resistance $R_c$ includes contributions from the system coupling to both terminals $1$ and $6$.
Each terminal offers a contact resistance to the established current $I$.
Thus, in the IQH regime the resistance $R_{16,16}=h/2Ne^2$ corresponds solely to the resistance generated by $2N$ propagating modes at the contacts \cite{Datta1995,Ihn2010}.

Let us use the B\"uttiker voltage probe model \cite{Buttiker1986b} to briefly discuss decoherence effects and why the longitudinal resistance vanishes in the IQH regime. The introduction of an extra terminal $\varphi$ to a $\Lambda$-terminal system can be trivially accounted for by writing Eq.~\eqref{current} as
\begin{align}
I_\alpha &= \mathcal G_{\alpha\varphi} \left(V_\alpha-V_\varphi \right) + \sum_{\beta=1}^{\Lambda} \mathcal G_{\alpha\beta} \left(V_\alpha-V_\beta \right),
\label{currentalpha}\\
I_\varphi &= \sum_{\beta=1}^{\Lambda} \mathcal G_{\varphi\beta} \left(V_\varphi-V_\beta \right).
\label{currentvarphi}
\end{align}
Assuming that terminal $\varphi$ is a voltage probe ($I_\varphi=0$), one can calculate $V_\varphi$ using Eq.~(\ref{currentvarphi}) and substitute the result into  Eq.~(\ref{currentalpha}) to obtain \cite{Datta1995}
\begin{align}
I_\alpha &= \sum_{\beta=1}^{\Lambda} { \mathcal G}_{\alpha\beta}^{\rm eff} \left(V_\alpha-V_\beta \right) ,
\end{align}
where
\begin{align}
{\mathcal G}_{\alpha\beta}^{\rm eff} \equiv \mathcal G_{\alpha\beta} + \frac{\mathcal G_{\alpha\varphi}\mathcal G_{\varphi\beta}}{\sum_{\beta=1}^\Lambda \mathcal G_{\varphi\beta}}.
\label{conductanceeff}
\end{align}

For simplicity, let us discuss the influence of a voltage probe on a two-terminal setting ($\Lambda=2$), with ``left'' and ``right'' leads.
In the IQH regime for $N=1$,
the conductance  reads  $\mathcal G_{LR}=\mathcal G_{RL}=2e^2/h$.
The addition of an extra terminal $\varphi$ interrupts the direct flow between terminals $L$ and $R$, leading to the non-vanishing conductance elements
$\mathcal G_{L\varphi}=\mathcal G_{\varphi R}=\mathcal G_{RL}=2e^2/h$.
The extra terminal changes $\mathcal G_{LR}$ from $2e^2/h$ to $0$ and preserves $\mathcal G_{RL}$.
On the other hand,
the effective conductance in Eq.~(\ref{conductanceeff}) reads ${\mathcal G}_{LR}^{\rm eff}={\mathcal G}_{RL}^{\rm eff}=2e^2/h$.
Thus, the introduction of an extra terminal does not change the effective conductance of the system.

Usually a reservoir incoherently populates the system channels producing additional resistance \cite{Protogenov2013, Datta1995}.
But in the IQH regime,
the chiral nature of the states forces the extra terminal to collect the entire electronic flow and inject it back.
Thus, decoherence processes do not introduce momentum relaxation here, leading to the absence of longitudinal resistance.

Let us now qualitatively discuss the local conductances and currents.
Figure~\ref{fig:hall_bar_QH_local_i_neq} sketches the nonequilibrium local electronic current, Eq.~(\ref{iklneq}) and \eqref{conductancelocal}, for  the voltages given by Eq.~(\ref{voltagesqh})
For $N=1$, the value of the local conductance at each edge is $2e^2/h$
while the voltages at terminals $1$, $2$ and $3$ correspond to $hI/2e^2$, leading to an electronic current $I$ flowing through the system upper edges.
In contrast, since $V_4=V_5=V_6=0$, there is no electronic current flowing at the bottom edge of the Hall bar.

\begin{figure}[htbp]
	\centering
		\includegraphics[width=0.75\columnwidth]{fig04.pdf
	\caption{
		Nonequilibrium component of the local electronic current map, $I_{\nu'\nu}^{\rm neq}$, estimated using Eq.~(\ref{iklneq}).
		Since the states are spin degenerate both spin components are sketched together for simplicity.
	}
		\label{fig:hall_bar_QH_local_i_neq}
\end{figure}

In turn, the local {\it equilibrium} current map is obtained from Eq.~(\ref{ikleq}).
As discussed in Sec.~\ref{sec:toymodel}, the magnetic field breaks time reversal symmetry causing a current imbalance at the central branch of the Hall bar, the sum of transmissions in Eq.~(\ref{ikleq}) does not vanish locally resulting in a clock-wise circulating electronic current.
Since Eq.~(\ref{ikleq}) involves an energy integration over all filled states,  $I_{\nu'\nu}^{\rm eq}$ contains contributions from both edge states and bulk ones.
We discuss these features in the next subsection for a given model system.

\subsection{Local currents in Quantum Hall systems: Numerical results}
\label{sec:QH_numerics}

Let us now present a quantitative analysis of the local currents.
For that purpose we consider, as an example, a graphene system modeled by a single-orbital nearest  neighbor tight-binding Hamiltonian, namely \cite{CastroNeto2009}
\begin{align}
	H = -\sum_{\left<i,j\right>,\sigma} t_{ij} c_{i\sigma}^\dagger c_{j\sigma},
	\label{hamiltonianqh}
\end{align}
where $c_{i\sigma}^\dagger$($c_{i\sigma}^{}$) is the operator that creates (annihilates) an electron at the $i$th site of a honeycomb lattice with spin projection $\sigma$, $t_{ij}=te^{i\varphi_{ij}}$ where $t= 2.7$ eV and $\varphi_{ij}$ is the Peierls phase discussed in Sec.~\ref{sec:toymodel} that accounts for a constant magnetic field ${\bf B}=B\hat{\bf z}$, and $\left< \cdots \right>$ restricts the summation to nearest neighbors sites. This simple Hamiltonian is very successful in describing the low energy properties of monolayer graphene samples \cite{CastroNeto2009}.

The Landau levels energies $E_{N_L}$ follow the analytical approximation $E_{N_L}=E_1\sqrt{N_L}$ \cite{CastroNeto2009}, where $E_1=\sqrt{3/2}a t/\ell_B$ and $\ell_B=\sqrt{\hbar/eB}=\sqrt{A_H/2\pi(\Phi/\Phi_0)}$ is the magnetic length.
Here, $\Phi \equiv BA_H$ is the magnetic flux enclosed by a single hexagon of the honeycomb lattice of area $A_H=a^2\sqrt{3}/2$, the magnetic quantum flux is $\Phi_0\equiv h/e$ and $a=2.46$ \AA.
The deviation of the numerical value of $E_N$ from its analytical prediction is small as long $|E_N|\ll t$ \cite{Lima2018}.
We set $\Phi=0.01$, so that $E_1=0.330t$.
The Hall bar has armchair (zigzag) edges along the horizontal (vertical) direction, which is about $900$ \AA ~($530$ \AA) long.

\begin{figure}[th]
	\centering
		\includegraphics[width=0.99\columnwidth]{fig05a.pdf}
		\includegraphics[width=0.99\columnwidth]{fig05b.pdf}
	\caption{
		Nonequilibrium local electronic charge current $\widetilde{\bf I}^{\rm neq}$ in the IQH regime calculated numerically using Eq.~(\ref{iklneq}) for (a) $E_F=0.333t$ (b) $E_F=0.4t$.
		The current is in units of $e^2V_{\rm bias}/h$ and includes both spin orientations.
		Positive values (blue) indicate current flowing towards the arrow, while negative values (red) indicate current flowing opposite to the arrow.
	}
		\label{fig:IQH_neq}
\end{figure}

Figure~\ref{fig:IQH_neq} shows that the {\it nonequlibrium} local charge current of a graphene Hall bar in the IQH regime for $E_F = 0.333 t$ and  $E_F = 0.4 t$ follows the behavior  predicted in Sec.~\ref{sec:QH_general}, see Fig.~\ref{fig:hall_bar_QH_local_i_neq}.
The local electronic current has contributions from terminals $1$, $2$ and $3$ that are equilibrated at the same voltage.
Here, terminals $4$, $5$ and $6$ play no role in Eq.~(\ref{iklneq}) since their voltages vanish.
The resulting current map shows an electronic flow from source (left) to drain (right) that runs through the upper edges, as expected from the qualitative discussion of the previous subsection \cite{Buttiker1988}.

The insets of Fig.~\ref{fig:IQH_neq}(a) show that for $E_F = 0.333 t$ the electronic current has two counter propagating components due to interference effects caused by scattering at the corners between zigzag and armchair edges in the middle of the Hall bar.
The component flowing towards the drain lies further away from the edges while the backscattered component lies in the edge vicinity.
This simple picture is ratified by noticing that the interference effect is energy dependent.
Changes on the electronic energy $E_F$ that lead to variations of the electron group velocity, modify the interference pattern.
We find (not shown here) that this effect is stronger for $E_F$ in the vicinity of the Landau levels energies $E_{N_L}$, where the group velocity displays an enhanced dependence on $E_F$. We show in Fig.~\ref{fig:IQH_neq}(b) that for $E_F = 0.4 t$ the counter propagating flow of electrons disappears because the ration between the Fermi wavelength and the distance between the corners changed.
Notice that the number of propagating channels is $3$ and the current integrated over a vertical cross section has the same value for both cases.

\begin{figure}[ht]
	\centering
		\includegraphics[width=0.99\columnwidth]{fig06.pdf}
	\caption{
		Equilibrium local electronic charge current $\widetilde{\bf I}^{\rm eq}$ in the IQH regime calculated numerically using Eq.~(\ref{ikleq}).
		The current is in units of $et/h$ and includes both spin orientations.
		The arrows indicate the color code for the electronic propagation direction. Propagation towards the right and top directions have positive values (blue) while propagation towards left and bottom have negative values (red).
	}
		\label{fig:IQH_eq}
\end{figure}

Figure~\ref{fig:IQH_eq} shows the local {\it equilibrium} current map, obtained by integrating Eq.~(\ref{ikleq}) from $E=0.300t$ up to $E_F=0.400t$ to capture the transition though the Landau level $N=1$, that is, the transition between $1$ to $3$ propagating edge modes per spin \cite{Lima2018}.
As predicted by the qualitative analysis of Sec.~\ref{sec:QH_general},
the equilibrium electronic current flows clockwise at the system edges.
The equilibrium local current has contributions from edge states with energies above and below $E_1$.
The superposition of different interference effects from different energies result in an additional counter-clockwise electronic current closer to the system center.
The resulting pattern of counter propagating components persists and follows the same behavior seen in Fig.~\ref{fig:IQH_neq}(a).

Figure~\ref{fig:IQH_eq} shows that the equilibrium current flows parallel to the system-leads interfaces.
This is a consequence of the mismatch between the system modes and the leads modes (doped at a high density of states energy $E_F=-t$).
The current injected by a given lead enters the system through all sites at the corresponding lead-system interface.
The magnetic field inside the central region forces the current to flow parallel to the system-lead interface towards the edge.
The superposition of the injection by all six terminals in the system renders the electronic flow that is parallel to the system-leads interfaces as seen in Fig.~\ref{fig:IQH_eq}.

We recall that the local {\sl equilibrium} electronic current integrated over any cross section of the system vanishes, since in the absence of bias the net current thought the system must be zero.
We have integrated $\widetilde{\bf I}^{\rm eq}$ over several cross sections and verified that, within numerical precision, this statement is correct.

We conclude this section by stressing that equilibrium charge currents cannot be assessed by nonequilibrium transport measurements, since they inevitably involve dissipation processes. However, since the equilibrium local currents in the IQH regime are chiral, they produce a magnetic field that can in principle be measured.
Having obtained $\widetilde{I}^{\rm eq}_{ij}$ for a given system, the change in magnetic field $\Delta B ({\bf r})$ can be calculated using the Biot-Savart law, as discussed in Ref.~\cite{Chen2020}. Experiments can follow the reverse path, for instance, use SQUIDs to measure $\Delta B ({\bf r})$
\cite{Uri2020} and infer the local equilibrium current.

\section{Local currents in Quantum spin Hall systems}
\label{sec:QSH}

In this section we study the electronic and spin transport of multi-terminal systems
in the QSH regime.
As above, for simplicity we consider the zero temperature limit.
For Fermi energies within the topological gap, the non-equilibrium electronic transport in a Hall bar geometry is ruled by helical states propagating at the system edges and the conductances ${\cal G}_{\alpha\beta} $ are quantized \cite{Hasan2010,Qi2011,Bernevig2013}.
We show that the equilibrium currents do not follow this simple picture.

\subsection{General discussion}
\label{sec:QSH_general}

The conductance matrix elements $\mathcal G_{\alpha \beta}$ of an electronic system in the QSH regime at $T=0$ in a Hall bar configuration (see Fig.~\ref{fig:hall_bars_setup_A}) reads
\begin{align}
	\boldsymbol{\mathcal{G}} =
	\frac{e^2}{h}
	\begin{pmatrix}
-2&1 &0 &1 &0 &0\\
1 &-2&1 &0 &0 &0\\
0 &1 &-2&0 &0 &1\\
1 &0 &0 &-2&1 &0\\
0 &0 &0 &1 &-2&1\\
0 &0 &1 &0 &1 &-2
\end{pmatrix}
\label{conductanceqsh}
\end{align}
The quantized entries reflect the fact that in the QSH regime each terminal $\alpha$ injects two modes that propagate towards opposite edges.

We use Eq.~(\ref{conductanceqsh}) to solve Eq.~(\ref{current}) for the unknown voltages to obtain
\begin{align}
\hspace{-0.04cm}
\left(
  V_1,
  V_2,
  V_3,
  V_4,
  V_5,
  V_6
\right)
=
\left(
 \frac{3}{2},
 1,
 \frac{1}{2},
 1,
 \frac{1}{2},
0
\right) \frac{h}{e^2}I,
\label{voltagesqsh}
\end{align}
where $I=(2e^2/3h)V_{\rm bias}$, in agreement with previous works \cite{Protogenov2013}.

The QSH transverse resistance is $R_{xy} =0$.  Here, since the voltage probes are not spin resolved, the transverse charge resistance vanishes
since there is a symmetry between the electronic propagation at the ``upper" and ``lower" edges of the Hall bar.
Thus, there is no charge imbalance between the transverse terminals.
In turn, the longitudinal resistance is  $R_{xx}=h/2e^2$ and the contact resistance is $R_c=R_{16,16}-R_{xx}=h/e^2$. As expected for two resistors in series, $R_c$ is twice the resistance corresponding to a terminal with two perfectly conducting propagating modes in a channel.

Here, in distinction to the IQH regime, dephasing processes can cause momentum relaxation and possibly explain why it has been so difficult to observe a perfect conductance quantization in QSH systems.
Although already addressed in the supplemental material of Ref. \cite{Roth2009} and by Refs. \cite{Protogenov2013, Mani2016}, this issue has been mostly overlooked.
As in Sec.~\ref{sec:QH_numerics}, we discuss dephasing effects using the phenomenological B\"uttiker voltage probe model.
In IQH systems, a voltage probe does not inject electrons back to the same edge state from which they have been drained, preventing interference effects between injected and drained electrons.
In distinction, in the QSH regime, the voltage probe injects electrons at both the clockwise and anti-clockwise helical propagating edge states, promoting interference between counter propagating helical edge modes with opposite spins orientations.

Using an analogous three-terminal system \cite{Buttiker2009} with a strongly coupled voltage probe, $\mathcal G_{\alpha\varphi}=2e^2/h$, we find that the QSH effective conductance, Eq.~(\ref{conductanceeff}), drops from $2e^2/h$ to $3e^2/2h$.
Here the voltage probe interferes with one of the edge states carrying half of the current between the original two terminals,
leading to an effective resistance of $2h/3e^2$ \cite{Protogenov2013,Mani2016}.
The supplemental material of Ref.~\cite{Roth2009}  studies weakly coupled voltage probes $\mathcal G_{\alpha\varphi} \ll e^2/h$ and shows the disappearance of perfect conductance, a hallmark of the QSH, with increasing dephasing.

Let us now qualitatively discuss the non-equilibrium currents.
We take $E_F$ within the topological gap.
In this case, only helical edge states are responsible for the electronic transport.

Using Eq.~\eqref{iklneq}, we build the local current map shown in Fig.~\ref{fig:hall_bar_QSH_local_i} by weighting each local conductance map by its respective voltage, see Eq.~(\ref{voltagesqsh}).
At the edge between terminals $1$ and $2$, only the local conductances from those terminals are not zero, with values $e^2/h$ per channel.
Since $V_1=(3/2)(hI/e^2)$ and $V_2=(hI/e^2)$ Eq.~(\ref{iklneq}) gives a net local current of $I/2$ propagating from terminal $1$ to terminal $2$ with spin up orientation.
Repeating this operation locally over the whole system, we find that (i) all edges have current $I/2$, (ii) terminals $1$ and $6$ are indeed injecting and receiving a total current $I$ and (iii) terminals $2$ through $5$ are voltage probes with zero net current.
Figure~\ref{fig:hall_bar_QSH_local_i} shows that there is no net spin current flowing from terminals $1$ to $6$, as expected in the absence of a spin polarized bias.

\begin{figure}[tbp]
	\centering
		\includegraphics[width=0.75\columnwidth]{fig07.pdf}
	\caption{Local electronic current map estimated using Eq.~\eqref{iklneq}. Here, the current $I=(2e^2/3h)V_{\rm bias}$ is established as a response of the bias $V_{\rm bias}$ for $E_F$ inside the topological gap.}
	\label{fig:hall_bar_QSH_local_i}
\end{figure}

Let us now consider equilibrium currents.
We base our discussion on
the local transmissions maps defined in Eq.~(\ref{localtcc}) for energies $E$ within the topological gap and postpone the quantitative analysis of the local currents, that involve, an integration of the whole band, Eq.~(\ref{ikleq}), to the next subsection.
Figure~\ref{fig:hall_bar_QSH_local_t} shows the local charge and spin transmission $\widetilde T_{i\sigma,j\sigma}(E)$.
The solid blue (dashed red) arrows indicate spin up (down) propagation.
Each terminal injects two modes with opposite spin orientations that propagate at opposite edges towards the neighboring terminal.

\begin{figure}[tbp]
	\centering
		\includegraphics[width=0.49\columnwidth]{fig08a.pdf}
		\includegraphics[width=0.49\columnwidth]{fig08b.pdf}
		\includegraphics[width=0.49\columnwidth]{fig08c.pdf}
		\includegraphics[width=0.49\columnwidth]{fig08d.pdf}
		\includegraphics[width=0.49\columnwidth]{fig08e.pdf}
		\includegraphics[width=0.49\columnwidth]{fig08f.pdf}
	\caption{ %
	Sketches of local transmissions $\widetilde T_{i\sigma,j\sigma}(E)$ for a system in the QSH regime.
	Panels (a) to (f) correspond to an electronic injection from the terminal $\alpha=1$ to $6$, respectively.
	}
		\label{fig:hall_bar_QSH_local_t}
\end{figure}

Notice that the edge between terminals $1$ and $2$ supports an electronic spin up propagation flowing from $\alpha=1$ to $2$ [Fig.~\ref{fig:hall_bar_QSH_local_t}(a)] and spin down propagation from $2$ to $1$ [Fig.~\ref{fig:hall_bar_QSH_local_t}(b)].
Thus, the net charge transmission vanishes.
The superposition of the transmission maps results in a zero net charge transmission, since at all edges of the Hall bar there are counter propagating electronic states that cancel each other.
Therefore, the equilibrium contribution to the local charge current in Eq.~(\ref{ikleq}) vanishes.
In turn, our analysis suggests the presence of finite equilibrium edge spin currents circulating in the system.

Our qualitative analysis shows that the local transmissions maps obtained considering a single terminal injections can be deceiving.
Only by taking into account all terminal injections one obtains correct results.
The pictures usually drawn in experimental studies of spin and valley Hall effects
\cite{Mihajlovic2009,Balakrishnan2014,Avsar2014,Shimazaki2015,Katosh2015,Park2017,Safeer2019,Herling2020}
can be misleading as they indicate that the electrons injected from the source terminal flow to neighboring terminals (depending on the spin or the valley degree of freedom) other than the drain terminal. This would imply breaking of charge conservation.
Once the bias establishes the source and drain terminals, the electrons enter the system from the source terminal, travel through a path determined by their interaction with the material, and leave the system at the drain terminal irrespective of their degrees of freedom.
One the other hand, the transmission between terminals it is indeed affected by electronic degrees of freedom.
The local transmission $\widetilde {\mathcal T}^{\beta}_{kl}$, often confused as the actual current path, indicates the local current path only in the case where the electrons are injected from a single terminal $\beta$ ($V_\beta \neq 0$) and all the other terminals $\alpha$ ($\alpha\neq\beta$) are grounded ($V_\alpha = 0$), see Eq.~(\ref{iklneq}).
In this case, the electrons can propagate to different terminals depending, for instance, on their spin, valley, and/or orbital quantum numbers.
Thus, one must not confuse the current exit terminal driven by the bias voltage, which does not depend on the electronic degrees of freedom, with the transmission exit terminal, which indeed depends on the electronic degrees of freedom.

\subsection{Numerical results}
\label{sec:QSH_numerics}

For the numerical analysis of the local currents in the QSH regime we use the Kane-Mele model  \cite{Kane2005},
that describes by a tight-binding Hamiltonian the low energy electrons in graphene
with strong spin-orbit interaction in the central region, namely
\begin{align}
	H = -\sum_{\left<i,j\right>,\sigma} t_{1} c_{i\sigma}^\dagger c_{j\sigma} + \sum_{\left<\left<i,j\right>\right>,\sigma,\sigma'} it_2\nu_{ij}S^z_{\sigma\sigma'} c_{i\sigma}^\dagger c_{i\sigma'}.
	\label{hamiltonianresults}
\end{align}
The first term accounts for nearest neighbors hopping processes where $t_{1}=t=2.7$ eV \cite{CastroNeto2009}.
The second term includes the spin dependent hopping amplitude between second neighbors where $t_2=0.065t$.
The presence of spin orbit coupling creates a topological gap $\Delta_{SO}=6\sqrt{3}t_2$ \cite{Kane2005}
such that the bulk bands onsets occur at the energies $\pm\Delta_{SO}/2=\pm 0.338t$.
Here $\mathbf S^z$ is the Pauli matrix $z$ component and $\nu_{ij}=+1$ ($\nu_{ij}=-1$) if the path from $i$ to $j$ follows the counterclockwise (clockwise) direction with respect to the hexagon centers of the honeycomb lattice.
The graphene Hall bar is attached to six semi-infinite leads, where the spin-orbit coupling is turned off.
The leads are doped at $E=-t$ to maximize the density of propagating modes.

In order to compute both non-equilibrium and equilibrium currents,
we first employ the multi-probe recursive Green's function method \cite{Lima2018} to calculate the local spin resolved transmission coefficients given by Eq.~(\ref{localtcc}).

Figure~\ref{fig:QSH_neq} shows the nonequilibrium component of the local current, Eq.~(\ref{iklneq}).
The electronic levels are assumed to be filled up to $E_F=-0.300t$.
We plot the local charge current driven by $V_{\rm bias}$ (including both spin orientations).
As expected from the discussion in Sec.~\ref{sec:QSH}, see Fig.~\ref{fig:hall_bar_QSH_local_i}, the electronic transport from source (left) to drain (right) occurs through both top and bottom edge states.
For the chosen model parameters, the top (bottom) edge carriers has $\sigma=1$ ($\sigma=-1$).

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.99\columnwidth]{fig09.pdf}
\caption{Nonequilibrium component of the local electronic charge current in the QSH regime calculated
numerically using Eq.~(\ref{iklneq}) {for $E_F=-0.3t$}.
The current is in units of $e^2V_{\rm bias}/h$ and includes both spin orientations.
The color code is the same as in {Fig.~\ref{fig:IQH_eq}}.
}
	\label{fig:QSH_neq}
\end{figure}

The computation of the equilibrium currents is more involved, since it requires an energy integration of $\widetilde T_{i\sigma,j\sigma'}^\alpha$ over all occupied states.
A recent paper \cite{Chen2020} has shown that the spin equilibrium currents vanish for Fermi energies at the Dirac point for the BHZ model, a feature that could simplify the integration in Eq.~\eqref{ikleq}. Unfortunately, our calculations do not show vanishing equilibrium bond spin currents at $E_F = 0$ for the Kane-Mele model.

According to the model parameters, for electronic energies $E<-0.338t$, the propagation is dominated by bulk states, while for $E>-0.338t$ the propagation occurs via edge states.
We chose to integrate Eq.~(\ref{ikleq}) from $E=-0.400t$ up to $E_F=-0.300t$ to capture information about from both edge and bulk propagation into the local equilibrium current.

Analogously to the absence of a net charge current due to counter propagating helical states, the bulk contribution to the equilibrium local charge current also vanishes due to the symmetric counter propagating states (not shown here).
On the other hand, a net spin current is present as we show in Fig.~\ref{fig:QSH_eq}.
The equilibrium spin current, Eq.~\eqref{eq:I_spin_local}, is defined as the difference between the spin up and the spin down components of the equilibrium local charge current, namely,
\begin{align}
\widetilde{I}_{ij}^{s,{\rm eq}} = \frac{\hbar}{2e}\left[\widetilde{I}_{i\uparrow,j\uparrow}^{\rm eq}-\widetilde{I}_{i\downarrow,j\downarrow}^{\rm eq}\right].
\end{align}
We find that even the equilibrium spin current in Fig.~\ref{fig:QSH_eq} has a counter propagating structure along the horizontal direction.
The counter propagating structure originates from the superposition between high intensity edge propagation ($E>-0.338t$) and bulk propagation ($E<-0.338t$) that has lower intensity near the edges since it spreads across the cross section.
Interference effects are also present as discussed using Figs.~\ref{fig:IQH_neq} and \ref{fig:IQH_eq} in the IQH regime.
But here, their dependence with the electronic energy is weak, since the dispersion of the helical edge states are linear, that is, the group velocity does not vary inside the topological gap.
Similarly to the IQH case, the equilibrium spin current shown in Fig.~\ref{fig:IQH_eq} flows clockwise at the system edges and has an additional counter-clockwise spin current closer to the system center.
The helical states lying inside the topological gap dominate the transport near the edges.

\begin{figure}[htbp]
	\centering
		\includegraphics[width=0.99\columnwidth]{fig10.pdf}
	\caption{
	Equilibrium component of the local electronic spin current in the QSH regime calculated numerically using Eq.~(\ref{ikleq}) in units of $t$.
	The color code is the same as in {Fig.~\ref{fig:IQH_eq}}.
	}
	\label{fig:QSH_eq}
\end{figure}

Note that the nonequilibrium current in both QSH and IQH regimes behaves differently at their terminal interfaces.
In the QSH regime, Fig.~\ref{fig:QSH_neq}, the current flowing through the central region crosses the system-terminal interface via the edges.
It enters the leads via one edge it is injected back again into the system via the opposite edge.
In this case, the interface has low resistance, so that all the electronic current crosses the interface at the point it arrives.
The change in its direction changes inside the leads where the QSH regime fades away and allows for bulk propagation connecting opposite edges.
That is not the case in the IQH regime.
In Fig.~\ref{fig:IQH_neq}, we clearly see current at the system-terminal interface connecting opposite edges.
Here, the interface has high resistance, so that only a fraction of the electronic current is absorbed at the point where it arrives and the remaining fraction continuous to travel in the system.
As a result, the change in its direction happens partially inside the system it self and inside the leads as well.
We stress that there is no conservation rule requiring the current to be entirely absorbed by any voltage lead.
It is only required that the current collected by a voltage lead must be injected back.
In an experiment, the current will mostly likely turn back to the system inside the leads where the metal constituting the leads most often does not share the same properties, such as QSH or IQH regimes, of the studied system.

We conclude by stressing that, like in the IQH regime, the local equilibrium spin currents $\widetilde{I}^{s,{\rm eq}}_{ii'}$ cannot be detected by nonequilibrium transport experiments. Yet, one can envisage strategies to indirectly measure equilibrium currents, for instance, $\widetilde{I}^{s,{\rm eq}}_{ii'}$ carries a magnetic moment that gives origin to an electric dipolar field \cite{Chen2020}, which can in principle be measured.

\section{Summary and conclusions}
\label{sec:conclusions}

We reviewed the Landauer-B\"uttiker formalism in terms of non-equilibrium Green's functions applied to multi-terminal systems and derived general expressions that allows one to compute local electronic currents in mesoscopic systems at both the equilibrium and non-equilibrium regimes.
We show that general expressions we put forward for the local electron current can be projected into a suitable basis allowing the calculation of the spin bond current. By analogy, this approach can be extended to obtain further transport properties, such as equilibrium and nonequilibrium valley and orbital local currents.

We demonstrate that in the absence of a bias voltage, the net electronic current through the system vanishes, but the local (equilibrium) current is not necessarily zero.
To illustrate this effect, we used a toy model with minimal elements, namely, a central region lattice with bonds forming a loop and a magnetic field breaking the time reversal symmetry.

Motivated by recent experiments \cite{Uri2020}, we further discuss the properties of equilibrium and non-equilibrium (linear response) local currents in the {IQH and the QSH} regime for realistic Hall bar setups. As a case to point, we consider graphene systems.
In the IQH regime, a non-vanishing local charge current is driven by the magnetic field induced broken time reversal symmetry.
The equilibrium local charge currents circulate without leaving the system.
Any transverse section of the system has net equilibrium charge current equal to zero, as expected.
We revisited the decoherence process introduced by a voltage probe to justify the absence  of momentum relaxation that leads to the absence of longitudinal resistance.

In the QSH regime, the counter-propagating helical edge states interact at the system-terminal interfaces, forcing momentum relaxation and the appearance of a finite longitudinal resistance.
At equilibrium, the local charge current vanishes due to time reversal symmetry, but the spin-orbit interaction gives rise to a finite local spin current that circulates the system with no net current flowing through any cross section of the system.

Our results are nicely interpreted by the B\"uttiker picture \cite{Datta1995} that consideres the injection of electrons by all terminals for electron energies ranging over the entire band, weighted by the Fermi distribution. In equilibrium, all
terminal attached to the system simultaneously inject the same current, leading to a zero {\it net} equilibrium electronic current.
The presence of an external magnetic field breaks time reversal symmetry and gives rise to non-vanishing dissipationless equilibrium {\it local} charge currents, while spin-orbit interactions {originate} non-vanishing dissipationless equilibrium local spin currents with zero local charge current.
Let us stress that despite not being detectable by nonequilibrium transport measurements, equilibrium currents can be assessed indirectly, as discussed in Sec. \ref{sec:IQH} and \ref{sec:QSH}.

Finally, we emphasize the distinction between transmission, conductance and current.
Thus, we learn that local transmissions maps obtained
considering single terminal injections can be deceiving
\cite{Mihajlovic2009,Balakrishnan2014,Avsar2014,Shimazaki2015,Katosh2015,Park2017,Safeer2019,Herling2020}.
Only by taking into account all terminal injections one obtains meaningful results.

\acknowledgments

We thank B. Nikoli\v{c} for enlightening discussions at the early stage of this work. We acknowledge financial support of the Brazilian Institute of Science and Technology (INCT) in Carbon Nanomaterials and the Brazilian agencies CAPES, CNPq, FAPEMIG, and FAPERJ.