# Robust nonequilibrium edge currents with and without band topology

Mark T. Mitchison<EMAIL_ADDRESS>School of Physics, Trinity College
Dublin, College Green, Dublin 2, Ireland  Ángel Rivas<EMAIL_ADDRESS>Miguel
A. Martin-Delgado<EMAIL_ADDRESS>Departamento de Física Teórica, Facultad de
Ciencias Físicas, Universidad Complutense, 28040 Madrid, Spain.  CCS-Center
for Computational Simulation, Campus de Montegancedo UPM, 28660 Boadilla del
Monte, Madrid, Spain.

###### Abstract

We study two-dimensional bosonic and fermionic lattice systems under
nonequilibrium conditions corresponding to a sharp gradient of temperature
imposed by two thermal baths. In particular, we consider a lattice model with
broken time-reversal symmetry that exhibits both topologically trivial and
nontrivial phases. Using a nonperturbative Green function approach, we
characterize the nonequilibrium current distribution in different parameter
regimes. For both bosonic and fermionic systems, we find chiral edge currents
that are robust against coupling to reservoirs and to the presence of defects
on the boundary or in the bulk. This robustness not only originates from
topological effects at zero temperature but, remarkably, also persists as a
result of dissipative symmetries in regimes where band topology plays no role.
Chirality of the edge currents implies that energy locally flows against the
temperature gradient without any external work input. In the fermionic case,
there is also a regime with topologically protected boundary currents, which
nonetheless do not circulate around all system edges.

The physics of boundary structures, whether they be dots, lines or surfaces,
has attracted a great deal of attention from various directions in the past.
The theoretical discovery, and its subsequent experimental verification, of
both insulating and superconducting topological materials [1, 2, 3, 4, 5, 6]
has further spurred the study of boundary physics. The main reason for this is
that nontrivial band topology endows edge phenomena with a remarkable
robustness. This feature opens up a myriad of possible applications that go
well beyond condensed matter physics. Historically, most of the mainstream
studies conducted on those topological materials shared two basic properties
in common: (i) the quantum system is considered as closed and thus isolated
from the detrimental effects of the surrounding environment; (ii) the
constituent particles are fermions. The combination of these two properties
underpins the stability of boundary effects in topological phases of matter.

In this work, we present a different paradigm of robust boundary physics in
which we depart from these two common features, entering less well-trodden
ground. The motivation is to study more demanding scenarios corresponding to a
quantum system coupled to thermal baths, which generate external noise
yielding fluctuations and dissipation [7, 8, 9, 10, 11]. In fact,
understanding noisy circumstances like these is crucial for the successful
development of scalable quantum technologies. In this context, we find that
robust edge currents can be generated without resorting to the standard band
topological mechanism, focusing our attention not only on fermionic systems,
but on bosonic particles as well.

To illustrate this unusual form of dissipative robustness, we consider a
bosonic variant of the model previously introduced by Qi, Wu and Zhang (QWZ)
for fermions on a square two-dimensional (2D) lattice [12]. This model has the
virtue of presenting two topologically different band structures depending on
the values of its coupling parameters. In the fermionic case at half filling,
these two correspond to different phases: one is a trivial insulator, and the
other a topological insulator. Of course, the single-particle band structure
of the QWZ model is the same for fermions and bosons. However, the statistics
of the particles determine how those bands are filled. The Pauli principle
forces fermions to fill the bands up to the Fermi level, thereby unveiling the
band topology. On the contrary, bosons at low temperatures tend to condense in
the single-particle ground state, making them largely insensitive to the
global band structure.

Figure 1: The QWZ model. (a) Phase diagram showing the topologically
nontrivial regimes with nonzero Chern number $\mathrm{n}=\pm 1$, a topological
invariant of the band structure under periodic boundary conditions [6]. (b)
Schematic of the system, featuring a 2D lattice hosting particles with two
internal states, coupled to thermal baths at different temperatures. An edge
current (red) emerges from an erasure effect where modes with the same current
circulation (black) cancel in the bulk but add constructively on the boundary.
(c) Single-particle eigenenergies $\omega_{\alpha}$ as a function of the mode
index $\alpha$ for a system with open boundary conditions, with
$L_{X}=L_{Y}=20$, $t_{X}=t_{Y}=t$.

Remarkably, for both bosonic and fermionic QWZ lattices, we find dissipatively
robust, chiral edge currents flowing between two thermal reservoirs in a
parameter regime leading to trivial band topology. Here, robustness is defined
by stability of the currents with respect to the introduction of defects on
the edges or in the bulk. This stability is not arbitrary, but is subject to
the fulfilment of discrete symmetries on the geometrical distribution of the
defects, which must be compatible with the underlying symmetries of the non-
equilibrium steady state (NESS). In the fermionic case, and when the
temperature of one of the two reservoirs is much smaller than the band gap, we
also observe topologically protected edge currents that are effectively
independent of the defect distribution but may not circulate around every edge
of the perturbed system.

The interplay between topology and symmetries in dissipative quantum systems
has recently been explored within a Markovian approximation [13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26], which generally requires weak system-
reservoir interactions (see also Ref. [27] for non-Markovian extensions in
fermionic systems). Here, we instead adopt a nonperturbative Green function
method to compute the exact NESS for arbitrary values of the system-bath
coupling, in order to monitor how the current distribution within the system
changes from strong to weak coupling regimes. Interestingly, edge currents
appear in the weak-coupling limit but are masked by bulk currents at strong
coupling. This contrasts with the common expectation that exotic thermodynamic
effects are more prone to arise in strong-coupling configurations [28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39]. Our exact analysis can also help to
identify the range of parameters in which these dissipative edge currents can
be experimentally realized. An appropriate way to achieve such realizations is
a setup with a high degree of control over microscopic degrees of freedom,
e.g., quantum simulation on platforms developed to deal with large systems
[40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51], whose fundamental
constituents may be bosons.

How can dissipatively robust edge currents arise without band topology? Let us
first explore this phenomenon and then we will explain its origin.

Model.—We consider the QWZ Hamiltonian [12] describing a collection of non-
interacting fermions or bosons with two internal “flavor” states. The
particles occupy a 2D square lattice with sites specified by the coordinates
$x=1,\ldots,L_{X}$ and $y=1,\ldots,L_{Y}$. The Hamiltonian is expressed in
terms of vectors of canonical ladder operators
${\bf\hat{a}}^{\dagger}_{x,y}=(\hat{a}^{\dagger}_{x,y,\uparrow},\hat{a}^{\dagger}_{x,y,\downarrow})$
for each lattice site, whose components create a particle with flavor
$\uparrow$ or $\downarrow$. Note that the flavors merely index distinct bands
and are unrelated to angular momentum. Explicitly, the Hamiltonian reads as
$\hat{H}=\hat{H}_{m}+\hat{H}_{X}+\hat{H}_{Y}$, with ($\hbar=1=k_{B}$)

$\displaystyle\hat{H}_{m}$
$\displaystyle=\sum_{x,y}{\bf\hat{a}}^{\dagger}_{x,y}\cdot(\omega_{0}\mathbb{1}+m\sigma_{z})\cdot{\bf\hat{a}}_{x,y},$
(1) $\displaystyle\hat{H}_{X}$
$\displaystyle=\frac{t_{X}}{2}\sum_{x,y}{\bf\hat{a}}^{\dagger}_{x+1,y}\cdot(\sigma_{z}+\mathrm{i}\sigma_{y})\cdot{\bf\hat{a}}_{x,y}+{\rm
h.c.},$ (2) $\displaystyle\hat{H}_{Y}$
$\displaystyle=\frac{t_{Y}}{2}\sum_{x,y}{\bf\hat{a}}^{\dagger}_{x,y+1}\cdot(\sigma_{z}+\mathrm{i}\sigma_{x})\cdot{\bf\hat{a}}_{x,y}+{\rm
h.c.},$ (3)

where $\sigma_{x,y,z}$ are Pauli matrices in flavor space, $\omega_{0}$ is the
on-site energy, $m$ is the flavor energy splitting, and $t_{X,Y}>0$ are the
tunnelling amplitudes in the $x$ and $y$ directions.

The QWZ Hamiltonian has two notable symmetries. The first is
$\hat{\Pi}\hat{R}_{\pi}$, a combined $\pi$-rotation about the $z$ axis in real
space,
$\hat{R}_{\pi}\hat{\mathbf{a}}_{x,y}\hat{R}_{\pi}^{\dagger}=\hat{\mathbf{a}}_{L_{X}+1-x,L_{Y}+1-y}$,
and flavor space,
$\hat{\Pi}\hat{\mathbf{a}}_{x,y}\hat{\Pi}^{\dagger}=\sigma_{z}\hat{\mathbf{a}}_{x,y}$.
The second symmetry is $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$, which combines
time reversal, $\hat{\Theta}\hat{H}\hat{\Theta}^{-1}=\hat{H}^{*}$, spatial
reflection about the $y$ axis,
$\hat{\Sigma}_{y}\hat{\mathbf{a}}_{x,y}\hat{\Sigma}_{y}^{\dagger}=\hat{\mathbf{a}}_{L_{X}+1-x,y}$,
and the flavor $\pi$-rotation $\hat{\Pi}$ defined above. We can already
anticipate that these discrete symmetries will play an important role in
stabilising edge currents out of equilibrium, as previously found in the
context of a bosonic Hofstadter model [19]. Unlike that model, however, the
QWZ Hamiltonian exhibits both topologically trivial and nontrivial phases
depending on the value of $m$ relative to $t_{X,Y}$. Nontrivial topology
manifests as a series of edge states with linear dispersion relation
connecting the two single-particle energy bands. Conversely, in the
topologically trivial regime there are no edge states and the bands are
separated by a finite energy gap. The phase diagram and corresponding band
structure is indicated in Fig. 1.

