Topology by Dissipation: Transport properties

TL;DR

This paper investigates how engineered dissipation can create topological 2D Chern insulator states, analyzing their transport properties and revealing differences from equilibrium systems, including conditions for quantized Hall conductance.

Contribution

It demonstrates the emergence of topological phases via dissipation without Hamiltonian tuning and explores their transport properties using open quantum systems tools.

Findings

Dissipative engineering can produce topological states efficiently.

The relation between Chern number and Hall conductance can be broken in non-equilibrium.

Conditions for quantized Hall conductance are identified.

Abstract

Topological phases of matter are the center of much current interest, with promising potential applications in, e.g., topologically-protected transport and quantum computing. Traditionally such states are prepared by tuning the system Hamiltonian while coupling it to a generic bath at very low temperatures; This approach is often ineffective, especially in cold-atom systems. It was recently shown that topological phases can emerge much more efficiently even in the absence of a Hamiltonian, by properly engineering the interaction of the system with its environment, to directly drive the system into the desired state. Here we concentrate on dissipatively-induced 2D Chern insulator (lattice quantum Hall) states. We employ open quantum systems tools to explore their transport properties, such as persistent currents and the conductance in the steady state, in the presence of various Hamiltonians. We find that, in contrast with equilibrium systems, the usual relation between the Chern topological number and the Hall conductance is broken. We explore the intriguing edge behaviors and elucidate under which conditions the Hall conductance is quantized.

Equations

$$\frac{d}{dt}{\rho }_D=0$$

$$\frac{d}{dt}\rho =-i\left[H,\rho \right]+\mathcal{D}\rho ,$$

$$H^{\mathrm{ref}}$$

$$H^{\mathrm{ref}}=\sum _{\mathbf{k},\lambda }{ϵ}_{\mathbf{k},\lambda }{a}_{\mathbf{k},\lambda }^{†}{a}_{\mathbf{k},\lambda },$$

$$H^{\mathrm{diss}}$$

$$-{\mu }^{\ast }\sum _{l_x,l_y}{b}_{l_x,l_y}^{†}{a}_{l_x,l_y}+h.c.,$$

$$H^{\mathrm{diss}}=\sum _{\mathbf{k},\lambda }\left({ϵ}_{\mathbf{k},\lambda }-{\mu }^{\ast }\right)b_{\mathbf{k},\lambda }^{†}{a}_{\mathbf{k},\lambda }+h.c.,$$

$${\mathcal{D}}^{\mathrm{out}}\rho$$

$${\gamma }_{\mathbf{k},\lambda }^{\mathrm{out}}=\frac{2\pi }{\hslash }{\nu }_0{\left|{ϵ}_{\mathbf{k},\lambda }-{\mu }^{\ast }\right|}^2,$$

$${\gamma }^0\equiv \frac{2\pi }{\hslash }{\nu }_0{\left|t_{\mathrm{diss}}\right|}^2.$$

$$\mathcal{D}\rho$$

$$+{\gamma }^{\mathrm{in}}\sum _{\mathbf{k},\lambda }\left(a_{\mathbf{k},\lambda }^{†}\rho a_{\mathbf{k},\lambda }-\frac{1}{2}\left\{a_{\mathbf{k},\lambda }{a}_{\mathbf{k},\lambda }^{†},\rho \right\}\right).$$

$$n_{\mathbf{k},\lambda }=\frac{{\gamma }^{\mathrm{in}}}{{\gamma }^{\mathrm{in}}+{\gamma }_{\mathbf{k},\lambda }^{\mathrm{out}}}.$$

$${\gamma }_{\mathrm{opt}}^{\mathrm{in}}=\sqrt{\max \left\{{\gamma }_{\mathbf{k},1}^{\mathrm{out}}\right\}\cdot \min \left\{{\gamma }_{\mathbf{k},2}^{\mathrm{out}}\right\}}.$$

$$H_{\mathrm{comp}}^{\mathrm{sys}}$$

$$H_{\mathrm{inc}}^{\mathrm{sys}}$$

$$H_S^{\mathrm{diss}}=U^{†}{H}_L^{\mathrm{diss}}U=-t_{\mathrm{diss}}\sum _{l_x,l_y,s=\pm 1}{b}_{l_x,l_y}^{†}\left(e^{is{l}_x\pi \alpha }{a}_{l_x,l_y+s}+e^{-is{l}_x\pi \alpha }{a}_{l_x+s,l_y}\right)-{\mu }^{\ast }\sum _{l_x,l_y}{b}_{l_x,l_y}^{†}{a}_{l_x,l_y}+h.c.,$$

$$H_{\mathrm{comp},S}^{\mathrm{sys}}=U^{†}{H}_{\mathrm{comp},L}^{\mathrm{diss}}U=-t_{\mathrm{comp}}\sum _{l_x,l_y}{a}_{l_x,l_y}^{†}\left(e^{i{l}_x\pi \alpha }{a}_{l_x,l_y+1}+e^{-i{l}_y\pi \alpha }{a}_{l_x+1,l_y}\right)+h.c.,$$

$${\tilde{H}}_{\mathrm{inc}}^{\mathrm{sys}}=U^{†}{H}_{\mathrm{inc}}^{\mathrm{diss}}U=-t_{\mathrm{inc}}\sum _{l_x,l_y}{a}_{l_x,l_y}^{†}\left(e^{-i{l}_x\pi \alpha }{a}_{l_x,l_y+1}+e^{-i{l}_y\pi \alpha }{a}_{l_x+1,l_y}\right)+h.c.$$

$$P_{l_x,l_y;l_x^{\prime },l_y^{\prime }}=\mathrm{Tr}\left\{\rho \left(t\to \infty \right)a_{l_x,l_y}^{†}{a}_{l_x^{\prime },l_y^{\prime }}\right\}.$$

$$\left[\frac{1}{2}\left({\Gamma }^{\mathrm{in}}+{\Gamma }^{\mathrm{out}}\right)+ih\right]P+P\left(\frac{1}{2}\left({\Gamma }^{\mathrm{in}}+{\Gamma }^{\mathrm{out}}\right)-ih\right)={\Gamma }^{\mathrm{in}},$$

$$\frac{d}{dt}{n}_{l_x,l_y}\left(t\right)$$

$${\dot{n}}_{l_x,l_y}^H\left(t\right)=-\left(j_{l_x+1,l_y}^x-j_{l_x,l_y}^x\right)-\left(j_{l_x,l_y+1}^y-j_{l_x,l_y}^y\right).$$

$${\dot{n}}_{l_x,l_y}^{\mathcal{D}}\equiv J_{l_x,l_y}^{\mathcal{D},\mathrm{in}}-J_{l_x,l_y}^{\mathcal{D},\mathrm{out}},$$

$$j_{l_x,l_y}^x=i{t}_{\mathrm{comp}}\left(P_{l_x+1,l_y;l_x,l_y}-P_{l_x,l_y;l_x+1,l_y}\right),$$

