Auxiliary master equation approach to non-equilibrium correlated impurities

We present a numerical method for the study of correlated quantum impurity problems out of equilibrium, which is particularly suited to address steady state properties within Dynamical Mean Field Theory. The approach, recently introduced in [Arrigoni et al., Phys. Rev. Lett. 110, 086403 (2013)], is based upon a mapping of the original impurity problem onto an auxiliary open quantum system, consisting of the interacting impurity coupled to bath sites as well as to a Markovian environment. The dynamics of the auxiliary system is governed by a Lindblad master equation whose parameters are used to optimize the mapping. The accuracy of the results can be readily estimated and systematically improved by increasing the number of auxiliary bath sites, or by introducing a linear correction. Here, we focus on a detailed discussion of the proposed approach including technical remarks. To solve for the Green's functions of the auxiliary impurity problem, a non-hermitian Lanczos diagonalization is applied. As a benchmark, results for the steady state current-voltage characteristics of the single impurity Anderson model are presented. Furthermore, the bias dependence of the single particle spectral function and the splitting of the Kondo resonance are discussed. In its present form the method is fast, efficient and features a controlled accuracy.


I. INTRODUCTION
Correlated systems out of equilibrium have recently attracted increasing interest due to the significant progress in a number of related experimental fields. Advances in microscopic control and manipulation of quantum mechanical many-body systems within quantum optics 1 and ultra cold quantum gases, for example in optical lattices, 2-6 have long reached high accuracy and versatility. Ultrafast laser spectroscopy 7,8 offers the possibility to explore and understand electronic dynamics in unprecedented detail. Experiments in condensed matter nanotechnology, 9 spintronics, 10 molecular junctions [11][12][13][14][15][16] and quantum wires or quantum dots, 17,18 are able to reveal effects of the interference of few microscopic quantum states. The non-equilibrium nature of such experiments does not only offer a new route to explore fundamental aspects of quantum physics, such as non-equilibrium quantum phase transitions, 19 , the interplay between quantum entanglement, dissipation and decoherence 20 , or the pathway to thermalization, 21,22 , but also suggests the possibility of exciting future applications. 11,23 Addressing the dynamics of correlated quantum systems poses a major challenge to theoretical endeavors. In this respect, quantum impurity models help improving our understanding of fermionic many-body systems.
In particular, the single-impurity Anderson model (SIAM), 24 which was originally devised to study magnetic impurities in metallic hosts, 25,26 has become an important tool in many areas of condensed matter physics. 27,28 Most prominently, it features nonperturbative many-body physics which manifest in the Kondo effect. 29 It provides the backbone for all calculations within dynamical mean field theory (DMFT), 28,30 a technique which allows to understand the properties of a broad range of correlated systems and becomes exact in the limit of infinite dimensions. 31 The basic physical properties of the SIAM in equilibrium are quite well understood 29 thanks to the pioneering work from Kondo,32 renormalization group 33 as well as perturbation theory (PT) [34][35][36][37] and the mapping to its low energy realization, the Kondo model. 38 The SIAM out of equilibrium provides a description for several physical processes such as, for example, nonlinear transport through quantum dots, 17,39 correlated molecules 13,14,[40][41][42] or the influence of adsorbed atoms on surfaces or bulk transport. 43 As in the equilibrium case, the solution of the SIAM constitutes the bottleneck of non-equilibrium DMFT [44][45][46][47][48][49][50][51] calculations. Therefore, accurate and efficient methods to obtain dynamical correlation functions of impurity models out of equilibrium are required in order to describe time resolved experiments on strongly correlated compounds 7,8 ant to understand their steady state transport characteristics. 23 However, nonequilibrium correlated impurity models still pose an exciting challenge to theory. Our work addresses this issue with special emphasis on the steady state. But before introducing the present work in Sec. I A, we briefly review previous approaches. In recent times a number of computational techniques have been devised to handle the SIAM out of equilibrium. Among them are scattering-state BA, 52 scattering-state NRG (SNRG), [53][54][55] non-crossing approximation studies, 56,57 fourth order Keldysh PT, 58 other perturbative methods 59,60 in combination with the renormalization group (RG), [61][62][63][64][65] iterative summation of real-time path integrals, 66 time dependent NRG, 67 flow equation techniques, 68,69 the time dependent density matrix RG (DMRG) [70][71][72][73][74][75] applied to the SIAM, 76,77 non-equilibrium cluster PT (CPT), 78 the non-equilibrium variational cluster approach (VCA), 79,80 dual fermions, 81 the func-tional RG (fRG), 82,83 diagrammatic QMC, 84,85 continuous time QMC (CT-QMC) calculations on an auxiliary system with an imaginary bias, [86][87][88][89][90] super operator techniques, 91,92 many-body PT and time-dependent density functional theory, 93 generalized slave-boson methods 94 , real-time RG (rtRG), 95 , time dependent Gutzwiller mean-field calculations 96 and generalized master equation approaches. 97 Comparisons of the results of some of these methods are available in literature 77,98,99 and time scales have been discussed in Ref. 100.
Despite this large number of approaches, only a limited number of them is applicable to non-equilibrium DMFT, and very few are still accurate for large times in steady state. Beyond the quadratic action for the Falicoff Kimball model, 46,101,102 iterated PT (IPT), 45 numerical renormalization group 48 , real time QMC, 48,103 the noncrossing approximation (NCA) 104,105 and recently Hamiltonian based impurity solvers 106 have been applied in the time dependent case. Some of the above approaches, such as QMC 49 and DMRG 73 are very accurate in addressing the short and medium-time dynamics, but in some cases the accuracy decreases at long times and a steady state cannot be reliably identified. Some other methods are perturbative and/or valid only in certain parameter regions or for restricted models. RG approaches (e.g., Ref. 61) are certainly more appropriate to identify the low-energy behavior.

