Describing many-body localized systems in thermal environments

In this work we formulate an efficient method for the description of many-body localized systems in weak contact with thermal environments at temperature $T$. For this purpose we exploit the representation of the system in terms of quasi-local integrals of motion ($l$-bits) to derive a quantum master equation using Born-Markov approximations. We show how this equation can be treated by using quantum-jump Monte-Carlo techniques as well as by deriving approximate kinetic equations of motion. As an example, we consider the one-dimensional Anderson model for spinless fermions including also nearest-neighbor interactions, which we diagonalize approximately by employing a recently proposed method valid in the limit of strong disorder and weak interactions. Coupling the system to a global thermal bath, we study the transport between two leads with different chemical potentials at both of its ends. We find that the temperature-dependent current is captured by an interaction-dependent version of Mott's law for variable range hopping, where transport is enhanced/lowered depending on whether the interactions are attractive or repulsive, respectively. We interpret these results in terms of spatio-energetic correlations between the $l$-bits.

In this work we develop an efficient formalism with which we can access the steady states of MBL systems weakly coupled to thermal environments. We apply this method to a spinless fermionic Hubbard chain with strong onsite disorder and weak interactions, where we achieve to solve for mesoscopic system sizes of the order of 100 lattice sites. Utilizing the l-bit representation of MBL systems in terms of quasi-local integrals of motion, we derive a quantum master equation by means of a standard Born-Markov approximation [34]. While deriving explicitly the lbit representation is in general a demanding task, here, we make use of a recently proposed efficient approximate method controlled in the limit of strong disorder and weak interactions [35]. We show how this derived quantum master equation can be solved by means of quantum-jump Monte-Carlo simulations [36]. In order to push our description to even larger systems we use a further simplification of the quantum master equation in terms of a kinetic theory, which as we show provides also an accurate description with the advantage that it allows us to reach system sizes of the order of 100 lattice sites.
The coupling of an MBL system to an environment induces relaxation washing out information about initial conditions and thereby destabilizes the nonergodic nature of the MBL phase. However, the dynamics in the resulting steady state can still be dominated by the localized character of the system Hamiltonian leading to slow and constrained dynamics. A prominent example constitutes variable-range hopping, which is well-studied for Anderson localized systems of noninteracting particles [37]. The coupling to a thermal bath induces nonvanishing transport coefficients, while the conductivity σ remains highly suppressed at low temperatures T following Mott's law, where d denotes the spatial dimension and T 0 the temperature scale at which Mott's law sets in. We take variable-range hopping as the starting point for our considerations. For that purpose we couple our studied Hubbard chain to a bosonic thermal bath at temperature T . Further, we couple the chain at both ends to fermionic reservoirs with a chemical potential difference such as to induce a current flowing through our system. This induces a nonequilibrium current-carrying steady state, which reduces to the conventional linear response regime in the limit of vanishing chemical potential difference.
We first study a noninteracting Anderson insulator in one spatial dimension. Starting from our microscopic model, we map out the full temperature dependence of transport. At asymptotically low temperatures we observe that the induced current becomes temperature-independent characterized by long-range hopping particles across the full chain, occurring with an amplitude which is exponentially suppressed as a function of system size. Upon increasing temperature beyond this finite-size dependent asymptotic regime, we recover a conductivity following Mott's law, which we take as a first indication that our approach captures the relevant physics. At even larger temperature, transport crosses over to a simple activated regime where σ ∝ exp[−T 1 /T ]. There, particles can overcome typical energy barriers imposed by the strong disorder potential on short distances as opposed to variable-range hopping which is characterized by longer-ranged tunneling processes.
As a next step, we then apply our approach to the weakly interacting case at strong disorder. In the temperature regime, where for the Anderson system we have observed variable-range hopping, we again find that transport follows Mott's law, however, the Mott temperature T 0 is modified by interactions. Depending on whether interactions are repulsive or attractive, T 0 increases or decreases with respect to the noninteracting Anderson insulating case. We interpret the modifications of T 0 in terms of interaction-induced changes in the spatio-energetic correlations among the local integrals of motion in the MBL system.
The paper is organized as follows. In Section 2, we present a general scheme to address MBL systems in the presence of thermal baths. This is followed by the description of the model under consideration (Section 3). In Section 4, we introduce the quantum-jump Monte-Carlo method and the kinetic theory that we utilize to solve the master equations that describe the dynamics of the open quantum system. We then discuss the dependence of transport in the noninteracting system in Section 5. In Section 6, we address transport through the system in the presence of weak interactions. Finally, a summary of the main results is given in Section 7 to conclude.