$$j_{l_x,l_y}^y=i{t}_{\mathrm{comp}}\left(e^{i{l}_{x}2\pi \alpha }{P}_{l_x,l_y+1;l_x,l_y}-e^{-i{l}_{x}2\pi \alpha }{P}_{l_x,l_y;l_x,l_y+1}\right).$$

$$j_{l_x,l_y}^x=i{t}_{\mathrm{inc}}\left(P_{l_x+1,l_y;l_x,l_y}-P_{l_x,l_y;l_x+1,l_y}\right),$$

$$j_{l_x,l_y}^y=i{t}_{\mathrm{inc}}\left(P_{l_x,l_y+1;l_x,l_y}-P_{l_x,l_y;l_x,l_y+1}\right).$$

$$\delta H_{\mathbf{E}}=-\sum _{l_x,l_y}\left(E_x{l}_x+E_y{l}_y\right)a_{l_x,l_y}^{†}{a}_{l_x,l_y},$$

$$\left[\frac{{\Gamma }^{\mathrm{in}}+{\Gamma }^{\mathrm{out}}}{2}+ih\right]\delta P+\delta P\left[\frac{{\Gamma }^{\mathrm{in}}+{\Gamma }^{\mathrm{out}}}{2}-ih\right]=i\left[P_0,\delta h\right],$$

$$G_{xx}=\underset{E_x\to 0}{\lim }\frac{1}{E_x}\sum _{l_y}\frac{1}{L_y}\delta j_{\frac{1}{2}{L}_x,l_y}^x,$$

$$G_{yx}=\underset{E_x\to 0}{\lim }\frac{1}{E_x}\sum _{l_x,l_y}\frac{1}{L_y{L}_x}\delta j_{l_x,l_y}^y,$$

$$G_{yy}=\underset{E_y\to 0}{\lim }\frac{1}{E_y}\sum _{l_x}\frac{1}{L_x}\delta j_{l_x,\frac{1}{2}{L}_y}^y,$$

Full text

I Introduction

The adverse affects of a dissipative environment or bath on a quantum system coupled to it are well-known: Processes of decay and decoherence lead to a degradation of the desired quantum state and shortening of its coherence time (Gardiner and Haken, 1991; Hirota et al., 2012). For these reasons, a prime goal of quantum device engineering is to mitigate those effects, mainly by trying to reduce the system-bath interaction strength. However, it has been realized that an engineered dissipative interaction to an external bath can possibly be exploited, in order to reliably prepare interesting non-trivial quantum states in open systems; This approach has recently been brought to many-body systems (Poyatos et al., 1996; Diehl et al., 2008; Kraus et al., 2008; Verstraete et al., 2009; Krauter et al., 2011; Otterbach and Lemeshko, 2014; Kapit et al., 2014a). This is typically done by tailoring the bath interaction such that the required quantum state would emerge as a unique steady state solution to the open quantum (Lindblad (Lindblad, 1976)) master equation, independent of the initial conditions. If the sought after state is represented by the density matrix $\rho_{D}$, one would aim to manipulate the system-bath coupling such that $\rho_{D}$ is a “dark state” of the dynamical evolution of the open system, i.e.,

[TABLE]

This idea of dissipative preparation contrasts the conventional Hamiltonian approach, that heavily relies on reaching sufficiently low temperatures to attain the interesting properties of the desired quantum ground state. This approach might prove ineffective in some cases. Of particular interest is the case of quantum simulators implemented using ultra-cold atomic gases (Georgescu et al., 2014; Bloch et al., 2008), where equilibrium at low temperatures compared to the trapping potential energy scale is experimentally challenging to achieve. A dissipative preparation scheme, albeit with its own challenges and complexities, may circumvent such issues, as the accuracy of the final state will be determined solely by the degree of control one has on the details of the engineered bath coupling.

Topologically non-trivial phases of matter (Wen, 1995) have been at the forefront of condensed matter physics during the last several decades, since the discovery of the integer and fractional quantum Hall effects (Klitzing et al., 1980; Tsui et al., 1982). The more recently discovered topological insulators (Kane and Mele, 2005; Bernevig et al., 2006; König et al., 2007), superconductors (Fu and Kane, 2008; Lutchyn et al., 2018), and semimetals (Armitage et al., 2018) have opened the door to novel exciting possibilities, e.g., topologically protected quantum computing (Nayak et al., 2008). The idea of employing the engineered bath interactions scheme to stabilize a topologically ordered ground state has already received some attention in recent years (Diehl et al., 2011; Bardyn et al., 2012, 2013; König and Pastawski, 2014; Kapit et al., 2014b, a; Budich et al., 2015; Iemini et al., 2016; Gong et al., 2017; Goldstein, 2019).

We focus on the 2D Chern insulator phase discussed in Ref. (Goldstein, 2019), stabilized by employing purely dissipative dynamics of Lindblad type, i.e., the master equation describing the evolution of the system density matrix $\rho$,

[TABLE]

has the Hamiltonian $H=0$, and where ${\cal D}$ is some general engineered Lindbladian dissipator super-operator acting on $\rho$. ${\cal D}$ is of the form $\mathcal{D}\rho=\sum_{j}\gamma_{j}\left(L_{j}\rho L_{j}^{\dagger}-\frac{1}{2}\left\{\rho,L_{j}^{\dagger}L_{j}\right\}\right)$, with $L_{j}$’s being a set of “quantum jump” operators, and $\gamma_{j}$’s the rates governing the dissipative dynamics (Gardiner and Haken, 1991). The jump operators and rates which we use in this work will be specified later on. This work addresses the effects of including $H\neq 0$ in the Lindblad dynamics (2). This step is crucial, as it allows one to define current operators in the system, such that transport properties of the dissipative state may be explored. These properties, e.g, the Hall conductance, are the hallmark of the well-known equilibrium counterpart of this phase, and as such are important in characterizing the engineered state. The interplay of Hamiltonian and Lindbladian dynamics in topological systems has also been the focus of some recent works (Dangel et al., 2018; van Caspel et al., 2019). We will show that although the dissipatively engineered states could be arbitrarily close to equilibrium quantum Hall states, they do not always present the same persistent current and transport features. These depend on the details of the coherent dynamics, the relative amplitude of those dynamics compared to the dissipation energy scale, and on the implementation of artificial fields within the dissipative scheme.

The rest of the paper is organized as follows. In Sec. II we briefly introduce the dissipative scheme developed in (Goldstein, 2019). We introduce the Hamiltonian dynamics, and the tools required to analyze the steady state in its presence in Sec. III. This step enables us to discuss the persistent steady state currents that develop for different classes of Hamiltonians in Sec. IV. Then, the electric conductance response is calculated in Sec. V for the different cases, and compared with the known results for the equilibrium scenario. Finally, we summarize our findings and conclusions in Sec. VI. Some technical details are relegated to the Appendices.

