Tensor-network method to simulate strongly interacting quantum thermal machines

We present a methodology to simulate the quantum thermodynamics of thermal machines which are built from an interacting working medium in contact with fermionic reservoirs at fixed temperature and chemical potential. Our method works at finite temperature, beyond linear response and weak system-reservoir coupling, and allows for non-quadratic interactions in the working medium. The method uses mesoscopic reservoirs, continuously damped towards thermal equilibrium, in order to represent continuum baths and a novel tensor network algorithm to simulate the steady-state thermodynamics. Using the example of a quantum-dot heat engine, we demonstrate that our technique replicates the well known Landauer-B\"{u}ttiker theory for efficiency and power. We then go beyond the quadratic limit to demonstrate the capability of our method by simulating a three-site machine with non-quadratic interactions. Furthermore, we demonstrate the capability of our method to tackle complex many-body systems by extracting the super-diffusive exponent for high-temperature transport in the isotropic Heisenberg model.


I. INTRODUCTION
The miniaturisation of technologies in combination with the exquisite control now available over nanoscale systems has motivated increasing interest in thermal machines that operate in the quantum regime [1][2][3][4][5]. While recent demonstrations with trapped ions [6][7][8][9], nanomechanical oscillators [10] and diamond colour centres [11] serve as impressive proofs of principle, practical applications such as thermoelectric power generation call for electronic devices. To that end, the focus of experiments in mesoscopic physics has expanded beyond traditional questions of charge transport to include the manipulation of heat currents in platforms such as semiconductor quantum dots [12], superconducting circuits [13] and molecular junctions [14]. Understanding the non-equilibrium thermodynamics of these systems is a formidable theoretical challenge, due to the simultaneous presence of strong system-reservoir coupling, interparticle interactions and finite temperatures.
Existing approaches to modelling energy transport in complex quantum systems typically depend on perturbative arguments, which require a clear separation of energy or time scales. For example, a quantum master equation can be derived under the assumption of weak system-reservoir coupling [15]. However, the approximations needed to ensure positivity of the density matrix may fail to capture quantum coherences far from equilibrium [16][17][18][19], while a first-principles derivation requires full diagonalisation of the system Hamiltonian and thus becomes infeasible for large open systems. A more tractable approach for many-body problems is a local master equation, where incoherent sinks and sources * Corresponding author: brenesnm@tcd.ie create and remove excitations at the system's boundaries. This method has been successfully applied to study infinite-temperature transport in strongly interacting systems [20], but its finite-temperature predictions may violate basic thermodynamic laws [21][22][23][24] unless a specific kind of periodically modulated system-bath interaction is assumed [25][26][27][28][29]. Alternatively, non-equilibrium Green functions [30] can be used to model energy transport under strong system-reservoir coupling, but at the cost of treating many-body interactions within the system perturbatively [31,32]. Another possibility is the numerical renormalisation group, which can handle strong interactions but is typically limited to near-equilibrium transport properties [33]. The related chain representation of unitary system-bath dynamics [34] is also capable of non-perturbative transport calculations [35] but its scalability to large system size remains unclear.
In this work, we put forward a general and efficiently scalable numerical approach to quantum thermodynamics that can deal with simultaneously strong intra-system and system-bath interactions and which works arbitrarily far from equilibrium. We focus on autonomous thermal machines, where macroscopic fermion reservoirs held at different temperatures and chemical potentials drive currents through a complex quantum working medium. We model the macroscopic reservoirs by a finite collection of fermionic modes that are continuously damped towards thermal equilibrium by an appropriate Lindblad master equation. We use a purification scheme based on auxiliary "superfermion" modes [36] to compute the nonequilibrium steady states of both non-interacting and interacting working media. For interacting systems, we develop a tensor-network algorithm to efficiently simulate the real-time dynamics of the entire configuration, working directly in the energy eigenbasis of the reservoirs. Our approach is well suited to far-from-equilibrium problems in which all energy scales are comparable, such that per-turbative or linear-response theories fail. To exemplify this, we demonstrate that the efficiency of a three-site quantum heat engine is enhanced by repulsive interactions and is further improved when the system-reservoir coupling is increased.
The concept of modelling infinite baths by a finite set of damped modes has been widely adopted and adapted since the seminal work of Imamoglu [37] and Garraway [38,39]. In the context of open quantum systems coupled to bosonic reservoirs, this representation has recently been placed on a mathematically rigorous footing [40,41], while its amenability to tensornetwork simulations has been demonstrated [42]. Related approaches have been used to study quantum heat engines [43,44] and thermalisation in few-level [45] and many-particle systems [46,47]. In the fermionic setting, a rigorous proof of equivalence between continuum baths and mesoscopic reservoirs comprising damped modes is not yet available. Nevertheless, such mesoscopic reservoirs have been used quite extensively for studying transport in non-interacting systems [36,[48][49][50][51][52], including under time-dependent driving fields [53]. For interacting systems, a mesoscopic-reservoir description was recently applied to study particle transport and Kondo phenomena in impurity models [54,55], while a related approach to simulating non-equilibrium many-body problems via an auxiliary master equation has been reported [56,57].
A key feature of our work that differs from previous approaches is a novel tensor-network algorithm that exploits the superfermion representation to simulate Lindblad dynamics directly in the energy eigenbasis of the baths (the so-called star geometry). This configuration is particularly favourable in fermionic systems, where only a limited energy window participates in the dynamics at finite temperature due to Pauli exclusion effects at low energies. Although we focus here on steady states of autonomous machines, our methods can be adapted to study transient dynamics or time-dependent Hamiltonians. Moreover, our tensor-network algorithm is inherently scalable to many-body problems, as we demonstrate by extracting the super-diffusive transport exponents of the isotropic Heisenberg model at high temperature. Our work thus paves the way for simulations of heat transport in strongly correlated systems that probe heretofore inaccessible regimes of temperature and system size.
In the remainder of the article, we build our methodology step by step. We begin with an introduction to autonomous thermal machines in Sec. II, where the problem to be solved is precisely defined. We then outline the mesoscopic-reservoir approach and demonstrate its connection to the infinite-bath scenario in Sec. III. Subsequently, in Sec. IV we detail the superfermion representation and use it to find an analytical expression for the non-equilibrium steady state of a non-interacting (quadratic) system. In Sec. V we explain how to compute particle and energy currents within our framework. Equipped with the exact solution for quadratic systems, in Sec. VI we study a non-interacting quantum-dot heat engine and compare the results with Landauer-Büttiker theory in order to identify the number and distribution of modes in the mesoscopic reservoirs needed to accurately reproduce the continuum limit. Next, in Sec. VII we detail our tensor-network algorithm for studying interacting problems. We then apply this algorithm in Sec. VIII to study a three-site interacting heat engine and a manybody Heisenberg spin model. Finally, we summarise and conclude in Sec. IX.

II. AUTONOMOUS QUANTUM THERMAL MACHINES
This work is concerned with autonomous thermal machines whose working medium is a quantum system S, which may be a complex entity comprising many interacting subsystems. The working medium is connected to multiple fermionic reservoirs labelled by the index α. These reservoirs are macroscopic systems described by equilibrium temperatures T α = 1/β α and chemical potentials µ α (we set k B = 1 = ). The total Hamiltonian of such a setup takes the form whereĤ S is the system Hamiltonian,Ĥ α is the Hamiltonian of bath α andĤ Sα describes its coupling to the system. We will consider exclusively HamiltoniansĤ tot that conserve fermion numberN =N S + αN α , wherê N S andN α are the total particle number operators for the system and each bath α, respectively. Crucially, the baths are taken to have an infinite volume and heat capacity, implying a diverging number of degrees of freedom, N → ∞. Moreover, it is typical to assume a factorised initial state of the form whereρ(0) is the initial system state andρ B = αρ α , withρ α = e −βα(Ĥα−µαNα) /Z α a thermal state and Z α the partition function of each reservoir. Evolving into the long-time limit the system S will generically relax to a steady state given bŷ where Tr B denotes the trace over all bath degrees of freedom. If the temperatures or chemical potentials of the reservoirs differ,ρ(∞) will be a non-equilibrium steady state (NESS) possessing currents of particles and energy. We focus especially on the simplest scenario depicted in Fig. 1, with two reservoirs labelled by α = L, R. The sustained fluxes of particles and energy in this setup can be exploited, for example by operating the device as an autonomous heat engine. In this case a temperature gradient, T L > T R , drives a current that performs work by FIG. 1. A simple thermal machine scenario in which the system S is coupled to two reservoirs L and R at temperatures TL > TR and possessing a chemical-potential difference µR − µL > 0. A particle J P and energy J E current is thus sustained through S. moving fermions against a chemical-potential difference V = µ R − µ L > 0. The power developed per unit time is given by where J P is the particle current, defined to be positive when flowing from left to right. The concomitant energy current J E (also from left to right) transfers heat out of the left lead and into the right lead at a rate [3] so that the first law of thermodynamics can be written as P =Q L −Q R . The second law of thermodynamics imposes the relation β RQR ≥ β LQL . The efficiency of heat-to-work conversion is thus given by where η C = 1 − T R /T L is the Carnot efficiency. Thus, the performance of an autonomous thermal machine depends on the currents and their relationship to the thermodynamic properties of the reservoirs. Evaluating the currents requires finding the NESS of the quantum system. In general, however, the computation of Eq. (3) is a difficult task. Analytical solutions are available only if the global Hamiltonian is noninteracting, while a direct numerical approximation with finite baths may require prohibitively large values of N in order to avoid Poincaré recurrences within the timescale of relaxation. On the other hand, perturbative schemes are limited to cases where either the internal interactions within S or its couplings to the reservoirs are weak. We thus take an alternative approach, in which the macroscopic reservoirs are replaced with mesoscopic leads comprising L sites, which are continuously damped towards thermal equilibrium by dissipative processes. As a consequence, convergence can be obtained with only moderate values of L, bringing the non-equilibrium thermodynamics of complex many-body quantum systems within reach.

III. FROM MACROSCOPIC RESERVOIRS TO MESOSCOPIC LEADS
In this section, we detail our approach to studying the problem described in Sec. II, where an infinite bath is replaced by a finite collection of damped modes. Here we outline the idea, leaving the mathematical details in Appendix A.
The system S is assumed to be a lattice of D sites, with arbitrary geometry and interactions, while the baths are modelled by infinite collections of non-interacting spinless fermionic modes. To illustrate the approach, we consider first the case of a single bath B, as shown in Fig. 2, described by the Hamiltonian whereb † m creates a fermion with energy ω m . Each site j of the system is described by a fermionic operatorĉ j . A particular site p of the system exchanges particles and energy with the bath via a tunnelling interaction where λ m is its coupling to bath mode m. The Heisenberg equation for the system operators reads as (9) Here, we have defined the noise operator and the memory kernel χ(t − t ) = {ξ(t),ξ † (t )} . The Gaussian statistics of the noise operator with respect to the initial product state Eq. (2) are defined by ξ (t) = 0 FIG. 2. The dynamics of a system coupled to a single thermal bath is determined by the bath's spectral density J (ω), with a bandwidth W , and the Fermi-Dirac distribution f (ω) corresponding to its chemical potential µ and temperature T . and where we have defined the spectral density as and introduced the Fermi-Dirac distribution f (ω) = (e β(ω−µ) + 1) −1 . The average system-bath coupling strength is typically quantified as where 2W denotes the reservoir bandwidth, namely the size of the energy range over which J (ω) has support [see Eq. (37), for example]. The state of S is completely determined by f (ω) and J (ω) via the noise statistics, since for an overall closed system the solution of Eq. (9) is sufficient to reconstruct all n-point correlation functions. Our approach is based on a key insight. Namely, that the open-system dynamics in Eq. (9), induced by an infinite bath with spectral function J (ω), can be accurately approximated by instead coupling the system to a finite collection of damped modes. Indeed, let us consider a lead of size L coupled to site p of the system, described by the Hamiltonian whereâ † k creates a fermion in the lead with energy ε k , and κ kp is the coupling strength. Each energy eigenmode k of the lead is coupled to an independent thermal bath modelled by an infinite non-interacting fermion reservoir B k , as illustrated in Fig. 3 (see Appendix A for details). These baths have identical temperatures and chemical potentials, but crucially they are characterised by a structureless frequency-independent spectral density J k (ω) = γ k , where γ k is a characteristic damping rate whose value may be different for each bath.
To analyse the steady-state physics it is sufficient to focus on long times, such that t γ −1 k , τ rel . Here τ rel represents the characteristic relaxation timescale of S due to its coupling with the bath [58]. In this limit, we find that the Heisenberg equations for the system variables in this configuration are identical to Eq. (9), but the statistics of the noise operator are now determined by an effective spectral density FIG. 3. (a) A Lorentzian spectral density J Lor (ω) is equivalent to coupling the system to a single auxiliary mode damped by a structureless reservoir. (b) A mesoscopic reservoir comprising many damped modes gives rise to an effective spectral density J eff (ω) that is a sum of Lorentzians. By tuning the damping of each mode and its coupling to the system J eff (ω) can approximate J (ω) of the infinite bath depicted in Fig. 2.
It follows that this damped mesoscopic lead configuration reproduces the correct steady state of S, so long as the true spectral density J (ω) can be well approximated by a sum of Lorentzians as above. This is depicted in Fig. 3.
In particular, Eq. (17) tends to Eq. (13) in the limit where L → ∞, γ k → 0, κ kp → λ k → 0 and ε k → ω k , i.e. an arbitrary continuous spectral density can be recovered in this limit. In particular, consider a given set of lead energies ε k that sample the spectral density and are arranged in ascending order, with energy spacing d k = ε k+1 − ε k . By taking κ kp = J (ε k )d k /2π and γ k ∝ d k , we obtain a controlled approximation of the bath spectral function as the lead size L increases and thus d k → 0.
In order to obtain a tractable description of the augmented system-lead configuration, we use the fact that both the damping rates γ k and the coupling constants κ kp are small in the large-L limit. Tracing out the baths, we derive a master equation describing the joint state of FIG. 4. In the limit L 1 modes in the lead each bath B k is sufficiently weakly coupled its corresponding lead mode that it can be accurately modelled by a Lindblad dissipator. The dissipator on a lead mode then that injects and ejects fermions at rates which in isolation damp the mode into a thermal state.
S and L, valid for times t γ −1 k , τ rel and up to second order in both the lead-bath and system-lead coupling (see Appendix A). We emphasise that the assumption that individual modes of the lead couple weakly to the system does not imply that the overall system-bath coupling Γ is weak. The quantum master equation is whereĤ =Ĥ S +Ĥ L +Ĥ SL denotes the Hamiltonian of the system and lead, while thermalisation of the lead is described by the Lindblad dissipator with f k = f (ε k ) denoting the sampling of the Fermi distribution by the lead modes. This master equation configuration is illustrated in Fig. 4. The above representation does not simplify the problem a priori, since it is strictly valid only in the large-L limit. However, a simplification may arise if the expectation values of operators converge with increasing L. We show numerically in later sections that this convergence occurs rapidly in several examples of interest for quantum thermodynamics. In such cases, a tractable number of lead sites L can be used to obtain a good approximation of an infinite bath with a continuous spectral density. For this, it is crucial that γ k remains the smallest energy scale in the physical configuration, to both model the spectral function correctly and accurately approximate the baths via the Lindblad equation [47,51].
So far we have considered a single bath coupled to a particular site of the system. However, the above results are easily generalised to describe the situation of several sites connected to multiple baths at different temperatures and chemical potentials. The steps of the above analysis are carried out independently for each bath, leading to additive contributions to the master equation.

IV. SUPERFERMION REPRESENTATION OF NON-EQUILIBRIUM DYNAMICS
In order to solve the dissipative dynamics under a master equation of the form in Eq. (18), we use the superfermion formalism introduced in Ref. [36]. For a non-interacting (quadratic) open system, this method provides numerically tractable analytical expressions for steady-state quantities. The superfermion representation is also central to our approach to simulating interacting systems, as discussed in Sec. VII. Here, we limit ourselves to a concise review of the formalism; for more details, see Appendix B.
The superfermion approach is akin to a purification or thermofield scheme for open systems. It doubles the system size by introducing a new fermionic ancilla mode for each of the modes present in the system and leads. To describe the formalism succinctly we stick for now to the single-lead setup of Eq. (18). In order to distinguish clearly between the ancillary modes and the physical modes of the system and lead, we introduce a unified notation for the latter. In this single-lead setup the total number of system and lead modes is M = D + L and so we define M fermion mode operatorŝ The ancillary modes are described by M additional canonical creation and annihilation operatorsŝ † k andŝ k . We use an interleaved ordering for the physical and ancillary operators, so that the Fock basis of the combined Hilbert space is defined by Here n are m are binary strings of length M that describe occupation numbers for the physical and ancillary modes, respectively. While the ordering used for the Fock basis is entirely arbitrary, we shall see shortly that interleaving has useful locality properties exploited later in Sec. VII. We now define a new (unnormalised) ket vector called the left vacuum as where the sum runs over all 2 M binary strings n. Using this ket, we can define a quantum state representing the system-lead density operator aŝ and the expectation values of any system or lead operator A as A key aspect of this formalism are the conjugation relations allowing physical creation (annihilation) operators to be swapped for ancillary annihilation (creation) operators. For the interleaved Fock ordering these conjugation relations are given bŷ Acting the master equation Eq. (18) on |I and using the conjugation relations yields a Schrödinger-type equation for the state, with the (non-Hermitian) generator of time evolution given byL whereĤ d⇔s is the same as the system-lead Hamiltonian H but with all physical operators replaced by their ancillary counterparts,d k →ŝ k . Crucially, dissipative processes are now described by non-Hermitian quadratic operators that, according to the interleaved mode ordering of Eq. (21), couple only nearest neighboursd k andŝ k . The formalism generalises straightforwardly to multiple leads by introducing an additional ancilla mode needed for each additional lead mode. So far the superfermion formalism is entirely general. In the special case where the system HamiltonianĤ S is non-interacting the formalism provides a compact expression for the exact solution of the NESS. In this case the system-lead Hamiltonian is quadratic with the form where H is an Hermitian M × M matrix. Next we define M × M diagonal matrices Γ + and Γ − containing the injection and ejection rates of fermions for each site.
Specifically, for the single-lead setup the first L follow the thermal damping rates contained in the dissipator Eq. (19), while the last D entries corresponding to the system modes are zero, giving Using these we define two additional diagonal matrices Λ = (Γ − + Γ + )/2 and Ω = (Γ − − Γ + )/2. Consequently, for the case of a non-interacting system the generatorL is quadratic with the form T is the full 2Mdimensional column vector of all physical and ancillary operators [59].
To determine the NESS we diagonaliseL by a similarity transformation, L = V V −1 , to find the complex eigenvalues = diag( 1 , . . . , 2M ) and the matrix of right eigenvectors V of L. As shown in Appendix B, the manybody NESS is a Fermi-sea-like state in which only modes with Im( µ ) > 0 are occupied, furnishing us with a complete solution of the problem. In particular, two-point correlation functions of physical modes in the NESS are found to be where D µν = δ µν Θ( Im{ µ }), with Θ(x) the Heaviside step function. This gives an efficient prescription to find steady state observables such as currents for noninteracting systems, while higher-order correlation functions follow from Wick's theorem.

V. NON-EQUILIBRIUM THERMODYNAMICS WITH MESOSCOPIC LEADS
The central focus of our work is autonomous thermal machines in the two-lead configuration illustrated in Fig. 5, with mesoscopic reservoirs labelled by α = L, R. These two leads of size L are described by Hamiltonians of the form Eq. (15) and Eq. (16), where the left lead couples to the first system site, p = 1, and the right lead to the last system site, p = D. Each lead is also acted on by a dissipator of the form given in Eq. (19). The master equation for this set-up thus reads as whereĤ =Ĥ S +Ĥ L +Ĥ R +Ĥ SL +Ĥ SR . To find expressions for the particle and energy currents, we need to consider the continuity equations for  Fig. 1 where some generic system S is coupled to two reservoirs with differing chemical potentials and temperatures.
the total particle-number operatorN =N S +N L +N R and total energy operatorĤ for the system and the leads.
where J P α and J E α are respectively the particle and energy currents flowing into the entire configuration via lead α, given by In the NESS, the time derivatives in Eqs. (32) vanish. Defining positive currents to flow across the system from left to right, we thus take J P = J P L = −J P R and similarly J E = J E L . Explicitly, we show in Appendix D that where the sum runs over only the modes of the left lead with f L (ε) = (e βL(ε−µL) + 1) −1 being its corresponding equilibrium distribution and f L, For sufficiently large systems with short-range interactions [60], it is possible to define current operatorsĴ P,E S supported only on S. In this case, we show in Appendix D that the expected values of these operators agree with the formulae given above, i.e. Ĵ P,E S = J P,E . Therefore, the currents determined from Eqs. (34) and (35) are expected to converge to the correct values in the continuum limit L → ∞, according to the discussion in Sec. III. However, in some cases, e.g. if S comprises just a single lattice site, no system operator for the currents can be defined, in which case the arguments of Sec. III do not strictly apply. Nevertheless, we show by example in Sec. VI that the currents computed from the mesoscopic-leads approach give an excellent approximation to the continuum limit in this case as well.

VI. NON-INTERACTING EXAMPLE: THE RESONANT-LEVEL HEAT ENGINE
In this section, we apply our methods to analyse the performance of an autonomous thermal machine with a non-interacting working medium. Since exact results are available here for the L → ∞ limit, this serves as a benchmark to evaluate the performance of the mesoscopicreservoir formalism which can also be solved numerically exactly using the superfermion formalism. Using this we estimate the number of lead modes needed to accurately reproduce the continuum limit of infinite baths. We take a single resonant level as our working medium, described by the HamiltonianĤ whereĉ † andĉ are the fermionic creation and annihilation operators in the system, respectively, and is the energy of the level. This models a single quantum dot in the spin-polarised regime running as a heat engine between two baths [61]. We note that a quantum-dot heat engine was recently realised experimentally [12].
In principle, our methods can handle structured spectral densities that are different for each bath. For simplicity, however, we take both reservoirs to be characterised by identical, flat spectral densities within a finite energy band, given by where Γ is the coupling strength between the system and the leads. In the continuum limit of macroscopic baths, the particle and energy currents for a non-interacting system can be computed from the Landauer-Büttiker (L-B) formulae where f α (ω) denotes the Fermi-Dirac distribution for lead α = L, R and τ (ω) is the transmission function. The latter is computed using the formalism described in Appendix C.
In the mesoscopic-reservoir approach, the spectral density is sampled by a finite number L of lead modes, as in Eq. (17). Taking the distribution of lead mode energies {ε k }, widths {γ k } and couplings {κ kp } to be identical for each lead, there remains significant freedom to choose these distributions in order to well approximate the continuum limit using moderate values of L.
System In particular, we use the logarithmic-linear discretisation scheme proposed in Refs. [54,62]. Here, L lin modes are placed in the energy window [−W * , W * ], with equally spaced frequencies, i.e. d k = ε k+1 − ε k = 2W * /L lin . Energies outside of this range are sampled by a smaller set of modes L log , with frequencies logarithmically spaced from W * (−W * ) to W (−W ), with energy intervals [ε n−1 , ε n ] = [±Λ −(n−1) , ±Λ −n ] for n = 1, · · · , L log and Λ −L log = W * . The dissipation rates are taken equal to these spacings, γ k = d k , while the coupling constants κ k are determined by the equation Γ = 2πκ 2 k /d k [36], in accordance with the considerations of Sec. III. For a given number of modes L = L log + L lin , this discretisation scheme gives better resolution within a smaller energy window [−W * , W * ] that includes the most relevant energy scales for the problem at hand. In our calculations, we henceforth set W = 8 and use this parameter as the overall energy scale, while W * = W/2. Moreover, we choose L log /L = 0.2.
Under these conditions, we show in Fig. 6 the behaviour of the particle current, where we have set equal temperatures in the leads T L = T R = W/8 but used different chemical potentials µ L = −µ R = W/16. In Fig. 6(a) we show the results for the particle current as a function of the system energy for different numbers of modes L in the leads and compare it with L-B theory. From both Fig. 6(a) and Fig. 6(b), it can be observed that a good agreement is obtained, the biggest difference observed as → 0, when the current reaches its maximum value. As expected, the agreement is improved with increasing L, although even moderate values of L ∼ O(10) approximately reproduce the continuum. Furthermore, in Fig. 6(c) we fix the energy of the level to study the behaviour with increasing L as a function of temperature T L = T R = T with system-lead coupling strength Γ fixed, and in Fig. 6(d) the behaviour with Γ for fixed T . For this particular choice of parameters we find the particle currents are robust to a wide range of T and Γ. Either low or high temperatures and weak or strong coupling yield similar results in both continuum or mesoscopic scenarios, even for a moderate number of modes in the mesoscopic leads.
In Fig. 7 we show the corresponding results for energy current. From Fig. 7(a) it can be observed that a better approximation is obtained when the number of modes in each lead is increased for a fixed set of parameters, with the absolute difference decreasing as a function of L, as can be concluded from Fig. 7(b). In Fig. 7(c) a key difference can be observed from the results obtained for particle current. The mesoscopic lead configuration is a good approximation as long as T is kept above a given threshold. This threshold is dictated by the smallest energy spacing in the leads d k and can be understood as follows. The effective spectral function of the mesoscopic leads is a sum of Lorentzian peaks, as in Eq. (17). When the temperature is smaller than the minimum energy spacing d k in the mesoscopic lead, these peaks are too far apart to properly resolve the variation of the Fermi-Dirac distribution. This significantly modifies the noise statistics given by Eq. (12) and the approximation is not reliable. It can be observed from Fig. 7(c) that the approximation at lower temperatures is much better for larger leads [63].
From Fig. 7(d) we observe that the approximation for energy current in the mesoscopic lead configuration is quite robust to a wide range of values of system-lead coupling Γ. The condition on the system-lead coupling is Γ γ k d min k [54], which is much less restrictive than the condition on temperature described above, and so a good approximation can be obtained for a wide set of Γ values.
Next we evaluate the power and efficiency given by Eqs. (4) and (6). In Fig. 8(a) we show the power output as a function of average chemical potential µ = (µ L + µ R )/2 and the potential difference V = µ R − µ L us- ing continuum leads, i.e., using Landauer-Büttiker predictions. In our calculations we set T L = 1.1W/8 and T R = W/8 and show the power output results only for the values of µ − and V for which the system acts as a power generator. It can be observed that the power output reaches a maximum value depending on bias and average chemical potential. In Fig. 8(b) we show the results for the same calculation, but instead we substitute the continuum leads with our mesoscopic lead configuration. The results are in good agreement up to the point where µ − reaches the boundary of linearly discretised and logarithmically discretised lead modes. Beyond W * and −W * , the spectral function is not sampled as finely and the power output results get distorted. We note that the window can be increased to resolve a bigger set of the parameter space, however, this would require more lead modes to resolve the maximum power output with the same accuracy. In Fig. 8(c) we show the maximum power output P max as a function of the system-lead coupling for both the L-B and mesoscopic lead predictions, which in turn reveals the value of Γ max for which P max reaches its maximum value. With our choice of parameters, Γ max lies very close in both configurations, as well as the overall behaviour as a function of system-lead coupling.
In Fig. 9(a) we show the efficiency obtained using continuum leads, normalised by the Carnot efficiency. It can be observed that the points of maximum efficiency lie close to the boundary where the system stops operating as an engine, i.e., where the potential difference becomes too large for the temperature gradient to drive electrons in the opposite direction of the bias. In Fig. 9(b) we present the results for the mesoscopic lead configuration. As before, we find that both predictions are quantitatively similar up to the point where the boundary of W * is reached. In Fig. 9(c) we show the efficiency at the point where the maximum power is obtained from the configuration as a function of Γ, where we observe that both the continuum and mesoscopic lead configurations predict very similar results. Not only is the strong system-lead coupling behaviour well-captured, but so is the Curzon-Ahlborn efficiency limit (approximately given by η C /2) at weak coupling [64].

VII. TENSOR NETWORK APPROACH
Having established that relatively modest sized mesoscopic leads can capture the continuum behaviour of a non-interacting system we now move on to consider the highly non-trivial problem of interacting systems. To do this we introduce in this section a tensor network based numerical method that can efficiently and accurately compute the interacting NESS of the the tworeservoir problem illustrated in Fig. 5. To describe the method we will return briefly to the single-lead configuration shown in Fig. 4 in which the first site p = 1 of the system S is coupled to the mesoscopic lead. Since we will exploit the superfermion formalism we continue to use the unified notation for modesd k given in Eq. (20).

A. Spin-1/2 representation
Our approach uses the matrix product state (MPS) decomposition that is a tensor network with a onedimensional chain-like geometry, as shown in Fig. 10(a). To apply this powerful ansatz to our setup we first map the lead and system modes into a one-dimensional chain. In doing so the coherent coupling between the lead modes and the system become long-ranged within this chain since they corresponding to a so-called "star geometry". Fundamentally this is because we use the energy eigenbasis of the lead.
Additionally, since MPS apply to systems built from a tensor product of local Hilbert spaces, to describe a spinless fermionic system requires that we transform it into a spin-1/2 representation. Our starting point is to introduce Fock states constructed from the unified physical modes with occupation-number vector n as which in the single-lead case has M = L + D and is ordered with lead modes first, as shown in Fig. 10(b). A spin-1/2 representation is then obtain via the well-known Jordan-Wigner (JW) transformation involving M spins [65,66] whereσ z q is the Pauli spin matrix in the z direction and σ ± q are the spin raising/lowering operators for the q-th spin. Correspondingly, the Fock states of Eq. (40) are equivalent to the spin states since each JW string vanishes on polarised spins it is applied to. Transforming the total HamiltonianĤ = H S +Ĥ L +Ĥ SL [from Eqs. (15) and (16)] to this representation giveŝ The star geometry, shown in Fig. 10(c), thus introduces JW strings to the lead-system coupling terms making them long-ranged multi-body spin operators. Similarly, the Lindblad dissipator of Eq. (19) becomes showing that the jump operators are now also non-local due to the JW strings.

B. Superfermion representation
By using the energy eigenbasis of the lead we have arrived at a master equation with a highly non-local multibody Hamiltonian and dissipator. The JW strings therefore appear to severely frustrate the use of MPS algorithms in this setup. Typically those arising from the star geometry of the Hamiltonian in Eq. (43) are dealt with by tridiagonalising the lead Hamiltonian, transforming it into a chain geometry and localising its coupling to the system. However, it is clear that this procedure profoundly complicates the dissipator in Eq. (44). The thermal damping of the lead induced by the dissipator is most naturally described in the lead's energy eigenbasis.
In the lead energy eigenbasis, the dissipator's JW strings can be eliminated by exploiting the superfermion representation of the open system introduced in Sec. IV. There, an interleaved physical and ancillary mode ordering was used, resulting in the dissipative processes becoming nearest-neighbour non-Hermitian Hamiltonian terms, as shown in Eq. (27). In this form, when moving to a spin-1/2 representation, the JW string of each system or lead site cancels with that of the corresponding ancillary mode, rendering the dissipator terms local.
To observe this explicitly, first note that the Fock basis of the combined Hilbert space of the physical and ancilla sites, namely Eq. (21), can be written in the spin-1/2 basis as The non-Hermitian generator of the superfermion time evolution thus becomeŝ showing that the dissipator contribution consists of onsite and nearest-neighbour terms.

C. Time evolving block decimation with swaps
To efficiently simulate the time evolution of the correlated system described by Eq. (46), we use one of the most well-known algorithms within the tensor network family, namely, the time-evolving block decimation (TEBD) [67,68]. Given some system governed by a HamiltonianĤ loc = iĥ i,i+1 , comprising a sum of 2site termsĥ i,i+1 along a chain of length M , the standard formulation of TEBD computes the MPS approximation of the propagation |ψ(t) = exp(−iĤ loc t) |ψ(0) . This is done by first breaking up the evolution into many small time-steps δt and then performing a second-order Trotter expansion as whereÛ i,i+1 = exp(− i 2ĥ i,i+1 δt). In this way, a time step of propagation is implemented by a staircase circuit FIG. 11. The sweeping sequence of two-site gatesÛ k,1 between the k-th lead mode and the first system site along with the fermionic SWAPsŜ f needed to implement a Trotter step for the star geometry couplings shown in Fig. 10(c).
of two-site gates sweeping right-to-left and then left-toright. Each two-site gate can be applied to the MPS and, via a singular value decomposition, the result can be re-factorised and truncated back into MPS form.
Here, we use a simple modification of TEBD that allows us to compute the time-evolution under fermionic star-geometry HamiltoniansĤ star = iĥ i,M , where all sites i < M interact with the last site M . The key ingredient is the fermionic SWAP gateŜ f , which is a conventional SWAP gate between spins j and j + 1 that exchanges their spin configurations, but also incorporates the application of the localσ z j operator from the JW string of Eq. (43). For two sites, the gate is given bŷ where the negative sign accounts for the anticommutation relation between two fermionic creation operators when both sites j and j+1 are occupied. By interspersing fermionic SWAP gates within the Trotter expansion, as shown in Fig. 11, distant sites are temporarily made adjacent, allowing the standard update nearest-neighbour two-gate gate update to be applied. Time-evolution under a long-ranged Hamiltonian is generally considered impractical for tensor network calculations, due to very fast growth of entanglement across the system. This conjecture has been challenged in recent studies of fermionic impurity models, where efficient tensor network calculations have been performed using a star-geometry setup [69,70]. The proliferation of correlations in these models is curtailed by Pauli blocking within the majority of the modes of the lead, limiting them to the range of modes around the Fermi energy. This favourable situation persists in the mesoscopic thermal lead setup considered here.

D. Non-equilibrium steady state solver
The TEBD algorithm works equally well for non-Hermitian Hamiltonians generating non-unitary propagation. Indeed, it has been widely used to study the NESS of incoherently driven chains where the coupling to the reservoirs is localised to one [71][72][73][74][75][76][77][78][79][80] or two sites [81][82][83][84] at the boundaries. We have now introduced all the elements required to extend the capabilities of TEBD to simulate the open-system governed by the Hamiltonian Eq. (43) and the dissipator Eq. (44).
First, we move to the superfermion representation where the generatorL is given by Eq. (46). We define dimer sites composed of a physical (system or lead) site and its corresponding ancilla, as shown in Fig. 12(a). This procedure squares the dimension of the local basis. The left vacuum state |I in this representation is a product state of dimers, with each dimer initialised in FIG. 12. (a) Ancilla modes are interleaved with the system and lead modes they are associated to. Computationally the system or lead site and its ancilla are bundled together as a dimer site. (b) A two-dimer site gateÛ dim,k is applied between the k-th lead mode dimer, and the first system site dimer. This is followed by four fermionic SWAPs to shuffle the system site and its ancilla through the lead and its ancilla making the next lead mode adjacent. This is repeated all the way along the chain and back to complete one time-step. the |↑↑ + |↓↓ state.
Next, we identify all the terms inL that correlate the dimers locate at lead site k and system site p = 1. Assuming these sites are adjacent to each other through SWAP operations, we expresŝ We identify spin 1 as the k-th lead eigenmode with spin 2 being its corresponding ancilla mode. On the other hand, spin 3 is the system site coupled to the lead with spin 4 as the corresponding ancilla mode. A JW string appears between interacting spins that are not adjacent, however, they remain local to the dimer pair. The exponential of this operator,Û dim,k = exp(−iL dim,k δt/2), defines a non-unitary gate for a half time step δt. This operator accounts for all the coherent interactions and the non-Hermitian terms, describing the dissipation between the lead mode and the system site. We have assumed a Hamiltonian of the form Eq. (36) in Eq. (49). Finally, the non-unitary gatesÛ dim,k are then applied along with fermionic SWAP gates that shuffle the system dimer along the chain, as shown in Fig. 12(b). Altogether, this computes the action of the propagator exp(−iLδt) and formally solves Eq. (26) for a single timestep. We take the initial state to be |ρ(0) = |I , and find the steady state |ρ(∞) by evolving towards the long-time limit. Expectation values and the trace of the density operator follow from the inner product with |I as given in Eq. (24).
The same simulation scheme can be readily extended to the two-lead configuration, as shown in Fig. 13(a), with the long-time limit now giving rise to a NESS. The approach to stationarity is assessed by evaluating the convergence of observables such as the particle and energy currents. In practice, we used a dynamicallyincreasing MPS bond dimension χ for different time-step parameters δt. We chose an initial χ and δt and evolved the system up to an intermediate time. The resulting state was then further evolved in time with a larger χ and an appropriately reduced δt. This procedure is repeated until the currents obtained converged up to a small tolerance of 1 − 2%. The largest bond dimension used in our calculations was χ max = 220, showing that a moderate computational effort was required to access the NESS. All MPS calculations in this work were performed using the open-source Tensor Network Theory (TNT) library [85,86].

VIII. INTERACTING EXAMPLES
In this section, we employ the tensor network algorithm from Sec. VII to study an autonomous thermal machine with an interacting working medium, as depicted in Fig. 13(b). Our methods enable us to consider the challenging problem of simultaneously strong interactions and system-bath coupling, far beyond the linearresponse regime.

A. Interacting three-site engine
Our first example is an autonomous quantum heat engine with a three-site interacting working medium, which is described by the Hamiltonian wheren j =ĉ † jĉ j is the density operator for site j and U is the interaction strength. The last term in the equation above corresponds to a density-density interaction of neighbouring particles. A small central system composed of D = 3 interacting fermionic sites can be interpreted as a three-site version of the interacting resonant level model [87]. We set the system hopping t S = W/8 and focus on the regime in which the temperature gradient and the difference in chemical potential between the mesoscopic reservoirs is strong. We set T L = 10t S , T R = t S , µ L = −t S /2, µ R = t S /2 and j = = t S . With these parameters, the system operates as a heat engine, i.e. particle current flows from the left reservoir to the right reservoir, driven by the temperature gradient against a chemical potential gradient. As in Sec. VI, both leads are assumed to have identical, flat spectral densities given by Eq. (37) and we use the logarithmic-linear discretisation scheme with W * = W/2 and L log /L = 0.2. This choice of parameters is a representative example useful in exposing the efficacy of the proposed methodology.
We first focus on the dependence of the currents on the system-lead coupling Γ, as shown in Fig. 14. In Fig. 14(a), the energy current for a particular value of the interaction strength U = 1.2t S is shown as a function of Γ. Remarkably, a density-density interaction yields a larger energy current flowing through the system compared to the non-interacting case in the regime described by our choice of parameters. The same observation holds for the particle current in Fig. 14(b), since for our choice of parameters the particle current and the power output are equivalent [see Eq. (4)]. The efficiency shown in Fig. 14(c), remains approximately constant as a function of system-lead coupling strength just like the noninteracting case. Future work will investigate a larger range of parameters to identify a maximum power output for a given interaction strength.
The insets in Figs. 14(a) and 14(b) show the error associated to employing a finite number of modes in each reservoir for a specific value of Γ = 3t S , where the currents in the interacting case reach the maximum value. The error is computed from an extrapolated value of the currents to the L → ∞ limit, based on the currents for finite L, for each respective case. It can be observed that for the specific choice of parameters, a good approximation can be obtained to a few percent accuracy using L = 50, compared to larger reservoirs. The energy current converges faster than the particle current (power) in this case. This behaviour is expected, as observing Figs. 16 and 17 for the non-interacting case in Appendix E, the largest deviation for the particle current occurs where the maximum value is obtained, while the largest deviation for the energy current is observed near the edges of the band.

B. High-temperature transport
The transport properties of spin chains has been studied extensively using standard open-system MPS approaches based on a boundary driving Lindblad master equation. This has been successful in accurately describing the high-temperature spin/particle transport behaviour of the integrable anisotropic XXZ Heisenberg model [72][73][74]88] as well as non-integrable versions of the model when integrability-breaking perturbations are introduced, such as magnetic impurities [80] or disorder [78,79,83,84]. However, driving on the boundary spins is formally equivalent to infinite temperature baths. Modelling energy currents therefore requires more elaborate multi-site boundary driving to mimic finite temperature differences. While this approach has proven successful for the very high temperature limit, its reliability as the temperature is lowered is questionable. The mesoscopic leads construction introduced here overcomes this deficiency.
The system Hamiltonian introduced in Eq. (50) is the spinless fermion equivalent of the anisotropic XXZ Heisenberg model. This model exhibits a range of distinct linear response particle and energy transport behaviour as a function of the anisotropy U . Specifically, these include ballistic transport which is characterised by a constant value of the current as a function of system size D, as well as diffusive transport where J P ∝ 1/D ν with ν = 1 [80]. Anomalous diffusion is signaled by 0 < ν < 1 and ν > 1, corresponding to superdiffusion and subdiffusion, respectively. A sharp transition in the system's transport properties is known to occur at the isotropic point U/t S = 2, with the system displaying ballistic transport for U/t S < 2, while for U/t S > 2 transport becomes diffusive. Furthermore, precisely at the isotropic point U/t S = 2, boundary driving calcu- lations have shown that transport is superdiffusive with ν = 1/2 [74]. These results are expected to hold only in the linear-response regime at high temperatures, where the structure of the thermal baths becomes irrelevant. We now corroborate these results using our mesoscopic reservoir formalism.
As before, we choose the same discretisation scheme and bath structure parameters. We focus on the isotropic point U/t S = 2 and set j /t S = /t S = 1. We set the temperature on each reservoir to a high value of T L = T R = 1000t S and choose a small chemical potential gradient µ L = −µ R = 0.025t S , where we expect the system to be in linear response regime. In Fig. 15 we show both the particle and energy currents as a function of system size D. We have used L = 20 modes for both left and right reservoirs. As can be observed, the currents fit a power law scaling with an exponent very close to ν = 1/2 in clear indication of super-diffusive behaviour. We remark that at high temperature, fewer reservoir modes can be used to obtain the correct transport exponent, as observed from boundary driving calculations [74]. We believe that our method is primed to accurately explore the finite temperature transport properties of strongly correlated many-body systems.