Master equation for many-body localized systems in the presence of thermal baths
When a quantum system is brought into weak contact with a thermal environment, the impact of the environment can be captured by quantum jumps between different eigenstates of the system. The rate at which such a process occurs depends on the energy difference of the corresponding states. For most interacting quantum systems, however, the full spectrum is hard to obtain, so already deriving the equations of motion of the system in the presence of a thermal environment is challenging. We approach this problem by exploiting the fact that MBL systems posses an extensive number of integrals of motion [38][39][40][41]. This makes them very special in the sense that they are interacting quantum systems whose full many-body spectrum can be accessed. We make use of this fact to find equations of motion for MBL systems in thermal environments.

Born-Markov master equation
We are interested in the typical setup of open quantum systems described by the total HamiltonianĤ Here,Ĥ denotes the Hamiltonian of the system under investigation,Ĥ B = αh ω αb † αbα describes the bath with bosonic or fermionic annihilation operatorb α for the mode with energyhω α , andĤ SB is the system-bath coupling operator, whose overall strength γ we assume to be small. The bath is assumed to be in thermal equilibrium with temperature T and chemical potential µ.
In the weak-coupling limit (γ → 0) and assuming large thermal baths, the dynamics of the reduced density matrixρ = Tr B {ρ tot } of the system after tracing out the bath, can be described using the Born-Markov and the rotating wave approximation. This gives rise to a Lindblad master equation [34] where {·, ·} denotes the anti-commutator. The first term describes the unitary evolution, and the second term describes the coupling of the system to the environment. Here the operatorsL α describe quantum jumps between the eignestates of the system. Let us assume a system with spectrum E k and eigenstates |k coupled to a phonon bath with bosonic operatorsb α via the system operatorv, Then, the master equation takes the form The bath induces quantum jumps from level q to k, mediated by the action of the jump operatorsL with the jump rates Here we have defined the bath correlation function with the Bose distribution n B (E, T ) = (e E/kBT + 1) −1 and the spectral density of the bath J(E) = α λ 2 α δ(E −hω α ). Note that the Born-Markov and rotating-wave approximation rely on the assumption of a separation of timescales [34]. For them to be valid, the timescales of the system dynamics τ S ∝ |E k − E q | −1 need to be fast when compared to the timescales of dissipation τ R ∝ R −1 kq (we seth = k B = 1 hereafter). This implies Under the dynamics of Eq. (5), the coherences, k|ρ|q with k = q, decay such that the steady state reached asymptotically in the long-time limit is diagonal in the energy basis,ρ →ρ ∞ = k p k |k k| [42][43][44]. This can be seen most easily by noting that in the limit γ → 0, the first term in Eq. (3) has to vanish for the steady state with ∂ tρ∞ = 0. After a diagonal state is reached, the asymptotic dynamics is governed by dissipation and described by the Pauli rate equation Thus, the dissipative part of the master equation determines the probability distribution p k of the steady state. When the system is coupled to a single bath of temperature T , the rates obey the equilibrium condition R kq /R qk = e −(E k −Eq)/T . In this case the system approaches thermal equilibrium in the long-time limit, p k ∝ exp(−E k /T ), and the system obeys detailed balance: the terms of the sum in Eq. (10), which correspond to the net probability currents from states q to k, vanish individually. If the equilibrium condition is broken, e.g. when the system is coupled to baths of different temperature or chemical potential, the system approaches a nonequilibrium steady state, for which the right-hand side of Eq. (10) vanishes only as a whole, so that probability currents break detailed balance.

