Fermionic reaction coordinates and their application to an autonomous Maxwell demon in the strong coupling regime

We establish a theoretical method which goes beyond the weak coupling and Markovian approximations while remaining intuitive, using a quantum master equation in a larger Hilbert space. The method is applicable to all impurity Hamiltonians tunnel-coupled to one (or multiple) baths of free fermions. The accuracy of the method is in principle not limited by the system-bath coupling strength, but rather by the shape of the spectral density and it is especially suited to study situations far away from the wide-band limit. In analogy to the bosonic case, we call it the fermionic reaction coordinate mapping. As an application we consider a thermoelectric device made of two Coulomb-coupled quantum dots. We pay particular attention to the regime where this device operates as an autonomous Maxwell demon shoveling electrons against the voltage bias thanks to information. Contrary to previous studies we do not rely on a Markovian weak coupling description. Our numerical findings reveal that in the regime of strong coupling and non-Markovianity, the Maxwell demon is often doomed to disappear except in a narrow parameter regime of small power output.


I. INTRODUCTION
Many problems in quantum transport are modeled by an impurity Hamiltonian H imp linearly coupled to a bath of free fermions described via Here, c ( †) k annihilates (creates) a fermion of energy k in the bath, which is tunnel coupled with complex amplitude t k to the system via a fermionic annihilation (creation) operator d ( †) of the system. The actual physical system under study is described by H imp , which could be one (or multiple) quantum dots, molecules in mechanically controllable break junctions or nanoelectromechanical systems, among other possible impurity systems.
To treat such systems theoretically, various approximate or formally exact techniques have been developed, such as quantum master equations (MEs) [1][2][3], the formalism of nonequilibrium Green's functions [4], or renormalization group techniques [5]. Whereas MEs easily allow to treat interactions in the impurity even under nonequilibrium situations, their use is limited to the weak-coupling, Markovian, and high-temperature ("sequential tunneling") regimes. Green's functions can overcome the latter problem, but have difficulties to treat interacting impurities (for instance, due to Coulomb forces). Whereas this problem can be tackled by using numerical renormalization group approaches, they in turn are hard to apply far away from equilibrium.
In the first half of this paper (Sec. II), we will introduce a technique which lies "in-between" these approaches. It allows to some extent to overcome the limitations of the standard perturbative approach commonly applied to obtain a ME while still retaining a description in terms of a ME such that interactions and nonequilibrium situations can be conveniently treated. The price to pay is an enlarged Hilbert space, in which a suitable redefined impurity HamiltonianH imp includes particularly chosen dominant degrees of freedom from the bath, which we will call fermionic reaction coordinates (RCs). In the context of linear bosonic reservoirs (Caldeira-Leggett or Brownian motion models), this technique has a longer tradition [6]. It has found various applications in the theory of open quantum systems [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and it is also closely related to the "time-evolving density matrix using orthogonal polynomials algorithm" (TEDOPA) [23][24][25][26][27][28]. We remark that, although it shares many similarities with the bosonic case, the RC mapping was not studied for fermionic reservoirs before. In addition, our paper adds additional insights to the recent attempts to find a meaningful thermodynamic description beyond the weakcoupling and Markovian regimes [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46]. In particular, our work shows that techniques based on a redefined system-bath partition using RC mappings [34,39,42,46] turn out to be useful for fermionic reservoirs, too.
In the second half of the paper (Secs. III and IV), we make use of our method to study two (spinless) Coulomb-coupled quantum dots in contact with three heat reservoirs. This setup is raising increasing attention within the context of quantum thermodynamics, as it provides a prototypical example of a thermoelectric device transporting electrons against a potential bias due to an energetic flow from a hot to a cold bath [47,48].
It is well studied in the weak-coupling and Markovian regimes, theoretically [49,50] as well as experimentally [51,52]. Moreover, Ref. [50] identified necessary conditions which guarantee that this device can be interpreted as an autonomous Maxwell demon (MD), i.e., a device which is capable of extracting work (in this case by charging a battery) due to a clever way of processing information. From this perspective, the device was further studied theoretically [53,54] and experimentally [55].
Unfortunately, the study in Ref. [50] has revealed that a proper operation of the device as a MD requires a strong coupling to a cold reservoir and structured (i.e., non-Markovian) spectral densities for the hot reservoirs. These three requirements challenge the usual range of validity of a ME description, but the question whether the transparent interpretation of a MD also holds under these conditions has not been answered yet. By using the method of fermionic RCs, we will indeed show that the device can still be interpreted as a MD, albeit in a very narrow parameter regime restricted to low-power output. Our results are also in qualitative agreement with a recent study [56] where cotunneling effects were taken into account. Another complementary paper studies the impact of strong-coupling effects for nonautonomous, i.e., measurement-based, feedback loops on the thermodynamic performance of a MD [57].

II. FERMIONIC REACTION COORDINATES
In this section, we develop the theory of fermionic reaction coordinates (RCs), which is useful to explore the physics of open systems beyond the weak-coupling, Markovian, and high-temperature assumptions. Technically speaking, this mapping is a unitary transformation applied to the bath Hamiltonian, which allows to separate out a particular degree of freedom in the bath, called the RC. The details of this transformation are reported in Sec. II A.
The hope is then that the original impurity together with the RC is only weakly coupled to a Markovian residual bath such that it is possible to apply standard master-equation (ME) procedures to this redefined system-bath partition. Whether this is possible is case dependent and will be discussed in Sec. II B.
To make the paper self-contained, we also give a short review in Appendix A of the ME approach including the definition of energy and particle currents, which we will use in the second part of this paper. Furthermore, in Appendix B we benchmark the fermionic RC method against the exact solution of the Fano-Anderson model (also known as single-electron transistor).

A. Mapping
The mapping will work whenever it is allowed to describe the interaction between an arbitrary impurity Hamiltonian H imp and the fermionic reservoir as in Eq. (1). An important quantity in the study of such open systems is the spectral density (SD) (also called hybridization function) of the bath, which is defined as It contains the complete information about the way in which the bath is coupled to the system and will be of central importance in the following.
For mathematical rigor one often demands that J (ω) is strictly greater than zero for ω ∈ [ω L ,ω R ] and zero outside this interval where ω L < ω R ∈ R are referred to as cutoff frequencies [17]. However, as long as all quantities converge, our equations also remain valid if the SD decays only exponentially or polynomially. Convergence problems for infinite cutoff frequencies arise only when the mapping is applied iteratively (see below). Furthermore, gapped SDs (having support only at disconnected intervals) can be treated by applying this mapping to each sub-SD separately.
For later reference, we start by considering the Heisenberg equation of motion for d and c k (suppressing the time dependence on all operators and settingh ≡ 1 throughout): We Fourier transform them according to the definitionf (z) ≡ ∞ −∞ dt e izt f (t) [with Im(z) > 0]. This yields (again dropping the explicit z dependence) to After some algebra we obtain a formally exact expression for d, which reads as Here, we introduced the Cauchy transform By the Sokhotski-Plemelj theorem this fulfills for ω ∈ R where P denotes the Cauchy principal value. We now perform the mapping by introducing a new set of fermionic creation and annihilation operators {C k } via C k = l kl c l where is a unitary matrix fulfilling † = 1 such that the fermionic anticommutation relations are preserved. In particular, we fix the first row of by the requirement where the parameter λ 0 is fixed by the requirement {C 1 ,C † 1 } = 1, which implies In analogy with the bosonic case, we will call C 1 a fermionic RC or collective coordinate. Furthermore, we demand for l = 1 and m = 1 that k k lk new fermions, the Hamiltonian (1) becomes with T k = m m t m km /λ 0 and The complex phase of the constant λ 0 is not fixed by this procedure, but it can be fully absorbed by redefining the operator C 1 , adding only an additional phase to the renormalized couplings T k . We will therefore consider λ 0 positive in our numerical investigations. The effect of the residual modes which we now finally determine as a functional of the original SD J (ω). For this purpose, we look again at the Heisenberg equations, but now within the transformed coordinates. After Fourier transformation, we obtain analogously to Eqs. (5) and (6) Again, we can formally solve for C k and C 1 to obtain an exact expression ford: Since both Eqs. (7) and (17) are exact, they must coincide. Thus, by comparison we obtain and due to relation (9), we get This finally completes the mapping: performing a specific normal-mode transformation on the Hamiltonian (1) yields a new Hamiltonian (12). The new parameters λ 0 , E 1 and the new SD J 1 (ω) are given by Eqs. (11), (13), and (19), respectively. They are all completely specified in terms of the initial SD J (ω). We remark that this mapping is formally exact; no approximation has been made. The mapping is summarized in Fig. 1.
Before proceeding, we remark that the impurity Hamiltonian H imp is allowed to be explicitly time dependent without invalidating any point in the derivation above as the explicit form of H imp never entered the derivation. In fact, even a global time dependence of the tunnel Hamiltonian H I of the form α(t)H I does not invalidate any point of our derivation above as we can include the time dependence α(t) in the definition of the operator d (t) ≡ α(t)d. The reason FIG. 1. Pictorial representation of the RC mapping (top) and the equations determining the transformation for the general case (middle column) and the example of a Lorentzian SD as used throughout the text (right column). The shaded gray area in the sketch indicates which part is treated as the impurity, whereas the remaining part is treated as the bath in spirit of the ME approach outlined in Appendix A. We used different colors for the free fermions in the original/residual bath to emphasize that they are not the same before and after the mapping.
for this generality comes from the fact that the mapping is a unitary transformation in the bath Hilbert space alone and does not touch the system Hilbert space.
Finally, we investigate what happens if we apply the mapping iteratively. In fact, if we define a new impurity Hamiltoniañ a new tunnel couplingH I = k (T * k C † 1 C k + T k C † k C 1 ) and a new residual reservoirH R = k E k C † k C k , we see that the Hamiltonian (12) has the same structure as (1). Thus, applying the mapping iteratively we obtain a chain of RCs with coupling constants λ 0 ,λ 1 , . . . , energies E 1 ,E 2 , . . . , and residual SDs J 1 (ω),J 2 (ω), . . . . These can be determined recursively in full analogy where we have additionally inserted Eq. (9). In analogy to the bosonic case [17], it is therefore natural to ask what is the limiting SDJ (ω) obtained after n → ∞ iterations. Assuming that the limit exists and denoting we obtain from Eq. (18) the condition with the two possible solutionsW ± (z) =Ē − z ± (Ē − z) 2 − 4|λ| 2 . The SD for (Ē − ω) 2 4|λ| 2 therefore becomesJ where the requirement of positivity fixes a unique solution of the square root. Thus, the limiting SD describes a semicircle with radius 2|λ| centered around ω = E. In fact,Ē and |λ| are fixed by the initial cutoff frequencies viaĒ = (ω L + ω R )/2 and 2|λ| = ω R − ω L . Indeed, if this were not the case, one would obtain a contradiction. This follows from Eq. (19) together with our initial assumption that the SD is strictly greater for ω ∈ (ω L ,ω R ) and zero outside this interval. We remark that our reasoning here does not allow to draw any conclusion about the behavior of convergence. In particular, it could happen that the sequence of SDs does not converge (e.g., for infinite cutoff frequencies). The necessary conditions for convergence were studied by Woods et al. [25].

B. Advantages of the RC mapping
The previous section described how to apply a unitary transformation to the bath such that it can be mapped to a new, redefined impurity coupled to a residual bath. It is different from the conventional bosonic mapping [17]. There, the final Hamiltonian is described in terms of position and momentum operators and has a different quadratic form, which, rewritten with creation and annihilation operators, displays counterrotating terms which are absent in the fermionic case. See also Ref. [25] for further details on this point.
We emphasize that this mapping is formally exact, it does not touch the system Hilbert space and, thus, it can be applied to arbitrary (even time-dependent) impurity Hamiltonians and also to tunnel Hamiltonians H I with a global time dependence. Furthermore, if the impurity is coupled to multiple reservoirs, the mapping can be applied to each reservoir separately such that it can be easily applied to the study of nonequilibrium scenarios (see Sec. IV).
In principle, the hope is that after the mapping, the problem is easier to tackle by using any preferred theoretical method. We here follow the idea to use a Markovian quantum ME for the extended system with redefined impurity HamiltonianH imp (as shaded in gray in Fig. 1) [19,22,34,39,42,46]. To get a feeling for this approach, let us consider an initial Lorentzian SD of the form where describes the overall coupling strength and the width of the Lorentzian centered around the resonance frequency ω 0 . We will indeed use this SD for our applications below. By applying Eqs. (11), (13), and (19), we obtain that the system couples to the RC with coupling strength λ 0 = √ /2 and energy E 1 = ω 0 , which is in turn coupled to a residual bath described by a SD of the form J 1 (ω) = 2 . Thus, we see that the residual SD is completely flat, which is commonly believed to describe Markovian behavior. In order to justify a ME approach, the redefined impurity should be additionally weakly coupled to the residual bath, which is exactly the case if the initial width of the Lorentzian is small. Therefore, our method should be especially suited for the study of very structured SDs (e.g., described by sharp peaks) opposite to the wide-band limit. The equivalence for the particular example of a single quantum dot coupled to a bath with Lorentzian SD and a double quantum dot with flat SD was already noticed before [58], but we emphasize that our mapping is in general valid for any impurity system and SD.
Even if the residual bath is not strictly Markovian or weakly coupled to the redefined system, one might hope that the ME including the RC provides nevertheless a convenient way to improve the accuracy of the results (compared to a conventional ME approach) as it takes a larger part of the model exactly into account. This is also supported by our benchmark in Appendix B where the coupling strength to the residual bath is rather moderate than weak. Moreover, if one is more interested in weak coupling than Markovianity, it is also possible to split the support of the SD into multiple intervals and to apply the RC mapping to each interval separately. This gives rise to more RCs, but also to a weaker coupling to the residual baths. Furthermore, it should be emphasized that the validity of any ME approach is likely to break down at very low temperatures. This, however, does not indicate a failure of the RC method itself.
A benefit of the present approach is that it straightforwardly allows for a consistent thermodynamic interpretation [34,39,42,46]. In fact, due to the ME approach, we have direct access to the internal energy and entropy of the redefined impurity as well as to the energy and matter currents I E and I M from the residual reservoir, which are defined in Eqs. (A13) and (A14). If the system is coupled to multiple reservoirs ν at inverse temperature β ν and with chemical potential μ ν , it is straightforward to establish the validity of the nonequilibrium first law (energy balance) and second law (positivity of entropy production rate˙ ) of thermodynamics: Here,Ẽ imp (t) = tr{H impρimp (t)} is the internal energy and S imp (t) = −tr{ρ imp (t) lnρ imp (t)} is the entropy of the redefined impurity and the heat flows are given byQ M . In case of an explicit time dependence,Ẇ mech (t) = tr{[d tHimp (t)]ρ imp (t)} denotes the mechanical work done on the system. Note that a clean derivation of the positivity of the entropy production rate, Eq. (29), requires a ME in Lindblad form although we have numerically observed no violation of positivity for our ME derived in Appendix A.
The strategy to reformulate the laws of thermodynamics for an extended system incorporating a part of the previous bath has been suggested in Refs. [34,39,42,46]. It has the advantage that it avoids the difficulties faced with other methods such as Green's functions [30][31][32]38,45] or hierarchy of equations of motions [36,37], where the interaction energy cannot be unambiguously assigned to the system or the bath when the system is not at steady state.
We remark that the RC method allows to treat a larger class of initial conditions for the redefined impurity. To compare it with the conventional ME method [1][2][3], one has to choose the initial state of the RC to be equilibrated and decorrelated from the impurity state. In comparison with Green's function techniques [30][31][32]38,45], one usually assumes initially a global equilibrium state instead. Finally, we remark that by applying counting field techniques to the residual baths, nonequilibrium fluctuation relations such as those derived in Ref. [59] also hold for our approach.
To close this section, let us compare our method with the TEDOPA algorithm [23][24][25][26][27][28], which can be straightforwardly extended to fermionic reservoirs [24,25]. The goal of this method is to provide an exact mapping of the whole reservoir onto a semi-infinite chain of coupled fermions by using the theory of orthogonal polynomials. Then, one usually solves the whole system exactly by using density matrix renormalization group (DMRG) methods [60]. In principle, our method also allows by iterative application to obtain a semi-infinite chain of coupled fermions, but we believe that the strength of our method lies in the possibility to apply it step by step and to treat the transformed system by an approximate, yet simple and intuitive ME approach.

III. AUTONOMOUS ELECTRONIC MAXWELL DEMON: A REVIEW
We here collect some recent results about autonomous Maxwell demons (MDs), especially within the context of electronic transport. We start by giving a simplified argument in Sec. III A to underline why it is important to extend the study of autonomous MDs beyond the weak-coupling regime. Section III B gives an overview over recent activities in the field with special attention paid to the electronic context. In Sec. III C we then consider a particular important model, which we will numerically study using the theory of fermionic RCs in Sec. IV. For comparison, Sec. III D finally reviews the thermodynamics of this model in the weak-coupling regime and shows under which conditions the device can be viewed as an autonomous MD. The last two subsections also fix most of the notation and parameters eventually used in Sec. IV to explore the non-Markovian regime.

A. A general argument
For readers unfamiliar with the working mechanism of an autonomous Maxwell demon, we here give a simplified description of it, which is in direct analogy with the device studied later on in this section.
Let us start by specifying what we mean by an autonomous MD. We consider bipartite systems where one part can be understood as the controlled system and the other part assumes the role of the detector and controller by the physical interaction with the controlled system. The complete device should be autonomous in the sense that there is no time dependence in the global Hamiltonian. Also, external interventions by means of measurement and feedback control are forbidden, very similar to the idea to use "coherent quantum control" to simulate measurement-based quantum control [61,62]. Thus, all physically relevant parts of the device are explicitly modeled. The device should also be a useful thermodynamic machine allowing, for instance, the extraction of work. These desiderata are also fulfilled by many other engines, but beyond that we especially require the following: (i) The way information is processed in the device should be particularly transparent and these "informational degrees of freedom" are called the demon part of the device (whereas the rest is simply called the system).
(ii) The energetics associated to the demon part should be negligible (though not necessarily strictly zero) compared to the energetics of the system itself, i.e., the device should be information dominated.
For the moment, let us denote by E S ( E D ) and τ S (τ D ) some characteristic energy scale and some typical relaxation time scale of the system (demon). Furthermore, we assume that the system (demon) has access to some heat reservoir at temperature T S (T D ). Whenever convenient, we will also use the notation of rates γ S,D = τ −1 S,D and inverse temperature β S,D = (k B T S,D ) −1 in our description. We emphasize that the relaxation rates γ S,D are within the conventional weak-coupling approach proportional to the coupling strength between the system/demon and their respective reservoir; compare also with Appendix A.
We start our analysis with the observation that an important feature of any MD is to extract work from fluctuations. Thus, in order to have non-negligible fluctuations in the system, we demand that Let us assume that the device delivers useful work and let us approximate the output power viȧ In the nonautonomous version of MD, the second law for feedback control predicts that the maximum amount of work in the case of an ideal classical feedback controller is bounded by [63,64] where I S:D denotes the mutual information between the system and the demon, Here, ρ SD denotes the density operator of the system and demon, ρ S,D = tr D,S {ρ SD } the marginal state, and S(ρ) ≡ −tr{ρ ln ρ} the von Neumann entropy. Roughly speaking, I S:D quantifies the amount of correlations shared between the system and the demon and it plays an essential role in the field of information thermodynamics [63]. Based on this insight, the goal of an autonomous MD will also be to establish a strong correlation between the system and the demon part in order to harness it via some intrinsic feedback loop; thus, the demon part has to act like a detector. This implies that the demon must be able to react quickly enough to changes in the system, which naturally leads to the requirement i.e., the demons typical timescale τ D is much smaller than τ S . This requirement inevitably tells us that in order to enhance the power output of the device by increasing γ S , this necessarily implies a corresponding increase in γ D , which in turn implies that the demon must be more strongly coupled to its own bath at temperature T D . However, requirement Ia is not enough to guarantee that the demon adapts with high probability to the correct state: given a certain state of the system, there should be a unique and stable state of the demon. Expressed differently, whereas we required the system to fluctuate relatively strongly, the demon should not fluctuate more than necessary. Therefore, for any system state the requirement of a reliable or precise demon translates into Requirements Ia and Ib are linked to our desideratum (i) mentioned initially. Point (ii) is fulfilled by the requirement From requirement Ib and II and condition (30) it also follows immediately that Hence, the heat reservoir of the demon must necessarily be much colder than the heat reservoir of the system. To conclude, the above argumentation shows us that an autonomous MD needs to be carefully tuned. The fact that the demon acts like a detector requires it to be fast and precise and, thus, to be much more strongly coupled to its own bath than the system. On top of that, a small energy consumption requires the demon to have access to a low-entropy (or lowtemperature) reservoir. Both requirements, strong coupling and low temperature, challenge the usual range of validity of commonly employed perturbative approaches. The rough estimates given here should be also compared with the analysis in Sec. III D.

B. Maxwell's demon in the electronic context
The central idea of an electronic MD is to find some feedback mechanism, which shovels particles against a chemical gradient instead of a thermal gradient as in the traditional thought experiment of Maxwell. The big advantage of this setup comes from recent technological advances, which allow to measure and manipulate single electrons in quantum dots (see, e.g., [51,52,55,[65][66][67][68][69]). Direct measurement and manipulation of individual phonons, which are predominantly responsible for thermal transport, is much harder instead. Thus, there have been many theoretical studies how to use charge fluctuation in mesoscopic conductors to extract work via some feedback loop, either in an autonomous way [50,54,[70][71][72][73][74] or not [57,64,[75][76][77][78][79][80][81].
In the following, we will focus on the autonomous MD introduced in Ref. [50]. Globally, it resembles a thermoelectric device based on two Coulomb-coupled quantum dots [47,49]. However, by a careful fine tuning of the parameters, it is possible to show that it resembles the phenomenological electronic MD introduced in Ref. [75] and, thus, it has a particularly transparent interpretation in terms of information flows [50,53,54,78]. Moreover, very similar experimental realizations were reported [51,52,55,69] although it remains unclear whether it is possible to reach the ideal informationdominated regime.

C. Model
The model is schematically depicted in Fig. 2. The Hamiltonian of the impurity (system and demon) reads as where d † s/d (d s/d ) is a fermionic creation (annihilation) operator for the system/demon dot. The Hamiltonian describes two dots with onsite energies s and d and Coulomb interaction U , which we will treat exactly with our method. Note that there are no electrons tunneling between the dots.
The noninteracting reservoirs ν ∈ {L,R,D} are modeled as free fermions as in Eq. (1) with Hamiltonian H (ν) R = k kν c † kν c kν . The reservoirs L and R are tunnel coupled to the system quantum dot via the Hamiltonian whereas the reservoir D is tunnel coupled to the demon dot via The total Hamiltonian then reads as H tot = H imp + ν∈{L,R,D} (H (ν) I + H (ν) R ). Within the weak-coupling approach, each reservoir is assumed to be well described by an equilibrium distribution according to a temperature T ν and chemical potential μ ν . Below, we will set for simplicity To complete the picture, we also sketched the Lorentzian SDs of each reservoir whose product with the Fermi function determines the transition rates (see Appendix A). The trajectory starts with a filled dot of the demon but an empty dot of the system (1). If an electron wants to enter the system, it needs an energy s + U and, although the Fermi factor is smaller for the right reservoir, this is overcompensated by the imbalance of the Lorentzian tunneling rates (2). After an electron has entered the system, it pushes the electron in the demon dot to a higher energy above the chemical potential such that it instantaneously jumps out of the dot due to requirement (44) (3). Finally, and again due to the imbalance of the Lorentizian tunneling rates, the electron in the system dot leaves to the left bath and the demon dot gets refilled (4). Overall, one electron was transferred against the bias by harnessing the correlation between the system and the demon dot and by transferring an energy amount U from the hot to the cold reservoir.
T L = T R ≡ T and μ L = s + V /2 and μ R = s − V /2 where V = μ L − μ R 0 denotes the bias voltage across the system.
The SD J (ν) (ω) of each reservoir ν is modeled by a Lorentzian as in Eq. (27) with coupling strength ν , peak width ν , and resonance frequency ω 0ν . In the following, we will set L = R ≡ S and L = R ≡ S , i.e., the system is coupled with equal strength to the left and to the right reservoir and only the position of the peaks will be different.

D. Thermodynamics in the idealized limit
In the idealized picture, the device is weakly coupled to three Markovian thermal reservoirs as depicted in Fig. 2. Within this approach, the dynamics of the system is governed by a master equation (ME) with rates obtained from Fermi's golden rule (see Ref. [50] for the full rate ME). Altogether, the system works as a thermoelectric device [49,50]. If we set T D < T , the device is able to convert an energy current I E ≡ −I (D) E > 0 flowing into the reservoir D into an electric current I M ≡ I (L) M < 0 flowing against the bias V provided that the SDs J (ν) (ω) and other parameters are chosen appropriately [for a definition of energy and matter currents, see Eqs. (A13) and (A14)]. Because the device is out of equilibrium, it produces entropy at a ratė It is also possible to give compact analytical expressions for I M and I E in terms of the steady state [49,50], but for our purposes it suffices to note the following proportionality relation: This means that the power output of the device is directly proportional to the coupling strength S in the weak-coupling regime and, hence, it is desirable to choose S as large as possible. Furthermore, we add that one can associate a second law to the local dynamics of the system only, which reads as [53] Here,İ S is the flow of mutual information between the system and the demon (for a precise definition, see Ref. [53]). Equation (43) shows that the ability to shovel electrons against the bias, βV I M < 0, can be influenced by energetic (βI E ) as well as entropic (İ S ) contributions. More importantly for our discussion, there is a clean limit in which the system can be viewed as an autonomous MD and in which the reduced dynamics of the system quantum dot coincides with the ideal, nonautonomous MD studied in Ref. [75]. In this regime, the stochastic trajectory shown in Fig. 3 becomes likely and transport against the bias possible. We will shortly review here this limit for completeness; more details can be found elsewhere [50,53]. The three different limits we need to take are the following (compare also with Sec. III A): 1a Fast demon. First, the demon needs to be relatively fast in order to adapt quickly enough to the system state such that it is guaranteed that the correct feedback loop is applied. This is ensured by demanding that D S , (44) implying that the demon is much more strongly coupled to its reservoir than the system [cf. Eq. (34)]. Within the formal limit D / S → ∞, the demon dot can be adiabatically eliminated and a closed effective description for the system dot alone emerges. In the strict limit, the flow of mutual information becomes exactly identical toİ S = −β U I E so that the effective second law coincides with the true second law˙ S =˙ .
1b Precise demon. In order to use the demon as a detector, it should also be precise in the sense that it is as correlated as possible with the system since this increases the mutual information [compare with Eq.
Then, in the strict limit of D / S → ∞ and β D U → ∞, the demon dot is occupied if and only if the system dot is empty and vice versa implying a perfect (anti)correlation between the system and the demon. It should be noted, however, that this strict limit also implies a diverging entropy production rate (41). Additional details concerning the use of such setups as detection devices can be found in Refs. [82,83].
2 Maxwell demon. Within the fast and precise demon limit, the demon can be reliably interpreted as a detector, yet it still disturbs the energetic balances of the system at the order of U . Not very surprisingly, the MD limit consists of demanding that which is in agreement with Eqs. (30) and (36) if we note that the energy current through the system is roughly proportional to S I M in the limit of negligible U . Thus, in this limit the difference in energy transferred from the right to the left reservoir (which is given by the flow of energy in the detector bath I E ) becomes immeasurably small compared to the energetic current S I M . In this regime, the term βI E becomes negligible and the second law reads aṡ i.e., the ability to shovel electrons against the bias is purely entropy (or information) dominated. The above three limits are, however, not sufficient to shovel electrons against the bias. In order to achieve this, we also need to break the left-right symmetry of the device even in the absence of any bias (V = 0). This is most easily done by choosing different SDs J (L) (ω) and J (R) (ω) fulfilling J (L) ( s ) > J (L) ( s + U ) and J (R) ( s ) < J (R) ( s + U ) as depicted in Fig. 3. More specifically, we choose for our Lorentzian SDs The imbalance in the SDs is then quantified by which was previously [50] described by the parameter e δ . Furthermore, we choose to fix the following parameters for all upcoming numerical calculations: We set the onsite energies equal s = d = , the inverse temperature to β = 1, and we choose a small, but finite positive bias V = 0.01 . In addition, we set U = 0.015 such that it is small compared to the system energy [see Eq.

IV. ELECTRONIC MAXWELL DEMON BEYOND WEAK COUPLING
In this section we report our main findings which are based on numerical comparison of the ideal (i.e., weakly coupled, Markovian, and high temperature) model (see Sec. III D), with the extended RC approach able to go beyond these limitations. 1 We will focus on the steady-state behavior only and compare three important quantities.
First, the electric current I M flowing through the system dot from the left to the right reservoir. It is directly proportional to (minus) the work output of the device and is therefore an important quantifier for the overall thermodynamic performance.
Second, we look at the relative energetic imbalance E denotes the energy flow into the demon reservoir and I (L) E the energy flow from reservoir L [for a definition of energy and particle currents, see Eqs. (A13) and (A14)]. If this quantity is large, the device works in an energy-dominated regime, whereas it goes to zero in the Maxwell demon limit (46).
Third, we will evaluate the mutual information I S:D between the system and the demon dot. It reveals key insights about the question whether the working mechanism of the device can be seen as some implicit measurement and feedback loop, which needs large correlations. Note that for our setup 0 I S:D 2 ln 2, but the maximum value is attained only for a pure and maximally entangled state. The classical limit is instead I cl S:D ln 2 such that values relatively close to this limit signify already strong correlations in our noisy setup. In our numerical studies we have found that it suffices to restrict the plots of the mutual information to the interval [0, ln 2], although this does not a priori imply that there are no quantum correlations present.

A. Extended models
The validity of the standard ME approach as sketched in Appendix A is limited by the coupling strength to the reservoirs, the degree of non-Markovianity, and the temperature of the reservoirs. In principle, all three limitations are challenged by the electronic MD model from Sec. III. We therefore compare it with three different extended models: (1) The working mechanism of our device requires a breaking of the left-right symmetry, which was achieved by choosing different SDs J (L) (ω) and J (R) (ω) peaked around s and s + U , respectively. In the ideal case, where the peaks are very narrow, the SDs act like perfect electron filters increasing the thermoelectric performance. Quite problematically, however, a strongly peaked SD is usually associated with strong non-Markovianity [3,84]. In order to capture the non-Markovian behavior of the left and right reservoirs, we introduce fermionic RCs C ( †) l and C ( †) r for the left and the right reservoirs in what we call "model 1." For the SDs centered at ω 0,L = s and ω 0,R = s + U , the resulting redefined impurity Hamiltonian 1 Numerical results were obtained using a Mathematica notebook, which is available upon request from the corresponding author. To compute marginal states, as needed for Eq. (33), we used the algorithm "Partial trace of a multiqubit system" by M. S. Tame. extended model (1) extended model (2) extended model (3) FIG . 4. Sketch of the geometry of the extended models studied in the text. Each circle represents a single fermionic site ("quantum dot"). The tunnel coupling is indicated by straight lines and the Coulomb interaction by a capacitor. The three reservoirs as well as the energy and electric current are also indicated as in Fig. 2. reads as (50) and the coupling to the residual left and right reservoirs is described by a flat SD J (ν) 1 (ω) = 2 S , ν ∈ {L,R}. A sketch of the setup is shown in Fig. 4. We remark that the inclusion of additional quantum dots to model peaked SDs acting as energy filters has been used elsewhere, too [49,85,86].
(2) Putting the problem of non-Markovianity aside, the most pressing problem is that the demon dot must be relatively strongly coupled to a low-temperature heat reservoir. This problem becomes more pronounced if we want to enhance the power output by increasing S . Model (2) therefore introduces a fermionic RC C ( †) d for the demon bath in order to study those effects. When we use ω 0,D = d + U/2, the resulting redefined impurity Hamiltonian reads as and the coupling to the residual demon reservoir is described by a flat SD J (D) 1 (ω) = 2 D . A sketch of the setup is again shown in Fig. 4.
(3) In order to capture both problems, model 3 finally uses a fermionic RC for all reservoirs (see Fig. 4). Within our framework, this will provide the ultimate test for the simplified model from Sec. III D. The redefined Hamiltonian is given by the sum of Eqs. (50) and (51) without counting H imp twice, The SDs of the residual reservoirs L, R, and D are given by 2 S , 2 S , and 2 D , respectively.

B. Numerical results
In what follows, we will compare the results of our extended models with the original MD treatment exposed in Sec. III D. We will refer to the latter as double quantum dot Maxwell demon (DQDMD) treatment and use red, dashed lines in the plots. The regime where the simple double quantum dot treatment and the MD interpretation are valid are shaded in gray in the plots.

RCs for the system reservoirs
We start with the extended model (1) to answer the question how narrowly peaked can we choose the SDs J (ν) (ω) before the Markovian description from Sec. III D breaks down? The answer is shown in Fig. 5. It clearly shows that there is an optimum value for the sharpness of the peak, which can be found by maximizing the power output −V I M while still retaining a valid Markovian description. For S → 0 the RC description indeed predicts I M → 0 because the coupling to the residual baths is directly proportional to S . It is worth to emphasize that the naive ME treatment of Sec. III D predicts a monotonically increasing power output in the limit S → 0 in strong contrast to the actual behavior. We conclude that one has to be very careful if one chooses strongly peaked SDs (sometimes called "spectral filters") to increase the power output of a thermodynamic device.
In the opposite limit, S → ∞, the engine also breaks down because the SDs become essentially flat and the left-right symmetry remains unbroken. Note that the amount of mutual information between the system and the demon is almost unaffected by the shape of the SD as expected.
Although the simple DQDMD description becomes inaccurate below a certain critical threshold * S / s ≈ 0.01, the plot also shows that a Markovian description is valid even if the SD is not completely flat but slightly structured [at the critical value, the imbalance (49) is for the chosen numerical parameters roughly 1.18]. Finally, we note that the regime, where the Markovian description is valid while at the same time the energetic imbalance is very small (less than 5%), is restricted to a very narrow window around * S / s demonstrating the careful fine tuning needed to ensure a working autonomous MD.

RC for the demon reservoir
Numerical results associated to model (2) concerning the question what happens to the demon in the strong-coupling and low-temperature regimes are shown in Figs. 6 and 7.
First, Fig. 6 shows our three relevant quantities as a function of S , which (in the weak-coupling regime) is directly proportional to the power output of the device [Eq. with the choice S = * S = 0.01 s . As expected, the plot of the electric current I M demonstrates that the DQDMD treatment is only valid for very small coupling strength S < 10 −4 s (shaded gray region), beyond that the demon fails to work. Also, the plot of the relative energy imbalance I E /I (L) E shows that we are leaving the MD regime of negligible energy consumption for larger S because more transport channels are opening up. More interestingly, however, is the plot of the mutual information I S:D , which reveals the physical reason why the demon fails. As explained in Sec. III, the working mechanism is based on a strong correlation between the demon and the system due to the Coulomb interaction and a careful tuning of the demon's reservoir. But, for stronger coupling D the electron in the demon dot gets more and more correlated with its reservoir than with the system, i.e., it becomes more and more delocalized. The only way to counterbalance this behavior is by increasing the Coulomb interaction U , but this increases the energy imbalance pushing us away from the MD regime. This is the ultimate reason why the MD is limited to the weak-coupling situation and, thus, due to the required timescale separation D S , to very low power output. To complement the analysis, Fig. 7 shows the same quantities for varying inverse temperature β D for relatively strong-coupling strength S = 10 −4 s and, shown as insets, for S = 10 −5 s , which was used in Fig. 5. As expected, for larger S we see stronger deviations from the ideal values confirming our previous observation. In addition, Fig. 7 demonstrates two more important features. First, the device works better for lower temperatures because the formation of correlations is hindered if the demon is subjected to more thermal noise.
Second, this figure also shows that a DQDMD description is usually only valid at high temperatures, whereas more sophisticated methods are needed for lower temperatures. We stress once more that, although the RC method allows to treat lower temperatures, its validity also does not extend down to zero temperature T D → 0.

RCs for all reservoirs
Finally, numerical results using model (3) are shown in Fig. 8. It is our most accurate analysis and combines models (1) and (2) and the plots demonstrate that the results, which we have drawn above in a separate analysis, hold true also in the complete picture. Furthermore, although the simple DQDMD picture from Sec. III D fails for most parameters as expected, it also agrees well if we pay careful attention to its range of validity. Thus, the analysis of the ideal MD [50], despite the various limits involved, remains qualitatively and quantitatively true for a narrow parameter window.

V. SUMMARY
In this paper, we have developed the theory of fermionic RCs, which provides a tool to extend the range of validity of the usual ME, especially for very structured, i.e., strongly non-Markovian, SDs. The benefit of our approach is that the ME approach still allows to treat interactions in the system exactly, it can be straightforwardly applied to nonequilibrium situations and has a transparent thermodynamic interpretation. One drawback of the method is that not every initial SD can be mapped to an effectively weakly coupled and Markovian situation such that the application of the ME has to be justified on a case-by-case study. Another drawback comes from the use of a ME itself, which becomes invalid at very low temperatures. Nevertheless, we believe that the fermionic RC mapping has the potential to find widespread application in quantum transport as a simple and transparent tool to treat structured SDs.
As a particular application we then considered a thermoelectric device, which can be interpreted as an autonomous MD for a specific range of parameters. Previous analyses in the field of information thermodynamics were done in the idealized weakcoupling and Markovian regimes, an exception being Ref. [56]. As we have argued in Sec. III A and as Ref. [50] explicitly shows, MD lives in a parameter regime where the use of these idealized assumptions becomes increasingly questionable. We here addressed systematically the question of what happens to the performance of an autonomous MD if we relax the weakcoupling and Markovian assumptions (also see Ref. [57] for the study of a nonautonomous MD in the strong-coupling regime).
Our numerical results clearly convey two messages: First, we have proven that it is indeed possible to reach the idealized MD regime, even if one uses a more sophisticated method, which is able to take into account strong-coupling and non-Markovian effects. Thus, MD is not a mere hypothetical being, but can be found in actual physical systems. On the other hand, our paper also indicates that the possible parameter regime of an ideal MD is very narrow and necessarily limited to lowpower output. It therefore still remains a challenge to find out to what extent MD will play a role in actual, practically useful devices, where already the use of a model Hamiltonian of the form (38) amounts to a strong assumption as it neglects, e.g., spin and vibrational degrees of freedom as well as multiple charges on a single quantum dot. Agreement No. 681456). G.S. and P.S. furthermore acknowledge support of the WE-Heraeus foundation for supporting the 640 WE-Heraeus seminar.

APPENDIX A: DERIVATION OF THE MASTER EQUATION
In this appendix we briefly state the essential steps to derive the quantum ME used in the main text. We start by focusing on an impurity (also often called the "system") coupled to a single reservoir with Hamiltonian H tot = H imp + H I + H R . Because within the weak-coupling approach there is no direct influence between the different reservoirs, we can simply add the contributions of them at the end. More detailed derivations of MEs can be found elsewhere [1][2][3]59]. Note that H imp , H I , and H R are arbitrary and left unspecified here. For the MEs used in Sec. IV, one would indeed need to replace H imp byH imp .
Our starting point is the second-order Liouville-von Neumann equation in the interaction picture after performing the Markovian approximation whereÃ(t) ≡ e i(H imp +H R )t Ae −i(H imp +H R )t denotes operators in the interaction picture and R 0 = e −β(H R −μN R ) /Z describes the equilibrium density operator of the reservoir with N R being the particle-number operator of the reservoir. For our purposes, it suffices to consider an interaction Hamiltonian of the form where d is an arbitrary fermionic system operator. We denote the eigensystem and transition frequencies of H imp as such that we can write the system coupling operator in the interaction picture as Then, after introducing the SD J (ω) = 2π k |t k | 2 δ(ω − k ) of the bath and after moving out of the interaction picture, we obtain terms like and their Hermitian conjugate. Here, To evaluate the integral over τ , we then use and neglect the imaginary principal value part in the following (cf. Appendix B). After this step, the full ME can be written as where we introduced the operators The ME (A8) is ready for numerical implementation. Note that we did not perform the commonly employed secular approximation, which guarantees a Lindblad form of the generator [1][2][3]59], but often predicts unphysical results especially for more complex systems (see also Appendix B). If we include coupling to multiple baths characterized by a distinct chemical potential or temperature, we can simply add up the contribution of each bath separately to the ME. The result is then with the dissipator The energy and particle currents into bath ν are then given by where N imp is the particle-number operator of the impurity.

APPENDIX B: EXAMPLE AND BENCHMARK: SINGLE-ELECTRON TRANSISTOR
To benchmark our approach, we consider the possibly simplest fermionic transport setup usually called a singleelectron transistor (SET): a spinless quantum dot with onsite energy coupled to two fermionic baths ν ∈ {L,R}. In this case, we have H imp = d † d and we model the contact to the baths again by a Lorentzian of the form (27), which we assume to be the same for the left and right reservoirs. Figures 9 and 10 now compare the energy and particle current through the SET using three different methods. The first (dotted green) is based on a standard rate equation for the probability to find the quantum dot empty or filled (i.e., the ME from Appendix A applied to H imp = d † d). The second method makes use of a single RC mapping individually applied to each reservoir and yields the Hamiltoniañ which now consists of three serially coupled quantum dots. The second method then treats the residual baths as weakly coupled and Markovian by using the ME treatment from Appendix A. In addition, it allows for different perturbative treatments by either including (solid dark-blue line) or neglecting (dashed light-blue line) Lamb shift terms or by performing an additional secular approximation on top (dashed-dotted black line). Finally, the model also admits an exact solution for the matter and energy current which reads as (solid red line; compare with, e.g., Ref. [87])  Here, = L + R denotes the Lamb shift As we can see, the RC method based on the ME (A8) (i.e., without secular approximation) gives an excellent agreement with the exact result for a wide parameter regime. Generally, we see that over all coupling strengths, including or neglecting the Lamb shift has little effect. For intermediate-coupling strengths, including the Lamb shift terms, it is even reducing the agreement with the exact solution, whereas for the ultrastrong-coupling regime it is improving the agreement, as is particularly visible in the energy current. More importantly, however, we see that the often employed secular approximation fails completely in the weak to intermediate parameter regimes. These reasons justify the use of the ME we derived in Appendix A.
It is nevertheless important to remark that the RC method is also limited. First of all, in order to justify using a weakcoupling ME forH imp , the width of the Lorentzian must be small enough because it is directly proportional to the coupling strength of the RC with the residual bath. Second, even for small our approach is still based on the use of a ME invalidating our results for very low temperature. For instance, the differential conductance of the SET lim V →0 dI M dV computed with the RC method completely fails to reproduce the exact results for T → 0 (not shown here for brevity).