II Dissipatively induced topological state

In this section, we recapitulate the dissipative recipe presented in Ref. (Goldstein, 2019) to realize a dissipative lattice integer quantum Hall state. The main components in the scheme are illustrated in Fig. 1. Consider the Harper-Hofstadter model (Harper, 1955; Hofstadter, 1976), describing nearest-neighbor hopping on a two-dimensional square lattice pierced by a magnetic flux, which we will use as a reference in the construction below,

[TABLE]

with $a_{l_{x},l_{y}}$ a fermionic annihilation operator on the site $\left(l_{x},l_{y}\right)$, $t_{{\rm ref}}$ the hopping strength, and $\alpha$ is the magnetic flux per plaquette in units of flux quantum, introduced using the Peierls substitution. $\alpha$ is assumed to be a rational fraction $\alpha=\frac{p}{q}$ (with $p,q,$ integers with no common factor). The Hamiltonian (3) is diagonalized by moving to two-dimensional momentum space,

[TABLE]

with ${\bf k}$ the two-dimensional momentum, and $\lambda=1,2,...,q$ the band index. The spectrum obtained for $\alpha=\frac{1}{7}$ with periodic boundary conditions along one direction is shown in Fig. 2(a), where the distinct bands are apparent, as well as the edge state spectrum.

The goal of the construction of (Goldstein, 2019) is to achieve a state as close as possible to the ground state of the reference Hamiltonian with chemical potential in the first gap (i.e., the lowest band completely filled and all the other bands completely empty) by purely dissipative dynamics. For that we do not implement the Hamiltonian (3)–(4), but rather introduce a different Hamiltonian, built using the matrix elements of the Hofstadter reference Hamiltonian (3),

[TABLE]

describing a coupling of the lattice fermions $a$ to an additional fermionic species $b$ (e.g., a different hyperfine state). A thorough discussion regarding the cold-atom implementation for this Hamiltonian is presented in (Goldstein, 2019), and is based on previous equilibrium constructions for the realization of artificial gauge fields (Jaksch and Zoller, 2003; Gerbier and Dalibard, 2010). In the eigenbasis of the reference Hamiltonian we obtain

[TABLE]

with $\epsilon_{\mathbf{k},\lambda}$ now proportional to $t_{{\rm diss}}$ instead of $t_{{\rm ref}}$. Note that now they do not represent energies in the $\mathbf{k},\lambda$ single particle basis (since we do not implement reference Hamiltonian (3)–(4), though we still use its eigenbasis to define the state). Rather, they represent the amplitudes of the coupling of the $a$ and $b$ particle in the $\mathbf{k},\lambda$ state.

We assume that the $a$-particles are optically trapped in the 2D system plane, whereas the $b$-particles are free to escape along the perpendicular $z$ direction. This dynamics of the $b$-particles is captured by supplementing $H^{{\rm diss}}$ by an Hamiltonian for the bath particles, $H^{{\rm bath}}=\sum_{l_{x},l_{y},q_{z}}\epsilon_{z}\left(q_{z}\right)b_{l_{x},l_{y},q_{z}}^{\dagger}b_{l_{x},l_{y},q_{z}}$, where $q_{z}$ is the momentum of the $b$-particles along the $z$ direction, and $\epsilon_{z}\left(q_{z}\right)$ the corresponding dispersion relation. The $a$-particles are coupled via Eq. (5) to the $b$ operators in the system plane $z=0$, i.e., to $b_{l_{x},l_{y}}=\sum_{q_{z}}\frac{1}{\sqrt{N_{z}}}b_{l_{x},l_{y},q_{z}}$, where $N_{z}$ is the system size in the perpendicular direction. We assume that, in the cold-atom implementation, $b$-particles which leave the 2D lattice escape to infinity, and as a consequence, $\left\langle b_{l_{x},l_{y},q_{z}}^{\dagger}b_{n_{x},n_{y},q_{z}^{\prime}}\right\rangle={\rm Tr}\left\{\rho_{B}b_{l_{x},l_{y},q_{z}}^{\dagger}b_{n_{x},n_{y},q_{z}^{\prime}}\right\}=0$, for all $q_{z}$, $q_{z}^{\prime}$, with $\rho_{B}$ the density matrix of the bath particles. Hence, the $b$-particles constitute a bath with infinitely negative chemical potential for the $a$-particles, and we can integrate them out.

Here the central point of this construction becomes clear: The Hamiltonian of the $a$ particles themselves is zero; they only experience a coupling to the $b$-particle bath. Hence, the dynamics is not set by the energetics of the $a$ system, but rather via controlling the matrix elements of the $a$-$b$ coupling by Eqs. (5)–(6). As detailed in Ref. (Goldstein, 2019), upon integrating out the $b$ bath, one arrives at a contribution to the Lindblad master equation for the density matrix of the $a$ particles,

[TABLE]

Calculation of the rates $\gamma_{\mathbf{k},\lambda}^{\mathrm{out}}$ can be done using Fermi’s golden rule (Dirac, 1927), leading to

[TABLE]

with $\nu_{0}$, the density of available $b$-states, taken as constant (Goldstein, 2019). We thus define a typical dissipative scale for the system,

[TABLE]

The implication of (8) is that given a flat band, i.e., $\epsilon_{\mathbf{k},\lambda}=\epsilon_{0}$ for one value of $\lambda$, it is possible to fine-tune $\mu^{*}$ such that the depletion rate goes to zero for that particular band alone. The system will exponentially decay into a state where the finely-tuned flat band is the only occupied band, with its occupancy depending on the initial state of the lattice filling. This is analogous to a low temperature equilibrium scenario where the chemical potential lies in a gap between the bands. Again, this is achieved by tuning the system-bath coupling matrix elements, not their energetics.

Since a topologically non-trivial (non-zero Chern number) exactly flat-band in a finite-range hopping Hamiltonian, separated by a finite gap from the other bands, is not possible (Chen et al., 2014; Read, 2017), we consider a band which is nearly flat, i.e., with a width much smaller compared to its minimal distance to the other bands. Then, a very small value of $\gamma_{\mathbf{k},\lambda}^{\mathrm{out}}$ for that particular band is attainable. We now introduce another ingredient to the dissipative scheme, a global filling reservoir, replenishing all bands at a rate $\gamma^{\mathrm{in}}$. This can be achieved by coupling the system to another reservoir of $c$-particles with infinitely positive chemical potential. For simplicity the corresponding reference Hamiltonian will be taken as proportional to the unit matrix in either real or $\mathrm{\boldsymbol{k}},\lambda$ space, leading to $\mathrm{\boldsymbol{k}},\lambda$-independent refilling rate (Goldstein, 2019). The dissipative master equation now reads

[TABLE]

The steady-state band occupation numbers, $n_{{\bf k},\lambda}={\rm Tr}\left\{\rho a_{\mathbf{k},\lambda}^{\dagger}a_{\mathbf{k},\lambda}\right\}$ can now be calculated,