Many-body localized systems
Note that in order to find the equations of motion for the system in the presence of a thermal bath, it is essential to diagonalize the full Hamiltonian to find its eigenstates and eigenenergies. For a quantum many-body system, this is generally difficult. Non-interacting systems however can be treated readily [45], since the many-body eigenstates are given by the Fock states of the single-particle Hamiltonian. As we discuss in the following, MBL systems provide another exception, due to an emergent form of integrability [38][39][40][41]. We consider a system of spinless lattice fermions. A typical example for such a fermionic system with an MBL phase is the one-dimensional Anderson model of spinless fermions with nearest-neighbor interactions V , which we introduce in Section 3. The generalization to spinful fermions or spin Hamiltonians is straightforward. We assume that the system Hamiltonian can be brought into the diagonal formĤ where the fermionic number operatorsn k describe quasi-local integrals of motion (l-bits) coupled by diagonal matrix elements that decay exponentially with the localization length ξ [46], U k1k2...kn ∼ e − mini,j |ki−kj |/ξ (assuming the l-bits to be ordered according to their spatial position). We can directly read off the full manyparticle eigenstates, which are the Fock states with respect to the l-bits |n = |{n k } , n k = 0, 1, and the corresponding many-particle spectrum E n = E({n k }).
We focus on the limit of weak interactions, where it is assumed that there exists a quasi-local transformation (adiabatic connection) between the l-bits and the local annihilation operatorsâ i for fermions on lattice sites i [4,40] such that the l-bit operators readn k =ĉ † kĉ k . In the MBL regime with strong disorder, the coefficients u i , u ijk essentially have a local support with some exponential tails decaying on length scale ξ. It can be shown [35] that the terms that lead to particlehole dressing u ijk and all higher order terms contribute at least in first order in the interaction strength V , so for weak interactions V and strong disorder the quasiparticles are in leading order given by a single-particle transformation.

MBL system in the presence of a phonon bath
Let us now discuss the coupling of an MBL system to a phonon bath with a coupling operatorv that is some single-particle operatorv = ij v ijâ † iâ j , The form ofĤ SB induces quantum jumps between eigenstates ofĤ and the corresponding dissipator is given by withL n n = √ R n n |n n| describing the jump from Fock state |n to the state |n . The associated jump rate is then given by with g(E) the bath correlation function [Eq. (8)]. Note that generally, in the full spectrum of a quantum many body system one expects that a large number of close degeneracies occur, which violates the validity condition in Eq. (9). However, for MBL systems in the limit of weak interactions and strong disorder that we are aiming at, we will find significant rates only for processes, where one or at most a few excitations are transferred. Moreover, typically the bath will transfer excitations only between nearby l-bits, since v ij decays on the correlation lengths of the bath. Both together leads to a much milder condition. For weak interactions and strong disorder, particle-hole dressing is suppressed, such that the transformation to the l-bits effectively becomes a single-particle In this case, the dissipator simplifies to withL qk (n) = R n qk ,n |n qk n| describing the jump from Fock state |n = |n 1 , . . . , n k , . . . , n q , . . . , n M to the state |n qk = |n 1 , . . . , n k − 1, . . . , n q + 1, . . . , n M by transferring a particle from the single-particle mode k to the mode q. The fact that we neglect particle-hole dressing leads to with the single-particle matrix element of the coupling operator For convenience, we introduce an effective single-particle ratẽ which in the case noninteracting case, with E n qk − E n = ε q − ε k , reduces to the singleparticle rateR qk (n) = R qk . This allows for the extraction of quantum statistical factors from the many-particle rate This expression resembles the expression found for the ideal gas [45], where the many-particle rate is the single-particle rate R qk multiplied by the occupation of the departure state n k and the Pauli blocking factor (1 − n q ) of the target state. The difference is that due to interactions, the transition rate depends on the whole configuration, rather than only on the two single-particle states involved in the transition.

MBL system in the presence of a particle reservoir
Let us now turn to the case where the system may also exchange particles with an external fermionic reservoir with temperature T and chemical potential µ. For simplicity, we again consider the regime where we can neglect particle-hole dressing.
Here the system-bath Hamiltonian is given bŷ whered = i ϕ(i)â i is an arbitrary single-particle mode in the system. The resulting dissipator reads whereL +,k (n) = Γ +,n |n n k↓ | andL −,k (n) = Γ −,n |n k↓ n|, with |n k↓ = |n 1 , . . . n k − 1, . . . being the Fock state obtained by removing one particle from mode k. The jump rates Γ ±,n are given by with the coupling matrix element and the bath-correlation function for the particle exchange with the fermionic bath where f (E, µ, T ) = (e (E−µ)/T + 1) −1 is the Fermi distribution of the bath. Figure 1: Schematic illustration of the model under consideration. An open tightbinding chain of interacting spinless fermions (green dots) in a random potential (brown curve) is coupled to a global heat bath at temperature T (blue box). This global bath induces phonon (blue dot)-assisted heat exchange with the system and two local baths (yellow boxes) at the ends with different chemical potentials (µ L and µ R ), which induce particle transport through the chain.