To study an out-of-equilibrium situation, we couple one side of the lattice
($x=1$) to a hot thermal bath at temperature $T_{h}=\beta^{-1}_{h}$ and the
other side ($x=L_{X}$) to a cold bath at temperature $T_{c}=\beta_{c}^{-1}$,
as depicted in Fig. 1(b). These baths are modelled by reservoirs of non-
interacting fermions or bosons, which can tunnel to and from the system via a
linear coupling. We assume that reservoir modes coupled to distinct sites of
the system are uncorrelated, and are initially populated according to the
distribution function
$\bar{n}_{h,c}(\omega)=(\mathrm{e}^{\beta_{h,c}(\omega-\mu)}\pm 1)^{-1}$,
where the plus (minus) sign pertains to fermions (bosons). In the fermionic
case, we include a chemical potential $\mu$ to fix the average density, while
in the bosonic case we set $\mu=0$. At long times, the system reaches a NESS,
which can be computed exactly [52]. The NESS is a Gaussian state and thus
fully characterized by its correlation matrix
$C_{jk}=\langle\hat{a}^{\dagger}_{k}\hat{a}_{j}\rangle$, where the indices
$j,k$ represent the coordinates $(x,y)$ as well as the flavor state.
Explicitly, we have

$\mathbf{C}=\int\frac{\mathrm{d}\omega}{2\pi}\,\mathbf{G}(\omega)\cdot\Big{[}\mathbf{\Gamma}_{h}\bar{n}_{h}(\omega)+\mathbf{\Gamma}_{c}\bar{n}_{c}(\omega)\Big{]}\cdot\mathbf{G}^{\dagger}(\omega),$
(4)

where $\mathbf{G}(\omega)$ is the retarded Green function of the system
obtained by tracing over the reservoirs and $\mathbf{\Gamma}_{c,h}$ are self-
energies describing the system-reservoir coupling [53, 54]. We work in the
wide-band limit, where the self-energies can be approximated by a frequency-
independent constant $\gamma$, which we assume to be equal for both hot and
cold reservoirs. This approximation, which is valid so long as the reservoir
spectral densities vary slowly in the relevant frequency range, significantly
simplifies the calculations but is not essential for our conclusions to hold.

The applied thermal gradient gives rise to particle currents flowing within
the system. We denote by $J^{X}_{x,y}$ the mean particle current flowing from
site $(x,y)$ to site $(x+1,y)$, while $J^{Y}_{x,y}$ denotes the current
flowing from $(x,y)$ to $(x,y+1)$. These are expectation values of one-body
observables and can be found from the NESS correlation matrix $\mathbf{C}$
[52].

Figure 2: Current distributions on a lattice of size $L=8$. Red arrows show
boson currents in a topologically trivial phase ($m=3t$), blue arrows show
fermion currents in a topologically nontrivial regime ($m=t$). Arrow sizes
indicate the magnitude of the currents scaled relative to the largest value
within each plot. (a,b) Nonequilibrium current profile with different
temperatures $T_{h}=t$ and $T_{c}=0.01t$, on-site energy $\omega_{0}=10t$,
chemical potential $\mu=\omega_{0}-0.01t$ for fermions (we always set $\mu=0$
for bosons), and coupling strength (a) $\gamma=0.5t$ and (b) $\gamma=0.005t$.
(c) Equilibrium case with temperatures $T_{c}=T_{h}=t$, weak coupling,
$\gamma=0.005t$, and other parameters identical to (a,b).

Results.—In the following examples, we focus on the symmetric case with
$L_{X}=L_{Y}=L$ and $t_{X}=t_{Y}=t$, so that $|m|<2t$ defines the
topologically non-trivial phase. We also fix $\omega_{0}=10t$ and consider
relatively low temperatures, $T_{c,h}\lesssim t$, to accentuate the role of
band topology and particle exchange statistics.

Figure 3: Average edge currents, $J_{\rm edge}$ (solid lines), and bulk
currents, $J_{\rm bulk}$ (dashed lines), as a function of the system-bath
coupling strength, $\gamma$, for bosons (left panel) and fermions (right
panel) in topologically non-trivial ($m/t=1$) and trivial ($m/t=2.1,3$)
regimes, and with $T_{h}=t$, $T_{c}=0.01t$, $\omega_{0}=10t$, and
$\mu=\omega_{0}+0.1t$.

Figs. 2(a,b) plot the nonequilibrium current distributions for two different
values of the system-reservoir coupling, $\gamma$. Red arrows show the
currents for a bosonic system with $m=3t$ (similar results are obtained for
$|m|<2t$). Remarkably, the currents become progressively localized on the
boundary of the system as $\gamma$ is reduced, even though the band topology
is trivial. These edge currents also arise in the equilibrium case,
$T_{h}=T_{c}$, as shown in Fig. 2(c). Qualitatively similar results are
obtained for fermions, as shown for a nontrivial phase ($m=t$) by the blue
arrows in Fig. 2. The direction of fermionic particle flow depends on the
chemical potential; for the parameters in Figs. 2(a,b), the thermoelectric
induced current flows in the opposite direction to the temperature gradient.
Conversely, when $\mu=\omega_{0}$, all currents vanish due to particle-hole
symmetry [52].

In order to quantify the emergence of boundary currents more precisely, we
define the average edge and bulk currents

$\displaystyle J_{\rm edge}$
$\displaystyle=\frac{1}{2L_{X}}\sum_{x}\left(J^{X}_{x,L_{Y}}-J^{X}_{x,1}\right),\qquad
J_{\rm bulk}=\frac{1}{L_{Y}}\sum_{y}J^{X}_{x,y}.$ (5)

The choice of $x$-coordinate in the definition of $J_{\rm bulk}$ is arbitrary
due to particle-number conservation; we take $x=\lfloor L_{X}/2\rfloor$. The
total current flowing between the two reservoirs is given by $J_{\rm
tot}=L_{Y}J_{\rm bulk}$.

We plot the edge and bulk currents in Fig. 3 for bosonic and fermionic systems
as a function of the system-reservoir coupling strength, $\gamma$. In the
bosonic case, both bulk and edge currents increase with $m$ and are thus
larger in the topologically trivial phase. The situation is reversed for
fermions, with the largest edge currents obtained in the topologically
nontrivial phase. In all cases, the bulk currents are proportional to $\gamma$
and thus vanish as $\gamma\to 0$, while the edge currents remains invariant
for a wide range of values of $\gamma$ and persist in the weak-coupling limit.
These quantitative results therefore confirm the qualitative picture of Fig.
2, i.e., edge currents arise for small coupling irrespectively of band
topology or particle statistics. In topologically trivial phases, edge
currents begin to dominate once $\gamma$ becomes comparable to the level
spacing of the single-particle Hamiltonian. However, in the topologically
nontrivial phase for fermions near half filling, currents remain localized on
the edges even for strong coupling, $\gamma\sim t$.

Figure 4: Effect of impurities on the current distribution in the weak-
coupling limit, with impurity positions marked by filled black circles. Red
arrows show boson currents and blue arrows show fermion currents. In all
plots, $m=t$, $T_{h}=t$, $T_{c}=0.01t$, $\omega_{0}=10t$ and
$\mu=\omega_{0}-0.1t$.

Finally, we investigate whether these effects withstand the presence of
imperfections, a property which defines the notion of dissipatively robust
currents. Specifically, we examine how the currents change when static
impurities are added to the lattice. Impurities are modelled by a large on-
site energy shift, i.e., a term
$\Delta{\bf\hat{a}}^{\dagger}_{x,y}\cdot{\bf\hat{a}}_{x,y}$ added to the
Hamitonian, where $\Delta\gg\omega_{0},t$ and $(x,y)$ are the coordinates of
the impurity site [19]. We focus hereafter on the weak-coupling limit where
the edge currents are most prominent in the absence of impurities. In this
limit, an analytical expression for the NESS [53] can be obtained, which is
identical to the solution of the Lindblad equation derived under the Born-
Markov and secular approximations [52].

