A pseudo-fermion method for the exact description of fermionic environments: from single-molecule electronics to Kondo resonance

We develop a discrete fermion approach for modelling the strong interaction of an arbitrary system interacting with continuum electronic reservoirs. The approach is based on a pseudo-fermion decomposition of the continuum bath correlation functions, and is only limited by the accuracy of this decomposition. We show that to obtain this decomposition one can allow for imaginary pseudo-fermion parameters, and strong damping in individual pseudo-fermions, without introducing unwanted approximations. For a non-interacting single-resonant level, we benchmark our approach against an analytical solution and an exact hierachical-equations-of-motion approach. We also show that, for the interacting case, this simple method can capture the strongly correlated low-temperature physics of Kondo resonance.

Introduction.-Theorthodox model of strongly correlated electron transport in mesoscopic physics [1][2][3] imagines one or more electronic reservoirs, parameterized by their spectral density, temperature, and chemical potential, coupled to an arbitrary system (typically fermionic impurities interacting with each other, or with additional, e.g., bosonic, environments).While this model can be often well described by weak -coupling theories, strong coupling between system and reservoirs plays an important role in molecular electronics [4], mesoscopic transport through quantum dots (artificial molecules) [5][6][7][8], and quantum thermodynamics [9][10][11][12].While some integrable limits exist [13][14][15], generally speaking one must resort to numerical approaches to capture the non-trivial correlations that can build up between system and reservoirs.
Recently several such approaches based on discrete fermion methods have been proposed and studied to capture this strong-coupling regime [16][17][18], partially inspired by similar techniques developed for strong coupling to bosonic environments [19][20][21][22][23][24][25][26].The hierarchical equations of motion (HEOM) [27][28][29][30][31][32][33] (which nonperturbatively accounts for system-bath entanglement by evolving a hierarchy of time-local equations obtained by repeatedly differentiating the system path-integral representation) have found great success for fermionic systems [34][35][36][37][38][39] and will be used as a benchmark of our results in this work (see also related stochastic methods in [40][41][42]).More recently, the reaction-coordinate method [43], which non-perturbatively models the most relevant degrees of freedom of the environment, was adapted to fermionic systems, and is simple and transparent, both conceptually and in terms of implementation.However, it is arguably limited as the approximations needed to treat the residual bath break down in the wide-band limit.Similarly, a recent approach [12] based on fitting the power spectrum with discrete Lorentzians shows promising results in terms of convergence, but, by relying on only physical modes, it inherits limitations on the width and positivity of the Lorentzian fitting functions.Finally, other methods [44][45][46][47][48][49][50] build up a continuum reservoir from a finite set of discrete damped physical modes, and can be combined with tensor network techniques (such as the process-tensor [51][52][53]) for efficient construction and time-evolution of system properties.
Here, we develop a "pseudo-fermion method" (akin to bosonic pseudomodes [20,23,24,26,54]) based on a discrete set of effective fermions which reproduce the key features of the correlation functions of the continuum bath.By employing unphysical modes with complex couplings, we model the reduced dynamics of the system by solving a Lindblad master equation in the augmented system+pseudo-fermion space.The use of a Lindblad form for the damping of the pseudo-fermions does not, by construction, introduce any approximation.The only approximation arises from how well the total correlation function of the pseudo-fermions matches that of the original bath.In contrast to other methods, [21,[55][56][57][58][59][60][61][62][63], these unphysical degrees of freedom do not have direct connection to the original physical environment allowing op-Figure 1: Illustration of the pseudo-fermion model.In (a) a Fermionic system S interacts with a continuum of environmental Fermionic modes at different frequencies.The spectral density function in (b) describes the distribution of the system-environment coupling.The dynamics of S can be equivalently computed by considering the interaction with a discrete set of pseudo-fermions (red/blue circles) whose dissipative properties are described by a Lindblad master equation (wiggly arrows).This equivalence holds as the pseudofermions are designed to reproduce the correlation functions of the original bath.Specifically, for the spectral density in (b), the red (blue) pseudo-fermions in (c) reproduce the resonant (Matsubara) contribution to the correlation function plotted in (d), see Eq. ( 9).timization of the modeling over an enlarged domain.In summary, we present a conceptually simple framework to simulate non-perturbative effects in Fermionic systems by using a Lindblad master equation.To show that this simplicity does not necessarily come at the expense of modeling power, we benchmark the accuracy of the method by reproducing non-equilibrium and Kondo-physics effects.We note that here the term pseudo-fermion has a different meaning than in non-Hermitian quantum mechanics [64][65][66][67].
Fermionic Open Quantum Systems.-Weconsider a Fermionic system S interacting with a Fermionic Gaussian environment E described by the Hamiltonian H = H S + H E + H I ( = 1), where H S and H E = k ω k c † k c k are the Hamiltonians for the system and environment (so that each environmental Fermion c k is characterized by the frequency ω k ).The interaction Hamiltonian is defined as , where g k are the coupling strengths and s is an odd-Fermion parity operator with support on the system.We assume the bath to be initially in a thermal equilibrium state with inverse temperature β and chemical potential µ.
In this context, the reduced system dynamics depends, in general, on the free-bath statistical properties encoded in correlation functions involving the fields B(t) = k g k c k e −iω k t in the interaction picture.When the free statistics is Gaussian, this dependence can be reduced down to two-point correlations, invoking Wick's theorem to write [68] ρ in terms of the Fermionic time ordering operator T [see, for example, [69], Eq. (5.84)] and the Fermionic influence superoperator F(t, s, C σ ) (see Eq. (A1) in [70]).This operator exactly describes the effects of the bath on the system through its dependence on the system coupling operator s and the translational-invariant two-point correlation functions Here we defined the spectral density , and σ = ±1 to denote the presence/absence of Hermitian conjugation, i.e., B σ=1 = B † and This analysis implies that all memory effects present in non-perturbative regimes can be encoded in the superoperator F(t, s, C σ ).Therefore, Gaussian open systems with identical environmental correlations C σ (t) (and same system-coupling operators s) must have identical influence superoperators and, consequently, an equivalent reduced system dynamics through Eq. ( 1).
Pseudo-Fermion method.-Wenow proceed as schematically shown in Fig. 1(c), i.e., instead of solving the original system + continuum environment Hamiltonian, we define a new model consisting of the original system (S) in contact with a set of discrete pseudo-fermions (pf), themselves in contact with residual environments (re).The key point is that our new environment (pf + re) is designed to produce the same correlation functions of the original environment, and hence produce the same system dynamics, as per Eq.(1).To satisfactorily mimic the original correlation functions with only a finite discrete set of artificial systems, we allow certain parameters defining these pseudo-fermions to be unphysical.Importantly, each residual environment is an idealized quantum white noise (defined by constant spectral densities and frequency-independent equilibrium distributions), thereby allowing the dynamics in the pseudofermion space to be exactly described by a simple Lindblad equation [wiggly arrows in Fig. 1(c)].Surprisingly, the use of Lindblads does not limit the accuracy of the result as the derivation does not rely on approximations as long as the correlation functions produced by said Lindblads closely match the original bath ones.
Explicitly, we consider N pf pseudo-fermions cj , j = 1, . . ., N pf interacting with the system as described by the Hamiltonian  where H pf = N pf j=1 Ω j c † j cj is the free pseudo-fermion Hamiltonian which depends on the (formal) energies Ω j ∈ C and the fields Bσ j = λ j cσ j , in terms of the interaction strengths λ j ∈ C. We further assume the pseudofermions and residual environments to be initially in their equilibrium state ρ eq pf+re .As mentioned, the residual environment associated to each pseudofermion j is modeled as quantum white noise characterized by a formal decay rate Γ j ∈ C and a formal Fermi distribution n j ∈ C. As shown in [70], the free correlations in the pseudo-fermion and residual-environment space, i.e., (4) are defined using operators in the interaction picture and they can be obtained by directly solving the Heisenberg equation of motion leading to This result shows a main feature of the pseudo-fermion method: unphysical properties (for example, the fields B † j not depending on the conjugate of the parameters λ j ) allow to model a more general set of correlation functions (for example, C σ pf,j (0) < 0 requires λ j to be imaginary).Thanks to the Gaussianity hypothesis, the reduced dynamics of a system embedded in the environment made of pseudo-fermions and their residual environments is equivalent to the reduced dynamics of the original model with Hamiltonian H as long as The practical advantage provided by this result is that, by choosing the residual environments to be idealized white noise, their effects on the system+pseudofermions space can be modeled (without approximations, see [68,70]) by the following Lindblad master equation where r = ±1 and j = 1 • • • N pf .Here, ρ r=± S+pf , are the projections of the density matrix into the space with even/odd Fermionic parity and the dissipators D [71]).This is the main result of this article: when the Fermionic environment of a Gaussian open quantum system has free correlations which are equivalent to (or can be approximated by) Eq. ( 6), the reduced system dynamics can be equivalently (approximately) computed by solving the master equation in Eq. (7).
Case study: Lorentzian spectral density.-Wenow provide an explicit pseudo-fermion model to approximate a Fermionic environment initially in an equilibrium state ρ eq E with inverse temperature β and chemical potential µ, and where the interaction with the system is characterized by a Lorentzian spectral density where the frequency parameters W and Γ specify the width and the overall interaction strength, respectively.By inserting this expression into Eq.( 2), it is possible [70] to write the decomposition: Here, the "resonant" and "Matsubara" contributions are defined through Our goal is to implement these contributions using a set of pseudo-environments characterized by the free correlations in Eq. ( 5).To describe this correspondence, we define one resonant and two Matsubara pseudoenvironments (for each Matsubara frequency) by identifying j → res and j → (k, r = ±) in Eq. ( 5).This leads to the following equivalences among correlations which hold [70] when the parameters for the pseudoenvironments are defined as where ∆ ∈ C with |∆| → ∞.We note that the limit |∆| → ∞ might introduce numerical instabilities which can be regularized using intermediate values, such that ∆ 1.The parameter ∆ allows to match the σ-dependence between the correlations in Eq. ( 9) and Eq. ( 5), see [70].However, it is also possible to build an alternative model which does not require any asymptotic limit, but uses four pseudo-fermions for each Matsubara contribution M σ k (t) [70].The mapping above is possible because of the enlarged parameter domain of the model which allows for complex (i.e., unphysical) values.We can appreciate this by noting that (i) M k has a non-trivial overall complex phase requiring complex couplings λ k,± between the pseudomodes and the system, (ii) M σ k contains terms proportional to e αt (α ∈ R) which neither represent a pure physical oscillation (α is not imaginary) nor a pure physical decay (t does not appear in absolute value) ultimately requiring the pseudomodes frequencies Ω k,± to acquire an imaginary part, and (iii) the unbounded nature of the average Fermion number n k,± which implies unphysical density matrices.
The relations in Eq. ( 10) lead to an equivalence between the full pseudo-environment and the original model, i.e., C σ L (t) = C σ pf (t).This shows that a Fermionic environment characterized by the Lorentzian spectral density J L (ω) can be approximated using N pf = 1 + 2K pf pseudo-fermions.Here K pf represents a cut-off in the evaluation of the Matsubara series in Eq. ( 8), i.e., k>0 → K pf k=1 .However, when several Matsubara frequencies effectively contribute to the original correlation C σ L (t) (i.e., K pf 1), it might be more efficient to alternatively finding a best-fit approximation of the full Matsubara series M σ (t) ≡ k>0 M σ k (t), following the ansatz M σ K fit j=1 M σ j,fit (t), where K fit ∈ N and and only then proceeding with the pseudo-fermion mapping, thereby optimizing the number N fit pf = 1 + 2K fit of pseudo-fermions.
Non-equilibrium Single Impurity Model.-Ourfirst example is a single spin-less fermion coupled to two environments with Hamiltonian H = H S + H E + H I , where with B α = k g α,k c α,k for α ∈ {L, R}, indexing the 'left' and 'right' leads.We choose both environments to be identical apart from their chemical potential, so that J L (ω) is generalized to have an α-dependent µ α , and we define ∆µ = µ L − µ R .We then apply the fitting procedure described in Eq. ( 12), so that each environment is described by N fit pf = 1 + 2 pseudo-fermions.As a benchmark we use the standard analytical result (see [14]) for the current, and the HEOM method using the BoFiN-HEOM package [38] with a Pade decomposition of the bath exponents.
For the pseudo-fermions, we can evaluate the current following the logic in [36,37,72], defining the occupation of bath α as N α = k c † α,k c α,k , and the current into the bath α as In the pseudo-fermion formalism we equate B α with a sum over all pseudo-fermions describing environment α.Equivalence of the results [73,74] for the analytical current, the HEOM method, and the pseudo-fermions method is demonstrated in Fig. 2.
Kondo resonance.-Todemonstrate that the pseudofermions can indeed capture non-trivial correlations between system and environment we turn to the example of Kondo resonance [75][76][77][78][79]. Previously, this was used to demonstrate the power of the HEOM method in dealing with Fermionic systems [35], and here we do the same for the pseudo-fermion method.We start with a system containing two interacting Fermions labelled by their spin coupled to a spin-labelled reservoir, where ν ∈ {↑, ↓}.The spectral density of the impurity with spin ν (an experimentally observable quantity [80,81] which can exhibit signatures of Kondo resonance) is defined as A ν (ω) = dte iωt {s ν (t), s † ν (0)} /2π.This can be evaluated using the pseudo-fermion equation of motion.Fig. 3 demonstrates the spectrum for a symmetric example where U = −2 .The pseudofermion method fits the predicted HEOM result (using a converged Padé decomposition of the bath correlation functions) remarkably well.The two side-peaks appear around the system energies, while the Kondo resonance appears at zero frequency, as expected.
Importantly, by considering the full Matsubara series M σ (t), we show in [70] that all terms associated with large (compared to the inverse time-scales of the system) Matsubara frequencies do not effectively contribute to the reduced system dynamics (while for Bosonic environments such terms introduce additional Markovian decay [28]).This implies that a finite truncation of the series is possible, and hence a finite number of pseudomodes can be employed even in the scaling limit.However, despite this, direct numerical solutions using standard Jordan-Wigner mappings become numerically challenging when more than a few pseudomodes are needed.Fortunately, by employing matrix product state techniques, akin to those used in [12,82], we have found that we can access this regime, as demonstrated in Fig. 3(b)(see also [70]).
Summary.-We presented a methodology to generate the reduced dynamics of an arbitrary system interacting with Fermionic environments beyond Born-Markov approximations.The model parameters can break certain physical constraints thereby allowing for greater possibility of optimization.This method allows to accurately In (a) we see that the HEOM result with N k = 2 has not converged, while all other results are (pseudo-fermion and HEOM, N k = 4 results overlap).In (b), which uses the Matsubara series in Eq. (B43), the HEOM results converge for N k = 5, while the standard pseudo-fermion results, relying on K fit = 1 terms in the series in Eq. ( 12) to fit the Matsubara contribution (shown in the inset), have not converged.To efficiently include more exponents, and hence more pseudo-fermions, we employ matrix product states, allowing us to include K fit = 3 term in the Matsubara series (yellow curve in inset) and achieve convergence (solid purple curve).
simulate challenging non-Markovian regimes by solving a Lindblad master equation without resorting to approximations as long as the correlation functions closely match the ones of the original bath.We demonstrated this balance between conceptual simplicity and numerical accuracy by reproducing the exact results of both nonequilibrium transport and Kondo resonance by benchmarking against the HEOM method.The latter is known [35] to perform as well as renormalization-group methods at finite temperatures, and out-perform continoustime quantum Monte Carlo, Green's function equations of motion, and slave-boson mean-field theory.It is important to iterate that, compared to the powerful HEOM method, the pseudofermion method produces identical results, with the advantage of a more transparent interpretation (discrete effective bath fermions) which, as we have demonstrated, are then more intuitively amenable to MPS techniques [83,84].In summary, the method can be applied to simulate the effects of any fermionic environment on a quantum system with a Lindblad master equation thereby allowing challenging regimes to become more accessible by a broader audience.This could allow to design control/environmental-engineering protocols, to study superconducting/Majorana leads or hybrid environments (involving both fermions and bosons).Another interesting future direction is a generalization to compute bath observables and correlations to study non-equilibrium physics within the environment itself, possibly leading to optimize transport properties.