IX. CONCLUSIONS AND OUTLOOK
In this work we introduced a novel methodology to simulate the heat and particle currents in thermal machines which comprise a complex working medium coupled to fermionic leads at fixed temperatures and chem-ical potentials. The method is based on the concept of mesoscopic reservoirs whose energy modes are damped in order the replicate the continuum. The method allows for calculations in highly non-equilibrium scenarios such as strong system-lead coupling and large biases. In order to cope with non-quadratic interactions in the working medium, we implemented a novel tensor network algorithm directly in the star geometry using auxiliary modes.
For the purpose of expounding the method, in this paper we considered only autonomous thermal machines where the working medium is time independent. In order to benchmark our technique we first focussed on replicating the steady-state thermodynamics of the resonantlevel heat engine. The simplicity of this quadratic model allows for direct comparison with the Landauer-Büttiker theory for quantum transport. We observed excellent agreement across a wide parameter regime. We then explored efficiency and power in a strongly interacting three-qubit machine in a parameter regime where other methods are known to struggle. In doing this we observed that, remarkably, the efficiency is enhanced as a function of the system-lead coupling in the presence of non-quadratic interactions. Furthermore, we demonstrated that our technique is suitable to perform highly non-trivial heat and particle transport calculations in strongly correlated many-body systems by performing a scaling analysis at the isotropic point of the paradigmatic Heisenberg model.
Due to the flexibility of our technique we expect that the method is extendable further in the direction of steady-state thermodynamics of complex interacting quantum systems. Beyond strong coupling and far-fromequilibrium scenarios, we expect that our technique may find useful applications in the study of time-dependent working media, bulk noise effects and non-trivial spectral densities; thus taking quantum thermodynamics to unexplored horizons.
Note added in proof. During the preparation of the manuscript we became aware of a recent, related work that uses a different tensor-network algorithm to study transport with mesoscopic reservoirs [89].
(A4) which is equivalent to Eq. (9). Here, the noise operator isξ(t) = −i m e −iωmt λ mbm (0) and the memory kernel is The solution of Eq. (A4) at time t depends in principle on the entire past history of the noise operatorξ(s) for s < t. Once found, the solution forĉ j (t) is sufficient to reconstruct all n-point correlation functions of S, which together uniquely specify the quantum state (amongst other information). Since the initial bath state is Gaussian, these correlation functions depend on the noise only via its two-time correlations