In Fig. 4 we show three examples of how impurities affect the nonequilibrium
current distribution in the topologically nontrivial phase. We observe
substantial differences between bosonic and fermionic systems when one of the
reservoirs is at low temperature. The bosonic edge currents are preserved only
when the impurity distribution is invariant under one of the symmetries
$\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$ or $\hat{\Pi}\hat{R}_{\pi}$, which
correspond to the purely spatial symmetries $\hat{\Sigma}_{y}$ and
$\hat{R}_{\pi}$ since the defects satisfy $\hat{\Pi}$ and $\hat{\Theta}$
automatically. In the case of impurities placed on the edge (and assuming the
symmetries are respected in the bosonic case), the currents simply detour
around the impurity sites [Fig. 4(a)]. If impurities are placed in the bulk of
a bosonic system in an appropriately symmetric way, counter-currents shield
the impurity sites by circulating in the opposite direction to the edge
currents [Fig. 4(b)]. In the absence of these symmetries, the bosonic edge
currents are strongly disrupted [Fig. 4(c)]. For fermions in a topological
phase, the situation is starkly different when at least one reservoir is at
low temperature: impurities placed in the bulk have no effect whatsoever on
the current distribution [Figs. 4(b,c)]. In the topologically trivial phase or
at high temperature, both bosonic and fermionic edge currents enjoy the same
symmetry-protected robustness as bosons in the nontrivial phase.

Discussion.—The question posed in the introduction on the nontopological
origin of the edge currents can now be answered. In the weak-coupling limit,
nonequilibrium coherences between energy eigenstates become small [53] and
$J_{\rm bulk}$, which is directly proportional to these coherences [55], is
negligible in comparison to the contributions from individual energy
eigenstates. Moreover, the single-particle eigenmodes in the QWZ model present
a nonzero Berry curvature $\bm{\mathcal{F}}_{\alpha}$ for $m\neq 0$ [12],
which breaks time-reversal symmetry and endows the eigenmodes with a
particular chirality. This can be understood in a semi-classical picture [56],
where a wavepacket propagates with velocity
$\mathbf{v}_{\alpha}=\partial_{\mathbf{k}}\omega_{\alpha}+\nabla
V\times\bm{\mathcal{F}}_{\alpha}$, with $\mathbf{k}$ the wavevector and
$V(x,y)$ a confining potential describing the edge of the system [57].
Assuming that the non-equilibrium distribution function only depends on
energy, the net current due to the derivative term
$\partial_{\mathbf{k}}\omega_{\alpha}$ can be shown to vanish in the bulk due
to cancellations between positive and negative wavevectors [52]. However, the
contribution of the axial force $\nabla V\times\bm{\mathcal{F}}_{\alpha}$ is
non-zero near the boundary, thus generating a net chiral edge current [52].
This erasure effect arises whenever there is broken time-reversal symmetry and
a smooth, energetically monotonic distribution of chiral eigenmodes [19]; it
is thus independent of exchange statistics or band topology. Intuitively, the
effect can be understood in terms of circulating currents that cancel in the
bulk but add constructively on the boundary [Fig. 1(b)].

A chiral current between two thermal reservoirs unavoidably entails a local
cross-current phenomenon: on one edge there are particles flowing from the
cold bath to the hot bath. In bosonic systems, this implies a “violation” of
the second law of thermodynamics within a one-dimensional subsystem, i.e.,
energy flows against the temperature gradient along one edge even though the
total rate of entropy production is positive. This effect was reported in
Refs. [19, 58] for bosonic lattices governed by the Hofstadter Hamiltonian.
Now, we see that the key requisite for this behavior is not the nontrivial
band topology of the lattice, but rather the nonzero Berry curvature in a
regime of weak system-bath coupling [Fig. 2(b,c)].

The edge currents are robust to impurities if the defect distribution
satisfies either of the symmetries $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$ or
$\hat{\Pi}\hat{R}_{\pi}$. Note that the relevant nonequilibrium symmetries are
determined both by the Hamiltonian and by the configuration of the baths.
These symmetries leave the NESS invariant and lead to a steady-state
correlation matrix that is independent of the spatial orientation of the
reservoirs. Under these conditions, the nonequilibrium distribution function
is simply the average of the reservoir distributions [52], i.e.,
$n(\omega_{\alpha})=\tfrac{1}{2}\left[\bar{n}_{h}(\omega_{\alpha})+\bar{n}_{c}(\omega_{\alpha})\right]$
is the population of the eigenmode with frequency $\omega_{\alpha}$. This
smooth distribution function yields boundary currents due to the erasure
effect in the bulk. In the absence of symmetry, however, the contribution of
each mode depends not only on energy but also on the spatial profile of the
corresponding wavefunction. This creates an erratically varying distribution
function that destroys the erasure effect [52].

For fermions at low temperature $T_{c}$, $n(\omega_{\alpha}$) undergoes a
sharp change near $\omega_{\alpha}\approx\mu$. In a topologically nontrivial
phase near half filling, this feature acts as a filter that populates only one
or two edge states [52], which are topologically protected against
perturbations in the bulk. The corresponding boundary current is thus
completely unaffected by impurities [Figs. 4(b,c)]. Remarkably, this holds for
arbitrarily large $T_{h}$, providing an instance of a topologically protected
property that is also stable under dissipation. The number and character of
edge modes that contribute to the current pattern is determined by the
chemical potential [52].

In summary, our nonperturbative analysis has shown that robust boundary
currents emerge at weak coupling in the QWZ lattice driven out of equilibrium
by a thermal gradient. Our results reveal a novel dissipative mechanism for
symmetry-protected edge transport, which arises in bosonic and fermionic
systems with and without nontrivial band topology. This makes quantum
simulators based on either fermionic or bosonic degrees of freedom attractive
candidates to observe this exotic boundary physics.

Acknowledgements.—M. T. M. acknowledges funding from the ERC Starting Grant
ODYSSEY (Grant Agreement No. 758403) and the EPSRC-SFI Joint Funding of
Research project QuamNESS. A. R. and M. A. M.-D. acknowledge financial support
from the Spanish MINECO grants MINECO/FEDER Projects FIS2017-91460-EXP,
PGC2018-099169-B-I00 FIS-2018 and from CAM/FEDER Project No. S2018/TCS-4342
(QUITEMAD-CM). The research of A. R. and M. A. M.-D. has been partially
supported by the U.S. Army Research Office through Grant No. W911NF-14-1-0103.
Calculations were performed on the Lonsdale cluster maintained by the Trinity
Centre for High Performance Computing. This cluster was funded through grants
from SFI.