Supplemental Material to:
A pseudo-fermion method for the exact description of fermionic environments: from single-molecule electronics to Kondo resonance Mauro Cirio, Neill Lambert, Po-Chen Kuo, Yueh-Nan Chen, Paul Menczel, Ken Funo, and Franco Nori In this supplemental material, we present technical details supporting the main text.
Appendix A: Fermionic Open Quantum System In this section, we present technical details about the Fermionic open quantum systems described by the original System-Environment Hamiltonian H = H S +H E +H I used in the main text.Under the hypothesis described in the main article, the reduced system dynamics depends solely on the influence superoperator presented in Eq. ( 1) which can be written as where Here we defined s σ=−1 ≡ s and s σ=1 ≡ s † .The parity operator is defined as where in terms of the Fermions d k S which can populate the system.The expression of the influence superoperator explicitly shows us that the effects of the environment are fully encoded in the correlation functions C σ 21 ≡ C σ (t 2 , t 1 ), where ) in terms of the interaction picture operators B(t) = k g k c k e −iω k t , the spectral density J(ω) = π k g 2 k δ(ω− ω k ), and the Fermi distribution n eq E (ω) = (1 + exp{β(ω − µ)}) −1 (in which µ represents the chemical potential).These definitions are the explicit version of the ones in Eq. ( 2) upon noticing that the invariance of ρ eq E under the free evolution of the bath implies time-translational invariance for the correlations.