Mesoscopic-lead configuration
We now turn to the mesoscopic-reservoir configuration, with total HamiltonianĤ =Ĥ S +Ĥ SL +Ĥ L +Ĥ LB +Ĥ B . HereĤ L andĤ SL describe the lead and its coupling to the system and are given explicitly by Eqs. (15) and (16).
Each mode of the lead is further coupled to an infinite reservoir according tô whereâ k describes mode k of the lead, while the ladder operatorsb kq describe the bath connected to mode k. Each bath is described by the flat spectral density We are interested in the evolution of the joint system-lead state ρ(t) starting from the initial product state Eq. (2), where all baths are initialised at the same temperature and chemical potential. As in Eq (A3), we formally solve the Heisenberg equation of motion for the bath variables to find (A10) Substituting this into the equation of motion forâ k (t), we obtain d dtâ Here, we defined the noise operatorŝ and the memory kernels χ k (t−t ) = {ξ k (t),ξ † k (t )} . For the flat spectral density in Eq. (A9), the noise correlations are given by Next we formally solve Eq. (A11) to find Considering long times, such that γ k t 1, the first term above is negligible and will be ignored in the following. Substituting this solution into the equations of motion for the system variables, we finally obtain an effective quantum Langevin equation This is of the same form as Eq. (A4), but with the noise operator and the memory kernel where the effective spectral density J eff (ω) is the sum of Lorentzians in Eq. (17). The second equality above follows via an identity which can be proved by contour integration: It remains to check the effective noise correlations. We have, using Eqs. (A13), (A14) and (A19), where we have neglected all terms proportional to e −γ k s or e −γ k s . This approximation is valid at long times, so long as the solution of Eq. (A16) depends only on the past history ofξ eff (s) within a time window that is essentially finite. This will generically be the case for any system that relaxes to a steady state when coupled to a bath, since any memory of environmental fluctuations in the far past is eventually lost. In particular, if τ rel is the (slowest) characteristic timescale of relaxation of S, then we need consider only arguments ofξ eff (s) in the range t − τ rel s < t. Hence, the approximations leading to Eqs. (A20) and (A21) are valid for all times such that If this holds, we have shown that the effective noise generated by the mesoscopic lead is equivalent to an infinite bath with a spectral density given by Eq. (17), giving rise to an identical equation of motion for the system, Eq. (A16).