## References

  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Oxford, 2002).
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer, Berlin, 2004).
  * Kamenev [2009] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, Cambridge, 2009).
  * Rivas and Huelga [2012] A. Rivas and S. F. Huelga,  _Open quantum systems. An introduction_ (Springer, Heidelberg, 2012).
  * Landi  _et al._ [2021] G. T. Landi, D. Poletti, and G. Schaller, Non-equilibrium boundary driven quantum systems: models, methods and properties (2021), arXiv:2104.14350 [quant-ph] .
  * Qi  _et al._ [2006] X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Phys. Rev. B 74, 085308 (2006).
  * Diehl  _et al._ [2011] S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, Nature Phys. 7, 971 (2011).
  * Viyuela  _et al._ [2012] O. Viyuela, A. Rivas, and M. A. Martin-Delgado, Phys. Rev. B 86, 155140 (2012).
  * Rivas  _et al._ [2013] A. Rivas, O. Viyuela, and M. A. Martin-Delgado, Phys. Rev. B 88, 155141 (2013).
  * Budich  _et al._ [2015] J. C. Budich, P. Zoller, and S. Diehl, Phys. Rev. A 91, 042117 (2015).
  * Iemini  _et al._ [2016] F. Iemini, D. Rossini, R. Fazio, S. Diehl, and L. Mazza, Phys. Rev. B 93, 115113 (2016).
  * Linzner  _et al._ [2016] D. Linzner, L. Wawer, F. Grusdt, and M. Fleischhauer, Phys. Rev. B 94, 201105 (2016).
  * Rivas and Martin-Delgado [2017] Á. Rivas and M. A. Martin-Delgado, Sci. Rep. 7, 6350 (2017).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Shavit and Goldstein [2020] G. Shavit and M. Goldstein, Phys. Rev. B 101, 125412 (2020).
  * Gau  _et al._ [2020] M. Gau, R. Egger, A. Zazunov, and Y. Gefen, Phys. Rev. Lett. 125, 147701 (2020).
  * Lieu  _et al._ [2020] S. Lieu, M. McGinley, and N. R. Cooper, Phys. Rev. Lett. 124, 040401 (2020).
  * McGinley and Cooper [2020] M. McGinley and N. R. Cooper, Nature Phys. 16, 1181 (2020).
  * Flynn  _et al._ [2021] V. P. Flynn, E. Cobanera, and L. Viola, Phys. Rev. Lett. 127, 245701 (2021).
  * Altland  _et al._ [2021] A. Altland, M. Fleischhauer, and S. Diehl, Phys. Rev. X 11, 021037 (2021).
  * Iles-Smith  _et al._ [2014] J. Iles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * Esposito  _et al._ [2015] M. Esposito, M. A. Ochoa, and M. Galperin, Phys. Rev. B 92, 235440 (2015).
  * Bruch  _et al._ [2016] A. Bruch, M. Thomas, S. Viola Kusminskiy, F. von Oppen, and A. Nitzan, Phys. Rev. B 93, 115318 (2016).
  * Carrega  _et al._ [2016] M. Carrega, P. Solinas, M. Sassetti, and U. Weiss, Phys. Rev. Lett. 116, 240403 (2016).
  * Strasberg  _et al._ [2016] P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, New J. Phys. 18, 073007 (2016).
  * Newman  _et al._ [2017] D. Newman, F. Mintert, and A. Nazir, Phys. Rev. E 95, 032139 (2017).
  * Perarnau-Llobet  _et al._ [2018] M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, and J. Eisert, Phys. Rev. Lett. 120, 120602 (2018).
  * Miller and Anders [2018] H. J. D. Miller and J. Anders, Nature Commun. 9, 2203 (2018).
  * Pancotti  _et al._ [2020] N. Pancotti, M. Scandi, M. T. Mitchison, and M. Perarnau-Llobet, Phys. Rev. X 10, 031015 (2020).
  * Rivas [2020] A. Rivas, Phys. Rev. Lett. 124, 160601 (2020).
  * Talkner and Hänggi [2020] P. Talkner and P. Hänggi, Rev. Mod. Phys. 92, 041002 (2020).
  * Popovic  _et al._ [2021] M. Popovic, M. T. Mitchison, A. Strathearn, B. W. Lovett, J. Goold, and P. R. Eastham, PRX Quantum 2, 020338 (2021).
  * Aidelsburger  _et al._ [2013] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Nature 496, 196 (2013).
  * Hafezi  _et al._ [2013] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Nature Photonics 7, 1001 (2013).
  * Jotzu  _et al._ [2014] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).
  * Aidelsburger  _et al._ [2015] M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbène, N. R. Cooper, I. Bloch, and N. Goldman, Nat. Phys. 11, 162 (2015).
  * Anderson  _et al._ [2016] B. M. Anderson, R. Ma, C. Owens, D. I. Schuster, and J. Simon, Phys. Rev. X 6, 041043 (2016).
  * Khanikaev and Shvets [2017] A. B. Khanikaev and G. Shvets, Nature Photonics 11, 763 (2017).
  * Viyuela  _et al._ [2018] O. Viyuela, A. Rivas, S. Gasparinetti, A. Wallraff, S. Filipp, and M. A. Martin-Delgado, npj Quantum Info. 4, 10 (2018).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * Chalopin  _et al._ [2020] T. Chalopin, T. Satoor, A. Evrard, V. Makhalov, J. Dalibard, R. Lopes, and S. Nascimbene, Nature Phys. 16, 1017 (2020).
  * Viebahn  _et al._ [2021] K. Viebahn, J. Minguzzi, K. Sandholzer, A.-S. Walter, M. Sajnani, F. Görg, and T. Esslinger, Phys. Rev. X 11, 011057 (2021).
  * Ferri  _et al._ [2021] F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, Phys. Rev. X 11, 041046 (2021).
  * [52] See the Supplemental Material for the derivation of the exact nonequilibrium steady state, the definitions of the currents, a discussion of particle-hole symmetry, a detailed analysis of the weak-coupling limit and a semiclassical picture for the edge circulation. This includes citations to Refs. [59, 60, 61, 62].
  * Dhar  _et al._ [2012] A. Dhar, K. Saito, and P. Hänggi, Phys. Rev. E 85, 011126 (2012).
  * Ryndyk [2016] D. Ryndyk,  _Theory of Quantum Transport at Nanoscale_ (Springer, Berlin, 2016).
  * Mitchison and Plenio [2018] M. T. Mitchison and M. B. Plenio, New J. Phys. 20, 033005 (2018).
  * Haldane [2004] F. D. M. Haldane, Phys. Rev. Lett. 93, 206602 (2004).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Xing  _et al._ [2020] B. Xing, X. Xu, V. Balachandran, and D. Poletti, Phys. Rev. B 102, 245433 (2020).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Martinez Alvarez  _et al._ [2018] V. M. Martinez Alvarez, J. E. Barrios Vargas, M. Berdakin, and L. E. Foa Torres, Eur. Phys. J. Spec. Top. 227, 1295 (2018).
  * Price and Cooper [2012] H. M. Price and N. R. Cooper, Phys. Rev. A 85, 033620 (2012).

## Supplemental Material

###  .1 Steady-state solution

In this section we briefly detail the solution for the non-equilibrium steady
state using the quantum Langevin formalism. We consider a general quadratic
Hamiltonian of the form $\hat{H}_{\rm
tot}=\hat{H}_{S}+\hat{H}_{SB}+\hat{H}_{B}$, with

$\displaystyle\hat{H}_{S}$
$\displaystyle=\sum_{j,k}H_{jk}\hat{a}^{\dagger}_{j}\hat{a}_{k},$ (S1)
$\displaystyle\hat{H}_{SB}$ $\displaystyle=\sum_{j\in\partial
S}\sum_{q}g_{qj}\left(\hat{a}^{\dagger}_{j}\hat{b}_{qj}+\hat{b}^{\dagger}_{qj}\hat{a}_{j}\right),$
(S2) $\displaystyle\hat{H}_{B}$ $\displaystyle=\sum_{j\in\partial
S}\sum_{q}\Omega_{qj}\hat{b}^{\dagger}_{qj}\hat{b}_{qj}.$ (S3)

In Eq. (S1), $\hat{a}^{\dagger}_{j}$ creates a boson or fermion on site $j$ of
the system and $H_{jk}$ are the elements of a hermitian single-particle
Hamiltonian matrix, $\mathbf{H}$. Similarly, the canonical operators
$\hat{b}^{\dagger}_{qj}$ create a particle in mode $q$ of the reservoir
connected to site $j$, with corresponding frequency $\Omega_{qj}$ and coupling
$g_{qj}$ to a site $j$ lying on the boundary of the system, $\partial S$.

Following the standard procedure, we formally solve the equations of motion
for $\hat{b}_{qj}$ in the Heisenberg picture and substitute the result into
the Heisenberg equation for $\hat{a}_{j}$, obtaining

$\mathrm{i}\partial_{t}\hat{a}_{j}(t)=\sum_{k}H_{jk}\hat{a}_{k}(t)+\int_{t_{0}}^{\infty}\mathrm{d}t^{\prime}\,\chi_{j}(t-t^{\prime})\hat{a}_{j}(t^{\prime})+\hat{\xi}_{j}(t).$
(S4)

Above, we have introduced the retarded memory kernel $\chi_{j}(t)$ and the
noise operator $\hat{\xi}_{j}(t)$, which for $j\in\partial S$ are given by

$\displaystyle\hat{\xi}_{j}(t)$
$\displaystyle=\sum_{q}g_{qj}\mathrm{e}^{-\mathrm{i}\Omega_{qj}(t-t_{0})}\hat{b}_{qj}(t_{0}),$
(S5) $\displaystyle\chi_{j}(t)$
$\displaystyle=-\mathrm{i}\theta(t)\sum_{q}g^{2}_{qj}\mathrm{e}^{-\mathrm{i}\Omega_{qj}t},$
(S6)

with $\theta(t)$ the unit step function, while for $j\notin\partial S$ we
define $\hat{\xi}_{j}(t)=0=\chi_{j}(t)$. Note that the memory kernel is
nothing but the retarded Green function of the noise operator, since
$\theta(t-t^{\prime})\langle[\hat{\xi}_{j}(t),\hat{\xi}^{\dagger}_{k}(t^{\prime})]_{\pm}\rangle=\mathrm{i}\delta_{jk}\chi_{j}(t-t^{\prime})$,
where the minus (plus) sign indicates the (anti-)commutator for bosons
(fermions).

Since we are interested in the steady state, we take the limit
$t_{0}\to-\infty$ in Eq. (S4) and solve in the Fourier domain. The solution
can be expressed in the compact form

$\tilde{\mathbf{a}}(\omega)=\mathbf{G}(\omega)\cdot\tilde{\bm{\xi}}(\omega),$
(S7)

by defining vectors
$\tilde{\mathbf{a}}=(\tilde{a}_{1},\tilde{a}_{2},\ldots)^{T}$ and
$\tilde{\bm{\xi}}(\omega)=(\tilde{\xi}_{1},\tilde{\xi}_{2},\ldots)^{T}$ of
Fourier-transformed operators and introducing the (matrix-valued) retarded
Green function

$\mathbf{G}(\omega)=\left[\omega\mathbf{1}-\mathbf{H}-\mathbf{\Sigma}(\omega)\right]^{-1}.$
(S8)

Here, $\mathbf{1}$ is the identity matrix and the retarded self-energy matrix
$\mathbf{\Sigma}(\omega)$ is diagonal with elements
$\Sigma_{jk}(\omega)=\delta_{jk}\tilde{\chi}_{j}(\omega)$ given by the Fourier
transform of the memory kernel. The latter decomposes as