A. Present work
In this paper we discuss a method, first proposed in Ref. 51, which addresses the correlated impurity problem out of equilibrium, and is particularly efficient for the steady state. The accuracy of the results is controlled as it can be directly estimated by analyzing the bath hybridization function (details below). Here, we extend, test and provide details of this approach and its implementation. The basic idea is to map the impurity problem onto an auxiliary open system, consisting of a small number of bath sites coupled to the interacting impurity and additionally, to a so-called Markovian environment. 107 The parameters of this auxiliary open quantum system are obtained by optimization in order to represent the original impurity problem as accurately as possible. The auxiliary system dynamics are governed by a Lindblad Master equation which is solved exactly with the non-hermitian Lanczos method. The crucial point is, that the overall accuracy of the method is thus solely determined by how well the auxiliary system reproduces the original one. This can be, in principle, improved by increasing the number of auxiliary bath sites.
In the present study we provide convincing benchmarks for the steady state properties of the SIAM coupled to two metallic leads under bias voltage. We include a discussion of convergence as a function of the number of bath sites and present a scheme to estimate the error and partially correct for it. In its presented form the method is fast, efficient and is directly applicable to steady state dynamical mean field theory 51 for which previously suggested methods are less reliable. Extending the method to treat time dependent properties and multi-orbital systems is possible, in principle, however with a much heavier computational effort.
The paper is organized as follows: In Sec. II A the SIAM under bias voltage is introduced. In Sec. II B we introduce non-equilibrium Green's functions and in Sec. II C and II D we outline the auxiliary master equation approach where we also focus on details of our particular implementation. Results for the steady state, including the equilibrium situation are presented in Sec. III. This includes the steady state current-voltage characteristics which we compare with exact results from Matrix Product State (MPS) time evolution 77 as well as data for the spectral function under bias which we compare with non-equilibrium NRG. 54 We conclude and give an outlook in Sec. IV.

II. AUXILIARY MASTER EQUATION APPROACH
As discussed above, the method is particularly suited to deal with non-equilibrium steady state properties caused by different temperatures and/or chemical potential in the leads of a correlated quantum impurity system. As such, it can be readily used as impurity solver for non-equilibrium DMFT. 46,51 Here we illustrate its application to the fermionic SIAM with two leads having different chemical potentials, and, in principle, different temperatures.
A. Non-equilibrium single impurity Anderson model We consider a single Anderson impurity coupled to electronic leads under bias voltage (see Fig. 1

(a))
The impurity orbital features charge as well as spin degrees of freedom and is subject to a local Coulomb repulsion UĤ Here f † σ /f σ denote fermionic creation/annihilation operators for the impurity orbital with spin σ ∈ {↑, ↓} respectively. The particle number operator is defined in the usual wayn f σ = f † σ f σ and the impurity on-site potential is f = (V G − U 2 ), with gate voltage V G = 0 at particle hole symmetry. The impurity is coupled to two noninteracting electronic leads λ ∈ {L, R} with dispersion λkĤ  1) consisting of an impurity with interaction U coupled via hybridizations t λ to noninteracting leads at chemical potential µ λ and temperature T λ , λ ∈ {L, R}. (b) Illustration of the auxiliary open quantum system Eq. (10a) with single particle parameters Eµν and Lindblad dissipators Γ κ µν consisting of the impurity at site f = 0, NB bath sites (NB = 4 in the plot) as well as a Markovian environment (shaded areas). When evaluating linear corrections (see Sec. C), an additional site NB + 1 is used.
The effect of a bias voltage φ is to shift the chemical potential and the on-site energies of the two leads by λ = ± φ 2 , respectively. For the energies λk of the leads we will consider two cases: (i) Two tight-binding semi-infinite chains with nearest-neighbor hopping t, corresponding to a semi-circular electronic density of states (DOS). In this case, the boundary retarded single particle Green's function of the two uncoupled leads is given by 108-110 with a bandwidth of D SC = 4 t. (ii) A constant DOS with a bandwidth D WB = π t, resulting in boundary Green's functions 109 The choice D WB = π t makes sure that the DOS at ω = 0 of both lead types coincide. The leads are coupled to the impurity orbital bŷ where we take the same hybridization t λ = −0.3162 t for both leads, and N k → ∞ is the number of k points.
Expressions presented below are valid for arbitrary temperatures, although we will show results for zero temperature only, which is numerically the most unfavorable case. 111 The setup chosen here represents by no means a limitation of the method and extensions to more complicated situations, such as non-symmetric couplings, off particle-hole symmetry, etc. are straightforward.
B. Steady state non-equilibrium Green's functions We are interested in the steady state behavior under bias voltage of the model described by Eq. (1). We assume that such a steady state exists and is unique. 112 We denote the single particle Green's function of the impurity in the non-equilibrium Green's function (Keldysh) formalism by [113][114][115][116][117] Fourier transformation to energy ω is possible since in the steady state the system becomes time translationally invariant. In that case, the memory of the initial condition has been fully washed away, so there is no contribution from the Matsubara branch. 118 We will use an underline · · · to denote two-point functions with the Keldysh matrix structure as in Eq. (4). The Green's function of the correlated impurity can be expressed via Dyson's equation where Σ(ω) is the impurity self-energy. The noninteracting impurity Green's function G 0 (ω) can be written in the form g 0 (ω) being the noninteracting Green's function of the disconnected impurity, 108 and is the hybridization function of the leads (a 2×2 Keldysh object, in contrast to the equilibrium case, where it is convenient to work in Matsubara space). We define an equilibrium Anderson width 29 for each lead ∆ 0 ≡ − 1 2 m (∆ R (ω = 0)) = t 2 λ t ≈ 0.1 t. Below, we will use ∆ 0 as a unit of energy and in addition we choose = e = 1.
The boundary Green's functions g λ of each disconnected lead is determined by (a) its retarded component g R λ (either (2) or (3)), (b) its advanced component g A λ = g R * λ , and (c) its Keldysh component, which satisfies the fluctuation dissipation theorem since the disconnected leads are in equilibrium. Here, p F (ω − µ λ ) is the Fermi distribution with chemical potential µ λ . For the noninteracting isolated impurity one can take g −1 0 R = ω − f , and g −1 0 K = 0, since infinitesimals 0 + can be neglected after coupling to the leads (unless there are bound states). As usual, the presence of the interaction U makes the solution of the problem impurity plus leads a major challenge both in equilibrium as well as out of equilibrium, which we plan to address in the present paper.
Similarly to the equilibrium case, the action of the leads on the impurity is completely determined by the hybridization function ∆(ω), independently of how the leads are represented in detail. In other words, if one constructs a different configuration of leads (e.g. with more leads with different temperatures, DOS, etc.), which has the same ∆(ω), i.e. the same ∆ R (ω) and ∆ K (ω) as Eq. (7), then the resulting local properties of the interacting impurity, e.g. the Green's function G(ω) are the same. This holds provided the leads contain noninteracting fermions only.
The approach we suggested in Ref. 51 precisely exploits this property. The idea is to replace the impurity plus leads system (Eq. (1)) by an auxiliary one which reproduces ∆(ω) as accurately as possible, and at the same time can be solved exactly by numerical methods, such as Lanczos exact diagonalization. Details on the construction of the auxiliary impurity system are given below.
The self energy Σ aux (ω) of the auxiliary system, obtained by exact diagonalization, is used in analogy to DMFT 28,119 as an approximation to the physical self energy of the original impurity system. Inserting Σ(ω) ≈ Σ aux (ω) into Eqs. (5), (6), together with the exact hybridization function ∆(ω) yields an approximation for the physical Green's function. From this, observables such as the current or the spectral function are then calculated. We emphasize that the accuracy of this approximation can be controlled by the difference between the ∆ aux (ω) of the auxiliary system and the physical one ∆(ω), and that this can be, in principle, systematically improved, as discussed below.

