Optimized auxiliary representation of a non-Markovian environment by a Lindblad equation

We present a general scheme to map correlated nonequilibrium quantum impurity problems onto an auxiliary open quantum system of small size. The infinite fermionic reservoirs of the original system are thereby replaced by a small number $N_B$ of noninteracting auxiliary bath sites whose dynamics is described by a Lindblad equation. Due to the presence of the intermediate bath sites, the overall dynamics acting on the impurity site is non-Markovian. With the help of an optimization scheme for the auxiliary Lindblad parameters, an accurate mapping is achieved, which becomes exponentially exact upon increasing $N_B$. The basic idea for this scheme was presented previously in the context of nonequilibrium dynamical mean field theory. In successive works on improved manybody solution strategies for the auxiliary Lindblad equation, such as Lanczos exact diagonalization or matrix product states, we applied the approach to study the nonequilibrium Kondo regime. In the present paper, we address in detail the mapping procedure itself, rather than the many-body solution. In particular, we investigate the effects of the geometry of the auxiliary system on the accuracy of the mapping for given $N_B$. Specifically, we present a detailed convergence study for five different geometries which, besides being of practical utility, reveals important insights into the underlying mechanisms of the mapping. For setups with onsite or nearest-neighbor Lindblad parameters we find that a representation adopting two separate bath chains is by far more accurate with respect to other choices based on a single chain or a commonly used star geometry. A significant improvement is obtained by allowing for long-ranged and complex Lindblad parameters. These results can be of great value when studying Lindblad-type approaches to correlated systems.