$\tilde{\chi}_{j}(\omega)=\fint\frac{\mathrm{d}\omega^{\prime}}{2\pi}\frac{\gamma_{j}(\omega^{\prime})}{\omega-\omega^{\prime}}-\mathrm{i}\frac{\gamma_{j}(\omega)}{2},$
(S9)

where $\fint$ denotes a principal-value integral and we introduced the
spectral density

$\gamma_{j}(\omega)=2\pi\sum_{q}g_{qj}^{2}\delta(\omega-\Omega_{qj}).$ (S10)

The non-equilibrium steady state of the system depends on the statistical
properties of the noise operator via Eq. (S7). Assuming that the initial state
of the reservoirs (at time $t=t_{0}$) is thermal and uncorrelated, we find the
noise spectrum

$\langle\tilde{\xi}^{\dagger}_{j}(\omega)\tilde{\xi}_{k}(\omega^{\prime})\rangle=2\pi\delta(\omega-\omega^{\prime})\delta_{jk}\gamma_{j}(\omega)f_{j}(\omega),$
(S11)

where $f_{j}(\omega)=(\mathrm{e}^{\beta_{j}(\omega-\mu_{j})}\pm 1)^{-1}$ is
the reservoir distribution function, characterized by an inverse temperature
$\beta_{j}$ and chemical potential $\mu_{j}$, and expectation values are taken
with respect to the state at time $t=t_{0}$. It is now straightforward to
write down the two-point correlation function of the system in its steady
state, $C_{jk}(t)=\langle\hat{a}^{\dagger}_{k}(t)\hat{a}_{j}(0)\rangle$, as

$\mathbf{C}(t)=\int\frac{\mathrm{d}\omega}{2\pi}\mathrm{e}^{\mathrm{i}\omega
t}\mathbf{G}\cdot(\mathbf{\Gamma}\circ\mathbf{F})\cdot\mathbf{G}^{\dagger}.$
(S12)

Here we have suppressed frequency arguments and defined diagonal matrices
$\mathbf{\Gamma}=\mathrm{i}(\mathbf{\Sigma}-\mathbf{\Sigma}^{\dagger})$ and
$F_{jk}=\delta_{jk}f_{j}$, while $\circ$ denotes the elementwise (Hadamard)
product. The steady-state correlation matrix $\mathbf{C}=\mathbf{C}(0)$ is
found simply by setting $t=0$ in Eq. (S12).

The setup considered in the main text features hot baths coupled to sites on
the boundary at $x=1$ and cold baths coupled to the other boundary at
$x=L_{X}$. Moreover, the wide-band limit approximation is made, which amounts
to replacing the spectral densities $\gamma_{j}(\omega)$ by frequency-
independent constants. This approximation is valid so long as
$\gamma_{j}(\omega)$ does not vary significantly over the frequency range
within which the integrand of Eq. (S12) differs appreciably from zero.
Assuming symmetric coupling, $\gamma_{j}=\gamma$, one then recovers Eq. (4).
In this case, $\mathbf{\Gamma}=\mathbf{\Gamma}_{h}+\mathbf{\Gamma}_{c}$, where
the self-energy $\mathbf{\Gamma}_{h}$ ($\mathbf{\Gamma}_{c}$) of the hot
(cold) bath is given by a diagonal matrix, with
$[\mathbf{\Gamma}_{h,c}]_{jk}=\gamma\delta_{jk}$ whenever $j$ corresponds to a
boundary site with $x=1$ ($x=L_{X}$) and $[\mathbf{\Gamma}_{h,c}]_{jk}=0$
otherwise. Likewise, we have $f_{j}=\bar{n}_{h}$ ($f_{j}=\bar{n}_{c}$) for $j$
corresponding to $x=1$ ($x=L_{X}$).

###  .2 Definition of currents

The current densities are defined via the continuity equation for the particle
density, $\hat{n}_{x,y}={\bf\hat{a}}^{\dagger}_{x,y}\cdot{\bf\hat{a}}_{x,y}$,
as

$\displaystyle\mathrm{i}[\hat{H}_{X},\hat{n}_{x,y}]$
$\displaystyle=\hat{J}^{X}_{x-1,y}-\hat{J}^{X}_{x,y}$
$\displaystyle\mathrm{i}[\hat{H}_{Y},\hat{n}_{x,y}]$
$\displaystyle=\hat{J}^{Y}_{x,y-1}-\hat{J}^{Y}_{x,y}.$ (S13)

Here $\hat{J}^{X}_{x,y}$ denotes the current density operator describing
particle flow from site $(x,y)$ to its nearest neighbour in the $x$ direction,
i.e. $(x,y)\to(x+1,y)$. Similarly, $\hat{J}^{Y}_{x,y}$ is the corresponding
current in the $y$ direction, $(x,y)\to(x,y+1)$. Using the canonical
(anti-)commutation relations, it is straightforward to derive

$\displaystyle\hat{J}^{X}_{x,y}$
$\displaystyle=\frac{t_{X}}{2\mathrm{i}}{\bf\hat{a}}^{\dagger}_{x+1,y}\cdot(\sigma_{z}+\mathrm{i}\sigma_{y})\cdot{\bf\hat{a}}_{x,y}+{\rm
h.c.},$ (S14) $\displaystyle\hat{J}^{Y}_{x,y}$
$\displaystyle=\frac{t_{Y}}{2\mathrm{i}}{\bf\hat{a}}^{\dagger}_{x,y+1}\cdot(\sigma_{z}+\mathrm{i}\sigma_{x})\cdot{\bf\hat{a}}_{x,y}+{\rm
h.c.}$ (S15)

The expectation values of these one-body observables can then be read off
directly from the elements of the steady-state correlation matrix.

The total particle current flowing, say, into the cold bath is represented by
the observable $\hat{J}_{\rm tot}=\mathrm{i}[\hat{H}_{\rm tot},\hat{N}_{c}]$,
where $\hat{N}_{c}=\sum_{j\in c}\sum_{q}\hat{b}_{qj}^{\dagger}\hat{b}_{qj}$ is
the particle-number operator of the cold bath. Following a similar calculation
to Sec. .1, we find the standard Landauer formula for the mean steady-state
current [54]

$J_{\rm
tot}=\int\frac{\mathrm{d}\omega}{2\pi}\,\mathcal{T}(\omega)\left[\bar{n}_{h}(\omega)-\bar{n}_{c}(\omega)\right],$
(S16)

where
$\mathcal{T}(\omega)=\mathrm{Tr}\left[\mathbf{\Gamma}_{c}\cdot\mathbf{G}\cdot\mathbf{\Gamma}_{h}\cdot\mathbf{G}^{\dagger}\right]$
is the transmission function. Particle conservation implies that
$\sum_{y}J^{X}_{x,y}=J_{\rm tot}$ in the NESS.

###  .3 Particle-hole symmetry

In the fermionic setting, the system obeys a generalized particle-hole (GPH)
symmetry [59] when $\mu=\omega_{0}$. It is simplest to first analyse this
situation assuming that $\omega_{0}=0$. In this case, the fermionic
Hamiltonian [Eqs. (1)–(3)] is invariant under the combined flavor swap and
particle hole transformation,
$\hat{\Upsilon}\mathbf{\hat{a}}_{x,y}^{\dagger}\hat{\Upsilon}^{\dagger}=(\sigma_{x}\mathbf{\hat{a}}_{x,y})^{\rm
T}$ or, explicitly,
$\hat{\Upsilon}(\hat{a}^{\dagger}_{x,y,\uparrow},\hat{a}^{\dagger}_{x,y,\downarrow})\hat{\Upsilon}^{\dagger}=(\hat{a}_{x,y,\downarrow},\hat{a}_{x,y,\uparrow})$.
Meanwhile, the current operators defined in Eqs. (S14) and (S15) are odd under
this transformation, i.e.
$\hat{\Upsilon}\hat{J}^{X,Y}_{x,y}\hat{\Upsilon}^{\dagger}=-\hat{J}^{X,Y}_{x,y}$.
It follows that any state $\hat{\rho}$ that is GPH-invariant, satisfying
$[\hat{\rho},\hat{\Upsilon}]=0$, has vanishing currents. Indeed, schematically
we have
$\braket{\hat{J}}=\mathrm{Tr}[\hat{J}\hat{\rho}]=-\mathrm{Tr}[\hat{\Upsilon}\hat{J}\hat{\Upsilon}^{\dagger}\hat{\rho}]=-\braket{\hat{J}}$,
where $\hat{J}$ stands for any current operator that is GPH-odd.

In general, the correlation matrix can be written in block form as

$\mathbf{C}=\left(\begin{matrix}\mathbf{C}_{\uparrow\uparrow}&\mathbf{C}_{\uparrow\downarrow}\\\
\mathbf{C}_{\downarrow\uparrow}&\mathbf{C}_{\downarrow\downarrow}\end{matrix}\right),$
(S17)