C. Auxiliary open quantum system
The idea presented here is strongly related to the ED approach for the DMFT impurity problem in equilibrium. 28,119 Here, the infinite leads are replaced by a small number of bath sites, whose parameters are optimized by fitting the hybridization function in Matsubara space. The reduced system of bath sites plus impurity is then solved by Lanczos ED. 120 This approach cannot be straightforwardly extended to the non-equilibrium steady state case for several reasons: (i) since the small bath is finite, its time dependence is (quasi) periodic, i.e. no steady state is reached, (ii) there is no Matsubara representation out of equilibrium, 121 thus, one is forced to use real energies but (iii) in this case m (∆ R aux (ω)) of the small bath consists of δ-peaks and can hardly be fitted to a smooth ∆ R (ω). The solution we suggested in Ref. 51 consists in additionally coupling the small bath to a Markovian environment, which makes it effectively "infinitely large", and solves problems (i) and (iii) above. Specifically, we replace the impurity plus leads model (Eq. (1)) by an auxiliary open quantum system consisting of the impurity plus a small number of bath sites, which in turn are coupled to a Markovian environment. The dynamics of the system (consisting of bath sites and impurity), including the effect of the Markovian environment is expressed in terms of the Lindblad quantum master equation which controls the time dependence of its reduced density operatorρ: 107,122 ρ =Lρ .
The Lindblad super-operator 123 consists of a unitary contribution as well as a non-unitary, dissipative term originating from the coupling to the Markovian environment where [Â,B] and {Â,B} denote the commutator and anti-commutator respectively. The unitary time evolution is generated by the Hamiltonian describing a fermionic "chain" (E µν is non-zero only for on-site and nearest neighbor terms). It is convenient to choose the interacting impurity at site f = 0 and N B auxiliary bath sites at µ, ν = 1, · · · , N B (see Fig. 1). As usual, d † µσ /d µσ create/annihilate the corresponding auxiliary particles. The quadratic form of the dissipator (Eq. (10b)) corresponds to a noninteracting Markovian environment. The dissipation matrices Γ (κ) µν , κ ∈ {1, 2} are hermitian and positive semidefinite. 122 The advantage of replacing the impurity problem by the auxiliary one described by Eqs. (9)- (11), is that for a small number of bath sites the dynamics of the interacting auxiliary system can be solved exactly by diagonalization of the super-operatorL in the space of many-body density operators (see Sec. II D 2).
Intuitively, one can consider the effective system as a truncation of the original chain described by Eq. (1), whereby the Markovian environment compensates for the missing "pieces". However, this would still be a crude approximation, and in addition, it would not be clear how to introduce the chemical potential in the Markovian environment (except for weak coupling). Our strategy, similarly to the equilibrium case, consists in simply using the parameters of the auxiliary system in order to provide an optimal fit to the bath spectral function ∆(ω). The parameters for the fit are, in principle, E µν and Γ (κ) µν . However, one should consider that there is a certain redundancy. In other words several combinations of parameters lead to the same ∆(ω). For example, it is well known in equilibrium that in the case of the E µν one can restrict to diagonal and nearest neighbor terms only. 124 The accuracy of the results will be directly related to the accuracy of the fit to ∆(ω), and this is expected to increase rapidly with the number of fit parameters, which obviously increases with N B . On the other hand, also the computational complexity necessary to exactly diagonalize the interacting auxiliary system increases exponentially with N B . The fit does not present a major numerical difficulty, as the determination of the hybridization functions of both the original model (Eq. (7)), as well as the one of the auxiliary system ∆ aux (ω) described by the Lindblad equation (10) require the evaluation of G 0 (cf. (6)), i.e. the solution of a noninteracting problem.
The fit is obtained by minimizing the cost function with respect to the parameters of the auxiliary system. The advanced component does not need to be considered as ∆ A = ∆ R * . Of course, like in ED based DMFT, there exists an ambiguity which is related to the choice of the weight function W α (ω), which also sets the integral boundaries. This uncertainty is clearly reduced upon increasing N B . Depending on the expected physics, it might be useful to adopt a energy-dependent weight function. This could be used for example to describe the physics around the chemical potentials more accurately.
Once the auxiliary system is defined in terms of E µν and Γ (κ) µν , the corresponding interacting nonequilibrium problem Eq. (10) can be solved by an exact diagonalization of the non-hermitian super-operatorL within the space of many-body density operators. The dimension of this space is equal to the square of the dimension of the many-body Hilbert space, and thus it grows exponentially as a function of N B . Therefore, for N B ≥ 4 a non-hermitian Lanczos treatment must be used. The solution of the noninteracting Lindblad problem is non standard (see e.g. Ref. 125), and a method particularly suited for the present approach is discussed in Sec. II D 1.