Appendix B: Pseudofermion Model
In this section we present details about the pseudofermion model presented in the main text.
The model is an approximation scheme which operates in an augmented Hilbert space involving ancillary pseudofermions.The prefix pseudo highlights the absence of physicality constraints other than the requirement to generate the correct reduced system dynamics.Explicitly, we consider N pf pseudo-fermions cj , j = 1, . . ., N pf , each interacting with its own residual bath of Fermions cjk defining a formal open quantum system with Hamiltonian Here, H S+pf = H S + H pf + H I S+pf is the free Hamiltonian in the augmented system+pseudo-fermions space, where H pf = N pf j=1 Ω j c † j cj is the free pseudo-fermion Hamiltonian which depends on the (formal) energies Ω j ∈ C and where (B1) are the system+pseudo-fermions and pseudo-fermions+residual-environment interaction Hamiltonians written in terms of the fields Bσ j (t) = λ j cσ j (t), the energies Ω jk ∈ R, and the interaction strengths λ j ∈ C, λ jk ∈ R. We further assume each auxiliary environment to act as an idealized white noise defined by constant spectral densities J j (ω) = π k λ 2 jk δ(ω − Ω jk ) = Γ j (in terms of rates Γ j ∈ C) and by constant (frequency independent) equilibrium distributions n eq jk = n j (n j ∈ C), see Eq. (B20).Following the ideas introduced in [25] for the Bosonic case, while some of the parameters in this formal model are allowed to take complex values, the dynamics is still defined using the same Schrödinger equation as if all parameters were physical, i.e., without introducing any extra complex-conjugation, hence taking the name pseudo-Schrödinger equation [87].This allows us to use an unphysical model to simulate a wider range of free correlation functions without adding extra complexity.To conclude its characterization, we assume the pseudo-environment to be initially in its equilibrium state ρ pf+re , see Eq. (B15).