[TABLE]

Given the values of the maximal rate for the bottom band, $\mathrm{max}\left\{\gamma_{\mathbf{k},1}^{\mathrm{out}}\right\}$, and the minimal one for the next (closest) band $\mathrm{min}\left\{\gamma_{\mathbf{k},2}^{\mathrm{out}}\right\}$, an optimal choice of the tunable refilling rate, such that $n_{{\bf k},1}\sim 1$ and $n_{{\bf k},\lambda>1}\sim 0$ is thus

[TABLE]

Fig. 2(b) shows an example of results for this scheme, where we have properly tuned $\mu^{*}$ and $\gamma^{\mathrm{in}}$, and achieved almost completely full or empty occupation of the lower band or upper bands, respectively. This mixed state is not only very close to the pure ground state of the reference Hamiltonian, but also shares with it the value -1 of the topological index, the Chern number (Goldstein, 2019).

III Introducing Hamiltonian dynamics

We now consider adding some coherent dynamics to the $a$ particles which live on the lattice. This will modify the master equation (2), since if $H\neq 0$, one must include the commutator term between the Hamiltonian and the density matrix. One may naturally combine the dissipative artificial gauge field engineering approach of Ref. (Goldstein, 2019) with the non-dissipative artificial gauge field engineering approach of Refs. (Jaksch and Zoller, 2003; Gerbier and Dalibard, 2010), to independently set the artificial gauge fields in both parts of Eq. (2). Of the many possibilities, in this work, we explore the consequences of two kinds of hopping Hamiltonians,

[TABLE]

[TABLE]

where $H_{{\rm comp}}^{{\rm sys}}$ is compatible with the dissipative interaction Hamiltonian, i.e., it has the dynamics of particles hopping on a square lattice with the same flux $\alpha$ as before, and $H_{{\rm inc}}$ is a simple nearest-neighbor hopping Hamiltonian which is incompatible. Understanding the different consequences for particles with a Hamiltonian, which is either compatible or incompatible with engineered dissipation, will enable us to determine how system observables change and the role of topology in the presence of both dissipative and coherent (Hamiltonian) dynamics.

The total Hamiltonian, Eq. (5) combined with either Eq. (13) or (14) [and hence the resulting Lindblad dynamics (2)] are invariant under a general gauge transformation $U=\exp\left[i\sum_{l_{x},l_{y}}f_{l_{x},l_{y}}(a_{l_{x},l_{y}}^{\dagger}a_{l_{x},l_{y}}+b_{l_{x},l_{y}}^{\dagger}b_{l_{x},l_{y}})\right]$. For example, choosing $f_{l_{x}l_{y}}=-\pi\alpha l_{x}l_{y}$ to go from the Landau to the symmetric gauge in Eq. (5) leads to

[TABLE]

The compatible system Hamiltonian then transforms to the symmetric form as well:

[TABLE]

while the incompatible Hamiltonian becomes

[TABLE]