D. Green's functions of the auxiliary Lindblad problem
In this section we present expressions for the Green's functions of the auxiliary system. Specifically, we will derive an analytic expression for the noninteracting Green's functions in Sec. II D 1, and illustrate the numerical procedure to determine the interacting ones in Sec. II D 2. The derivations make largely use of the formalism of Ref. 126, see also Ref. 127. For an alternative appealing approach to the noninteracting case, see also Ref. 125. All Green's functions discussed in Sec. II D are the ones of the auxiliary system, which are different from the physical ones for N B < ∞.
The dynamics of the auxiliary open quantum system described by the super-operatorL (Eq. (10)) can be recast in an elegant way as a standard operator problem in an augmented fermion Fock space with twice as many sites. [125][126][127][128] Specifically, one introduces "tilde" operators where |S are many-body states of the original Fock space, |S the corresponding ones of the tilde space 126 and N S the number of particles in S. The nonequilibrium density operator can be written as a state vector in this augmented space The Lindblad equation is rewritten in a Schrödinger like fashion 123,126 where nowL is an ordinary operator in the augmented space.L =L 0 +L I is conveniently represented in terms of the operators of the augmented space in a vector notation: 129 Its noninteracting part L 0 reads in the augmented space 123,126 where Tr denotes the matrix trace and the matrix h is given by with Λ = Γ (2) + Γ (1) , Its interacting part has the form 126 In this auxiliary open system, dynamic two-time correlation functions for two operatorsÂ andB of the system can be expressed as whereρ U is the density operator of the "universe" U composed of the system and Markovian environment, tr is the trace over the system degrees of freedom, tr E the one over the environment, tr U = tr ⊗ tr E the one over the universe,Ô U (· · · ) denotes the unitary time evolution of an operatorÔ according to the Hamiltonian of the universeĤ U . Here, 107 Notice that the time evolution ofρ U (t), as well as the one in Eq. (19) are opposite with respect to the Heisenberg time evolution of operators. This is the convention for density operators. For t = t 2 − t 1 > 0 one can use the quantum regression theorem 107 which holds under the same assumptions as for Eq. (9). It states that In the augmented space, in the same way as for (14)- (15), one can associate the operator (19) with the state vector |A t1,t =Â t1,t |I . For this vector, (20) translates into Considering its initial value (time t = 0) the solution of (21) reads Therefore, we have for the correlation function (18) for is the non-hermitian time evolution of the operatorB, and we have exploited the relation 126 I|L = 0. For the steady state correlation function, which depends on t = t 2 − t 1 , we have whereρ ∞ is the steady state density operator. Since the quantum regression theorem only propagates forward in time, for t < 0 one has to take the complex conjugate of Eq. (18), which gives for the t < 0 steady state correlation function denoted as G − BA : Using (24), the steady state greater Green's function for times t > 0 reads 130 We can use (24) also for the lesser Green's function, however for 130 t < 0 For the opposite sign of t, we can use (25), so that for both Green's functions one has 123,130 For the Fourier transformed Green's function, defined, with abuse of notation as relation (27) translates into We need the retarded and the Keldysh Green's functions whereby both relations hold for the time-dependent as well as for the Fourier transformed ones.

Noninteracting case
To solve the noninteracting Lindblad problem described by (16), one first diagonalizes the non-hermitian matrix 126 h in (17): where ε is a diagonal matrix of eigenvalues ε µ . The noninteracting Lindbladian Eq. (16) can then be written as in terms of the normal modes and a constant η. The normal modes still obey canonical anticommutation rules but are not mutually hermitian conjugate. The steady state |ρ ∞ obeys the equation Let us now consider the time evolution (22) of a state initially consisting of the normal mode operators applied to the steady state density matrix: If m (ε µ ) < 0 this term diverges exponentially in the long-time limit, which would be in contradiction to the fact that |ρ ∞ is a steady state, unless the state created by ξ µ is zero. Therefore, we must have Similarly, we must havē These equations, thus, define the steady state as a kind of "Fermi sea". In addition, by requiring that expectation values of the form I| ξ µ (t)ξ ν |ρ , do not diverge for large t, we obtain that From (34d) it follows that an expectation value of the form I|ξ µ ξ ν |ρ ∞ vanishes for the case m (ε µ ) < 0. For m (ε µ ) > 0 we make use of the anticommutation rules (33) together with (34b) and the fact that 126 I|ρ ∞ = tr ρ ∞ = 1 and arrive at where the matrix Similarly, The expression for the steady state correlation functions of the eigenmodes ξ ofL 0 can be now evaluated by considering that, due to the anticommutation rules the Heisenberg time evolution (23) gives Thus, In this way, the greater Green's function for t > 0 becomes where we have used (32). The Green's functions are defined with operators d µ /d † µ in the original Fock space, so that it is sufficient to know the first N B + 1 rows (columns) of V (V −1 ). For this purpose we introduce With this, the Fourier transform (28) of (35) is given by 131 and G >− 0µν (ω) is obtained with the help of (29). Similarly, the lesser Green's function for t < 0 with the Fourier transform and G <− 0 (ω) is obtained from (29). Using (30) together with (36) and (37), we get and for the Keldysh Green's function using also (29) In principle, one could just carry out the diagonalization (31) and then evaluate (38) and (39) numerically, which is a rather lightweight task. However, it is possible to obtain a (partially) analytical expression for the Green's functions. Indeed, a lengthy but straightforward calculation yields for the retarded one Similarly, for the Keldysh component of the inverse Green's function we obtain To sum up, (40) and (41) are the main results of this subsection. To evaluate ∆ aux (ω), one then uses (6), whereby one should consider that the matrix G 0 in Keldysh space is just the local one, i.e., in terms of the components local at the impurity G R 0f f and G K 0f f : In turn, G K 0f f , the f f component of G K 0 has to be obtained from (41) by the well-known expression 116