Interacting Anderson insulator as quantum wire
We benchmark the explained method studying a disordered extended Hubbard chain of spinless fermions (Fig. 1). The system Hamiltonian [47][48][49] readŝ with the single-particle term given bŷ whereâ i (â † i ) is the fermionic annihilation (creation) operator at site i in a chain of length M . Moreover, w i are random fields uniformly distributed in the interval [−W, W ], J and V are respectively the hopping and the interaction strength.Ĥ 0 describes the non-interacting Anderson-model [50] (see appendix 8.3) and all its singleparticle eigenstates are exponentially localized for any amount of disorder W [51]. For V = 0 the model has an MBL transition [3][4][5][6]. At large disorder W V , the system is in the MBL phase and it is a perfect insulator, while for weak disorder W V the system is ergodic, describing a thermal phase. We couple the system to two leads at the ends of the chain. The leads have the same temperature T but different chemical potentials, µ L and µ R , which exchange particles with the system, inducing thus particle transport through the chain. Moreover, the system is also coupled to a global thermal bath also at temperature T (unless stated otherwise), with which it exchanges energy, as sketched in Fig. 1.
We are interested in the limit of weak interactions at strong disorder, where in leading order l-bits can be approximated by the non-interacting Anderson operators [35], so that we arrive at the approximate Hamiltonian Heren k =ĉ † kĉ k withĉ † k = i ψ k (i)â † i , and ε k , ψ k are the single-particle eigenenergies and eigenmodes, respectively. The interaction coefficients U kq are given by Since the orbitals are exponentially localized, ψ k (i) ∝ e −|i−i k |/ξ loc (k) , with localization center i k and the single-particle localization length ξ loc [52], after a relabeling of the indexes (k, q) according to their spatial positions, we have U kq ∼ V e −|k−q|/ξ loc (here ξ loc is the localization length in the middle of the single-particle spectrum). This approximation [Eq. (30)] is equivalent to discarding the off-diagonal elements of the full Hamiltonian [Eq. (28)] in the non-interacting basis. Its reliability for weak interactions has been shown in Ref. [35]. Its validity is strongly supported also by the spatio-energetic anti-correlations between the single-particle Anderson orbitals. If two single-particle eigenstates are close to each other in space, their energy difference is more likely to be large, suppressing scattering events between these states. Figure. 2 gives evidence of the existence of this anti-correlation, showing the distribution of the spatial distance between pairs of eigenstates ∆R = | i i(|ψ k (i)| 2 − |ψ q (i)| 2 )| as a function of their energetic distance ∆E = |ε k − ε q |. The distribution is mainly concentrated at the left bottom corner of Figures (a) It is important to point out that a better approximation for the l-bits can be obtained systematically by treating the omitted matrix elements perturbatively, giving rise to higher order corrections in interaction strength [28]. Such a procedure is beyond the scope of this paper but it would not render our approach inefficient. We couple this system to a heat reservoir via the system-bath Hamiltonian which couples each site to an independent bath of bath-correlation length zero. C α and γ are the coupling strengths. Furthermore, we assume an Ohmic bath, with spectral density J(E) proportional to E, J(E) ∝ E. The total rate that enters in Eq. (17) is therefore the sum of the rates for individual heat baths that couple to a single site i only.
We investigate the transport properties of the steady state, when the system is coupled to the two leads and the thermal bath. The lead coupling is described bŷ The resulting coupling rates η L (k) = |ψ k (1)| 2 , η R (k) = |ψ k (M )| 2 that enter in Eq. (25) are proportional to the value of the wavefunctions at the ends of the chain. For the leads we assume a constant density of states. The effective single-particle energies are shifted due to interactions. Within the scope of the approximate Hamiltonian in Eq. (30), it is convenient to define an interaction-shifted energy operator for each l-bit The energy difference occurring in the rate for the heat bath, Eq. (20), can then be expressed as For the particle reservoir, the energy difference in Eq. (24) is given by Note that writing energy differences in terms ofε k does not constitute an additional approximation.

Solving the Master equation
We employ two methods to compute the non-equilibrium steady state of the master equation. The first method is a quantum-jump Monte Carlo technique, which makes use of the fact that the equations of motion for the Fock-space occupation probabilities p n can be mapped to a classical random walk. This method gives accurate results after sufficient statistical sampling. The second method consists in deriving kinetic equations of motion by employing a mean-field decomposition of density-density correlations. Importantly, both approaches are found to agree very well.