where $\mathbf{C}_{\uparrow\uparrow}$ is the correlation matrix for the
$\uparrow$ flavor,
$\mathbf{C}_{\uparrow\downarrow}=\mathbf{C}_{\downarrow\uparrow}^{\dagger}$
describes coherences between the different flavors, etc. The condition for
GPH-invariance, $[\hat{\rho},\hat{\Upsilon}]=0$, implies that the fermionic
correlation matrix satisfies

$\mathbf{C}=\mathbf{1}-\mathbf{S}\cdot\mathbf{C}^{*}\cdot\mathbf{S},$ (S18)

where $\mathbf{S}$ is a matrix that swaps the flavors, and $\mathbf{C}^{*}$ is
the complex conjugate of $\mathbf{C}$. With the basis order implicit in Eq.
(S17), we have

$\mathbf{S}=\left(\begin{matrix}\mathbf{0}&\mathbf{1}\\\
\mathbf{1}&\mathbf{0}\end{matrix}\right),$ (S19)

where $\mathbf{0}$ is the zero matrix. The GPH symmetry of the many-body
Hamiltonian translates to the condition

$\mathbf{H}=-\mathbf{S}\cdot\mathbf{H}^{*}\cdot\mathbf{S}.$ (S20)

In turn, this implies the following constraint on the retarded Green function:

$\mathbf{G}(-\omega)=-\mathbf{S}\cdot\left[\mathbf{G}(\omega)\right]^{*}\cdot\mathbf{S}.$
(S21)

Note that here we assume the self-energy is also GPH-symmetric, i.e.
$\mathbf{\Sigma}(\omega)=-\mathbf{S}\cdot[\mathbf{\Sigma}(-\omega)]^{*}\cdot\mathbf{S}$.
In particular, this applies when the system-bath coupling is identical for
both flavor states, so that $[\mathbf{\Sigma},\mathbf{S}]=0$, and when the
spectral densities are symmetric, $\gamma_{j}(\omega)=\gamma_{j}(-\omega)$.

At the GPH-symmetric point, $\mu=0$, the distribution functions satisfy
$\bar{n}_{h,c}(\omega)=1-\bar{n}_{h,c}(-\omega)$. Thus, reversing the dummy
integration variable as $\omega\to-\omega$ in Eq. (4), we obtain

$\displaystyle\mathbf{C}$
$\displaystyle=\int\frac{\mathrm{d}\omega}{2\pi}\,\mathbf{G}(\omega)\cdot\Big{[}\mathbf{\Gamma}_{h}\bar{n}_{h}(\omega)+\mathbf{\Gamma}_{c}\bar{n}_{c}(\omega)\Big{]}\cdot\mathbf{G}^{\dagger}(\omega)$
$\displaystyle=\int\frac{\mathrm{d}\omega}{2\pi}\,\mathbf{G}(\omega)\cdot\Big{[}\mathbf{\Gamma}-\mathbf{\Gamma}_{h}\bar{n}_{h}(-\omega)-\mathbf{\Gamma}_{c}\bar{n}_{c}(-\omega)\Big{]}\cdot\mathbf{G}^{\dagger}(\omega)$
$\displaystyle=\mathbf{1}-\int\frac{\mathrm{d}\omega}{2\pi}\,\mathbf{G}(-\omega)\cdot\Big{[}\mathbf{\Gamma}_{h}\bar{n}_{h}(\omega)+\mathbf{\Gamma}_{c}\bar{n}_{c}(\omega)\Big{]}\cdot\mathbf{G}^{\dagger}(-\omega)$
$\displaystyle=\mathbf{1}-\mathbf{S}\cdot\mathbf{C}^{*}\cdot\mathbf{S}.$ (S22)

Above, we have made use of Eq. (S21) as well as the identity
$\mathbf{G}\cdot\mathbf{\Gamma}\cdot\mathbf{G}^{\dagger}=\mathrm{i}\left(\textbf{G}-\mathbf{G}^{\dagger}\right)$
and the property

$\int\frac{\mathrm{d}\omega}{2\pi}\mathrm{i}\left(\textbf{G}-\mathbf{G}^{\dagger}\right)=\mathbf{1},$
(S23)

which follows from the (anti-)commutation relations [9]. Since the NESS is the
unique Gaussian state with correlation matrix given by Eq. (.3), it follows
that the state is GPH-symmetric when $\mu=\omega_{0}=0$. This implies that all
currents vanish. Consistent with this observation, GPH symmetry implies that
the transmission function is even, $\mathcal{T}(\omega)=\mathcal{T}(-\omega)$,
and therefore the integral in Eq. (S16) is odd and equates to zero.

The same GPH symmetries and associated vanishing currents are obtained when
$\mu=\omega_{0}\neq 0$. This can be shown simply by shifting the frequency
variable $\omega\to\omega+\mu$ in Eq. (.3). The final result then follows from
the GPH symmetry of the shifted Hamiltonian
$\mathbf{H}^{\prime}=\mathbf{H}-\omega_{0}\mathbf{1}$. We have numerically
confirmed the GPH symmetry of the NESS in the wide-band limit when
$\mu=\omega_{0}$, both in weak- and strong-coupling regimes.

###  .4 Weak-coupling limit

Following Ref. [53], we show that in the weak-coupling limit the NESS reduces
to an energy-diagonal density matrix that coincides with the solution of the
secular-Born-Markov master equation [7, 10, 19]. The system Hamiltonian can be
written in its spectral representation

$\mathbf{H}=\mathbf{U}\cdot\mathbf{\Omega}\cdot\mathbf{U}^{\dagger},$ (S24)

where $\mathbf{\Omega}={\rm diag}\\{\omega_{\alpha}\\}$ and $\mathbf{U}$
comprise the eigenvalues and eigenvectors of $\mathbf{H}$. The many-body
Hamiltonian in this basis reads as

$\hat{H}=\sum_{\alpha}\omega_{\alpha}\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha},\qquad\hat{c}_{\alpha}=\sum_{j}U^{*}_{j\alpha}\hat{a}_{j}.$
(S25)

Now consider the retarded Green function (S8) in the wide-band limit,
$\mathbf{\Sigma}=-\mathrm{i}\mathbf{\Gamma}/2\propto\gamma$. As $\gamma\to 0$,
$\mathbf{G}(\omega)$ becomes approximately diagonal in the energy
representation. Following identical arguments to Ref. [53], we obtain the
steady-state correlation matrix in the energy eigenbasis as

$\braket{\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha^{\prime}}}=\delta_{\alpha\alpha^{\prime}}\frac{s_{\alpha}\bar{n}_{h}(\omega_{\alpha})+r_{\alpha}\bar{n}_{c}(\omega_{\alpha})}{s_{\alpha}+r_{\alpha}},$
(S26)

where $s_{\alpha}$ and $r_{\alpha}$ describe the dimensionless coupling
strength of eigenmode $\alpha$ to the hot and cold baths, respectively:

$\displaystyle
s_{\alpha}=\gamma^{-1}\sum_{j,k}U^{*}_{j\alpha}[\Gamma_{h}]_{jk}U_{k\alpha}=\sum_{y=1}^{L_{Y}}\sum_{\sigma=\uparrow,\downarrow}|U_{(1,y,\sigma),\alpha}|^{2},$
(S27) $\displaystyle
r_{\alpha}=\gamma^{-1}\sum_{j,k}U^{*}_{j\alpha}[\Gamma_{c}]_{jk}U_{k\alpha}=\sum_{y=1}^{L_{Y}}\sum_{\sigma=\uparrow,\downarrow}|U_{(L_{X},y,\sigma),\alpha}|^{2}.$
(S28)

We now consider the Lindblad master equation obtained within the Born-Markov
and secular approximations. Following the standard derivation [7, 10] yields
the master equation

$\dot{\hat{\rho}}=-\mathrm{i}[\hat{H},\hat{\rho}]+\mathcal{L}_{h}\hat{\rho}+\mathcal{L}_{c}\hat{\rho},$
(S29)

with dissipators for the hot and cold baths given respectively by

$\displaystyle\mathcal{L}_{h}=\gamma\sum_{\alpha}s_{\alpha}\left(\bar{n}_{h}(\omega_{\alpha})\mathcal{D}[\hat{c}^{\dagger}_{\alpha}]+[1\mp\bar{n}_{h}(\omega_{\alpha})]\mathcal{D}[\hat{c}_{\alpha}]\right),$
(S30)
$\displaystyle\mathcal{L}_{c}=\gamma\sum_{\alpha}r_{\alpha}\left(\bar{n}_{c}(\omega_{\alpha})\mathcal{D}[\hat{c}^{\dagger}_{\alpha}]+[1\mp\bar{n}_{c}(\omega_{\alpha})]\mathcal{D}[\hat{c}_{\alpha}]\right),$
(S31)

where
$\mathcal{D}[\hat{L}]\bullet=\hat{L}\bullet\hat{L}^{\dagger}-\tfrac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\bullet\\}$
is a Lindblad superoperator and the plus (minus) sign is for bosons
(fermions). It is easy to check that the steady-state solution of Eq. (S29),
$\dot{\hat{\rho}}=0$, is a Gaussian state with correlation matrix given
exactly by Eq. (S26).