Interacting case
The next step consists in solving the interacting auxiliary Lindblad problem described by (10a) in order to determine the Green's function and the self energy at the impurity site. This is done by Lanczos exact diagonalization within the many-body augmented Fock space.
First, the steady state |ρ ∞ has to be determined as the right-sided eigenstate of the Lindblad operatorL with eigenvalue l 0 = 0. For convenience we introducê which is a kind of non-hermitian Hamiltonian with complex eigenvalues . The dimension of the Hilbert space can be reduced by exploiting symmetries similar to the equilibrium case. The conservation of the particle number per spinN σ is replaced here by the conservation ofN σ −N σ . 51 The steady state lies in the sector N σ −Ñ σ = 0. Starting from Eq. (26), the steady state greater Green's function of the impurity reads in a non-hermitian Lehmann representation, for t > 0 where the identity n |R The analogous expression for the lesser Green's function G < µν (ω) is obtained by inserting a complete set of eigenstates in the N σ −Ñ σ = −1 sector and exchanging the elementary operators accordingly. G K µν (ω) and G R µν (ω) are obtained using Eq. (30), see also (29).
For a small number of bath sites N B ≤ 3 the dimension of the augmented Fock space is still moderate, and eigenvalues and eigenvectors can be determined by full diagonalization. For N B ≥ 4 a non-hermitian Lanczos procedure has to be carried out. Especially extracting the steady state is not an easy task, since it lies in the center of the spectrum. Details of our numerical procedure are given in App. A.
Once the interacting and non-interacting Green's functions of the auxiliary system at the impurity site G(ω) and G 0 (ω), respectively, are determined, the corresponding self energy is obtained via Dyson's equation in Keldysh space Eq. (5). The individual components are explicitly 51 As discussed in Sec. II B, this is used in the Dyson equation (5) for the physical Green's function.

III. RESULTS
In this section, results for the steady state properties of a symmetric, correlated Anderson impurity coupled to two metallic leads under bias voltage are provided. We assess the validity of the proposed method by discussing the fit of the hybridization function and outline how uncertainties are estimated. Results for the current voltage characteristics and the non equilibrium spectral function are presented and compared with data from TEBD 77 and SNRG 54 calculations, respectively. The effect of a linear correction of the calculated Green's functions is illustrated.

A. Hybridization functions
The optimal representation of the exact bath ∆(ω) by the auxiliary one ∆ aux (ω) is obtained by minimizing the cost function Eq. (12). In practice this is done by employing a quasi-Newton line search. 132,133 In particular, we chose an equal weighting of the retarded and the Keldysh component W R (ω) = W K (ω) = Θ(ω c − |ω|). After finding our results to be robust upon different values for the cut-off ω c , as well as upon using different norms (n = 1, 2) in Eq. (12), we finally choose ω c = 50 ∆ 0 and consider imaginary parts ( m (∆ α (ω) − ∆ α aux (ω))) 2 in the cost function only. This is justified since ∆ K aux (ω) is purely imaginary and the real part of ∆ R aux (ω) is connected to its imaginary part via the Kramers-Kronig relations. 134 The asymptotic behavior of ∆ R aux (ω) is determined by Λ f f whereas the one of ∆ K aux (ω) by Ω f f . Therefore, the correct asymptotic limit lim ω→±∞ ∆ aux (ω) = 0 is guaranteed by taking Γ (1) f µ = 0 due to the requirement of semi-positive definiteness of Γ (κ) µν . Particle-hole symmetry allows for a further reduction of the auxiliary system parameters. 135 In this work we use an even number of auxiliary bath sites N B = 2, 4, and 6 in a linear setup (see Fig. 1 (b)) with an equal number to the left and to the right of the impurity (only Fig. 6 displays one calculation for an odd number of bath sites). In Fig. 2, the obtained auxiliary hybridization functions are compared with the exact ones for various bias voltages. We find a quick convergence as a function of N B , which degrades for large bias voltage φ. The Fermi steps at the chemical potentials in ∆ K aux (ω) cannot be properly resolved in the case of N B = 2. Espe-cially in the case of φ = 10 ∆ 0 the auxiliary hybridization functions for N B = 6 as well as for N B = 4 agree fairly well with the exact one and capture all essential features, in particular the Fermi steps. The auxiliary bath develops spurious oscillations in ∆ R aux (ω) at the energies of the Fermi levels of the contacts. Here the discrepancy with ∆ R (ω) is considerable in magnitude, but extends over small ω-intervals, thus inducing only small errors in the self energies.
When following the absolute minimum of the cost function Eq. (12) as a function of some external parameter, such as, e.g., the bias voltage φ, spurious discontinuities appear due to the fact that local minima cross each other. This occurs for large bias voltages and large U , and/or small N B , for which the approach is more challenging. An example for such a situation is shown in Fig. 2 for the case N B = 4, when comparing the hybridization functions just before and after such a crossing, i.e. for φ = 27 ∆ 0 and φ = 28 ∆ 0 . Even though the changes in the exact hybridization function are only minor, ∆ aux (ω) displays a considerable difference. The influence of this spurious effect on observable quantities is shown in Fig. 3 (right panel, orange circles) for a different parameter set of N B = 6 at around φ c = 33 ∆ 0 . The artificial discontinuity in the current is caused by the shift of spectral weight in ∆ aux (ω).
To deal with these discontinuities, we adopt a scheme which is suitable for obtaining a continuous dependence  (1) with tight-binding leads and on-site interaction U = 12 ∆0 (left) and U = 20 ∆0 (right). Results for three different auxiliary systems with NB ∈ {2, 4, 6} are displayed and compared with reference data from TEBD (magenta dotted and ×). 77 We plot the averaged mean values connected by lines together with error bars determined according to Sec. III A and App. B. The additional data marks for U = 20 ∆0 are as follows: The circles for NB = 6 display j(φ) when considering the absolute minimum in the fit (12). NB = 2, lc and NB = 4, lc present the results of a linear correction of the current values of the absolute minima as described in App. C. The inset displays the difference ∆j of the calculated currents to the TEBD results.
of observables on external parameters and in addition, allows to estimate their uncertainties (see Fig. 3). We first identify a set of local minima of the cost function Eq. (12), obtained by a series of minimum searches starting with random initial values. These local minima are then used to calculate an average and variance of physical quantities, such as the current. We consider the distribution of local minima with a Boltzmann weight associated with an artificial "temperature", whereby the value of the cost function Eq. (12) is the associated "energy". This artificial temperature for the Boltzmann weight is chosen in such a way, that the averaged spectral weight of the hybridization function as a function of φ is as smooth as possible. Details are outlined in App. B. A possible pitfall however is, that physical discontinuities, i.e., real phase transitions could be overlooked. It is thus compulsory to additionally investigate the results for the absolute minima and for different bath setups carefully. This approach has a certain degree of arbitrariness. However, we point out that it only affects regions with large error bars in Fig. 3, i.e. large φ and large U for which also other techniques are less accurate.