A Pseudo-Environment
The Pseudo-Environment described by the Hamiltonian is supposed to be, initially, in the following equilibrium state ) in which Z pf+re ensures that the trace is 1 and where each β j , β jk ∈ R can be found by imposing on top of the already stated white-noise assumption

This leads to
which concludes the characterization of this "pseudoenvironment".We note that the condition in Eq. (B4) allows to write the exponent defining ρ pf+re as a weighted sum over the total number of Fermions in the different pseudo-environments, which is a constant of motion under the free evolution induced by H pf+re = H pf +H I S+pf + H re (as it only contains interaction terms which are written in a "rotating wave" style, i.e., they preserve the total number of bare excitations).

Computing the correlation functions through the Heisemberg equation of motion
The most direct way to compute the correlations in Eq. ( 4) is through the Heisemberg equation of motion.For clarity, in this section we will be omitting the label j used throughout the text to describe independent pseudo-environments.In fact, we will be focusing on a single pseudo-fermion c and its residual environment made of fermions ck .With this notation, the free Hamiltonian of the environment, i.e., the part of the Hamiltonian in Eq. (B2) which has no support on the system, reads which, in Laplace space, becomes into the first equation we get (B8) This equation can be made more explicit by computing To proceed, we consider Which of these two alternative should we use in Eq. (B8)?
The choice depends on whether we are interested in evaluating the dynamics of c(t) for t > 0 or t < 0. For t > 0, we should chose the condition Re(s) > 0 as it is always compatible with the integration path needed to define the inverse Laplace transform.For t < 0, the opposite choice must be made.This means that where the ± depends on whether t > 0 or t < 0. A similar result holds for the conjugate Note that, even if some of the original parameters Ω and Γ might be complex, no complex conjugation is introduced on them, as a pseudo-Schrödinger equation is used to generate the dynamics.
Our goal is now to use Eq.(B12) and Eq.(B13) to compute the correlations in Eq. ( 4) which read where c(t) = exp{iH pf E t}c exp{−iH pf E t} and where (B16) Note that these constraints are essential to ensure invariance under time translation of the correlation.In fact, using them, the equilibrium state takes the form as a function of the total number of Fermions in the structured environment This number is conserved by the free dynamics, i.e., [H pf , N pf ] = 0 which implies that where t = t 2 − t 1 .Therefore, using Eq.(B12) and Eq.(B13), we can immediately compute where L −1 t is the inverse Laplace transform which can be computed to obtain which, reintroducing the indexes j used in the main text to label independent pseudo-environments, leads to Eq. ( 5) by linearity.