Method I: Quantum-jump Monte Carlo simulation
As we discuss in Section 2.1, the dynamics of the many-body occupation probabilities p n = n|ρ|n decouples from the off-diagonal elements, which decay over time. For our model, the equation of motion for the probabilities is given by whereε k = n|ε k |n , and we use the short notation f α (E) = f (E, µ α , T α ) with α = L, R. The first line in Eq. (38) describes the heat exchange processes, with the first term denoting the increase of the probability p n by jumping from |n qk to the state |n , and the second term denoting the inverse process. The second line in Eq. (38) describes particle exchange with the leads, with the first term denoting particle gain, and the second term the particle loss.
We generalize the quantum-jump Monte-Carlo method described in Ref. [36] to the case of number-dependent single-particle ratesR qk (n). The dynamics of the occupation probability p n is simulated by taking random walks in the classical space composed by the Fock states |n = |n 1 , . . . , n k , . . . , n M (not their superpositions) [45] corresponding to the given jump rates. In our case, we have two types of jumps, the first one due to the global thermal bath and the second one due to the leads. For the heat exchange processes, a jump transfers one particle from one single-particle eigenstate to another, which conserves the total particle number. For the particle exchange process, a jump adds (or removes) a particle to (or from) one eigenstate. We perform these simulations by using a Gillespie-type algorithm (see appendix 8.1). We can then compute steady state expectation values (e.g. n k , n knq ), by averaging over the long-time dynamics of many trajectories.

Method II: Kinetic theory
In order to treat larger systems and also to gain some intuitive understanding of the non-equilibrium steady state (NESS), we will now derive the kinetic equations of motion for the mean occupation numbers n l . The time evolution of n l , again for the special case of our model system, is given by with (see appendix 8.2) and In Eq. (40), we have defined an effective rate operator Using Eq. (36), we findR The effect of interactions is to modify the energy difference between the states (ε l −ε k − U lk ), implying that the transition rate depends on the whole configuration, rather than only on the two single-particles states involved in the transition (ε l − ε k ).
In order to obtain a closed set of equations in terms of mean occupation values, we employ the mean-field approximation In this way, we get a set of nonlinear kinetic equations of motion In the following, we will always consider the steady state and, therefore, drop the subscript NESS. Using the single-particle wavefunctions ψ l we can find the mean occupation number (density profile) in real space The nonlinear term in Eq. (45) prevents us from finding an analytic solution. However, in the absence of the global thermal bath (γ = 0), the solution of Eq. (46) is given by Let us compare the methods described in Sec. 4.1 and 4.2 by computing the density profile n i , which also determines the current to be discussed later. Figure 3 shows the density profile n i obtained by numerically solving the kinetic equations [Eq. (46), solid lines] and the one calculated using the quantum-jump Monte-Carlo technique [Eq. (38), markers]. The two results are almost indistinguishable both for the non-interacting and the interacting case. The chemical potentials in the leads are taken equal to µ L,R =μ ± δµ/2 with relatively large chemical potential difference δµ = 5J, andμ is fixed so that the system is at half-filling ( N = M/2, wherê N counts the total number of particles). Later on, we will use only small chemical potential offsets δµ, in order to compute conductivities. The large difference between the chemical potentials produces a density gradient in the absence of the global thermal bath (γ = 0), as shown in Fig. 3 (black dashed lines). The presence of the global thermal bath, which induces phonon-assisted tunneling, erases this gradient for weak disorder, as shown in (a) and (c). Furthermore, by increasing the disorder strength W , which prompts the localization of the wavefunctions and thus the suppression of particle tunneling, the density gradient is restored, as shown in Fig. 3(b) and (d).

Results I: Observation of variable-range hopping for noninteracting Fermions in a microscopic model
We investigate the particle current in the NESS, which is defined as the net particle flow at the ends of the chain Having shown that the kinetic equations, Eq. (45), are reliable, we now use this approach to study the dependence of the current I, Eq. (49), on disorder strength W , on the temperature of the bath T , and on interaction strength V . As shown in Appendix 8.4, for the current, the results of the kinetic equations also agree very well with those from Monte-Carlo simulation. In the following, we will only show the results of the kinetic theory with disorder averaged.  50)]. The parameters are system size M = 100, dissipation rate κ = 0.1γ 2 , interaction strength V = 0, and µ L =μ + δµ/2, µ R =μ − δµ/2, with chemical potential imbalance δµ = 0.1J, andμ is set to make N = M/2 = 50. Except when otherwise stated, we will take the same parameters for the rest of the work.