B. Current voltage characteristics
After evaluating the interacting impurity Green's function of the physical system according to (5) with the self energy evaluated in Sec. II D, we are able to determine the steady-state current. This is done with the help of the Meir-Wingreen expression 116,136,137 in its symmetrized form, where we have already summed over spin γ λ (ω) = −2|t λ | 2 m (g R λ (ω)) are the "lead self-energies" and p F,λ (ω) = p F (ω − µ λ ) denotes the Fermi distribution of lead λ with chemical potential µ λ .
To quantify the accuracy of the method we compare the results for the current voltage characteristics with quasi exact reference data from TEBD. 77 We find very good agreement for interaction strength U < 12 ∆ 0 . Since in this paper we want to benchmark the approach in "difficult" parameter regimes, in the following, we will discuss U 12∆ 0 only. In Fig. 3 we display data for U = 12 ∆ 0 and U = 20 ∆ 0 . The data points and error bars shown are obtained by using the averaging scheme as described in App. B. For the universal physics at small and medium bias voltages φ 20 ∆ 0 , the current as a function of the auxiliary system size (N B ∈ {2, 4, 6}) converges rapidly to the expected result. The convergence is even monotonic in a broad region of the parameter space. The zero bias response is linear for all N B and approaches the results expected from the Friedel sum rule 29 j(φ = 0 + ) = 2 e 2 h φ quickly for increasing N B . For U = 12 ∆ 0 already the N B 4 results yield a good reproduction of the current in this bias regime. For U = 20 ∆ 0 and φ 20 ∆ 0 a larger difference between the N B = 4 and N B = 6 results is observed. Notice that also other available methods do not yield a satisfactory result in this parameter regime. In the lead dependent high bias regime the fit becomes more challenging and large variances appear in the calculated quantities. This indicates the presence of many competing local minima with similar values for the cost function whose value tends to increase with increasing φ. For φ ≥ 40 ∆ 0 the densities of states of the left and the right contact do not overlap anymore and the current has to vanish. This limit cannot be exactly reproduced by the proposed approach due to spurious long-range Lorentzian tails present in the auxiliary Markovian environment. Nevertheless, j(φ = 40 ∆ 0 ) approaches zero as one increases the number of bath sites. This holds true for quantities obtained at the absolute minimum of the cost function as well as for averaged ones.
To extrapolate our results to larger N B , a scheme for linear corrections is discussed in App. C. Data for N B = 2, lc and N B = 4, lc, whereby "lc" denotes "linear correction", is shown in Fig. 3. For large U = 20 ∆ 0 and small to medium bias voltages φ 20 ∆ 0 , a solid improvement towards the TEBD reference values is observed (see inset Fig. 3). Correction ratios r (see App. C) close to one indicate a good applicability of the linear correction scheme. We find on average r ≈ 0.75 for φ 20 ∆ 0 (N B = 2, lc and N B = 4, lc). In the high bias regime, however, the linear correction cannot be applied with large magnitude and r drops below 0.5 for N B = 2, lc. Nevertheless, the calculation of the effective, auxiliary hybridization function ∆ aux,r (ω) as described in App. C, successfully avoids an "over-correction" of the current values and automatically allows one to estimate the reliability of the results.
Judging from the larger uncertainty from the averaging procedure and the strong effects of the linear corrections, we conclude that the high bias regime is more sensitive to the details of the fitted, auxiliary hybridization function. The universal low and medium bias regime are however very well reproduced even with a small number of auxiliary bath sites.

C. Non-equilibrium spectral function
The bias-dependent single particle spectral function is evaluated from the physical steady state Green's function of the impurity A(ω) = − 1 π m (G R (ω)). Results obtained using N B = 6 for U = 12 ∆ 0 and U = 20 ∆ 0 are presented for the whole bias range of interest in Fig. 4. Data for N B = 4 are similar but here the Kondo physics cannot be reproduced as accurately as in the case of N B = 6. Our approach does preserve the local charge 2π m (G K (ω)) = 1 and magnetization m f = 0 as well as the spectral sum rule. 138 The presented method reproduces qualitatively correctly also the equilibrium physics at φ = 0, since A(ω) displays a Kondo resonance at ω = 0 and two Hubbard satellites at the approximate positions ω ≈ ±U/2. This renders the application to equilibrium DMFT problems an interesting perspective. The width and magnitude of the Kondo resonance is discussed in comparison with (S)NRG data in Sec. III C 1.
Upon increasing the bias voltage, the Kondo resonance splits up and two excitations are observed at the energies of the Fermi levels of the leads. 78,139,140 For U = 12 ∆ 0 , the splitted resonances merge into the Hubbard bands at approximately φ ≈ 15 ∆ 0 and cannot be clearly identified thereafter. In contrast, in the case of U = 20 ∆ 0 the resonances overlap with the Hubbard satellites and can still be observed in the spectrum A(ω) at higher voltages. Calculations with increasing U in the high bias regime φ ≈ 40 ∆ 0 have shown the consistency of this effect and that a minimum value of U ≈ 15 ∆ 0 is needed in order for the resonances at the Fermi energies to be perceptible after having crossed the Hubbard bands.

Comparison with scattering states numerical renormalization group
We compare the computed spectral functions with results obtained by means of SNRG. 53 For this purpose, we use a flat DOS Eq. (3) for the leads, as in Ref. 53. Focusing on the low bias regime and N B = 6, the obtained spectral functions are depicted in Fig. 5. Compared with SNRG, our results do not achieve the same accuracy in the low energy domain, i.e. in the vicinity of ω ≈ 0. However, our data provides a better resolution at higher energies. When inspecting the Kondo peak in the equilibrium case φ = 0, our results do not fully fulfill the Friedel sum rule. 29,141,142 Depending on parameters the height of the Kondo resonance is underestimated. This is due to the fact that the imaginary part of the self-energy at ω = 0 has a small finite value which is due to the Lorentzian tails of the Markovian environment.
The resolution does not suffice to tell whether a two or a three peak structure is present for very low bias voltages φ 2 ∆ 0 . Nevertheless, one can say that the higher bias regime φ > 4 ∆ 0 is resolved more accurately and one is able to clearly distinguish the excitations at the Fermi energies of the contacts from the Hubbard satellites. The observed linear splitting is consistent with experiments on nanodevices. 139,140 Within second-order Keldysh PT 58 and QMC results 143 the resonance does not split but is suppressed only. In fourth-order and in NCA it splits into two, which are located near the chemical potentials of the two leads. 58 Other methods yield a splitting with features slightly different in details: real-time diagrammatics, 144 VCA, 78 imaginary potential QMC 145 or scaling methods. 146 Overall, a good qualitative agreement with the SNRG results is achieved which underlines the reliability of the calculated spectral functions.  Results are obtained for NB = 6 and at the absolute minimum of Eq. (12). For a comparison with SNRG 53 (Fig. 2a therein), note that their Γ = 2 ∆0.