As is well known, the Lindblad equation (S29) can be written as

$\dot{\hat{\rho}}=-\mathrm{i}\hat{H}_{\rm
eff}\hat{\rho}+\mathrm{i}\hat{\rho}\hat{H}_{\rm
eff}^{\dagger}+\mathcal{J}_{h}\hat{\rho}+\mathcal{J}_{c}\hat{\rho},$ (S32)

where the non-Hermitian Hamiltonian is

$\displaystyle\hat{H}_{\rm eff}$
$\displaystyle=\hat{H}-\frac{\mathrm{i}}{2}\sum_{\alpha}\gamma
s_{\alpha}\left[\bar{n}_{h}(\omega_{\alpha})+(1\mp
2\bar{n}_{h}(\omega_{\alpha}))\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha}\right]$
$\displaystyle\qquad-\frac{\mathrm{i}}{2}\sum_{\alpha}\gamma
r_{\alpha}\left[\bar{n}_{c}(\omega_{\alpha})+(1\mp
2\bar{n}_{c}(\omega_{\alpha}))\hat{c}^{\dagger}_{\alpha}\hat{c}_{\alpha}\right],$
(S33)

and the superoperators $\mathcal{J}_{c,h}$ represent “quantum jumps” due to
exchange of particles with the baths [60],

$\displaystyle\mathcal{J}_{h}\hat{\rho}$ $\displaystyle=\sum_{\alpha}\gamma
s_{\alpha}\left[\bar{n}_{h}(\omega_{\alpha})\hat{c}_{\alpha}^{\dagger}\hat{\rho}\hat{c}_{\alpha}+\left(1+\bar{n}_{h}(\omega_{\alpha})\right)\hat{c}_{\alpha}\hat{\rho}\hat{c}_{\alpha}^{\dagger}\right],$
(S34) $\displaystyle\mathcal{J}_{c}\hat{\rho}$
$\displaystyle=\sum_{\alpha}\gamma
r_{\alpha}\left[\bar{n}_{c}(\omega_{\alpha})\hat{c}_{\alpha}^{\dagger}\hat{\rho}\hat{c}_{\alpha}+\left(1+\bar{n}_{c}(\omega_{\alpha})\right)\hat{c}_{\alpha}\hat{\rho}\hat{c}_{\alpha}^{\dagger}\right].$
(S35)

By including both the non-Hermitian Hamiltonian and the jump terms, we ensure
a completely positive and physically consistent description of the full many-
body state of the non-equilibrium lattice system, valid at arbitrary times.
This is in contrast to non-Hermitian models that neglect the jumps (valid for
rare post-selected trajectories of continuously measured systems) or the
effective non-Hermitian Hamiltonians appearing in single-particle Green
functions, the topological properties of which have recently been of intense
interest; e.g., see Ref. [61] and references therein.

Figure S1: Current patterns for the same asymmetric impurity distribution
shown in Fig. 4(c) but with different parameters. (a,b) Fermionic currents in
the topologically nontrivial regime, $m=t$, with on-site energy
$\omega_{0}=10t$, chemical potential $\mu=\omega-0.5t$, and (a) low
temperatures, $T_{h}=t$ and $T_{c}=0.01t$, or (b) higher temperatures,
$T_{h}=100t$ and $T_{c}=t$. (c) Same as (a) but showing bosons (red) and
fermions (blue) in a topologically trivial regime, $m=3t$.

###  .5 Weak-coupling distribution functions

Figure S2: Distribution functions in the NESS for bosons (red) and fermions
(blue). (a) Symmetric impurity configuration shown in Fig. 4(b), with
$T_{h}=t$, $T_{c}=0.01t$, $\omega_{0}=10t$ and $\mu=\omega_{0}-0.1t$. (b)
Asymmetric impurity configuration shown in Fig. 4(c), but with
$\mu=\omega_{0}-0.5t$ as in Fig. S1. Circles show the data for $T_{h}=t$ and
$T_{c}=0.01t$ [Fig. S1(a)], while light blue squares show the result for
fermions at high temperature, $T_{h}=100t$ and $T_{c}=t$ [Fig. S1(b)]. The
values of $n(\omega_{\alpha})$ for bosons are multiplied by $10^{3}$ to be
visible on the same axes.

We now investigate in more detail the form of the weak-coupling distribution
function,
$n(\omega_{\alpha})=\braket{\hat{c}_{\alpha}^{\dagger}\hat{c}_{\alpha}}$,
defined by Eq. (S26). When the reservoirs are at the same temperature,
$\bar{n}_{h}=\bar{n}_{c}=\bar{n}$, we obtain simply
$n(\omega_{\alpha})=\bar{n}(\omega_{\alpha})$. This is a smooth, monotonically
decreasing function of frequency that is independent of the mode-specific
couplings $s_{\alpha}$ and $r_{\alpha}$. As a consequence, the equilibrium
currents—as seen, for example, in Fig. 2(c)—are unaffected by the addition of
defects in any configuration [19].

Out of equilibrium, an analogous result is obtained in the clean QWZ model
described by Eqs. (1)–(3) due to the symmetries
$\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$ and $\hat{\Pi}\hat{R}_{\pi}$. In
particular, the symmetry $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$ implies that
$|U_{(1,y,\sigma),\alpha}|^{2}=|U_{(L_{X},y,\sigma),\alpha}|^{2}$ and
therefore $s_{\alpha}=r_{\alpha}$ identically. Similarly, the symmetry under
$\hat{\Pi}\hat{R}_{\pi}$ means that
$|U_{(1,y,\sigma),\alpha}|^{2}=|U_{(L_{X},L_{Y}+1-y,\sigma),\alpha}|^{2}$, and
again $s_{\alpha}=r_{\alpha}$ upon summing over $y$ in Eqs. (S27) and (S28).
As a consequence, the mode occupations (S26) become independent of
$s_{\alpha}$ and $r_{\alpha}$, being given simply by the average of the
reservoir distribution functions:

$n(\omega_{\alpha})=\frac{1}{2}\left[\bar{n}_{c}(\omega_{\alpha})+\bar{n}_{h}(\omega_{\alpha})\right].$
(S36)

This is again a function of the frequency $\omega_{\alpha}$ only, and
therefore leads to stable current patterns as in the equilibrium case [19].
Any distribution of defects that respects one of the above discrete symmetries
will also yield $s_{\alpha}=r_{\alpha}$, leading to similar behaviour.

Note that the relevant symmetries are determined both by the Hamiltonian and
by the configuration of the baths driving the system out of equilibrium. For
example, the transformation $\hat{\Theta}\hat{\Sigma}_{x}$ is a symmetry of
the QWZ Hamiltonian, where $\hat{\Sigma}_{x}=\hat{R}_{\pi}\hat{\Sigma}_{y}$ is
a reflection about the $x$ axis. However, because this reflection does not
interchange the hot and cold baths it is irrelevant for the NESS. Therefore,
defect distributions that are only symmetric under $\hat{\Sigma}_{x}$ destroy
the bosonic edge currents whenever $T_{c}\neq T_{h}$.

Examples of the nonequilibrium distribution functions obtained for a symmetric
defect configuration are shown in Fig. S2(a), corresponding to the system
depicted in Fig. 4(b). As expected, the distributions are well behaved
functions of frequency. Conversely, Fig. S2(b) shows that when the defect
configuration does not respect either of the discrete symmetries, as in Figs.
4(c) or S1, the mode occupations are erratically varying functions of
frequency. This disrupts the erasure effect and leads to the destruction of
the boundary current patterns, in general.

Fig. S2 also shows that the fermionic distributions exhibit a sudden change at
the chemical potential whenever $T_{c}$ is sufficiently small. This sharp
feature explains the absence of counter-currents around the impurities in
Figs. 4(b,c). In this case, the chemical potential is just below the particle-
hole symmetric point $\mu=\omega_{0}$. This leads to a single “hole” that
occupies an edge state, visible in Fig. S2 as the isolated point just below
$\omega_{\alpha}=\omega_{0}$. The edge current pattern observed in Fig. 4(b)
is primarily due to this hole state which, occupying an edge mode, is robust
against any bulk perturbation. Conversely, when the chemical potential is
shifted down to $\mu=\omega_{0}-0.5t$, several holes are created in other edge
states [Fig. S2(b), blue circles], including edge states localized around the
impurities. Consequently, we see the appearance of counter-currents shielding
the impurities in Fig. S1(a).

At higher temperature, the sharp feature in $n(\omega_{\alpha})$ disappears
[Fig. S2(b), light blue squares]. The current distribution is then determined
mainly by bulk states and thus depends intimately on the symmetries of the
impurity configuration. If the configuration is asymmetric, as in Fig. S1(b),
the boundary currents are strongly disrupted.

In the topologically trivial regime, there are no edge states and both bosons
and fermions behave very similarly with respect to the introduction of
defects. In particular, asymmetric impurity distributions tend to inhibit the
formation of clean boundary current patterns. An example is shown in Fig.
S1(c).