Transport as a function of disorder strength W
In the absence of the thermal bath (γ = 0), the current in the NESS is given by which is obtained by substituting the solution of the occupation numbers Eq. (48) into Eq. (49). The terms of I γ=0 are proportional to the product η L (l)η R (l) ∝ |ψ l (1)| 2 |ψ l (M )| 2 ∼ exp(−4M/ξ loc ), where for strong disorder ξ loc ∼ 1/ log (W ). This implies that without coupling to a thermal environment the current will dramatically decay with disorder strength. This is confirmed in Fig. 4 (dashed lines), where the behavior of the current is shown as a function of W for the non-interacting chain of M = 100 lattice sites. Here and in the following, we consider the weak chemical potential difference δµ = 0.1J. Figure 4 also reports the case in which the thermal bath is present (γ = 0), showing that the presence of a global thermal bath enhances transport for sufficiently large T .

Transport as a function of temperature T
In this Section we study the dependence of the current I on temperature T and disorder strength W for the noninteracting case. Figure 5 (VRH) [37], I ∝ exp(− T 0 /T ). Finally, at 'high' temperatures (T J), I decreases as a function of T .
Note that the locations of the three regions depend on various parameters, such as the disorder strength W (as shown in Fig 5(a)), the coupling to the leads κ, and so on. The independence of the current on temperature in the 'low' temperature regime is due to the suppression of heat exchange, which can be inferred from the good agreement of the results with that without coupling to the global thermal bath [Eq. (50), dotted lines]. Nevertheless, the remaining current (T → 0) is a finite system size effect [53]. Indeed, as T → 0, I γ=0 decays exponentially with system size M , as shown in Fig. 6(a). For a given system size, I γ=0 depends on the disorder strength, and as expected it is smaller for stronger disorder [ Fig. 5(a)].
The intermediate-temperature VRH regime can be understood using the following argument due to Mott [37]. Let's consider the hop between states with localization centers separated by the distance ∆R and with energy difference ∆E. On the one hand, the probability to hop is proportional to the envelop overlap between the two states, thus it decays exponentially with the distance ∆R [∼ exp(−2∆R/ξ loc )]. On the other hand, the probability to produce excitations of order ∆E due to the presence of the heat bath is given by the Boltzmann factor exp(−∆E/T ). This leads us to assume that the current (conductivity) to leading order is given by As already mentioned, the spatial distance ∆R and the energy separation ∆E, are not independent, but show clear anticorrelations (Fig. 2), which shall be captured by ∆R ∼ (∆Eν) −1 , where ν is the density of states. Thus, the current is determined by the competition between the overlap term exp(−2∆R/ξ loc ) which favors short hops and the energy activation exp(−∆E/T ), which favors long hops. Maximizing this probability over ∆E, one finds in one spatial dimension with Mott's hopping length ∆R 0 = ξ loc /(2T ν) and Mott's temperature T 0 = 2/(ξ loc ν). Moreover, Mott's hopping energy ∆E 0 = √ T 0 T is the energy scale that defines the width of the energy interval of the activated eigenenergies.
By fitting the current in the intermediate temperature regime to Mott's law I = I s exp(− T 0 /T ), we can extract T 0 . Figure 6(b) shows T 0 as a function of ξ loc ν, where ξ loc is the single-particle localization length in the middle of the energyband (ε k ∼ 0). We find that at strong disorder, T 0 ∝ (ξ loc ν) −1.29±0.25 , agreeing rather well with the Mott's prediction ∝ (ξ loc ν) −1 .
Finally, in the 'high' temperature regime, we attribute the decrease of the current with respect to temperature to the fact that the difference between the Fermi distributions of the leads (f L ( ε l ) − f R ( ε l )) is washed out. To support this idea, we fix the temperature of the leads at T L = T R = 0.01J. In this case, the current I increases with the temperature of the global thermal bath T as exp(−T 1 /T ) in the 'high' temperature regime, as shown in Fig. 7(a). The reason is that the thermal energy T is so high that, despite of their large energy separation, already activated hopping between nearest neighboring localized orbitals dominates over the variablerange hopping behavior. From the current, we can deduce the conductance G = I/δµ, as well as the conductivity Figure 7(b) shows the conductivity σ at temperature T = 0.5J as a function of the system size M for W = 5J. We can observe that σ converges to a finite value in the limit M → ∞.