Linear correction of Green's functions
Here we consider the effect of a linear correction of the Green's functions, as outlined in App. C. In the left panels (right panels) of Fig. 6, we show data for N B = 2 (N B = 4) including linear corrections (r = 1) for a high interaction strength in the low bias regime. We benchmark to data obtained using N B = 6 without corrections.
For N B = 2 without linear corrections, the spectral function of the auxiliary system does not feature excitations at the Fermi energies of the contacts (ω = ±2 ∆ 0 ), which are present in the N B = 6 data. Also the spectra appear washed out. The linearly corrected result, however, features not only the two resonances at the appropriate energies but also the shoulders present in the reference data. Again in the Keldysh Green's function a large correction towards the more accurate N B = 6 results is observed. To highlight the fact that the improvement of the linear correction is not only due to the inclusion of one additional bath site, also a calculation for an auxiliary system with N B = 3 is shown. Evidently, the N B = 3 spectral function exhibits a large weight at low frequencies, but, the resolution is rather low and only a single, smeared out peak at ω = 0 ∆ 0 is observed. It clearly does not account for the splitting of the Kondo resonance.
For N B = 4, a similar enhancement is found. Clearly the size of the corrections is much smaller. Especially in the Keldysh component, the Green's function for N B = 6 and for the corrected N B = 4 system nearly coincide. In general, the difference between the N B = 6 and the N B = 4 calculations (raw and corrected) are quite small, so that the presented spectral functions in Fig. 5 for larger values of φ 12 ∆ 0 can be assumed to be quite accurate.
Overall, the linear correction enables a vast improvement in the universal low and medium bias regime for all U , which becomes especially important for large U . For large bias voltages, when lead band effects become prominent, the linear correction is more challenging (see also Sec. III B).

IV. CONCLUSIONS
We have presented a numerical approach to study correlated quantum impurity problems out of equilibrium. 51 The auxiliary master equation approach presented here is based on a mapping of the original Hamiltonian to an auxiliary open quantum system consisting of the interacting impurity coupled to bath sites as well as to a Markovian environment. The dynamics of the auxiliary open system are controlled by a Lindblad master equation. Its parameters are determined by a fit to the impurity-environment hybridization function. This has many similarities to the procedure used for the exactdiagonalization dynamical mean field theory impurity solver, but has the advantage that one can work directly with real frequencies, which is mandatory for nonequilibrium systems.
We have illustrated how the accuracy of the results can be estimated, and systematically improved by increasing the number of auxiliary bath sites. A scheme to introduce linear corrections has been devised. We presented in detail how the non-equilibrium Green's functions of the correlated open quantum system are obtained by making use of non-hermitian Lanczos diagonalization in a super-operator space. These techniques make the whole method fast and efficient as well as particularly suited as an impurity solver for steady state dynamical mean field theory. 51 In this work we have applied the approach to the single impurity Anderson model, which is one of the paradigmatic quantum impurity models. We have analyzed in detail the systematic improvement of the current-voltage characteristics as a function of the number of auxiliary bath sites. Already for four auxiliary bath sites, results show a rather good agreement with quasi exact data from time evolving block decimation 77 in the low and medium bias regime. In the high bias regime, the current deviates from the expected result with increasing interaction strength. However, we have shown how to estimate the reliability of the data from the deviation of the hybridization functions and how results can be corrected to linear order in this deviation. The impurity spectral function obtained in our calculation features a linear splitting of the Kondo resonance as a function of bias voltage. Good agreement with data from scattering state numerical renormalization group 53 was found.
Applications of the presented method to multi-orbital correlated impurities or correlated clusters is in principle straightforward, although numerically more demanding. Such systems are themselves of interest as models for transport through molecular or nanoscopic objects and as solvers for non-equilibrium cluster dynamical mean field theory. In this case, a larger number of auxiliary sites might be necessary to obtain a good representation of the various hybridization functions. For this situation, one should use numerically more efficient methods to solve for larger correlated open quantum systems, such as matrix product states and density matrix renormalization group, possibly combined with stochastic wave-function approaches, [147][148][149] , sparse polynomial space, 150,151 or configuration interaction approaches. 152 A more accurate determination of low-energy, and possibly critical properties might be achieved by a combination with renormalization group iteration schemes, similar to the numerical renormalization group. Work along these lines is in progress.
Although we have presented results for the steady state, where the method is most efficient, also extensions to time dependent phenomena provide an interesting and feasible perspective. While other approaches, such as time-dependent density matrix renormalization group, 73 or quantum Monte Carlo, 49 are certainly more accurate at short times, the present approach could be used to estimate directly slowly-decaying modes by inspecting the behavior of the low-lying spectrum of the Lindblad operator.
index for the sake of clarity. The biorthogonal Lanczos vectors are determined by the three term recurrence and a normalization constant c n such that φ n L |φ n R = 1. One has a certain degree of freedom in the choice of c n and k n due to the relation k * n = c n , which is fulfilled for example by k n = k * n = c n . In the Krylov basis,L takes on a tridiagonal form T nm = φ n L |L |φ m R with the matrix elements T nn = e n , T n−1n = k n and T nn−1 = k * n . When n + 1 becomes as large as the degree of the minimal polynomial ofL, the eigenvalues and eigenvectors of T represent those of L. 153,161 If one truncates the Krylov basis, this statement holds still approximately true, especially for the largest eigenvalues in magnitude. Analogous to the hermitian case, 165 an exponential convergence of the eigenspectrum of T towards the one ofL is observed, which is of particular importance for the calculation of Green's functions. A peculiarity of the two-sided Lanczos scheme is, that not every Krylov subspace guarantees that m ( n ) < 0 for all eigenvalues n of T . In order to obtain the appropriate pole structure for the estimated Green's functions, one has to check m ( n ) < 0 together with convergence criteria. In cases in which m ( n ) < 0 cannot be fulfilled exactly, it has to be ensured at least that the corresponding weights of these eigenvalues are negligible.
For the calculation of the Green's functions needed here it is convenient to choose appropriate initial vectors, which are in the case of the greater Green's function (43) When denoting by n and U k,n the eigenvalues and rightsided eigenvectors of T , respectively, Eq. (43) can be cast into the form This section contains details on the approach we used to determine the artificial "temperature" for the Boltzmann weights as described in Sec. III A. We consider the situation that a set of local minima for which Eq. (12) becomes stationary is known. Let us specify by a y (φ) the vector of parameters {E µν , Γ (κ) µν } y corresponding to one certain local minimum for a set of model parameters, labeled by y. In order to quantify the spectral weight distribution of the corresponding hybridization function ∆ aux (ω; a y (φ)), we define which are similar to the second and third moment of ∆ R aux and ∆ K aux , respectively. For the Keldysh component a definition analogous to the first moment would yield the desired information as well but the choice above has been found to be more sensitive to details of ∆ K aux . The value of the corresponding cost function χ(a y (φ)) of the y-th minimum is used as an artificial "energy" and enables one to define weights when making use of Boltzmann's statistic where we introduced an artificial "temperature" β −1 . For each bias voltage separately, we are then able to calculate averaged quantities as well as m K 3 (φ, β) and χ(φ, β) in an analogous manner. The quantities m R 2 (φ, β) and m K 3 (φ, β) provide an estimate of the center of the spectral weight for the averaged set of hybridization functions for each bias voltage φ.
Our goal is that these quantities vary in a smooth way when changing the bias voltage. To achieve this, we employ a minimum curvature scheme, 133 meaning that we optimize the function with respect to β. This determines the optimal artificial "temperature", which ensures that the averaged cost function as well as the averaged spectral weight are as smooth functions of φ as possible, given the set of calculated minima {a y (φ)}. As in many optimization problems, an arbitrariness exists in the definition of the quantities m R 2 (φ, β) and m K 3 (φ, β), as well as in choosing the values of the weights w R , w K and w χ . In our case all of the weights were chosen to be equal to one in units of t.
An improvement of the results, to a certain degree at least, could be expected when making use of extensions like a bias-dependent β(φ). This has not been considered in the present work since already a single variable β provided quite smooth observables. As mentioned in the main text, in any case, it is obligatory to examine besides the averaged results also the ones for the absolute minima and/or for different averaging schemes, in order to avoid that physical discontinuities are averaged out. We stress that this approach has to be taken with due care, since it is in some aspects arbitrary. However, it is useful to give an estimate of the error of the calculation, and can certainly identify regions in parameter space where the error is negligible small.