Quantum master equation
Finally, we briefly discuss the derivation of the quantum master equation. In the limit of large lead size, L → ∞, both the lead-bath couplings κ kp and the system-lead coupling γ k must tend to zero in order to recover the continuum spectral density J (ω) (see the discussion below Eq. (17)). In this limit, we may derive a quantum master equation using perturbation theory up to second order in the lead-bath coupling. Following the standard procedure [15], and working in an interaction picture with respect to the free Hamiltonian (A23) In the interaction picture, the free evolution of the lead operators is given bŷ In this Appendix we give further details the superfermion [36] steady state solution of the master equation in Eq. (18) for a non-interacting system of size D coupled a single mesoscopic lead of size L.
This open system has a quadratic generatorL = f † Lf − η defined by the 2M × 2M non-Hermitian matrix L where M = D + L. To compute its NESS we proceed to diagonalise this matrix as L = V ε V −1 to giving a diagonal matrix ε of complex eigenvalues ε µ . These eigenvalues come in conjugate pairs and we shall denote the half with Im{ µ } > 0 as set Ξ + and the other half with Im{ µ } < 0 as Ξ − .
We identify the corresponding normal mode operators asξ † =f † V andχ = V −1f . Althoughχ µ andξ µ mix physicald k and ancillary modesŝ k via a similarity transformation, and so are not Hermitian conjugates of one another, they still obey canonical anticonmmutation relations [90], e.g.
The equations of motion for the normal mode operators follow from the commutator withL giving so in vector form the time-evolved mode operators arê A defining property of the NESS isL |ρ(∞) = 0. Using this we compute the time-evolution of the NESS when acted upon by a normal mode operator to obtain and also For these time-evolved states not to diverge in time we require thatξ † µ |ρ(∞) = 0 when µ ∈ Ξ + andχ ν |ρ(∞) = 0 when ν ∈ Ξ − . This pair of constraints is analogous to those of a Fermi sea state |FS whereĉ † j |FS = 0 when mode j is occupied, andĉ j |FS = 0 when it is empty. Similarly for the left vacuum state |I we get implying the complementary constraints I|ξ † µ = 0 when µ ∈ Ξ − and I|χ ν = 0 when ν ∈ Ξ + . Together these relations fully define the 2M × 2M matrix D of normal mode two-point correlations of the NESS with elements We immediately see that D µν = 0 whenever µ ∈ Ξ − and/or ν ∈ Ξ − . The case µ, ν ∈ Ξ + is then determined using Eq. (B1) to find that D µν = δ µν . Hence in general we have indicating that the set Ξ + of normal modes are the unit filled Fermi sea of the NESS. Using this result we can evaluate physical quantities such as the single-particle Green function G ij (t, t ) = ĉ † i (t)ĉ j (t ) = I|ĉ † i (t)ĉ j (t )|ρ(∞) for the system S. Transforming back from the normal modes we havê and thus the Green function follows as where we have used that D is diagonal and the indices i, j = (L + 1), . . . , M give the physical system S modes.
This reduces to the NESS expectation value in Eq. (30) once t = t = 0. The Fermi sea structure of the NESS allows Wick's theorem to be applied to breakup expectation values for high-order correlations into two-point ones, for example leaving products of terms that can be readily evaluated using the NESS normal mode constraints determined above. When the central system is a single-level withĤ S from Eq. (36), the transmission function can be proven to be of Lorentzian form and equivalent to while a central system composed of D fermionic sites witĥ H S from Eq. (E1) has a transmission function which corresponds to a convolution of Lorentzian functions whose form depends on the site energies and hopping amplitudes t S , as observed from Eq. (C5). With the previous expressions for τ (ε), Eqs. (38) and (39) can then be evaluated numerically to obtain particle and energy currents for a given system.