Results II: Variable-range hopping for interacting Fermions
Let us now discuss transport in the presence of weak interactions, for which the system is in the MBL phase. Here, we exploit the advantages of our method, allowing for the treatment of interacting systems in the presence of a thermal environment. All of our results are obtained in the limit of weak interactions and strong disorder in which the approximation [Eq. (30)] is justified.
6.1. Transport as a function of temperature T Figure 8(a) shows the temperature dependence of the current for several interaction strengths V . Interestingly, also for the interacting case (V = 0), we find a regime of temperatures, where the current is explained by Mott's law (I ∼ exp(− T 0 /T )). We find that in this VRH regime, attractive interactions (blue dashed line) enhance transport, while for repulsive interactions (orange dash-dotted line) the current decreases. This means that the interaction changes Mott's temperature T 0 . Fig. 8(b) shows the dependence of Mott's temperature T 0 on the interaction strength V . We can see how attractive interactions decrease T 0 , while repulsive interactions increase T 0 .
Can we understand this behavior? From our previous discussion we know that T 0 depends on the localization length ξ loc and on the density of states ν. However, a change of the l-bits and their localization length is expected to be of second order with respect to interactions (and is also not taken into account in the approximate Hamiltonian (30)). We should, therefore, be able to explain the interaction-induced shift of Mott's temperature in terms of interaction-induced shifts of the density of states.

Interaction-shifted density of states
In Fig. 9(a), we show the energies ε k in the steady state for V = −0.4J, V = 0 and V = 0.4J all at T = 0.5J. As expected from Eq. (35), repulsive interactions (V > 0, orange dash-dotted line) shift the single-particle energies ε k up (black solid line) and attractive interactions (V < 0, blue dashed line) decrease them. However, this does not imply a change in the average level splitting ε k −ε k−1 . What is crucial is rather that the interaction-dependent energy shift depends on k. The absolute value of the energy shift δE k ≡ ε k − ε k is shown in Fig. 9(b). For the eigenstates in the middle of the spectrum, which are the ones that contribute predominantly to the transport, the energy shift increases with energy.
This behavior can be understood from considering Eq. (35): The energy shift δE k = q U kq n q depends via U kq on the overlap of the involved single-particle wave functions. Therefore, the main contribution to the shift originates from eigenstates that are close by in space. Now, due to the anti-correlation property between energetic and spatial distance, these eigenstates q will have a large energy difference with respect to k. Thus a state k slightly above the Fermi energy will have more likely neighboring states below the Fermi energy with a large occupation probability, while a state slightly below the Fermi energy will more likely have neighboring states above the Fermi energy, with a small occupation probability. In this way, the positive/negative energy shift of repulsive/attractive interactions, will be larger for states above the Fermi level than for states below it. This implies that the level spacing between neighboring energy levels ( ε k −ε k−1 ) is increased for repulsive interactions and is decreased for attractive interactions. In other words, repulsive interactions reduce the density of states ν k = 1/( ε k −ε k−1 ) and attractive interactions enhance it.

Explanation of the change of Mott's temperature
The suspected behavior is confirmed in Fig. 10(a), where we show the averaged density of states (over states with k between 40 and 60),ν, as a function of interaction strength at temperature T = 0.5J. We can now check whether the dependence of Mott's temperature T 0 on interactions can be explained with their effect on the density of states. When we plot T 0 versusν, as shown in Fig. 10(b), we observe a behavior T 0 ∝ν −1.18±0.15 , which agrees well with the predicted ν −1 dependency in Mott's temperature.
We would like to emphasize that the origin of Mott's law we find is different from the case of Coulomb interactions, which has been shown to lead to a variable-range hopping conductivity σ ∝ exp[−(T 0 /T ) 1/2 ] independent of dimensionality due to a nonanalytic modification of the density of states near the Fermi energy [54]. In our case of local interactions of nearest-neighbor type the density of states becomes modified only smoothly which continuously connects the interacting with the noninteracting result.