I. INTRODUCTION
Strongly correlated systems out of equilibrium have recently attracted considerable interest due to progress in several experimental fields, such as ultrafast pump-probe spectroscopy 1,2 , ultracold quantum gases 3-7 , solid-state nanotechnology [8][9][10] . These advances have also prompted the interest in related theoretical questions concerning thermalisation [11][12][13] , dissipation and decoherence 14 , and nonequilibrium quantum phase transitions 15 . An interesting aspect is the interplay between correlation and dissipation in systems in which the latter is not included phenomenologically but is part of the microscopic model. The challenge lies in the fact that the Hilbert space for correlated fermionic systems increases exponentially with system size. For a finite system, on the other hand, the spectrum remains discrete and dissipation does not occur. When considering purely fermionic correlated systems, dissipation is usually modeled by infinite reservoirs of noninteracting fermions. These reservoirs are in contact with a correlated central region of interest. A paradigmatic example of such a system is the single-site Kondo or Anderson impurity model 16 . If there is just one reservoir with a single chemical potential µ and temperature T , then the whole system (typically) reaches thermodynamic equilibrium. Alternatively, one can consider a nonequilibrium situation in which several reservoirs with different µ and T are in contact with the central region.
Since the reservoirs are infinite they act as dissipators and the system in most cases reaches a nonequilibrium steady state in which a particle and/or heat current flows across the central region. 67 There are several approaches to treat such systems numerically. Some of them start out from the situation in which the central region and the reservoirs are decoupled which allows the individual systems to be treated exactly. There are different schemes to include the missing coupling between the reservoirs and the central region. First of all, one could carry out a perturbative expansion in terms of the reservoir-central region coupling. Low energy properties are better addressed within a renormalisation-group treatment of the perturbation (see, e.g. Ref. 17). Alternatively, one can try and compute the self-energy (most nonequilibrium quantities of interest follow from Dyson's equation) for the correlated sites based on finite clusters consisting of the central region plus a small number N r of reservoir sites. This is done in nonequilibrium cluster perturbation theory 18,19 , whose accuracy increases with increasing N r . A generalization of this idea is the nonequilibrium variational cluster approach, [20][21][22] , where single-particle parameters of the model are optimized self-consistently, which allows for the adjustment of the self-energy to the nonequilibrium situation.
In a different type of approach one tries to "eliminate" the degrees of freedom of the reservoir and take into account its effects on the dynamics of the interacting central region. Formally, this can be expressed in terms of a functional integral whereby the part of the action describing the reservoirs, which is quadratic, is integrated out and one obtains an effective action restricted to the central region only, whereby the effects of the reservoir introduce couplings with retardation effects. These physically describe processes in which particles jump from the central region to the reservoir and then come back after a certain delay. This retarded action can be treated, e.g., via continuous time Monte Carlo approaches 23,24 , which, however are plagued by the minus sign problem. Due to retardation, exact diagonalization approaches are not appropriate. There are several other ways to achieve this elimination of the reservoir degrees of freedom. Renormalisation-group approaches are certainly convenient whenever one is interested in the low-energy sector. 25,26 The numerical renormalisation group has proven extremely powerful for quantum impurity models. 16

Markovian approximations and beyond
Another approach consists in treating the coupling to the reservoir within the Born-Markov approximation. In this way, the effect of the reservoir is to introduce nonunitary dynamics in the time dependence of the reduced density operator ρ f of the central region leading to the Lindblad equation 27 , which is a linear, time-local equation for ρ f preserving its hermiticity, trace, and positivity. One important precondition for the validity of this mapping, however, is the Markovian assumption that the decay of correlations in the reservoir is much faster than typical time scales of the central region. 27 As pointed out, e.g. in Refs. 27,28 the approximations leading to the Markovian Lindblad master equation are justified provided the typical energy scale Ω of the reservoir is much larger than the reservoir-central region coupling. However, for a fermionic system, Ω can be estimated as min(W, max(|µ − ε|, T ), where W is the reservoir's bandwidth, and ε is a typical single-particle energy of the central region. Therefore, even in the wide-band limit W → ∞, the validity of the Markov approximation is limited either to high temperatures or to chemical potentials far away from the characteristic energies of the central region. As a matter of facts, the effect of a noninteracting reservoir with W, |µ| → ∞ (or T → ∞ with finite µ/T ) can be exactly written in terms of a Lindblad equation. This can be easily deduced from the "singular coupling" derivation of the Lindblad equation 27 . This is valid independently of the strength of the coupling between central region and reservoir. A nontrivial situation is obtained by introducing different reservoirs with different particle densities. The pleasant aspect of this limit is that the Lindblad parameters depend on the properties of the reservoir and of its coupling with the central region only, but not on the ones of the central region.
This is in contrast to the more standard weak-coupling Born-Markov version in which the Lindblad couplings (see, e.g. 27,29 ) depend on the central region's properties. To illustrate this, consider a central region consisting of a single site with energy ε f , i.e. with Hamiltonian (omitting spin) and reduced density matrix ρ f . The part of the Lindblad operator L b describing the coupling to a noninteracting reservoir is given by with Here, Γ is proportional to the reservoir's density of states at the energy ε f , and f F is the Fermi function which obviously contains the information on the chemical potential and temperature of the reservoir but also on the onsite energy ε f in the central region. This could be unsatisfactory since one would like to describe the effect of the reservoir in a form which is independent of the properties of the central region, especially when the latter consists of many coupled sites. One possible way to eliminate the dependence of the Lindblad couplings on the parameters of the central region is to use an intermediate auxiliary buffer zone (mesoreservoir) between the Lindblad couplings and the central region (see, e.g. [30][31][32] The buffer zone consists of isolated discrete sites (levels) each one coupled to a Markovian environment described by Lindblad operators with the same T and µ as given in Eq. (2),Eq. (3). If the buffer zone is sufficiently large, i.e. if its levels are dense enough, then one can show that the buffer zone including Lindblad operators yields an accurate representation of the reservoir, which becomes exact in the limit of an infinite number of levels. Importantly, the parameters of this buffer zone do not depend on the central region's properties. The disadvantage of this approach is that one needs quite a large number of buffer levels, especially at low temperatures where the Fermi function is sharp. Consequently, the many-body Hilbert space is too large and the treatment of a correlated problem becomes prohibitive.

This work
In this paper, we show that the accuracy of the bufferzone idea can be improved significantly even with a moderate number of auxiliary buffer levels (sites) by allowing for more general Lindblad couplings, which are adjusted to optimize the representation of the physical reservoirs independently of the parameters of the central region. In particular, we show that allowing for long-ranged, and even complex Lindblad terms (see below) dramatically improves the accuracy of the reservoir's description for a fixed number of auxiliary sites. In the case of a singleimpurity model, already a small number of sites (4 to 6) is enough in order to reach a very good accuracy [33][34][35] sufficient to resolve the splitting of the Kondo peak at finite bias. This is crucial, since the Hilbert space of such a small system can still be treated by Krylov-space methods 68 Larger systems can be tackled by matrix-productstate approaches for open quantum systems [36][37][38][39] . Due to its rapid convergence this scheme can be used as an accurate nonequilibrium solver for correlated impurity problems, which even in equilibrium becomes competitive with other established approaches. Here we discuss in particular several way to optimally represent a physical ("ph") reservoir by means of an auxiliary ("aux") one consisting of a small number of noninteracting fermionic sites and arbitrary Lindblad terms. Presenting results of the application to interacting systems is not the main goal of this work, so there will be only a short discussion in Sec. II E, and we refer to previous publications 33, 34,39 for details. Here, we are interested in a systematic analysis of the performance of different geometries of the auxiliary reservoir, see Fig. 1, including a scaling analysis of the accuracy as a function of the number of bath sites N B and a discussion of the importance of long-ranged Lindblad terms. This paper is organized as follows: in Sec. II A, we introduce the models we are interested in and define the basic notation. In Sec. II B, we illustrate the most important aspect of this work, namely the mapping of the physical Hamiltonian problem onto an auxiliary open quantum system described by a Lindblad equation. In Sec. II C, we present the expressions for the noninteracting Green's function of the auxiliary system, and in Sec. II D we illustrate the fit procedure. In Sec. II E, we briefly discuss the relation with the interacting case. In Sec. III, we present in detail the convergence of the fit as a function of N B for the different geometries presented in Sec. II F and for different temperatures, and discuss the advantages and disadvantages of these setups. Finally, in Sec. IV we summarize our results and discuss possible improvements and open issues. In three appendices we present technical details of the minimization procedure (Sec. A), show the explicit form of the matrices for the different geometries (Sec. B), and discuss certain redundancies of the auxiliary system (Sec. C).

A. Model
We begin with a general discussion, which we eventually apply to the single-site Anderson impurity model. In the general case the central region may represent a small cluster or molecule. The Hamiltonian of the physical system at study is written as where H f is the Hamiltonian of the central region describing a small cluster of interacting fermions, H α is the Hamiltonian of the reservoir α describing an infinite lattice of noninteracting particles, and H αf is the coupling between central region and reservoirs.
consists of a noninteracting part and an interaction term H U . The fermions in the reservoirs can be described by in usual notation. For simplicity, spin indices are not explicitly mentioned here. Quite generally, a suitable single particle basis "star representation" can be chosen such that ε αp,p ∼ δ p,p . H αf is taken to be quadratic in the fermion operators: and d i (f i ) are the fermionic destruction operators on the reservoir's (central region's) sites i. We are interested in a steady state situation, although the present approach can be easily extended to include time dependence, especially if this comes from a change of the central region's parameters only. In the steady state we can Fourier transform with respect to the time variables so that the Green's functions depend on a real frequency ω only, which here is kept implicit. We assume that initially the hybridization H αf is zero and the reservoirs are separately in equilibrium with chemical potentials µ α and temperatures T α . Then the H αf are switched on and after a certain time a steady state is reached. We use the non-equilibrium (Keldysh) formalism [40][41][42][43][44] whereby the Green's function can be represented as a 2 × 2 block matrix where the retarded G R , advanced G A , and Keldysh G K components are matrices in the site indices (i, j) of the central region. We will adopt the above underline notation in order to denote this 2 × 2 structure. We use lowercase g (g α ) to denote Green's function of the decoupled central region (reservoir α), while uppercase G represent the full noninteracting Green's function of the central region. For simplicity we omit the subscript 0 , since in this paper we deal mainly with noninteracting Green's functions anyway. We use the subscript int for interacting ones. G is easily obtained from the Dyson equation as where is the reservoir hybridization function (commonly called bath hybridization function) in the Keldysh representation. The retarded Green's functions g R α for reservoirs with non-interacting fermions in equilibrium can be determined easily by standard tools, and the Keldysh components g K α can be obtained from the retarded ones by exploiting the fluctuation-dissipation theorem: which is valid since the uncoupled reservoirs are in equilibrium. Here, and f F (ω, µ α , T α ) is the Fermi function at chemical potential µ α and temperature T α . From now on, for simplicity of presentation, we restrict to the Anderson impurity model (SIAM) in which the central region, described by Eq. (5), consists of a single site, i.e. there is only one value for the index i, which we drop, and The idea we are going to present in Sec. II B can be immediately extended to the case of a central region consisting of many sites in which each site is connected to separate reservoirs. In the most direct fashion, this can be done with exactly the same approach as formulated here for the SIAM by just mapping each reservoir independently onto auxiliary Lindblad bath sites. An interesting application is for example, the case of an interacting chain coupled on both sides to reservoirs with different chemical potentials 45 . Also the extension to the case of arbitrary (quadratic) couplings with the reservoirs that intermix the central region sites, relevant, e.g., for cluster-DMFT, is conceptually straightforward, although more complex.

B. Mapping onto an auxiliary master equation
A crucial point in the following considerartions is the fact that, even in the interacting case, the influence of the reservoirs upon the central region is completely determined by the bath hybridization function ∆(ω) only. In other words, any interacting correlation function of the central region solely depends on the central region Hamiltonian H f and on ∆. This result is well known, at least in equilibrium, and can be easily proven, for example diagrammatically.(see footnote 69 ) The argument holds independently on whether one works with equilibrium or nonequilibrium Green's functions. Moreover, it crucially depends on the fact that the reservoir is noninteracting. This can be exploited to choose different representations for the reservoir depending on convenience. In equilibrium, especially in connection with numerical renormalisation group (NRG), one uses either the diagonal ("star") representation in which the ε αp,p ∼ δ p,p are diagonal, as in Eq. (7), or the "chain" representation in which they describe a nearest-neighbor chain (see, e.g. Ref. 46). While for a continuous density of states one needs, in principle, an infinite number of sites for the reservoirs, one can approximate the physical 70 ∆ ph Eq. (11) by an auxiliary ∆ aux corresponding to a bath with a finite number of sites and optimize their parameters ε αpp and v αp via a best fit. Notice that for this Hamiltonian representation the space of parameters is redundant, so that one can restrict, for example to diagonal ε αpp ∼ δ p,p and real v αpj . This is the "star" representation mentioned above. The "chain" representation is given by a single nonvanishing v αpj and local or nearest-neighbor ε αpp and is obtained from the star representation via a unitary transformation. Therefore, for N B sites of the auxiliary system one has 2N B parameters available for the fit. This approach is used, for example, for exact-diagonalisation-based Dynamical-Mean-Field Theory (ED-DMFT) 47,48 . Here, the parameters are optimized by fitting the bath hybridization function in Matsubara space. The auxiliary system of bath sites plus impurity is then solved by Lanczos exact diagonalisation (ED). 49 .
Clearly, the same fit strategy is inconvenient out of equilibrium for several reasons. First of all, the auxiliary system cannot dissipate, since it is finite, and a steady state cannot be reached. In Refs. 33,50 we have suggested a different approach (Auxiliary Master Equation Approach, AMEA), which adopts an auxiliary reservoir, consisting of a certain number N B of bath sites which are additionally coupled to Markovian environments described by a Lindblad equation Here, the Hamiltonian for the auxiliary system is given by (we reintroduce spin) and enters the unitary part of the Lindblad operator The dissipator L D describes the coupling of the auxiliary sites to the Markovian environment and is given by The indices i, j in Eq. (17), Eq. (18) run over the impurity i = f (we identify c f σ = f σ ) and over the N B bath sites. 71 Similarly to the case of the ED-DMFT impurity solver mentioned above, the idea is to optimize the parameters of the auxiliary reservoir in order to achieve a best fit to the physical bath hybridization function Eq. (11), i.e., for a given N B , ∆ aux (ω) should be as close as possible to ∆ ph (ω): As for the ED-DMFT case, one can choose a singleparticle basis for the auxiliary bath such that the matrix E is sparse, 72 i.e. it has a "star" or a "chain" form, and is real valued. However, there is no reason why the matrices Γ (1) and Γ (2) should be sparse and real in the same basis as well, and, in fact, as discussed below, for an ED treatment of the Lindblad problem it is convenient to allow for a general form in order to optimize the fit. This larger number of parameters allows one to fulfill Eq. (19) to a very good approximation. The introduction of dissipators (18) additionally allows to carry out the fit directly for real ω, see Sec. II D below, since ∆ aux (ω) is a continuous function. This makes this approach competitive with ED-DMFT for the equilibrium case as well.
Notice that Eq. (18) is not the most general form of the dissipator, and one could think of including Lindblad terms that contain four or more fermionic operators, or also anomalous and spin-flip terms. This would increase the number of parameters available for the fit. However, the latter would violate conserved quantities and the former would describe an interacting bath, so that the argument of Sec. II B (footnote 69 ) does not apply. As a matter of fact, the exact equivalence to a noninteracting bath 71 only holds for a quadratic form of the Lindblad operator as in Eq. (18).
Once the optimal values of the matrices E, Γ (1) and Γ (2) for a given physical model are determined for the non-interacting system, one could solve for the dynamics of the correlated auxiliary system defined by Eq. (15), which amounts to a linear equation for the reduced manybody density matrix. If the number of sites of this system is small, one can solve exactly for the steady state and the dynamics of this interacting system by methods such as Lanczos exact diagonalization or Matrix Product States (MPS) [36][37][38][39] .

C. Computation of the Auxiliary bath hybridization function
In order to carry out the fit Eq. (19), we need to compute the auxiliary reservoir hybridization function ∆ aux (ω) for many values of the bath and Lindblad parameters. This can be done in an efficient manner since only noninteracting Green's functions are needed, see also Eq. (10) and the discussion above. Computing the singleparticle Green's function matrix G of Eq. (15) amounts to solving a noninteracting fermion problem, which scales polynomially with respect to the single-particle Hilbert space N B + 1. A method to deal with quadratic fermions with linear dissipation based on a so-called "third quantisation" has been introduced in Ref. 51. We adopt the approach of Ref. 31 in which the authors recast an open quantum problem like Eq. (15) into a standard operator problem in an augmented fermion Fock space with twice as many sites and with a non-Hermitian Hamiltonian. 31,52,53 This so-called super-fermionic representation is convenient for our purposes, not only to solve for the noninteracting Green's functions but also to treat the many-body problem in an analogous framework to Hermitian problems. An analytic expression for the the noninteracting steady-state retarded and Keldysh auxiliary Green's functions was derived in Ref. 33. An alternative derivation, which does not rely on super-fermions is given in Ref. 54. For the retarded component we get 72 and the Keldysh component of the inverse Green's function reads yielding the Keldysh Green's function The f f component of G is the auxiliary impurity Green's function From this one can determine the retarded component of For the Keldysh component, one has to carry out two inversions of Keldysh matrices (see, e.g. Ref. 43) yielding where the contribution from g K is infinitesimal.

D. Fit procedure
From the equations above we can efficiently compute ∆ aux (ω) for a given set of parameters of the auxiliary reservoir. The numerical effort for a single evaluation is low and scales only at most as O(N 3 B ). We introduce a vector of parameters x which yields a unique set of matrices E, Γ (1) and Γ (2) within a chosen subset (see, e.g. Fig: 1 and App. B), quantify the deviation from Eq. (19) through a cost function and minimize χ(x) with respect to x. The normalization χ 0 is hereby chosen such that χ(x) = 1 when ∆ aux (ω) ≡ 0. It is important to note that both, the retarded and the Keldysh component must be fitted. Due to Kramers-Kronig relations, the real part of ∆ R ph (ω) is fully determined by its imaginary part, provided the asymptotic behavior is fixed. Therefore, we can restrict to fit its imaginary part, while ∆ K ph (ω) is purely imaginary. Furthermore, in Eq. (26) we introduced a cut-off frequency ω c and a weighting function W (ω). In this work we take W (ω) = 1 and ω c = 1.5 D, with D the half-bandwidth of ∆ ph (ω). Different forms of W (ω) can be used in order to increase for instance the accuracy of the fit near the chemical potentials. The minimization of Eq. (26) constitutes a multi-dimensional optimization problem and appropriate numerical methods for it are discussed in Sec. A.
As asymptotic limit we require here ∆ aux (ω) → 0 for ω → ±∞, which is obtained when setting Γ For simplicity, we restrict here to the particle-hole symmetric case. This reduces the number of free parameters in E, Γ (1) and Γ (2) . For the case that the impurity site f is located in the center and that one has an even number of bath sites N B , particle-hole symmetry in the auxiliary system is obtained when for i, j ∈ {1, . . . , N B + 1}. More details for the particular form of E, Γ (1) and Γ (2) are given below in App. B.

E. Interacting case
Despite the fact that the solution of the interacting impurity problem is not the main topic of the present work, it is the main purpose of the overall approach. We thus briefly discuss here some relevant issues, in connection to the evaluation of particular observables of the physical system from results of the auxiliary system. More details can be found in Refs. 33,39 As already discussed, by mapping onto an auxiliary interacting open quantum system of finite size described by the Lindblad equation Eq. (15), we obtain a manybody problem which can be solved exactly or at least with high numerical precision, provided N B is not too large. In Ref. 33 we presented a solution strategy based on exact diagonalization (ED) with Krylov space methods, and in Ref. 39 one based on matrix product states (MPS). In the end both techniques allow us to determine the interacting impurity Green's function G aux,int (ω) of the auxiliary system. As discussed above, in the limit ∆ aux (ω) → ∆ ph (ω) (i.e. for large N B ) this becomes equivalent to the physical one G ph,int (ω). However, this equivalence only holds for impurity correlation functions, and, for example, it does not apply for the current flowing from a left (α = l) to a right (α = r) reservoir across the impurity. Therefore, the current evaluated within the auxiliary Lindblad system does not necessarily correspond to the physical current even for large N B , unless one fits the bath hybridisation functions ∆ ph,α (ω) for the left and right reservoirs separately. Such a separate fit, however, is not necessary and would simply worsen the overall accuracy for a given N B . Once the approximate G ph,int (ω) ≈ G aux,int (ω) is known, the current of the physical system can be evaluated by means of the wellknown Meir-Wingreen expression 43,55,56 , however, by using the Fermi functions and density of states (hybridisation functions) of the two physical reservoirs separately. Therefore, the knowledge of G aux,int (ω) enables one to compute most quantities of interest.
An additional step consists in extracting just the selfenergy from the solution of the auxiliary impurity system and inserting it into the Dyson equation for the physical system with the exact physical noninteracting Green's function Clearly, this step is only useful when the relation Eq. (19) is approximate, since for ∆ aux (ω) → ∆ ph (ω) also the noninteracting Green's functions G ph (ω) and G aux (ω) would coincide, i.e. in the hypothetical N B → ∞ case, and one could just set G ph,int (ω) → G aux,int (ω). For finite N B this substitution has the advantage that in (28) the noninteracting part G ph (ω) is exact, and the approximation Eq. (19) only affects the self energy.

F. Different geometries for the auxiliary system
With the goal in mind of providing the best approximation to the full interacting impurity problem described by the Hamiltonian Eq. (4), we would like to approximate ∆ ph (ω) by ∆ aux (ω) as accurately as possible for a given number of bath sites N B . In principle, one has the freedom to choose different geometries for the auxiliary system and a generic set of five different setups is depicted in Fig. 1. (An explicit form of the corresponding matrices for N B = 4 is given in App. B.) For large N B they all converge to the exact solution ∆ aux (ω) → ∆ ph (ω), the "1 chain n.n." "star" "2 chains onsite" "2 chains n.n." "full" FIG. 1: Sketch of the five geometries (setups) for the auxiliary system Eq. (15). An explicit form of the corresponding matrices for NB = 4 is given in App. B. The impurity is represented by a red circle while the bath sites are filled green ones. The hoppings described by the matrix E are represented by thick black lines. The couplings to the Markovian environments given by Γ (1/2) are expressed by grey lines connected to empty (Γ (1) ) or full (Γ (2) ) reservoirs. On-site terms in the Γ (1/2) -matrices are illustrated as a double grey line. The setup "full" represents the most general case with dense Γ (1) and Γ (2) matrices, which couple each bath site with every other one via the Γ . For simplicity, we don't depict all terms for this "full" case. For the other (sparse) cases all couplings are drawn, and n.n. denotes nearest neighbor terms in Γ (1/2) . question is how fast. In Sec. III we want to elaborate on this point in detail and present results obtained with those geometries, which we briefly discuss and motivate here.
In all cases one can restrict the geometries to a sparse (e.g. tridiagonal) and real-valued matrix E. As commonly true for impurity problems, the physics on the impurity site is invariant under unitary transformations among bath sites only. For an arbitrary unitary tranformation U with U if = U f i = δ if to new fermionic operators, one obtains an analogous auxiliary system with modified bath parameters E = U † EU , Γ (1) = U † Γ (1) U and Γ (2) = U † Γ (2) U . It is easy to check that the f f -component of the Green's functions Eqs. (20) and (22) is not affected by this transformation. Therefore, we choose without loss of generality E to be sparse as well as real, and for Γ (1/2) in the most general case dense matrices with O(N 2 B ) parameters. The particular form of E is irrelevant, i.e. whether it is diagonal for bath sites (star) or tridiagonal (chain), as long as the Γ (1/2) matrices are transformed accordingly.
Such a general geometry with sparse E and dense Γ (1/2) is referred to as "full" setup in the following.
Here, we will further distinguish between the case in which the Γ (1/2) are real or they have complex elements ("full complex"). In addition we consider the four sparse cases "2 chains n.n.", "2 chains onsite", "star", and "1 chain n.n.", in which also the Γ (1/2) are sparse. The meaning of these abbreviations is given in Fig. 1, see also App. B. These sparse geometries are however not linked to each other by unitary transformations and represent inequivalent subsets of the "full" setup. Which one of these is advantageous in practice is not obvious a priori, and discussed in the next section. 73 The "full" geometry comprises all other ones and thus, obviously, gives the best possible fit for a given N B . In addition, one can allow for the off diagonal matrix elements of the Γ (1/2) to be complex, thus extending the set of fit parameters. Nevertheless, the sparse setups may be of great value for sophisticated manybody solution strategies for the interacting Lindblad equation, such as MPS. We made use of the "full" setup (with real parameters) in the ED treatment Ref. 33, which is applicable to dense Γ (1) and Γ (2) matrices, and could consider up to N B = 6. Larger systems are prohibitive due to the exponentially increasing Hilbert space. 68 In the recent MPS implementation Ref. 39, on the contrary, we could consider as many as N B = 16 bath sites. However, in favour of the applicability of MPS methods one should avoid long-ranged hoppings and we thus employed the "2 chains n.n." geometry. As becomes evident also from the results below, the gain in N B hereby outweighs the restriction of the fit setup, so that the MPS approach is clearly superior. Also the other sparse setups investigated below are possible candidates for MPS, see also Ref. 57. Besides this, approaches such as the above men-tioned buffer zone scheme and variations of it, 30-32 which are often applied concepts in Lindblad-type representation of noninteracting environments, are related to the "star" geometry, see also the discussion below.

III. RESULTS
As discussed above, while the "full" geometry is the most efficient one, for the purpose of employing efficient many-body eigenvalue solvers such as MPS, it is of great relevance to consider setups which feature only sparse E, Γ (1) and Γ (2) matrices. Furthermore, it is also of general interest to investigate the importance of long-range terms in the Γ (1/2) -matrices, and why they are crucial in order to improve the fit. These are the questions that are addressed in this section. Moreover, we will analyze the rate of convergence as a function of N B for the different setups shown in Fig. 1, and for different temperatures of the physical system. The detailed knowledge of the convergence properties is important in order to be able to estimate whether certain systems can be accurately treated or not.
We consider a physical system consisting of an impurity site coupled to two reservoirs (leads) at different chemical potentials, corresponding to a bias voltage φ across the impurity, and with a flat density of states as plotted in Figs. (2-4). Typical results for a given φ and temperature T are shown in Figs. (2-4). For the different setups the quality of the fit is measured by the minimum of the cost function Eq. (26). As discussed above, the "full" setups give the best results. Already for a rather small number of bath sites N B 4, a good agreement between ∆ aux and ∆ ph is achieved, and the convergence is fast as a function of N B . Allowing for complex matrix elements produces a significant improvement. The accuracy obtained with N B = 8 for the real case is essentially achieved already with N B = 6 in the complex case (see also We now consider the sparse geometries. In contrast to the "full" setups, no improvement is obtained by allowing the matrix elements to be complex in this case. Among the sparse geometries, the ones with two chains are the most accurate. Both setups perform quite well. Again, a good agreement for small N B is obtained and a quick improvement shows up when increasing N B . "2 chains n.n." has off-diagonal Γ (1/2) -terms in contrast to "2 chains onsite", which leads to a faster convergence as seen e.g. for N B = 12. The "star" and most notably the "1 chain n.n." geometry are clearly worse. Both exhibit a rather poor convergence as a function of N B . For the "star" setup, this is due to the fact that the fitted auxiliary hybridization function consists of a sum of Lorentzian peaks. These enter in the Keldysh component with either positive or negative weights and can thus cancel each other. However, the rather broad Lorentzians with long 1/ω 2 tails make it apparently difficult to resolve the Fermi edges properly. The problem with slow convergence is most severe for the "1 chain n.n." geometry. Here, the single chain is clearly inadequate to represent at the same time the desired density of states and the sudden changes in the occupation number, see also the discussion below. While the Keldysh component is roughly reproduced, this comes at the price of large oscillations in the retarded one. In addition, the improvements with increasing N B are minor and the results for N B = 4 and N B = 12 are very close to each other.
The behavior just discussed is even more visible in the convergence study presented in Fig. 5 and Fig. 6. In Fig. 5 the minimal values of the cost function χ, Eq. (26), for various values of N B and the different setups are shown. Four different temperatures and each of them with φ = 0 and φ = 3 Γ are considered. As expected, the "full complex" setup gives the lowest values of χ in all cases, and, moreover, the fastest rate of convergence as a function of N B . The "full" setup, without complex terms also performs quite well. The sparse geometries "2 chains n.n." and "2 chains onsite" perform not as well, which is not surprising since only restricted subsets of the full available fit parameters are used in this case. Nevertheless, these setups achieve a rather high rate of convergence. This shows that of all possible geometries, "2 chains" ones apparently contain the most relevant contributions. In most cases studied here, the off-diagonal Γ (1/2) -terms in "2 chains n.n." result in a significant improvement compared to "2 chains onsite", which is the reason why we favored the former in our MPS many-body calculations performed in Ref. 39. In that work we found that an accuracy of at least χ ≈ 10 −2 was necessary in order to properly account for Kondo physics. This could be reached already for N B ≈ 12.
We now discuss the "star" setup. In order to present a fair comparison with the other geometries we optimize all available parameters within this geometry, namely all . In this way we obtain an exponential convergence as for the other setups, although with a significantly smaller rate. One should note that in standard buffer zone approaches 30-32 an equidistant energy spacing ∆ i ≈ 2D/N B with equal onsite Γ (1/2) -terms is often assumed for the bath sites. Clearly, such a discretization approach cannot converge exponentially and it is only first-order accurate in the spacing ∆ i . Therefore, the value of the cost function presented here for the "star" setup can be seen as a lower bound for the buffer zone approach. Despite of the exponential conver-"full" "full complex" FIG. 2: Fit to the bath hybridization functions for the "full" setups (real and complex) (see Fig. 1). The physical ∆ ph (ω) (black lines) describes a reservoir with a flat density of states with hybridization strength Γ and a half bandwidth of D = 10 Γ which is smeared at the edges. An applied bias voltage φ = 3 Γ shifts the chemical potentials of the two reservoirs (leads) anti-symmetrically and a temperature of T = 0.1 Γ is considered here. "2 chains onsite" "2 chains n.n." Fig. 2 for the "two-chains n.n." and "two-chains onsite" setups. . In order to resolve the scaling with temperature more reliably, we exclude the two data points with the smallest NB from each of the linear fits, which have not enough structures to resolve low-energy scales. Dotted lines represent results of linear fits in these semi-logarithmic plots. The temperature dependence of the convergence rates (as a function of NB) obtained in this way are illustrated in Fig. 7.

FIG. 3: Same as
FIG. 7: Estimated convergence rates obtained from the data in Fig. 5 plotted as a function of temperature. The rates for the sparse setups are obtained by assuming χ ∝ exp[−r(T )NB]. For the "full" setups, the exponent is quadratic in NB, therefore we plotted the differential rate, defined as − d log χ dN B evaluated at NB = 6. gence of the "star" geometry, it becomes apparent from Fig. 5 that a very slow rate of convergence is achieved. To reach an accuracy χ ≈ 10 −2 for the case T = 0.05 Γ and φ = 0 for instance, much larger auxiliary systems with N B ≈ 40 would be needed. For the MPS-solver used in Ref. 39 such large auxiliary systems are clearly out of reach. Therefore, the present analysis clearly demonstrates the huge advantage of optimizing the bath parameters of the auxiliary system, and furthermore, of choosing an appropriate geometry when considering only a restricted subset of the "full" setup.
Let us now turn to the results for the "1 chain n.n." setup in Fig. 5. Despite of the poor performance and the strongly limited practical use, the observed behavior is interesting from a fundamental point of view. As becomes evident from the results, a single chain with local dissipators is a particularly bad choice in order to represent a partially filled bath. The convergence is very slow and an extremely long chain would be needed in order to achieve results comparable to the other geometries. As shown above, a drastic improvement is obtained when using two chains instead. This would be more or less intuitive for the nonequilibrium case, in which the physical system also consists of two baths. However, the advantage of the "2 chains" geometry over the "1 chain" case is even more pronounced in the equilibrium case (see Φ = 0). Another important observation to better understand this is the following: In Ref. 39 we found nearly identical accuracies when considering the "2 chains" geometry as used here, or a filled/empty restriction of it. In the latter case one chain has the purpose of representing the filled spectrum and the other chain the empty spectrum of the physical hybridization function 74 , and not necessarily the two physical reservoirs. This shows that a single chain of small size is very well-suited to reproduce a certain density of states, but not simultaneously a Fermi edge or other sharp changes in the occupation number. Furthermore, a "2 chains" filled/empty setup seems to be a rather natural representation which contains the most relevant subset of the "full" geometry. Here, the resolution of sharp features in ∆ ph (ω), which either correspond to band edges or to sudden occupation changes at the Fermi edges, are resolved by appropriate Hermitian couplings E and correspondings broadenings/couplings stemming from Γ (1/2) . In this way, the filled and empty chain together can well reproduce sharp features in ∆ R ph (ω) and ∆ K ph (ω). 75 Additionally to the convergence as a function of N B we depict in Fig. 6 the cost function versus the number of available fit parameters C(N B ). As can be seen, the trends in the semi-logarithmic plot are well described by straight lines in all cases, which clearly shows the achieved exponential convergence with respect to C(N B ). For the sparse setups this means that χ ∝ exp[−O(N B )] whereas for the "full" setups even χ ∝ exp[−O(N 2 B )]. Due to this, the "full" geometries converge much quicker, as observed in the results above. With respect to the number of fit parameters, however, the "2 chains" se-tups perform best. Again, this signifies that these setups contain the most relevant subset of all possible fit parameters.
Another important aspect is the dependence of the convergence rate r(T ) on temperature. The estimated rates r(T ) for each setup are depicted in Fig. 7. Of course, the superior scaling of the "full" and the "2 chains" setups is also apparent in the magnitude of r(T ). Furthermore, in all cases one observes the trend that the higher the temperature the faster the convergence. This can be understood from the fact that at high T the Keldysh component ∆ K ph (ω) is weakly ω-dependent so that less bath sites are necessary for a reliable fit. Eventually, in the T → ∞ and wide-band limit the Markov approximation becomes even exact. In the other extreme T → 0 limit, discontinuous functions are present in ∆ K ph (ω), produced by the abrupt Fermi edges. However, each of the frequency dependent functions in the effective set given by Eqs. (20)(21)(22)(23)(24)(25) is continuous. Therefore, T → 0 can only be reproduced in the limit N B → ∞. This explains the observed trend that, for a given N B , the high-temperature regime is generally better represented than the low-temperature one. Furthermore, a nonzero φ tends to result in larger values for the cost function, see also Fig. 5. 76 A. Discussion of further aspects The present approach is equally suitable to describe a system in equilibrium as well as out of equilibrium. In the first case it becomes competitive with conventional ED-and MPS-impurity solvers for DMFT based on a bath without Lindblad terms. The distinction between the equilibrium or nonequilibrium situation shows up in the properties of the bath hybridisation function ∆ ph (ω). In the equilibrium case its Keldysh and retarded part will fullfill the fluctuation-dissipation theorem. Interestingly, the equilibrium problem is mapped onto an auxiliary nonequilibrium one, since a current will typically flow from Γ (2) to Γ (1) dissipators. An example is the case discussed in Sec. III of a two-chain geometry with a completely empty and a full one. Such geometry can be used to describe an equilibrium situation at the impurity as well, as long as ∆ K aux (ω) and ∆ R aux (ω) are chosen to fullfill the fluctuation-dissipation theorem. Nevertheless, a current will flow from one chain to the other across the impurity, which, however, will be in equilibrium. Notice that since the mapping will be approximate for finite N B , there will be small deviations from the fluctuationdissipation theorem.
Another interesting aspect is the role of chemical potential(s) µ α and temperature(s) T α of the different physical reservoirs. These determine only indirectly the values of the parameters of the auxiliary system E, Γ (1) , Γ (2) . More specifically µ α and T α first determine ∆ K ph (ω) via Eq. (11) and Eq. (12). In a second step, via the requirement Eq. (19) and the corresponding fit procedure, they finally determine the auxiliary parameters. For low temperatures, the chemical potentials µ α will then appear as sharp changes in ∆ K ph (ω). This is in contrast to more direct approaches, such as buffer-zone based ones, in which the parameters are directly determined form µ α and T α , in equations such as Eq. (3) for each bath level. Both methods have their advantages: Direct approaches can be more convenient, for example in NRG 54 . On the other hand, a fit procedure like the present one produces a much faster, exponential, convergence.

IV. SUMMARY AND CONCLUSIONS
In this work, a scheme for mapping the hybridization function of correlated quantum impurity problems with non-Markovian fermionic reservoirs onto an auxiliary open quantum system was developed and presented in a general framework and discussed in detail. The approach as outlined here can be used to model transport through interacting impurities, Hubbard chains or small clusters and molecules. The key aspect is to replace the infinite fermionic reservoirs of the original problem by a combination of a small number N B of bath levels plus Markovian terms. By this we arrive at a finite open quantum system described by a Lindblad equation, whose manybody problem can be solved with high accuracy by numerical techniques. However, despite of the Markovian environments for the bath levels the thereby approximated hybridization function is clearly non-Markovian at the impurity site in the sense that it has a frequency dependence, which is a consequence of the memory effects of the environment. While this idea is not new, the key point of our work is the formulation of an optimization procedure in order to determine the parameters of the auxiliary bath levels. This allows us to achieve an exponential convergence of the mapping, as clearly demonstrated in this work.
In the mapping one has certain degrees of freedom and different geometries for the auxiliary system are possible. In this work we discussed a variety of choices in detail and compared their performance.
When using Krylov-based many-body approaches it is convenient to take advantage of as many fit parameters as possible. For these cases the "full" setups are the best choice. Here we also showed that further allowing for complex matrix elements ("full complex") drastically improved the accuracy of the mapping with respect to the plain real "full" setup, which we used in Ref. 33. With efficient manybody solution techniques, such as MPS, in mind, it is of advantage to restrict the auxiliary quantum system to a sparse form. For this we analyzed four sparse setups. The results revealed the most relevant degrees of freedom in the auxiliary system and demonstrated clearly that the performance of different sparse setups may differ by orders of magnitude. In particular, the well-known "star" geometry turned out to exhibit a very slow rate of convergence when increasing N B , and also a geometry with "1 chain" and local Lindblad drivings performed much worse than the other cases. In contrast, setups with "2 chains" and local Lindblad drivings yielded very good results, with an accuracy orders of magnitude below the other two sparse cases. With this knowledge it is possible, on the one hand, to employ efficient sparse setups which yield already for modest values of N B very accurate results, and on the other hand, one can better understand the underlying mechanisms of the mapping. Together with the findings mentioned in Ref. 39, we can conclude that a so-called filled/empty geometry with "2 chains" is essentially a natural representation of a non-Markovian reservoir by auxiliary Lindblad levels. In this geometry one chain has the purpose of reproducing the filled spectrum of the original reservoir whereas the other chain the empty spectrum. This is achieved in each chain separately by an optimal combination of hoppings between the bath levels and couplings to one Markovian environment, which is then either completely filled or empty. By this separation it is possible to resolve sharp features in the original hybridization function in great detail, which may correspond to sudden occupation changes at the Fermi edges or band borders. A single chain coupled to filled and empty Markovian environments, on the contrary, cannot simultaneously represent a particular density of states and a partially filled spectrum appropriately, as evident from the "1 chain" setup.
Besides comparing different auxiliary setups to each other, we also analyzed the general convergence properties in detail. As mentioned above, a clear exponential convergence was found in all cases, which can be accounted to the optimization strategy for the bath parameters. Furthermore, a generic set of equilibrium and nonequilibrium reservoirs with various temperatures was chosen for the original system. From this we found the expected trend for all setups that the high-temperature regime is better represented by the auxiliary system than the low-temperature one, i.e. the rate of convergence of the mapping increases with temperature. Therefore, to achieve a given accuracy it is more challenging to resolve low temperatures and larger auxiliary systems must be considered. The plain exponential convergence shown here yields a simple tool to extrapolate results for low N B to higher values, and by this to judge the feasability of treating certain physical situations.
While we did not discuss this in the present work, it would be probably useful to exploit the freedom in the choice of the cost function Eq. (26), and in particular, of the corresponding weight function W (ω). For example, one could increase the weight around the chemical potential in order, possibly, to achieve a better resolution at low energies. An even more appropriate strategy in this sense would be to combine the present approach to NRG ideas such as the logarithmic discretisation, work along these lines is in progress, see also Ref. 54. Eq. (26) strongly disfavours delta-function-like or strongly oscillating spectra. In some cases, such oscillations may be unimportant, especially if they occur at high energies. Therefore, it would be useful to adopt a cost function which does not penalize them. This could be achieved, for example, by convoluting ∆ ph (ω) and ∆ aux (ω) with a suitable "smoothing" function before evaluating their difference Eq. (26). Alternatively, one could use cost functions based on differences in spectral moments up to a certain order.
In this work we considered the simplest case of a single impurity Anderson model (SIAM) in order to focus on the mapping itself. To treat more extended cluster or multi-level problems essentially the same approach can be used. In the simple case of diagonal cluster hybridization functions exactly the same equations are applicable to model each reservoir separately by auxiliary Lindblad levels. But also non-diagonal hybridization functions can be treated, of course. For this purpose the approach was presented here in a more general framework before focusing particularly on the SIAM.
Besides the technical aspects, we believe that the presented study contains relevant information to the general question of the representability of non-Markovian fermionic reservoirs by open quantum systems, and in particular by Lindblad-type equations. We expect that the insights gained in this work may contribute also to other closely related fields on Markovian and non-Markovian quantum master equations.
In this section we provide detailed informations for readers interested in an actual implementation. Furthermore, a working code is available on request. To obtain it just contact us via e-mail. Much of the information below is contained in standard textbooks and reviews. However, for completeness we outline here the standard algorithm in detail and point out choices we made, which turned out to be convenient for the specific problem.
As stated above, a single evaluation of Eqs. (20-25) is rather cheap numericaly since it involves only one matrix inversion and multiplications of matrix of size (N B + 1) × (N B + 1) . Thus, the increase in computation time with N B is rather moderate. However, the multi-dimensional optimization problem itself is demanding and it strongly depends on the particular behavior of χ(x) when varying the set of parameters x. In the worst case scenario, when χ(x) is a rough potential landscape with many local minima and short-scaled variations, one could imagine that it becomes necessary to nearly explore the whole parameter space. However, x is a continuous vector and even when assuming a fixed number of discrete values for each component in x, one faces a number of points in parameter space that grows exponentially with dim(x). In the other extreme, for the case that χ(x) is quadratic in x it is well-known that a conjugate gradient scheme leads to the exact minimum in dim(x) iterations. What we found in practice, when performing the minimization within AMEA, 33,39 is that we have an intermediate situation which exhibits local minima, but gradient-based methods still work fine especially for smaller values of N B . In the first work on the EDsolver, Ref. 33, we employed a quasi-Newton line search with many random starting points. This is particularly useful for N B < 6. However, the necessary number of starting points increases rapidly with N B . Therefore, in the course of the work on the MPS-solver, Ref. 39, we looked for more efficient solution strategies. In the end we implemented a parallel tempering (PT) approach with feedback optimization, which is a Monte Carlo scheme that is able to overcome local minima. We describe it in the following in more detail. In this way, the minimization problem for the ED-solver with N B = 6 and for the MPS-solver with up to N B = 16 can be solved in reasonable time. This amounts to minimizing in a space of ≈ 30 − 60 parameters in both geometries, depending on whether one has particle-hole symmetry or not. 77

Markov chain Monte Carlo
The PT algorithm is outlined in detail below. For completeness, let us first briefly recap the basic ideas of the underlying Markov chain Monte Carlo (MCMC), and of the related simulated annealing algorithm.
MCMC techniques were originally developed to evaluate thermodynamic properties of classical systems which exhibit a very large phase space where simple sampling strategies fail. For our purposes here, we are interested in minimizing the cost function χ(x) as defined in Eq. (26) with respect to the parameter vector x. For such high-dimensional minimization problems one can adapt MCMC schemes by viewing χ(x) as an artificial energy and by introducing an artificial inverse temperature β. In the so-called simulated annealing one samples from the Boltzmann distribution p(x) = 1/Z exp (−χ(x)β) at a certain β, and then successively cools down the artificial temperature. Motivated by the behavior of true physical systems one expects to end up in the low-energy state when letting the system equilibrate and when cooling sufficiently slowly. Analogous to thermodynamics one can calculate the specific heat C H = β 2 ∆χ(x) 2 and by this locate regions with large changes, i.e. phase transitions, where a slow cooling is critical. However, in practice it may be time consuming to realize the equilibration and sufficiently slow cooling, and for tests within AMEA we often ended up in local minima. In order to obtain a robust algorithm, which can also start from previous solutions as needed for instance within DMFT, we sought for a method which is able to efficiently overcome local minima and still systematically targets the low-energy states. For this a multicanonical and a PT algorithm were tested, whereby the latter turned out to be more convenient. In the following we briefly outline the PT scheme used within AMEA, and refer to Ref. 58-62 for a thorough introduction to MCMC, simulated annealing, multicanonical sampling and PT.
As just stated, in a MCMC scheme one typically samples from the Boltzmann distribution p(x) = 1/Z exp (−χ(x)β) at some chosen inverse temperature β. This is done through an iteratively created chain of states {x l }, whereby one avoids the explicit calculation of the partition function Z. An effective and well-known scheme for this is the Metropolis-Hastings algorithm 58,59 . One starts out with some state x l and proposes a new configuration x k , whereby it has to be ensured that every state of the system can be reached in order to achieve ergodicity. The proposed state x k is accepted with probability 7858,59 p l,k pacc. = min 1, (A1) If the proposed configuration is accepted, then the next element x l+1 in the chain is x k , otherwise x l again. From Eq. (A1) it is obvious that p l,k pacc. = 1 when p(x k ) > p(x l ), so that an importance sampling towards regions where p(x) is large is achieved. One can show that the algorithm fulfills detailed balance and draws a set of samples {x l } that follow the desired distribution p(x). However, stemming from the iterative construction, correlations in the chain are present which require a careful analysis for the purpose of statistical physics 58,59 . For optimization problems, on the other hand, the situation is much simpler and one is just interested in the element in {x l } which minimizes χ(x). Since a proposed step with χ(x k ) < χ(x l ) is always accepted the algorithm targets minima, however, also uphill moves in configuration space are allowed with a probability depending exponentially on the barrier height ∆χ k,l = χ(x k ) − χ(x l ) and β. Effectively, uphill moves take only place when ∆χ k,l β O(1). For small values of β large moves in configuration space with large ∆χ k,l are likely to be accepted, whereas for large β the distribution p(x) is peaked at minima in χ(x), so that especially those regions are sampled. For the latter case configurations in the chain {x l } are generally more correlated and once a x l corresponds to a local minimum the algorithm may stay there for very long times.
One has great freedom in defining a proposal distribution from which the new state x k is drawn given the current configuration x l . 79 Common choices are, for instance, a Gaussian or a Lorentzian distribution with the vector difference x k − x l as argument. We favored the former and updated each component i with a probability according to 58 Hereby, a different step size σ i for each component is expedient since the potential landscape χ(x) around x l is typically highly anisotropic. Ideally, one should make use of the covariance matrix Σ l of χ(x l ) and consider as argument for the Gaussian instead ( 58 . However, we encountered the problem that the estimation of the covariance matrix at run time was strongly affected by noise and thus not feasible. The adjustment of the step sizes σ i , on the contrary, can be done after a short number of updates by demanding that a value of p l,k pacc. ≈ 0.5 is reached on average when modifying the component i. For this we implemented a check at every single proposal, that increases σ i → 1.1 σ i when p l,k pacc. > 0.6 and decreases σ i → 0.9 σ i when p l,k pacc. < 0.4. Analogous to the treatment of spin systems, we define one sweep as a single update of all the components of x. 80

Parallel tempering
In a PT algorithm one considers instead of sampling at one certain temperature a set of different temperatures β −1 m and corresponding replicas x m l , each of which is evolved through a Markov chain. The largest β m thereby target local minima whereas low β m values allow for large moves in configuration space. The key idea of the PT approach is to let the individual replicas evolve dynamically in the set of β m . By this one achieves that a replica at high β m values systematically targets local minima but can overcome potential barriers again when its inverse temperature is changed to lower values. As a result, the time scales to reach an absolute minimum are drastically reduced and an efficient sampling of the low-energy states is achieved. For the purpose of calculating thermodynamic properties one usually chooses a Metropolis-Hastings probability to swap two replicas with adjacent temperatures 61,62 p m,m+1 swap,l = min 1, with the Boltzmann distribution for each β m given by p m (x) = 1/Z m exp (−χ(x)β m ). Such swap moves are conveniently proposed after a certain number of sweeps, which satisfies the sufficient condition of balance for thermodynamics 62 . In practice, we chose 10 sweeps before swapping replicas. For the exchange to effectively take place the underlying requirement is that the adjacent β m and β m+1 values are close enough to each other, so that the two energy distributions Ω[χ(x)]p m (x) and Ω[χ(x)]p m+1 (x) overlap, with Ω[χ 0 ] = dx δ(χ 0 − χ(x)) the density of states of the cost function. This means that a replica at one temperature must represent a likely configuration for the neighboring temperature 62,63 . In order to achieve this, a crucial point in the PT algorithm is to adjust the distribution of the inverse temperatures properly to the considered situation. Various criteria for this have been devised, see e.g. Ref. 62. A common choice is to demand that the swapping probability Eq. (A3) becomes constant as a function of temperature 64,65 , and in Ref. 66 a feedback strategy was presented which optimizes the round trip times of replicas. We tested the latter within AMEA but favored the simpler former criterion in the end, since it allows for a rapid feedback and quick adjustment to large changes in χ(x m l ). In the simple situation of a constant specific heat C H with respect to energy χ for instance, an optimal strategy is known since a geometric progression β m /β m+1 = const. of temperatures yields a constant swapping probability 62,63 . For interesting cases in practice this is rarely fulfilled, but within AMEA it served as a good starting point. The set of inverse temperatures is then optimized by averaging p m,m+1 swap,l over a couple of swappings to obtain the mean probabilitȳ p m,m+1 swap and adjusting the β m thereafter. For this we chose a fixed lowest and highest β m value and changed the spacings in between according to with ∆β m = β m+1 − β m and c adjusted properly so that max(β m )−min(β m ) = max(β m )−min(β m ). In the works by Ref. 64,65 it was shown that a constant swapping probability of 20% − 23% seems to be optimal. We determined the highest and lowest β m values by the changes in χ(x) we want to resolve or allow for, and the number of inverse temperatures β m was then set accordingly in order to roughly obtainp m,m+1 swap ≈ 0.25. Fixing the smallest and largest β m is, for our purposes, the most convenient choice among the many possibilities.
However, despite of the feedback optimization of temperatures as just described above, we encountered in practice the unwanted behavior that the set of parallel replicas effectively decoupled into several clusters. In order to suppress this we found it advantageous to introduce the following simple modification to Eq. (A3) 81 p m,m+1 swap,l = max p m,m+1 swap,l , p th.
with a certain threshold probability p th. swap , e.g. p th. swap = 0.1 or 0.05. In this way one avoids that the β m are shifted unnecessarily close to each other and avoids very long time scales, in which replicas oscillate only between two neighboring inverse temperatures. For the sake of clarity we present here for the different setups of Fig. 1 the form of the (hermitian) matrices E and Γ (1) for the case N B = 4 in the particle-hole symmetric case, i.e. under the constraint (27) which also fixes Γ (2) . In addition, we quote the number of available fit parameters C(N B ) for each setup. The fit parameters are denoted below as x i for i = 1, C(N B ), with the only constraint that Γ (i) should be semipositive definite. This, together with the requirement that ∆ aux vanishes for ω → ∞ further requires Γ In the first four setups, the impurity is in the center (i = 3). In the "1 chain n.n." it is on the first site (i = 1).
"full" geometry The parameters x 9 to x 14 can be complex. Therefore, it is straightforward to see that, for general N B the number of independent (real) parameters is C(N B ) = N B 2 * (N B +3) for the real case and C(N B ) = N B * (N B + 1) for the complex case. Appendix C: Reduction of bath to a "star" form In principle, one can represent a noninteracting dissipative bath consisting of N B sites (i = 1, . . . N B ) coupled to an impurity (i = f , we take f = 0) by specifying the single-particle parameters E ij , Γ (1) ij , and Γ (2) ij (i, j = 0, . . . N B ), with corresponding hermitian, and in the case of Γ (1) , Γ (2) semipositive definite matrices. We show here that for the sake of fitting the retarded component of a given bath spectral function ∆ R aux , these parameters are redundant.
We are interested in G R aux , which is the 00 component of G R . By a well known result of matrix inversion, this is given by which identifies ∆ R aux = T T (ω − F ) −1 T + δF 0 , where δF 0 ≡ F 0 − ε f , which, for simplicity, we set to zero. The first term can be rewritten by introducing the matrix V which diagonalizes F , 83 i.e.
This gives We can thus replace in Eq. (C1) F with a diagonal, complex matrix F diag and T ( T ) with T ( T ), and we get (C4) Here, G R has the same 00 element as G R from Eq. (C1), i.e. the same G R aux and ∆ R aux . In this way, by the requirement that E and Γ (+) must be hermitian, we can construct new E = (F +F † )/2 and Γ (+) = (F −F † )/(2i), i.e. a new auxiliary system which yield the same ∆ R aux and have the "star" geometry (cf. 1). 84 This means that, concerning the retarded part, one can restrict to the case of diagonal bath energies and Γ (+) , i.e., as in the nondissipative case, ∆ R aux is fixed by only O(N B ) independent bath parameters, the rest being redundant. This is also the case when the bath hybridisation function is represented by a completely empty and a completely full chain, as discussed in Sec. III, since in that case one simply fits the retarded components of the two chains sep-arately. On the other hand, for the most generic case, Γ (1) and Γ (2) will not commute and cannot be simultaneously diagonalized, so that the Keldysh component ∆ K aux (ω) appears to still depend on O(N 2 B ) bath parameters (cf. Fig. 5). Further investigations should be carried out in order to clarify this issue.