Markovian Regime
In this section, we review the limit in which the environment is Markovian.This regime is defined when all the memory effects described by the Fermionic influence superoperator F(t, s, C σ ) in Eq. ( 1) are negligible.The Markovian limit can be recovered [17,68] in the quantum white-noise limit characterized by a constant spectral density and a constant equilibrium distribution, i.e., In fact, these conditions, once inserted in Eq. ( 2) lead to which, arguably, more commonly defines the Markov approximation.As shown in [68], in the case of a deltacorrelated environment, the expression for the functional superoperator in Eq. ( 1) drastically simplifies, leading to an effective generalized Lindblad equation of motion In non-Markovian regimes, the influence superoperator F is able to model more complex physical scenarios.For example, in the limit where the system is a singlelevel impurity one recovers the single-impurity Anderson model, which is integrable in the non-interacting limit.Beyond this case, with e.g., interacting or non-linear impurities, one must resort to numerical methods.This analysis can be easily adapted to the problem of tracing out the residual environments of the pseudofermion model defined in Eq. (B2).In fact, the residual environment of each pseudo-fermion is modeled as idealized white noise and, as a consequence, it can be traced out by simply considering the composite sys-tem+pseudofermions in place of S in Eq. (B21) leading to Eq. ( 7).

Correlations for a Lorentzian spectral density
We now consider the Lorentzian spectral density where a = ω 0 + iW and ā = ω 0 − iW .With this spectral density, the correlation functions in Eq. ( 2) take the form and can be evaluated by noticing that the poles of the integrand are located at a, ā, and ω k = µ + ix k , where x k = (2k − 1)iπ/β, for k ∈ N.For positive (negative) time arguments, we can close the contour in the upper (lower) complex plane to derive, for t > 0 B24) Noticing that ωk = 2µ − ω k , the above expressions are compatible with the general relation C σ (−t) = Cσ (t).
(B25) Using the Matsubara expansion (see [1], Eq. (3.5), pag.110) we can alternatively write, for t > 0, where we noticed that 1/(a − 2µ + ω k ) − 1/(ω k − ā) = 0, ∀k.Using C σ (−t) = Cσ (t), together with Eq. (B25), we can write, for t > 0, (B28) Using C σ (−t) = Cσ (t) again, we can extend these results to t < 0 to write, for any t ∈ R where ) where sg(t) = t/|t| for t = 0, is the sign function.One interesting feature of this decomposition is that the pole at x k = W is explicitly removed from the notation (because of the presence of a corresponding zero in the numerator).The first line of Eq. (B30) reproduces the first line of Eq. ( 9) in the main text.In order to reproduce the second line of Eq. ( 9) there is a little more work to do.To achieve this, we use the identity in Eq. (B48) to write which is the second line in Eq. (9).
We finish noting that, in the zero-temperature limit (β → ∞), the Matsubara frequencies x k approach a continuum so that ) . (B34)