Appendix C: Linear corrections
In this section, we present a scheme to correct physical quantities up to linear order in the difference 166 with r = 1. We evaluate the functional derivative numerically in the following way. One first evaluates O[∆ aux ] at the optimum ∆ aux (ω). Then O is evaluated at a "shifted" m (∆ α (ω)), obtained by adding a delta function peaked around a certain energy ω 0 δ ω0 (ω) ≡ δ(ω − ω 0 ) , multiplied by a small coefficient . The functional derivatives are then approximated linearly, by making use of the equations which become exact in the → 0 limit. A (quasi) delta-peak correction δ ω0 to ∆ α (ω) can be obtained by attaching an additional bath site (N B + 1) with on-site energy E N B +1,N B +1 = ω 0 directly to the impurity site with a hopping E N B +1,f = /π. The sum of Γ (1) N B +1,N B +1 and Γ (2) N B +1,N B +1 is proportional to the width of δ ω0 and thus, should be taken as small as possible. In practice, one uses a discretization of the integration over ω 0 in Eq. (C1) and the width of the delta peaks has to be adjusted accordingly. Setting one of the components Γ κ N B +1,N B +1 to zero yields a peak in the Keldysh component with a coefficient ±2 , respectively, as used in Eq. (C2).
Notice that the functional derivative Eq. (C2) amounts to carrying out two many-body calculations for each point ω 0 on a system with N B + 1 bath sites. However, it is not necessary to repeat the calculation for each physical quantity of interest. In the linearly corrected current values presented in Sec. III B, a ω 0 mesh of 200 points was used, whereby this number is likely to be reduced when optimizing the method.
Strictly speaking the coefficient r in Eq. (C1) should be 1. However, for cases in which the linear correction is not small, this could produce an "over-correction". In order to avoid this, we introduce a smaller ratio r which is determined as follows: We evaluate the corrected self energy at each ω via Eq. (C1) and O = Σ(ω) with some value of r < 1 and denote it Σ r (ω). We do the same for the Green's function of the auxiliary system and denote it G r (ω). Using Equations (5) and (6) we now have an estimate of an effective r-dependent auxiliary hybridization function of the linearly corrected system via ∆ aux,r (ω) ≡ g −1 0 (ω) − G −1 r (ω) − Σ r (ω) .
In principle, for r = 1 this gives ∆ ex (ω) up to O(D 2 ). In practice, for finite D(ω), one can introduce a cost function χ(r) analogous to Eq. (12) to minimize the difference |∆ aux,r (ω) − ∆ ex (ω)| as a function of r. We checked that for the case in which the linear correction is a good approximation, the minimum occurs at r = 1. If the minimum of χ(r) is situated at some value r min < 1 then one corrects also other physical quantities according to Eq. (C1) with the same r = r min . Alternatively to the correction Eq. (C1) discussed above, one can use the numerical functional derivative evaluated via Eq. (C2) in order to estimate the sensi-tivity of the value of O with respect to variations of m (∆ α aux (ω)) as a function of ω and α. This is of use, in a second step, to adjust the weight function W α (ω) in Eq. (12), so that more sensitive ω regions acquire a larger weight.