Conclusion
In this work we have introduced a method which allows for an efficient description of MBL systems at strong disorder and weak interactions, when weakly coupled to thermal environments. We have benchmarked our method for noninteracting systems, where we showed that our method recovers Mott's law from variable-range hopping starting from a microscopic model. Upon adding weak interactions we have found that Mott's law persists while leading to perturbative corrections. We explain our observations by an interaction-induced modification of the density of states due to spatio-energetic correlations.
Already in recent years the study of MBL systems coupled to environments has received substantial interest not only in experiments [24,32]. On the theoretical side these previous studies [26-31, 33, 55-62] have been either limited in the accessible system sizes [27,29,33,55] or make rather specific assumptions regarding the properties of the environment, e.g. it is described by dephasing noise corresponding to an infinite temperature environment [27-29, 31, 56], that it is described by classical noise [30,57], or that it is given by a small bath of size comparable to the system [26,[58][59][60][61]. In the present work we provided an alternative approach, which allows for an efficient solution for mesoscopic system sizes in contact with conventional thermal environments. This might also be of particular importance in view of experiments in quantum simulators, where the influence of environments has been studied systematically recently [24,32].
Within the presented formulation the method is limited to MBL systems at strong disorder and weak interactions. However, this limitation originates solely from the method used to diagonalize the system Hamiltonian following the recent work [35]. This can be improved upon finding a more accurate diagonalization of the MBL Hamiltonian, i.e., a more accurate l-bit description, following, e.g., some recent ideas [63,64]. This would not influence substantially the structure of the quantum master equation, but would allow to study the steady states in broader parameter range, e.g., also closer to the MBL transition towards ergodicity.
Concluding, our work provides a framework to study open system dynamics at mesoscopic scales for various scenarios involving MBL systems at strong disorder and not too strong interactions. This includes a wide range of phenomena such as MBLspin glasses, MBL topological phases, or time crystals.
(t h ), gaining (t g ) and losing (t l ) a particle given bȳ According to the choice made, the corresponding jump operation is performed. Specifically, if τ = τ h , which means heat exchange process is chosen, then a particle is transferred from a randomly drawn departure mode k to the randomly drawn target mode q. This single-particle jump has the probability P (k → q, n) =t h R n qk ,n . If pumping is chosen, with τ = τ g , a particle is added to a mode randomly drawn with probability P (k) =t p α=L,R η α (k)f α (ε k ). Likewise, the particle loss process is performed when τ = τ l . These two steps are repeated until the desired evolution time is reached. An ensemble of trajectories is calculated individually, from which the wavefunctions |n (λ) (t) obtained are then used to compute the expectation value of an observable O as

Derivation of Eq. (40) and (41) in the main text
The time evolution of the mean occupation due to heat exchange is governed by with 2L † qk (n)n lLqk (n) −n lL † qk (n)L qk (n) −L † qk (n)L qk (n)n l = 2|n n qk |n l |n qk n| −n l |n n| − |n n|n l (57) Let us assume that Fock state |n has n l particle in l mode, i.e., n l |n = n l |n .
Then it is easy to verify that for state |n qk , which is obtained from state |n by transferring a particle from k mode to q mode, there is the following propertŷ n l |n qk = (n l + δ q,l − δ k,l )|n qk .
It then reduces to Eq. (40) by using Eq. (21) and (42) in the main text. Likewise, by making use of 2L knlL † k −n lLkL † k −L kL † kn l = 2δ k,l |n k↓ n k↓ |, 2L † kn lLk −n lL † kL k −L † kL knl = −2δ k,l |n n|, we can obtain Eq. (41) in the main text. In this section we show some further data concerning the non-interacting system. A straightforward basis transformation recastsĤ 0 into the diagonalized formĤ 0 = k ε knk withn k =ĉ † kĉ k andĉ † k = i ψ k (i)â † i , which creates a particle in the singleparticle eigenstate |k = i ψ k (i) |i with energy ε k that are distributed between −2J − W and 2J + W , as shown in Fig. 11 (e). Figure 11(a), (b) show some selected eigenstates for two different disorder strengths. It is clear that as disorder becomes stronger, the wavefunctions become more localized. To characterize the localization of the wavefunction, we define the localization length using the inverse participation ratio ξ loc ≡ 1/ i |ψ k (i)| 4 . As shown in Fig. 11(c), ξ loc decays rapidly as disorder increases. Moreover, it also depends on the eigenenergy, as shown in Fig. 11(d), and ξ loc has its maximum in the center of the energy-band [ Fig. 11(e)]. The density of states ν k ≡ 1/(ε k − ε k−1 ) (inverse of the energy gap between neighboring eigenstates) also depends on the disorder strength W . From Fig. 11(f), it is clear that larger values of W lead to a smaller density of states. Figure 12(a) shows the current as a function of disorder strength W . The kinetic theory results (solid lines) agree well with the results obtained from quantum-jump Monte-Carlo simulation (markers). The agreement is better for lower temperature and weaker disorder. The reason is that the error of the Monte-Carlo simulation scales as ∆I/ √ N M C , with N M C being the number of trajectories for the Monte-Carlo simulation and ∆I the fluctuation of the current. From (b) we can see that, the fluctuation of the current ∆I increases with disorder strength W and temperature T . In addition, as shown in (a), the mean value of the current decreases with increasing W . These imply that to maintain a small relative error, much more trajectories are needed for stronger disorder and higher temperature. In other words, for a given number of trajectories, which is 1000 in our numerical calculation, the error of Monte-Carlo simulation will be larger for stronger disorder and higher temperature.