Correspondence to pseudo-environments
Here we provide details about modeling the correlations C σ res (t) and M σ k (t) in Eq. ( 9) using Fermionic pseudo-environments.We start from the resonant contribution C σ res (t).We want to find the parameters of a pseudo-environment such that its free correlation function C σ pf,res (t), obtained using the identification j → res in Eq. ( 5), fulfills Using Eq. ( 9) and Eq. ( 5) the equation above translates to The equivalence in Eq. (B36) can be imposed by defining which fully characterize the resonant pseudoenvironment.Similarly, for each k, the Matsubara contribution M σ k (t) in Eq. ( 9) can be reproduced using two pseudo-environments.Explicitly, identifying j → (k, r), r = ± in Eq. ( 5), we want to impose which, using Eq. ( 9) and Eq. ( 5) is equivalent to B39) Imposing the equation above requires a bit of attention as the frequency r(W − x k )/2 appearing in the expression of M σ k (t) is not multiplied by the parameter σ as in the correlation for the pseudo-environment.On the contrary, the coefficients multiplying the exponential in the correlation for the pseudo-environment do depend on σ while M k does not.Luckily, the Matsubara contributions M σ k (t) are written in terms of a difference between exponentials with opposite frequencies which offers the opportunity to define the parameters characterizing the Matsubara pseudo-environments as where ∆ ∈ C in the limit |∆| → ∞ so that (B41) which, inserted in Eq. (B38), leads to B42) where, in the last step, the limit |∆| → ∞ was taken.The above equation is equivalent to the expression for M σ k (t) given in Eq. ( 9), thereby completing the proof.In numerical applications, the limit ∆ → ∞ might introduce some numerical instabilities which can be regularized using intermediate values such that ∆ 1.
Interestingly, it is also possible to build a pseudoenvironment which does not resort to any asymptotic parameter (∆ above).To achieve this, we need to introduce four pseudo-environments to model each M σ k (t) to impose through the identification j → (k, r, σ ), r, σ = ± in Eq. ( 5).In order to fulfill Eq. (B43), we can choose the parameters which correspond to a β → σ × ∞ limit, so that, see Eq. ( 5), which, inserted in Eq. (B43), leads to ) which, similarly as in the previous analysis, is equivalent to the expression for M σ k (t) given in Eq. ( 9), thereby completing the proof.