Figure S3: The emergence of edge currents in the weak coupling limit can be
understood by a semiclassical argument. (a) The factor
$n[\omega_{-}(\mathbf{k})]-n[\omega_{+}(\mathbf{k})]$ in Eq. (S48) for the
bosonic case with $\omega_{0}=10t$, $T_{h}=t$, $T_{c}=0.01t$ and $m=2.1t$. (b)
The Berry curvature ${\mathcal{F}\\!\\!}_{-}(\mathbf{k})$, which has a
balanced positive and negative contributions over B.Z. (the zero value plane
is depicted in light blue color) as for $m>2t$ the model is in a trivial phase
and ${\mathcal{F}\\!\\!}_{-}(\mathbf{k})$ integrates out to zero. (c) The
product
$\\{n[\omega_{-}(\mathbf{k})]-n[\omega_{+}(\mathbf{k})]\\}{\mathcal{F}\\!\\!}_{-}(\mathbf{k})$,
which breaks such a balance, leads to a nonvanishing integral in Eq. (S48).
(d) The semiclassical current pattern generated according (S48) for the
potential $V(x,y)=(x/8)^{10}+(y/8)^{10}$.

###  .6 Semiclassical picture

The behaviour in the weak coupling can be understood by using a semiclassical
picture [57, 62]. Under periodic boundary conditions, the QWZ Hamiltonian can
be written in the form

$\hat{H}=\sum_{\mathbf{k}\in\mathrm{B.Z.}}\mathbf{a}_{\mathbf{k}}^{\dagger}\cdot\mathbf{H}(\mathbf{k})\cdot\mathbf{a}_{\mathbf{k}},$
(S37)

with $\mathbf{k}=(k_{X},k_{Y})$ the quasimomentum and

$\mathbf{H}(\mathbf{k})=\omega_{0}\mathbf{1}+\mathbf{h}(\mathbf{k})\cdot\bm{\sigma},$
(S38)

where

$\mathbf{h}(\mathbf{k})=[t_{Y}\sin(k_{Y}),t_{X}\sin(k_{X}),m+t_{X}\cos(k_{X})+t_{Y}\cos(k_{Y})].$
(S39)

Note that we have taken a slightly different Hamiltonian with respect to the
original reference [12]. The two energy bands of the model are then given by

$\omega_{\pm}(\mathbf{k})=\omega_{0}\pm|\mathbf{h}(\mathbf{k})|,$ (S40)

and for the isotropic case considered throughout this work
($t_{X}=t_{Y}=t>0$),

$|\mathbf{h}(\mathbf{k})|=t\sqrt{\sin^{2}(k_{X})+\sin^{2}(k_{Y})+[m/t+\cos(k_{X})+\cos(k_{Y})]^{2}}.$
(S41)

The action of a force $\mathbf{F}$ modifies the semiclassical equation for the
velocity $\mathbf{v}_{\alpha}$ of a Bloch wave-packet in the band
$\omega_{\alpha}(\mathbf{k})$ (assuming the single-particle Hamiltonian (S38)
is gapped, i.e. $m\neq 0,\pm 2t$, and $m$ is small enough to justify the
semiclassical approximation) as

$\mathbf{v}_{\alpha}=\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}}-\mathbf{F}\times\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k}).$
(S42)

Here,
$\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k})=[0,0,{\mathcal{F}\\!\\!}_{\alpha}(\mathbf{k})]$,
with ${\mathcal{F}\\!\\!}_{\alpha}(\mathbf{k})$ the Berry curvature of the
band, which, for a two-band model, can be written in the form

$\displaystyle{\mathcal{F}\\!\\!}_{\pm}(\mathbf{k})$
$\displaystyle=\mp\frac{1}{2}\mathbf{n}(\mathbf{k})\cdot\left[\frac{\partial\mathbf{n}(\mathbf{k})}{\partial
k_{X}}\times\frac{\partial\mathbf{n}(\mathbf{k})}{\partial k_{Y}}\right],$
(S43)

with
$\mathbf{n}(\mathbf{k}):=\mathbf{h}(\mathbf{k})/|\mathbf{h}(\mathbf{k})|$. For
the QWZ Hamiltonian we find

${\mathcal{F}\\!\\!}_{\pm}(\mathbf{k})=\frac{\pm\big{[}\cos(k_{X})+\cos(k_{Y})+\tfrac{m}{t}\cos(k_{X})\cos(k_{Y})\big{]}}{2\big{\\{}\\!\sin^{2}(k_{X})+\sin^{2}(k_{Y})+\big{[}\tfrac{m}{t}+\cos(k_{X})+\cos(k_{Y})\big{]}^{2}\big{\\}}^{\frac{3}{2}}}.$
(S44)

We can now argue for the presence of an edge circulation under open boundary
conditions by introducing a confining potential $V(x,y)$ in the semiclassical
equation for $\mathbf{v}_{\pm}$ (S42). Assuming that the potential varies
slowly on the scale of the lattice, the energy of the wavepacket in the semi-
classical approximation is simply $\omega_{\alpha}(\mathbf{k})+V(x,y)$ [57].
Moreoever, the potential induces a confining force
$\mathbf{F}=-\bm{\nabla}V(x,y)$, which approximately vanishes in the bulk and
becomes very large in the boundary of the confined region,

$\displaystyle\mathbf{v}_{\alpha}^{\mathrm{(bulk)}}\simeq\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}},\quad\mathbf{v}_{\alpha}^{\mathrm{(edge)}}\simeq\bm{\nabla}V(x,y)\times\bm{\mathcal{F}}\\!\\!_{\alpha}(\mathbf{k}).$
(S45)

Since $\mathbf{F}$ is normal and points inward from the confined boundary, a
circulating current is induced along this boundary with direction given the
sign of the Berry curvature, which depends on the band $\alpha$, the sign of
$m$, and the specific value of $\mathbf{k}$: clockwise for
${\mathcal{F}\\!\\!}_{\alpha}(\mathbf{k})>0$ and counterclockwise for
${\mathcal{F}\\!\\!}_{\alpha}(\mathbf{k})<0$. Thus, for the NESS at the weak
coupling regime, we may estimate the current density (current per unit width)
by [57]

$\displaystyle\mathbf{I}(x,y)=\sum_{\alpha=\pm}\int_{\mathrm{B.Z.}}d^{2}\mathbf{k}\,n[\omega_{\alpha}(\mathbf{k})+V(x,y)]\mathbf{v}_{\alpha}(x,y,\mathbf{k}),$
(S46)

where the integral extends over the first Brillouin zone. Here we have assumed
that the NESS occupation number $n$ only depends on energy. As discussed in
Sec. .5, this is the case under a $\hat{\Pi}\hat{\Theta}\hat{\Sigma}_{y}$ or a
$\hat{\Pi}\hat{R}_{\pi}$ symmetric configuration.

Now, if we consider regions far from the boundary, so that
$V(x,y)\ll\omega_{\alpha}(\mathbf{k})$, we may approximate Eq. (S46) by

$\displaystyle\mathbf{I}(x,y)$
$\displaystyle\simeq\int_{\mathrm{B.Z.}}d^{2}\mathbf{k}\,\\{n[\omega_{-}(\mathbf{k})]-n[\omega_{+}(\mathbf{k})]\\}\mathbf{v}_{-},$
(S47)

using the fact that $\mathbf{v}_{+}=-\mathbf{v}_{-}$. Since there is no
privileged direction of $\mathbf{v}_{-}^{\mathrm{(bulk)}}$ on each energy
shell: $\omega_{\alpha}(\mathbf{k})=\omega_{\alpha}(-\mathbf{k})$ and hence
$\frac{\partial\omega_{\alpha}(\mathbf{k})}{\partial\mathbf{k}}=-\frac{\partial\omega_{\alpha}(-\mathbf{k})}{\partial\mathbf{k}}$.
Therefore, the contribution of $\mathbf{v}_{\alpha}^{(\rm bulk)}$ to the
integrand is odd and integrates to zero. We conclude that

$\displaystyle\mathbf{I}(x,y)$
$\displaystyle\simeq\int_{\mathrm{B.Z.}}d^{2}\mathbf{k}\,\left\\{n[\omega_{-}(\mathbf{k})]-n[\omega_{+}(\mathbf{k})]\right\\}\mathbf{v}_{-}^{(\rm
edge)},$ (S48)

which clearly vanishes in the bulk of the system where $V(x,y)=0$ (see Eq.
(S45)). Nevertheless, since $n(\omega)$ is a monotonically decreasing function
of energy [Eq. (S36)], there is an imbalance of the contributions to
$\mathbf{I}(x,y)$ within the Brillouin zone in regions where
$\mathbf{v}_{\alpha}^{\rm edge}\propto\nabla V(x,y)\neq 0$. Thus, Eq. (S48)
predicts a net chiral edge current, see Fig. S3. This erasure effect sustains
the appearance of the chiral current in the symmetric situation, but it is
spoilt in the nonsymmetric case, as the occupation number is no longer only a
function of the energy alone.