Appendix D: Definitions of currents
In this Appendix, we discuss the energy and particle current in more detail. In the mesoscopic-lead configuration, the currents are found from the continuity equation and given by Eq. (33) Since this superoperator acts only on the lead degrees of freedom, we find the explicit expressions quoted in Eqs. (34) and (35) with straightforward algebra. In sufficiently large central systems, an alternative definition of the currents can be derived from the continuity equations within the system itself. Let us focus on 1D systems with two-body interactions coupled to two baths at the first and final sites j = 1, D, as considered in the examples of Secs. VI and VIII. In this case, the fermion number and Hamiltonian can be written aŝ wheren j =ĉ † jĉ j is the local fermion density on site j andĥ j,j+1 denotes a local energy density operator. Sincê h j,j+1 has support only on sites j and j +1, we derive the continuity equation for number density from the Heisenberg equation forn j : where we defined the particle current operator which clearly depends only on system variables. In the steady state, the time derivatives of all expectation values vanish and we find that the current is homogeneous, i.e. Ĵ P j−1→j = Ĵ P j→j+1 . Eq. (D3) holds only for j = 1, D.
Here we used the fact that [Ĥ SL ,N L +n 1 ] = 0, which merely reflects the overall conservation of fermion number and the fact that L couples only to site j = 1. Combining Eqs. (D5) and (D6) and assuming steady-state conditions we deduce that Therefore, so long as the system comprises D ≥ 2 sites, the current computed via Eq. (34) coincides with the expectation value of a system operator. For the energy current, one similarly finds in the bulk of the system Now, considering the Heisenberg equations for bothĤ SL andĤ L and assuming steady-state conditions, we conclude that Therefore, the energy current computed from Eq. (35) also coincides with the expected value of a system operator, so long as D ≥ 3. The above arguments, although developed for the specific case of two-body interactions in one dimension, are based only on conservation laws and the locality of interactions, which are general principles. Similar arguments can thus be developed for more general n-body interacting systems in higher-dimensional geometries, so long as a sufficiently large region of the central system is not directly connected to the baths. In panels (c) and (d) we fix every parameter and study the particle current as a function of temperature and system-lead coupling, respectively. FIG. 17. Energy current from L-B and mesoscopic reservoir predictions flowing from the left lead and into the system (a) as a function of the on-site energy (same parameter for every site) for a central system with L = 100 sites and a fixed number of modes in the leads N = 50 and (b) as a function of the hopping amplitude tS (same parameter for every site). In panels (c) and (d) we fix every parameter and study the energy current as a function of temperature and system-lead coupling, respectively.
The Landauer-Büttiker calculations are done by evaluating Eq. (38) using the transmission function obtained as described in Appendix C. It can be observed that for a fixed number of modes in the leads L and a fixed number of sites in the central system D the approximation to the continuum limit using mesoscopic reservoirs is robust to a wide range of on-site energies. The small oscillations that can be observed near the band edges at | | |W * | are due to the logarithmic spacing of modes. Furthermore, from Fig. 16(b), the same can be said when is fixed and t S is changed to different values. Given that the energies in the central system are bounded by −2t S and 2t S , the oscillations due to logarithmic discretisation are observed close to t S ≈ W/2. The same observations hold for energy current in Figs. 17(a) and 17(b) As a function of temperature, a similar behaviour as for the single-level system can be observed. In particular, for particle current and energy current in Figs. 16(c) and 17(c), respectively, the continuum is properly approximated with the exception of the values of temperature that are lower than the minimum energy spacing of the modes in the leads. For these small temperatures, the Fermi-Dirac distributions of the leads resemble a Heaviside step function and the discontinuity can no longer be well-captured by discrete and broadened energy modes. Following from our previous discussion for the singlelevel system, to obtain a better approximation at lower temperatures one can either increase the number of total energy modes or decrease the width of the window [−W * , W * ]. The former choice comes with the cost of a larger computational complexity, while with the latter one can then only provide a good approximation of the continuum for a smaller range in the parameter space of , t S , µ L and µ R . If these values are fixed, a good choice of [−W * , W * ] can be used to obtain better approximations at lower temperatures with its limit, as discussed for the single-level system, related to the minimum value of d k in the linearly-discretised region. As a function of the system-lead coupling, the results are very robust to a wide range of values as observed from Figs. 16(d) and 17(d). Due to the ballistic (coherent) nature of transport in the central system, currents become independent of D in the asymptotic regime.