Proof of an identity for the sign function
Let us define sg(x) = x/|x| for x = 0. We have for any ω ∈ C. In a more compact notation, the previous equation can be written as e ω(t+|t|) − e ω(|t|−t) = sg(t)(e 2ω|t| − 1) .

(B48)
Appendix C: Effective model for fast decaying terms in the Matsubara correlation For Bosonic environments, when a term in the Matsubara series has a decay rate which is larger than the highest frequency Ω S which can be associated to the system, it is possible to approximate its effect on the system by adding an extra dissipator to the master equation [28].In this section, we analyze the same limit in the case of a Fermionic environment interacting with the system through the Lorentzian spectral density J L (ω).In this case, the Matsubara series takes the form described in Eq. ( 9), i.e., , (C1) where x k = (2k − 1)π/β.As in the Bosonic case, our starting point is the following limit-representation of the Dirac delta [28] Among the exponentials present in the Matsubara series above, the ones which can be modeled with a delta contribution are those for which either x k Ω S or W Ω S .We note that the former possibility corresponds to a "high-temperature" limit (as it corresponds to βΩ S (2k − 1)π); while the latter corresponds to a "broad spectral density" limit.Whenever we are in one of these regimes, the corresponding contribution in the Matsubara series can be approximated as where A eff = iΓ eff in terms of an effective decay rate Γ eff .We can now compute the corresponding contribution to the influence superoperator, i.e., we can compute Eq. ( 1) with the replacement C σ (t 2 , t 1 ) → M eff (t 2 − t 1 ).We find where where t = t 2 − t 1 ≥ 0, which justifies the last step.We then obtain This means that, in the even/odd sector, in the Schrödinger picture, we have which, if s is such to satisfy the Fermionic anticommutation rules, is zero.
In conclusion, when a term in the Matsubara correlation function in Eq. ( 9) can be modeled as a delta function, it does not bring any effect on the system dynamics.For example, this implies that the Matsubara correlation function can be neglected when both the following conditions are satisfied: the high temperature limit (i.e., 1/β much bigger than the highest frequency associated with the system) and a wide Lorentzian spectral density (i.e., with a width W much bigger than the highest frequency associated with the system).Here we describe how to simulate Eq. ( 7) of the main text with tensor networks.In the algorithm we have used the superfermion representation of the master equation and the Swap-gate technique to treat long-ranged hopping which arises due to the use of energy eigenbasis for the environments.Since we closely follow Ref. [12], we only elaborate some key aspects in designing the algorithm instead of presenting all details here.
The superfermion representation is a way to map the fermionic master equation onto a non-Hermitian Schrödinger equation, similiar to the purification technique of simulating master equations [85].In the superfermion representation, an extra auxiliary fermion is introduced to each physical fermion, doubling the dimension of the Hilbert space.The arrangement of physical and auxiliary fermions is in principle arbitrary, but for the purpose of tensor network simulation, it is beneficial to arrange them in a intertwined order, see Fig. 1.
A key object in the construction is the so-called leftvacuum state defined by Note that Eq. (D3) has no explicit dependence on the fermion parity, which is one of the main advantages in taking the superfermion representation.
For the tensor network simulation, it is more efficient to combine a physical fermion and its partner into one MPS site, so that the physical index has dimension four in our simulation.One immediately sees in this way that the Lindblad terms in L are all local, which act on a single MPS site; Hence easy to be handled in the simulation.
The hopping terms in the Hamiltonian are between the pseudofermions and the system, which constitutes the so-called star-geometry, needs some special care.In our implementation we use the Swap gate Therefore any long-ranged two-site gate can be decomposed into a sequence of Swap gates and nearestneighbour gates, which can be implemented in tensor networks efficiently.We have used a second-order Trotter decomposition of the propagator and by carefully choosing the order of gates, many Swap gates arising due to long-range hopping can be merged into the identity operator, thus avoiding numerical overheads.