Of course, gauge invariance dictates that all physical observables calculated using Eq. (5) and (13) are the same as those calculated with (15) and (16) in the compatible case, and similarly for (5) and (14) vs. (15) and (17) in the incompatible case. However, another possibility in the latter incompatible case is to work with Eq. (15) together with the unmodified incompatible Hamiltonian (14). This latter possibility (to be referred to as "incompatible-S” in the following) is not related by a gauge transformation, and hence physically different from working with (5) and (14), or equivalently (15) and (17) (to be referred to as “incompatible-L” in the following). In the following we will mostly work with the incompatible-L scheme, but also show the difference with the incompatible-S scheme when relevant. We would like to stress that both these schemes are equally natural from the point of view of cold atom implementation.

First, we want to understand how the presence of either (13) or (14) affects the steady state we arrived at using the engineered dissipative scheme. We notice right away that the steady state solution for $\rho$ with $H=0$ commutes with $H_{{\rm comp}}^{{\rm sys}}$, since it is diagonal in the same ${\bf k},\lambda$ basis as $\gamma^{{\rm out}}$. As a result, Eq. (11) holds exactly even for a non-zero compatible Hamiltonian $H=H_{{\rm comp}}^{{\rm sys}}$.

Conversely, if $H=H_{{\rm inc}}^{{\rm sys}}$ the steady state can be much different than (11), depending on the ratio $\frac{t_{{\rm inc}}}{\gamma^{0}}$. We use the full master equation to calculate the steady state single particle density matrix, whose eigenvalues are the occupation numbers $n_{\mathrm{{\bf k}},\lambda}$,

[TABLE]

This requires us to solve a Sylvester-type matrix equation (see Appendix A for details)

[TABLE]

with $h$, $\Gamma^{\mathrm{in}}$, and $\Gamma^{\mathrm{out}}$ the single-particle matrices corresponding to $H$, $\gamma^{\mathrm{in}}$, and $\gamma^{\mathrm{out}}$, respectively. Eq. (19) is similar in essence to the Lyapunov equations derived for the covariance matrix in Refs. (Kosov et al., 2011; Koga and Yamamoto, 2012; Nicacio et al., 2016; Schwarz et al., 2016). We find that the engineered nearly-pure steady state deteriorates with increasing $t_{{\rm inc}}$. As shown in Fig. 3, the incompatibility of the Hamiltonian requires one to use a faster refilling rate $\gamma^{{\rm in}}$ to maintain the high occupation of the bottom band. This in turn leads to non-negligible occupation of the higher energy levels, deviating from the desired state. This phenomenon becomes significant at about $t_{{\rm inc}}\sim 0.1\gamma^{0}$, i.e., when the Hamiltonian is no longer negligible compared to the dissipative energy scale. We note that the gap in the spectrum of $P$ (purity gap) is still finite, and the associated Chern number retains its value of -1 (Goldstein, 2019).

IV Persistent currents

Including an Hamiltonian for the $a$ species finally allows us to define a sensible current operator in the system, by looking at the time evolution of the local particle density expectation value $n_{l_{x},l_{y}}\left(t\right)\equiv{\rm Tr}\left\{\rho\left(t\right)a_{l_{x},l_{y}}^{\dagger}a_{l_{x},l_{y}}\right\}$ (see Appendix B for the full derivation). Using the master equation (2) for the evolution of the density matrix, one can separate the local change in particle number into a coherent contribution $\dot{n}^{H}$ and a dissipative part $\dot{n}^{{\cal D}}$,

[TABLE]

The coherent part, controlled by the $\left[H,\rho\right]$ term in the master equation, must obey a continuity equation involving the particle current, $\dot{n}^{H}+\nabla\cdot{\bf j}=0$, where $\left(\nabla\cdot\right)$ is a lattice version of the divergence operator, or more explicitly,

[TABLE]

By examining the expression for $\dot{n}^{H}$, which depends on the details of the Hamiltonian, one can extract and define a proper current operator. The current operator which we investigate indeed describes the local $a$-particle flow rate in the system, and is thus a physical measurable quantity. As shown in Appendix B, the steady state expectation values of the current are fully given in terms of elements of $P$, so no further complicated calculations are required in order to obtain them. Since in the steady state $\frac{d}{dt}n=0$, the dissipative steady state contribution can be calculated from the divergence of the current, $\dot{n}^{{\cal D}}=-\dot{n}^{H}=\nabla\cdot{\bf j}$, and can be decomposed into incoming and outgoing terms,

[TABLE]

which are proportional to the single particle matrix elements of $\Gamma^{{\rm in}}$ and $\Gamma^{{\rm out}}$, respectively (see Appendix B). An alternative derivation for these current operators, which relies on coupling the Hamiltonian to a probe gauge field, is given at the end of Appendix B. We present our results below for the different classes of Hamiltonian we have considered.

IV.1 Compatible hopping

The expectation value of the electric current operator in the different directions is given in terms of elements of the $P$ matrix by (Appendix B)

[TABLE]

Since $P$ itself is not affected by the presence of the Hamiltonian, the normalized current $\frac{j_{l_{x},l_{y}}^{\mu}}{t_{{\rm comp}}}$ is completely independent of $t_{{\rm comp}}$. As shown in Fig. 4, we find chiral currents strongly localized near the edges, much like in the ground state of the equilibrium quantum Hall effect. We note that the current in this compatible case is divergence-free, as expected for a chiral edge mode. As a consequence, the dissipative current $\dot{n}^{{\cal D}}\propto\nabla\cdot{\bf j}$ is zero, signaling a zero local net flux of particles into or out of the baths.

IV.2 Incompatible hopping

In the case of the hopping Hamiltonian (14), the current is similarly given by

[TABLE]

We find that the steady state currents which develop in the system in such a scenario are completely different than those of the compatible Hamiltonian, as we show in Fig. 5. The current is not localized near the edges, but rather flows along the $y$ direction with a structure periodic in $l_{x}$, determined by the magnetic unit cell whose size is $q$ sites in this particular dissipative configuration. The current trajectories terminate at the system edges where they reverse their flow direction. This reversal appears to be mediated by a local interchanging flux of particles from or into the bath, as indicated by the non-vanishing dissipative current $\dot{n}_{\mathcal{D}}$ near the edges.

This result seems to be somewhat peculiar at first, since the current in the system appears to prefer one direction ($y$) over the other ($x$). This is a result of using, in the incompatible L-scheme, the dissipative Hamiltonian (5), with its anisotropic phase factors, in combination with the isotropic incompatible system Hamiltonian (14). This is not the case in the gauge-inequivalent incompatible S-scheme, Eq. (15) and (14). And indeed, in that case, as evident in Fig. 6, the current pattern which emerges has circulating currents around plaquettes with $q\times q$ sites (recall that $\alpha=\frac{p}{q}$, so for, e.g., $\alpha=\frac{1}{7}$ one has $q=7$).

Qualitatively, the current patterns do not change appreciably as $t_{{\rm inc}}$ increases beyond the $t_{{\rm inc}}\ll\gamma^{0}$ limit. What we observe is a steep decline in the amplitude of the normalized current $\frac{j}{t_{{\rm inc}}}$ as $t_{{\rm inc}}$ becomes comparable to the dissipative scale $\gamma^{0}$. This is understood by the depletion of current carriers in the system as the steady state occupation deteriorates, and is compensated by sufficiently increasing the value of the refilling rate $\gamma^{{\rm in}}$.

V conductance

One of the most remarkable properties of the quantum Hall ground state is its exactly quantized transverse electrical conductance. In order to reveal the response and transport properties of the dissipatively prepared state, we introduce a small perturbation to the Hamiltonian, playing the role of an electric field along the $x$ direction. Using the relation between the electric field ${\bf E}$ and the electrostatic potential $\phi$, $\mathbf{E}=-\nabla\phi$, we incorporate a perturbative term synonymous with $\int d\mathbf{x}\rho\left(x\right)\phi\left(x\right)$,

[TABLE]

with the electric field ${\bf E}=(E_{x},E_{y},0)$. As we are only interested in the linear response regime, we calculate the first-order change in the occupation numbers matrix $\delta P$ due to a finite value of $\mathbf{E}$, which amounts to solving another matrix Sylvester equation,

[TABLE]

with $P_{0}$ the solution of Eq. (19) without an applied field, and $\delta h$ is the single particle matrix corresponding to $\delta H$. To the best of our knowledge, a Sylvester/Lyapunov equation of the kind of Eq. (26) has not been considered before as a method for calculating linear response properties. We calculate the modification of the values of the current due to this perturbation, $\delta j_{l_{x},l_{y}}^{\mu}=j_{l_{x},l_{y}}^{\mu}\left(P\rightarrow P+\delta P\right)-j_{l_{x},l_{y}}^{\mu}\left(P\right)$, and find the conductance tensor, which we define by

[TABLE]

with $E_{x}=0$. Notice that in all terms of $G$ we average over the dimension perpendicular to the electric field applied. In all calculations presented below we always have periodic boundary conditions in that dimension. For more details and the justification of Eq. (27), see Appendix (C). Similarly to the previous section, we discuss our findings separately for the two possible Hamiltonian classes.

V.1 Compatible hopping

Unlike our discussion regarding the persistent currents, it is clear from Eq. (26) that even in the compatible case, the relative amplitude of the hopping Hamiltonian, as compared to dissipation rate, may play an important role. We indeed find two regimes for the conductance in our system, corresponding to $\frac{t_{{\rm comp}}}{\gamma^{0}}\gg 1$ and $\frac{t_{{\rm comp}}}{\gamma^{0}}\ll 1$ with a smooth crossover around $t_{{\rm comp}}\approx\gamma^{0}$, see Fig. 7a. Interestingly, in both regimes the transverse currents flows near system edges in a similar fashion (Fig. 7b,c).

In the large $\frac{t_{{\rm comp}}}{\gamma^{0}}$ case, which we refer to as the Hamiltonian regime, we observe quantum-Hall-like behavior, namely the quantization of the transverse response in accordance with the topological Chern number associated with the steady-state density matrix, which is equal to the one associated with the lowest band of $H^{\mathrm{diss}}$ (Goldstein, 2019). The value at the conductance plateau stabilizes roughly when $t_{{\rm comp}}$is larger than $10\gamma^{0}$. For the parameters used in Fig. 7, the value of the Hall conductance at the plateau is $G_{yx}=-0.984\frac{e^{2}}{h}$, and we have checked numerically that the deviation from ideal quantization vanishes at the thermodynamic limit and scales as $\propto\frac{1}{L_{x}}$. Alongside this approximate quantization, we observe a vanishing longitudinal response with increasing $t_{{\rm comp}}$, as expected from the limit $t_{{\rm comp}}\gg\gamma^{0}$, where an insulating behavior is predicted. The correction to this asymptotic limit seems to scale approximately with $\frac{\gamma^{0}}{t_{{\rm comp}}}$, indicating that the first order correction in the dissipation to Eq. (26) is dominant and does not cancel out.

In the opposite limit, the dissipative regime, both transverse and longitudinal responses are small, and increase as $t_{{\rm comp}}$ becomes larger. While $G_{xx}$ grows linearly with $\frac{t_{{\rm comp}}}{\gamma^{0}}$, the transverse conductance $G_{yx}$ is proportional to $\left(\frac{t_{{\rm comp}}}{\gamma^{0}}\right)^{2}$. To understand why the linear contribution to $G_{yx}$ vanishes, although the current itself is proportional to $t_{{\rm comp}}$, one must examine the $h=0$ limit of Eq. (26). Performing a Fourier transform of all the fermionic operators in the system with respect to the $y$ direction, $a_{l_{x},l_{y}}\sim\sum_{k_{y}}e^{ik_{y}l_{y}}a_{l_{x},k_{y}}$, one finds that in this basis the matrices $\left(\Gamma^{\mathrm{in}}+\Gamma^{\mathrm{out}}\right)$, $P_{0}$, and $\delta h$ are purely real, which would make $\delta P$ in this limit purely imaginary. Since $\delta P$ is hermitian, its diagonal elements $\left\langle a_{l_{x},k_{y}}^{\dagger}a_{l_{x},k_{y}}\right\rangle$ vanish. But $j^{y}$ is composed of exactly these vanishing averages. Hence $G_{yx}$ is zero to first order in $\frac{t_{{\rm comp}}}{\gamma^{0}}$.

We have also investigated how the quantized conductance (in the Hamiltonian regime) changes when the parameters controlling the dissipative scheme, i.e., $\mu^{*}$ and $\gamma^{{\rm in}}$, are modified, and do not assume their optimal values. We find that a precise tuning of $\mu^{*}$ is required for the nearly precise quantization (Fig. 8a), and that a small deviation from the optimal value deteriorates $G_{yx}$ significantly. This is because $\mu^{*}$ should be tuned to minimize $\gamma_{{\bf k},1}^{{\rm out}}$, such that the $\lambda=1$ band would be nearly filled in the steady state. Once $\mu^{*}$ is detuned away from the middle of the first band, this band will always be partially filled (unless we increase $\gamma^{{\rm in}}$ to the point that higher bands are no longer nearly empty). In contrast, the transverse response is somewhat less sensitive to changes in $\gamma^{{\rm in}}$, and we find that only changes of order of magnitude have an appreciable effect (Fig. 8b): A small change of the refilling rate does little to change the occupation of the different bands, leading to small deviations from the quantized value.

V.2 Incompatible hopping

Similarly to what we have seen for the steady state persistent currents, things change when one has an Hamiltonian incompatible with the magnetic flux in (5). Whereas in the dissipative regime with $t_{{\rm inc}}\ll\gamma^{0}$ the conductance features similar dependence on $\frac{t_{{\rm inc}}}{\gamma^{0}}$ as in the coherent case, at higher $\frac{t_{{\rm inc}}}{\gamma^{0}}$ the conductance begins to decline in amplitude, never reaching the topological quantized value for $G_{yx}$, see Fig. 9. Again, this should not be surprising, as we have already seen that a large $t_{{\rm inc}}$ negates our ability to manipulate the particles into the desired dissipatively engineered steady state, which possesses some QHE-like features. Moreover, increasing $t_{{\rm inc}}$ beyond the dissipative regime without altering $\gamma^{{\rm in}}$, we in fact get a much less populated steady state, which in turn can carry significantly less current. We note that in the dissipative regime, the current distribution is also indistinguishable from the behavior for the compatible Hamiltonian, e.g., Fig. 7b.

Lastly, we find that the conductance matrix is anisotropic, namely that $G_{xx}\neq G_{yy}$ and $G_{yx}\neq-G_{xy}$ in this incompatible dissipative regime, which is not the case for the compatible Hamiltonian. Once more, this is nothing but a consequence of the anisotropy of the dissipation Hamiltonian (5), which explicitly has a preferable axis, when combined with the incompatible system Hamiltonian (14). Using instead the incompatible-S scheme, Eqs. (15) and (14), removes this anisotropy, as one would expect. In Fig. 10 we show calculation of different elements in the conductance matrix for both cases. For the longitudinal conductance, i.e., $G_{xx}$ and $G_{yy}$, we find not only different values, but also different functional form, which depends on both the bath coupling convention and the direction. In the incompatible-L scheme, $G_{yy}$ appears to vanish (up to numerical noise). This is due to the oscillatory nature of the underlying steady state along the $\hat{x}$ direction, as implied by the persistent currents in Fig. 5. There is nearly perfect cancellation of counterflowing current within the magnetic unit cell, giving rise to a vanishing total response to the electric field perturbation. Similar oscillatory behavior apparent from Fig. 6 also explains the much weaker longitudinal response that scales as $\propto\left(\frac{t_{{\rm inc}}}{\gamma^{0}}\right)^{3}$in the more symmetric incompatible-S scheme. It appears that the finite response in this case stems from the violation of perfect periodicity of the the persistent current flow by the edges of the cylinder, see again Fig. 6. This violation does not exist in the asymmetric incompatible-L scheme since according to Fig. 5, the edges of the cylinder are parallel to the persistent current flow and do not affect it. The transverse response however, shows a more universal trend of $\propto\left(\frac{t_{{\rm inc}}}{\gamma^{0}}\right)^{2}$, albeit with different amplitudes for $G_{xy}$, $G_{yx}$ in the asymmetric case and for $G_{yx}=-G_{xy}$ in the symmetric one. Again, it appears that the reason for the L/S difference, also in the transverse conductance, is the susceptibility of the underlying steady state to current flow, which can be inferred from the persistent current patterns, Figs. 5 and 6.

VI Conclusions

We investigate the transport properties of the purely dissipative theoretical scheme, presented in Ref. (Goldstein, 2019), which reproduces a state synonymous with a very low temperature equilibrium quantum Hall state in a 2D lattice exploiting engineered dissipation. These properties can not be probed without introducing some coherent Hamiltonian dynamics for the lattice particles. We have demonstrated that a departure from the completely dissipative scheme, crucial for the transport study, as well as being more experimentally realistic, does not necessarily harm the engineered steady state, provided that the Hamiltonian is compatible with the evaporative part, or alternatively, that the incompatible part of that Hamiltonian is sufficiently small in magnitude compared to the dissipative energy scale.

Having introduced the Hamiltonian dynamics, we could explore the persistent currents that flow once the dissipative steady state is reached in the system. In the case where the lattice Hamiltonian has the same magnetic field present in the evaporative dynamics the main observation is the well-known chiral edge-modes, characteristic of equilibrium quantum Hall and Chern insulator systems, and the lack of any other currents in the system. On the other hand, a hopping Hamiltonian lacking the compatible magnetic field, induces bulk currents in a pattern and a direction determined by the chosen pattern of phases of the artificial gauge field (the gauge inequivalent L vs. S schemes). Near the edges the net flow to/from the reservoirs becomes finite locally (though it sums up to zero globally), facilitating current backscattering at the edges via the particle reservoir.

The electrical conductance of the dissipative steady state was also examined. Quantization of the transverse conductance, consistent with the theoretical Chern number, was shown to arise for the compatible Hamiltonian, provided it is sufficiently larger as compared to the dissipative rates. The quality of this quantization is affected by the accuracy of the assigned dissipative parameters, namely the refilling rate $\gamma^{\mathrm{in}}$ and $\mu^{*}$, which allows to tune the evaporation rate of the filled band. We find the quantization is much less sensitive to the former compared to the latter. A weaker compatible Hamiltonian featured mainly currents in the vicinity of the edges mediated by a dissipative current, together with a weak second order transverse response, also located near the edges of the system. The conductance matrix in the incompatible regime is slightly more complicated in structure, as the choice of system-bath coupling details (incompatible-L vs. incompatible-S schemes) may render it anisotropic: we observe longitudinal response in the $x$ direction identical to the one obtained with a weak compatible Hamiltonian, and different power law trends depending on the scheme. The transverse response seems to have a universal $\propto t_{{\rm inc}}^{2}$ behavior, but with anisotropic and scheme dependent absolute values.

Under the circumstances of a compatible hopping Hamiltonian, sufficiently larger in amplitude than the dissipative energy scale, the potential for obtaining some of the equilibrium quantum Hall properties, such as chiral edge-states and quantized conductance, may be put to the test by checking its robustness to some finite amount of disorder, both in the bath interaction and lattice Hamiltonians, which will be the subject of a future study. The way the transverse current is carried through the system in its dissipative steady state, under such perturbations, can be probed using the tools we have developed, and may shed further insights into the evaporative processes occurring in the system. This also opens up the possibility for engineering increasingly complex coupling to the bath, allowing exploration of more interesting possibilities and topological traits through the engineering of dissipative two-particle interactions. This in turn would allow the exploration of the dissipative analogues of exotic fractional states (Jain, 1989; Levin and Stern, 2009), and perhaps the equivalents of anyons (Schmidt et al., 2008; Bauer et al., 2014; Nayak et al., 2008).

VII Acknowledgments

We would like to thank B. A. Bernevig, B. Bradlyn, J. Budich, J. I. Cirac, N. Cooper, S. Diehl, M. Hafezi, Y. Gefen, A. Gorshkov, A. Kemenev, N. Schuch, H. H. Tu, and P. Zoller for useful discussions. This work has been supported by the Israel Science Foundation (Grant No. 227/15), the German Israeli Foundation (Grant No. I-1259-303.10), the US-Israel Binational Science Foundation (Grants No. 2014262 and 2016224), and the Israel Ministry of Science and Technology (Contract No. 3-12419).

Appendix A Steady state occupation numbers

In this appendix, we re-derive some of the results of Ref. (Schwarz et al., 2016), and bring the parts relevant to our discussion for the reader’s convenience. We begin by introducing a quadratic Hamiltonian for the lattice particles $a$ (which in this work will be either the compatible Hamiltonian (13), or the incompatible Hamiltonian (14)),

[TABLE]

with generalized indexes $A,B,$ which may each represent for example two spatial indexes (e.g., $A=l_{x},l_{y}$ and $B=l_{x}+1,l_{y}-1$).

The dissipator is also assumed to be quadratic, so that the Lindbladian master equation for the density matrix is

[TABLE]

with the $\Gamma$ matrices encapsulating the dissipative processes. In this work, we have the matrices

[TABLE]

[TABLE]

We define the following matrix of expectation values $P_{AB}\left(t\right)=\left\langle a_{A}^{\dagger}a_{B}\right\rangle_{t}$, with the shorthand notation $\left\langle M\right\rangle_{t}=\mathrm{Tr}\left\{\rho\left(t\right)M\right\}$. According to Eq. (29), the dynamics of this matrix is given by

[TABLE]

Using the definition for $P_{AB}\left(t\right)$, Wick’s theorem, and the fermionic anti-commutation relations, we find the matrix equation

[TABLE]

The steady state version of Eq. (33) can be manipulated into a Sylvester equation for the matrix $P$, Eq. (19) of the main text. If $\Gamma^{\left(1\right)}$ and $\Gamma^{\left(2\right)}$ can be simultaneously diagonalized, such as in our case, where $\Gamma^{\left(2\right)}$ is a constant times unity matrix, the steady solution for (19) in the purely dissipative regime ($h\rightarrow 0$) is given (in the basis where the $\Gamma$ matrices are diagonal) by

[TABLE]

This reduces to Eq. (11) in the main text. Notice however, that this result remains intact even in the presence of an Hamiltonian $h$ which is diagonal in the same basis as the $\Gamma$ matrices, due to the cancellation of the commutator term in (33). This is the reason for the distinction between the compatible and the incompatible Hamiltonians: In the former case (34) holds, but not in the latter.

Appendix B The current operators

Consider the expectation value of the particle density in the site $\left(l_{x},l_{y}\right)$ of the lattice, $n_{l_{x},l_{y}}\left(t\right)\equiv\left\langle\rho\left(t\right)a_{l_{x},l_{y}}^{\dagger}a_{l_{x},l_{y}}\right\rangle=P_{l_{x},l_{y};l_{x},l_{y}}\left(t\right)$. Its time-evolution is given by

[TABLE]

with

[TABLE]

[TABLE]

Plugging in $H_{\mathrm{comp}}$, we get

[TABLE]

The r.h.s. can be written as a discrete divergence,

[TABLE]

which allows us to define the currents as in Eq. (23). A transformation $t_{{\rm comp}}\rightarrow t_{{\rm inc}}$, $\alpha\rightarrow 0$, then gives the current in the presence of $H_{{\rm inc}}$, Eq. (24).

A closer look at $\dot{n}^{\mathcal{D}}$ reveals that it can be separated into two terms, corresponding, respectively, to flow of particles from the refilling reservoir into the system, and particles leaving the system into the evaporative reservoir,

[TABLE]

Also note that for a compatible Hamiltonian, plugging in the steady state $P$, given by Eq. (34), leads to the steady state $\dot{n}^{\mathcal{D}}$ vanishing exactly. More generally, in the steady state $\frac{d}{dt}n=0$, hence the “dissipative current” is immediately found (for any choice of Hamiltonian) to be

[TABLE]

with the r.h.s. given by Eq. (39).

For completeness, we also list here the expressions for the current operators for all choices of $H^{{\rm sys}}$. With the system Hamiltonian $H_{{\rm comp},S}^{{\rm sys}}$ given by Eq. (16), we find

[TABLE]

[TABLE]

For the choice of $\tilde{H}_{{\rm inc}}^{{\rm sys}}$ in Eq. (17), the current operators are given by

[TABLE]

[TABLE]

Of course, for choices related by a physical gauge transformation, although the expression for the current may look different, its expectation value is gauge independent. On the other hand, we note that whether we work with the gauge-inequivalent incompatible-L or incompatible-S choices, the expressions for the current operators, Eqs. (24a)–(24b) are unaltered. However, since $H^{{\rm diss}}$ differs, the steady-state $P$ matrix will vary, and thus the calculated currents.

An alternative way to derive the form of the different current operators on the lattice is outlined below. One begins with the total Hamiltonian, comprise of $a$-particles hopping terms, $a$-$b$ interactions, and the bath particles Hamiltonian. One introduces a probe gauge field ${\cal A}_{l_{x},l_{y};n_{x},n_{y}}$ (in addition to the physical one), acting on links between the different sites. Taking as a concrete but representative example the compatible-L case, Eqs. (5) and (13), the Hamiltonian is given by

[TABLE]

The particle currents $J^{i}$ are then given by taking the derivative w.r.t. ${\cal A}$, and taking ${\cal A}\to 0$,

[TABLE]

As can be easily seen, $J^{i}$ is the sum of an $a$-$a$ current, whose definition exactly coincides with our definitions for the particle current in the main text, Eqs. (23a),(23b), and of an $a$-$b$ contribution, describing exchange of particles between the bath and the lattice $a$-particles. The same procedure may be repeated for any choice of $H^{{\rm sys}}$, arriving at a similar outcome.

Let us now turn to evaluate the $a$-$b$ part, and show how it gives rise to the dissipative current, Eq. (40). For this purpose let us use the notation $H^{{\rm diss}}=\sum_{\alpha\beta}b_{\beta}^{\dagger}h_{\beta\alpha}^{{\rm ref}}a_{\alpha}+{\rm h.c.}$, where $\alpha,\beta$ are generalized site indices. According to our present derivation of the $J^{i}$ operators, the operator corresponding to the total outgoing particle flux from a site $\alpha$ is

[TABLE]

Eq. (48) is linear in the dissipative coupling, so in order to be consistent with the derivation of the Lindblad equation itself, which is second order in the system-bath coupling (Goldstein, 2019), we should evaluate it to first-order in $h^{{\rm ref}}$ (remembering that the zeroth-order term vanishes). Using the Kubo formula,

[TABLE]

where in inside the integral the averaging is performed in the $h^{{\rm ref}}\to 0$ state. Expanding the commutator term inside of the integral, and using the separability of the $\left\langle a^{\dagger}a\right\rangle$ and $\left\langle b^{\dagger}b\right\rangle$ correlations for $h^{\mathrm{ref}}=0$, we find

[TABLE]

Now, we employ the fact that b-particles quickly evaporate after escaping the 2D lattice, so that we may set $\left\langle b_{\alpha}^{\dagger}\left(t\right)b_{\beta}\left(t^{\prime}\right)\right\rangle=0$. In addition, by the Markov approximation used in deriving the Lindblad equation, which implies that the bath dynamics is much faster than the system dynamics, one may set $\left\langle b_{\alpha}\left(t\right)b_{\beta}^{\dagger}\left(t^{\prime}\right)\right\rangle=\frac{2\pi}{\hbar}\nu_{0}\delta_{\alpha\beta}\delta\left(t-t^{\prime}\right)$, where $\nu_{0}$ is the density of states of the $b$-particles, as introduced in Sec. II. Plugging these identities in, we obtain

[TABLE]

and the time argument is suppressed, as all the operators are evaluated at the same time. Performing the summation over the dummy index $\beta^{\prime}$, and using the definition in Eq. (30),

[TABLE]

This result reproduces exactly the second half of Eq. (40). The first half, which originates in a constant refilling process captured by $\Gamma^{\left(2\right)}$, is obtained in a similar way, by introducing a refilling reservoir of $c$ atoms with infinitely positive chemical potential. This would give us the remaining $\left\{\Gamma^{\left(2\right)},1-P\right\}$ term. Thus, we have verified that the current operator we calculate is indeed the physical current of particles within the lattice. (Despite the seeming ambiguity in its derivation from the continuity equation, adding a divergence-free field to $j^{x/y}$ would not correspond to any measurable transport current of either particles.)

Appendix C Definition of the conductance tensor

We elucidate in this appendix the definition and physical meaning of the conductance we calculate in this work, as they appear in Eq. (27). We consider a cylindrical geometry, with periodic boundary conditions applied perpendicular to the direction of applied voltage (or more specifically electric field gradient in our case), as demonstrated in Fig. 11a.

For the transverse response, e.g., $G_{yx}$, we consider the total current going around the circumference of the cylinder along a longitudinal cut, averaged around the cylinder. As $G_{yx}$ is the ratio of total current through this cross section to the voltage difference applied, we divide the total current by the length times the applied field, e.g., for $G_{yx}$ we divide by $L_{x}E_{x}$.

The issue of longitudinal current is a bit more subtle, as the notion of a particle current flowing between two leads does not exist in the system we consider. Moreover, in the steady state the total particle number in the system is conserved. We quantify our response by examining the total flux of particles entering the system from one side and leaving it in the other. As can be seen in Fig. 11b–c, the electric field gradient induces an incoming flux of bath particles (negative $j_{{\cal D}}$) on the right side and an outgoing one (positive $j_{{\cal D}}$) on the left. It is then sensible to define this response by the longitudinal current exactly in the middle of the cylinder, since it is equal in amplitude to the total dissipative current integrated over one half of the system.

Lastly, let us comment on a different method of applying a voltage difference between the two edges. Namely, imposing “voltage leads” on the system, by imposing a finite positive chemical potential on a narrow region near one edge, and a negative on near the other edge. The perturbation to the Hamiltonian can be written as $\delta H_{{\rm leads}}=\frac{V}{2}\sum_{l_{y}}\delta H_{l_{y}}$, with $V$ the applied voltage,

[TABLE]

and $L_{{\rm lead}}$ is the size of the leads. As illustrated in Fig. 12, the longitudinal current in this case flows only at the vicinity of the perturbation. Moreover, the balanced flux of particles in and out of the bath (captured by the dissipative current) is localized as well, as bath particles enter and leave the system in a very narrow region. In this case, there is no sense of particle flow between the two sides of the system or bath, and it appears $\delta H_{{\rm leads}}$ is not a good choice for studying the longitudinal current. We note that for the transverse conductance however, one recovers qualitatively the results of using $\delta H_{\boldsymbol{E}}$, including the conductance quantization at the appropriate parameter regime.