Figure 2 :
Figure 2: Current into the right reservoir from a single impurity (coupled to two reservoirs) as a function of chemical potential difference ∆µ = µL − µR.For both reservoirs we choose = Γ, W = 2.5Γ, and β = 1/(0.2Γ).The truncation parameter for the HEOM results is nmax = 2.

3 Figure 3 :
Figure 3: Spectral density A(ω) for a single impurity containing two interacting fermions coupled to a spin-labelled reservoir.Here µ = 0, U = 3πΓ, = −U/2, and β = 1/(0.2Γ).(a) shows the results for width W = 2.5Γ, while (b) is for W = 10Γ.In (a) we see that the HEOM result with N k = 2 has not converged, while all other results are (pseudo-fermion and HEOM, N k = 4 results overlap).In (b), which uses the Matsubara series in Eq. (B43), the HEOM results converge for N k = 5, while the standard pseudo-fermion results, relying on K fit = 1 terms in the series in Eq. (12) to fit the Matsubara contribution (shown in the inset), have not converged.To efficiently include more exponents, and hence more pseudo-fermions, we employ matrix product states, allowing us to include K fit = 3 term in the Matsubara series (yellow curve in inset) and achieve convergence (solid purple curve).

Figure 1 :
Figure 1: Schematics for the superfermion representation.Blue circles correspond to physical fermions and orange circles auxiliary ones.They are combined in pairs to give one MPS site, shown as dashed ellipses.The first one comes from the system while others from the bath.Solid lines below indicate the long-ranged hopping terms in the Hamiltonian.
|I = |00, 00, • • • , 00 , (D1) where 0 (1) stands for empty (filled) fermionic site and we have combined the physical fermion and the corresponding auxiliary one (marked by underlines) into pairs.By applying the density matrix and master equation to |I from the left and making use of the following conjugation rules d † |I = −d|I , d|I = d † |I (D2) where d indicates system fermions (s) and pseudofermions (c) and d its auxiliary counterpart, we arrive at a non-Hermitian Schrödinger equation d|ρ /dt = −iL|ρ (D3) with |ρ = ρ|I and j (1 − n j ) cj cj − 1 2 (c † j cj + c † j cj ) .

S
= (I 2 ⊗ S ⊗ I 2 )(S ⊗ S)(I 2 ⊗ S ⊗ I 2 ) (D4)to exchange the local states of nearest-neighbour MPS sites, here I 2 is the 2 × 2 identity matrix and S can be explicitly written as in the basis made by {|00 , |01 , |10 , |11 }.Intuitively, S interchanges the states of two adjacent fermionic sites, while S swaps the states of two adjacent MPS sites (each consisting of two fermionic sites).For example, by labeling an MPS state as |n 1 n 1 , n 2 n 2 , the operator I 2 ⊗S ⊗I 2 swaps n 1 ↔ n 2 , the operator S ⊗ S swaps n 1 ↔ n 2 and n 1 ↔ n 2 , and the operator I 2 ⊗ S ⊗ I 2 swaps n 1 ↔ n 2 , resulting in the expression presented in Eq. (D4).