Renormalized Lindblad Driving: A Numerically-Exact Nonequilibrium Quantum Impurity Solver

The accurate characterization of nonequilibrium strongly-correlated quantum systems has been a longstanding challenge in many-body physics. Notable among them are quantum impurity models, which appear in various nanoelectronic and quantum computing applications. Despite their seeming simplicity, they feature correlated phenomena, including emergent energy scales and non-Fermi-liquid physics, requiring renormalization group treatment. This has typically been at odds with the description of their nonequilibrium steady-state under finite bias, which exposes their nature as open quantum systems. We present a novel numerically-exact method for obtaining the nonequilibrium state of a general quantum impurity coupled to metallic leads at arbitrary voltage or temperature bias, which we call"RL-NESS"(Renormalized Lindblad-driven NonEquilibrium Steady-State). It is based on coherently coupling the impurity to discretized leads which are treated exactly. These leads are furthermore weakly coupled to reservoirs described by Lindblad dynamics which impose voltage or temperature bias. Going beyond previous attempts, we exploit a hybrid discretization scheme for the leads together with Wilson's numerical renormalization group, in order to probe exponentially small energy scales. The steady-state is then found by evolving a matrix-product density operator via real-time Lindblad dynamics, employing a dissipative generalization of the time-dependent density matrix renormalization group. In the long-time limit, this procedure converges to the steady-state at finite bond dimension due to the introduced dissipation, which bounds the growth of entanglement. We test the method against the exact solution of the noninteracting resonant level model. We demonstrate its power using an interacting two-level model, for which it correctly reproduces the known limits, and gives the full $I$-$V$ curve between them.


I. INTRODUCTION
Quantum impurity models have fascinated theoreticians for several decades. These models seem extremely simple -they describe a small, typically interacting, quantum system, i.e., the impurity, coupled to a noninteracting environment. The quantum impurity consists of only a few degrees of freedom, so that its spectrum can be obtained exactly. However once this interacting impurity is coupled to the seemingly innocent quadratic environment, it gives rise to highly correlated behavior and exotic phenomena which cannot be explained solely in terms of the bare impurity, such as the Kondo effect (including its non-Fermi-liquid multichannel varieties) [1,2]. Quantum impurities can thus be seen as the basic building blocks of higher-dimensional strongly-interacting systems. The most striking feature of these arising phenomena is that they can occur at emergent energy scales which, a priori, do not appear in the Hamiltonian of either the bare impurity or the environment. An example is the Kondo temperature, which can be smaller by several orders of magnitude than any bare energy scale. Thus, in order to expose the physics of these models, they must be analyzed in a renormalization group (RG) framework. As of today, the most successful method for treating such problems in or close to equilibrium, is Wilson's numerical renormalization group (NRG) [3,4], a numerically exact RG procedure for integrating out high-energy modes and probing arbitrarily small energy scales.
A wide range of devices with various nanoelectronic and quantum computing applications, including semiconductor quantum dots [5,6], carbon nanotubes coupled to metallic leads [7,8], and molecular junctions [9,10], can be described as quantum impurity models, with the environment corresponding to two macroscopic leads. Most of the their applications involve imposing a voltage (chemical potential) or temperature bias between the leads will result in a nonequilibrium steadystate (NESS), with a tunneling current flowing through the impurity. Experimental results for such systems have successfully been explained in different limits, e.g., by linear response theory together with equilibrium NRG for small bias, or by solving a master equation at large temperature or voltage bias [11]. However, for arbitrary bias, a quantitative theoretical description of the NESS properties is still an open challenge. Any complete solution to this problem must (i) capture interaction induced many-body correlations, (ii) resolve a wide range of energy scales, and (iii) deal with an open system at its steady-state.
Attempts to generally tackle this problem analytically, e.g., in an RG framework [12][13][14] or by Keldysh field integral formulation [15] are so far restricted to uncontrolled approximations. Bethe ansatz approaches have also been tried [16,17], but are typically case-specific. Therefore much focus has been placed on finding a general numerical solution. A class of such attempts is based on capturing the many-body correlations by modeling the environment as large (but finite) leads, and evolving the many-body state of this finite system towards a finitetime quasi-steady-state, e.g., using the time-dependent density matrix renormalization group (tDMRG) method [18][19][20]. This approach has further been extended by treating the finite leads as open systems, governed by Lindblad dynamics, and similarly evolving in time towards a well defined steady-state [21,22]. The Lindblad approach has also been recently investigated in the context of density functional theory [23]. However, these attempts are typically limited in terms of the range of energy scales explored by the finite number of energy levels in the leads, with no RG procedure exploited in order to integrate out high-energy modes. Other numerically exact approaches applied to this problem are reported in [24][25][26][27], but are also not designed to explore the wide range of energy scales. Two attempts to leverage the unrivaled success of NRG in equilibrium and extend it out of equilibrium, are the so-called scattering-states NRG [28], and the NRG-tDMRG scheme [29], with the latter a predecessor of the method presented in this paper. These attempts have been quite successful at resolving a wide range of energy scales, while also capturing the many-body correlations. Yet while the former is plagued by logarithmic discretization artifacts within the dynamical energy window, the latter is based on non-dissipative time evolution of a finite, and thus closed system, which results in a quasi-steady-state in a limited time interval, making it challenging to extract steady-state observables. In this work we present a novel algorithm combining the full power of NRG and tDMRG for capturing manybody correlations at a wide range of energy scales, together with open system dynamics, in order to obtain an actual nonequilibrium steady-state. In what follows, we will refer to this approach as the Renormalized Lindbladdriven NESS (RL-NESS) method. The starting point of the presented method is a general impurity coupled to continuous leads (i.e., leads with a continuous spectrum). As shown in Fig. 1, each continuous lead is separated into a finite set of representative discrete energy levels, which in turn are coupled to the remaining continuous modes. The impurity together with this finite set of energy levels (large enough to allow the coherent formation of, e.g., the Kondo screening cloud, and the emergence of energy scales such as the Kondo temperature), is considered as an open system, coupled to an environment consisting of the remaining continuous modes, which are traced out. The latter is performed under the Born and Markov approximations, i.e., that the environment is memory-less and indifferent to the state of the system. As a result the dynamics of the system is governed by a Lindblad [30] master equation: The Liouvillian super-operator L can be separated into a von Neumann term consisting of the discrete system Hamiltonian H, and a dissipative super-operator D, describing a suitably-modeled dissipation into the environment. The two key elements of our method are: (i) the specific choice of discrete energy levels, such that highenergy modes can be integrated out, and (ii) the numerical solution of the Lindblad equation in the low-energy dynamical regime, formulated as a tensor-network algorithm.
With these requirements in mind, the Lindblad equation is obtained as follows: The Hamiltonian of the discrete system is derived by employing a mixed discretization scheme that crosses over from logarithmic to linear level spacing at the bias scale [29]. This permits integrating out modes whose energies are high compared to the bias voltage or temperature by means of NRG, with the logarithmic RG flow eventually cut off at this scale. Instead of formally deriving the dissipators, they are chosen based on two criteria: (i) the solution of the Lindblad equation reproduces the continuum limit, and (ii) Eq. (1) can be numerically solved efficiently. An important property of the chosen dissipators is that they are local in the basis in which the leads are diagonal. A set of exact transformations, dubbed the Lindblad driven discretized leads (LDDL) scheme [31], is then applied to the Lindblad equation, mapping it to a so-called chain geometry, which, due the short-rangedness (or locality) of interactions is more favorable for treatment in the tensornetwork framework, e.g., by tDMRG [32,33]. At this stage, high-energy modes (far above the bias voltage and temperature scales) are integrated out using equilibrium NRG, arriving at a local Lindblad equation in an effective low-energy basis. The state of the system is represented as a matrix-product density operator (MPDO), and is evolved in real time by a dissipative variant of tDMRG in Liouvillian space until convergence to a steady-state is obtained. Due to the dissipation induced by the environment, the entanglement entropy of the system saturates as function of time, rather than diverging, as is the case in the absence of dissipation. Hence the long-time limit steady-state can be obtained with finite MPDO bond dimension. A full description of the method will be presented in Sec. II.
By repeating the simulations for different bias voltages, a full I-V curve can be obtained. When numerically differentiated, one obtains the differential conductance. The method is demonstrated on two spinless fermionic models: the non-interacting resonant level model (RLM), and an interacting two-level model (I2LM). The RLM, discussed in Sec. III can be solved exactly in the single particle basis (in and out of equilibrium). It will therefore serve as a benchmark for the presented method. The I2LM, discussed in Sec. IV, contains two interacting dot levels with level spacing ∆ and interaction energy U . Our method recovers known results in the limits of small and large bias, yet goes beyond them by giving the full I-V curve. Conclusions and future directions are discussed in Sec. V, followed by a series of appendixes covering technical details.

II. METHOD
In this section the RL-NESS method is outlined in detail. The initial part follows much of the strategy in Ref. [29]. We start by presenting a general impurity model with continuous leads in Sec. II.A. The leads are then discretized in Sec. II.B, resulting in a Lindblad equation for a discrete system. In Sec. II.C we follow by a short overview of the LDDL scheme, used to bring this equation into a local form, both in the Hamiltonian and in the dissipators. In Sec. II.D we integrate out highenergy modes by NRG in order to obtain a renormalized impurity. In Sec. II.E we describe a matrix-product density operator procedure for real-time evolution towards the steady-state. Finally, in Sec. II.F we discuss the extraction of observables from the obtained steady-state.
Steps B-D are described schematically in Fig. 1, and step E is described in Fig. 2. Throughout this section, superoperators acting on the density matrix will be represented in calligraphic script, while regular operators will be represented in Roman script. Tensor-network calculations (NRG, MPDO evolution) were implemented using the QSpace tensor library, which can exploit both abelian and non-abelian symmetries on a generic footing [34,35].

A. Model
The total Hamiltonian of an impurity system can be generically separated into three parts: the (interacting) impurity, the non-interacting leads with a continuum density of states, and the coupling between them: The dot Hamiltonian is given by m (here spinless) levels λ with onsite Coulomb repulsion U : with fermionic creation operators d † λ , and total impurity occupation n d = λ n dλ , where n dλ = d † λ d λ . More com-plicated local interactions, such as exchange interactions or spin Hund's coupling, may also be incorporated. The lead Hamiltonian in this work is described by two metallic, i.e., non-interacting leads located left and right of the impurity. They are assumed to be featureless, with constant hybridization Γ αλ of lead α ∈ {L, R} with the impurity level λ ∈ {1, . . . , m} over a bandwidth ε ∈ [−D, +D], resulting in the total hybridization strength v αλ ≡ 2D·Γ αλ π . The lead and coupling Hamiltonians can therefore be written in the diagonal bath basis as: where c † αε creates an electron in lead α at energy ε. As indicated, c † α0 defines the normalized bath level that the impurity couples to, i.e., at the location of the impurity, obeying {c α0 , c † α0 } = 1. The generalization to spinfull and multi-channel leads is straightforward, while the generalization to a featured hybridization function is conceptually also possible. Throughout, we assume the limit of large bandwidth, i.e., that all energy scales and parameters are much smaller than D. Without loss of generality then, the voltage bias is chosen symmetric with respect to the Fermi energy, so that the chemical potentials of the leads are µ L = −µ R = − V 2 (taking unit of charge e = 1, throughout). For concreteness we will mostly concentrate on the case of zero temperature in both leads, but the described procedure also applies to finite and non-equal temperatures.

B. Lindblad Equation
The Lindblad equation is a first-order linear differential equation. Its general solution, given some initial condition ρ 0 , can thus be written by exponentiating the Liouvillian super-operator: The dynamics in our case is designed to have a unique nonequilibrium steady-state defined by i.e., either as a solution of a linear equation (l.h.s.), or as the state to which an arbitrary initial states decays to in the long-time limit (r.h.s). The first stage in the RL-NESS method is obtaining a Lindblad equation for a discrete system from the original  The system of interest is a general impurity coupled to two macroscopic and thus continuous leads at chemical potential difference V . This system can be mapped exactly onto (b), where the bath has been coarse-grained into distinct intervals. Each of these is written as a representative level that the impurity couples to, and a continuous bath consisting of the remainder of states in that interval. The discrete set of energy levels together with the impurity then forms a finite system. The remaining lead modes are integrated out, resulting in a Lindblad bath coupled to each discrete energy level. (b') The width of the intervals is chosen according to a logarithmic-linear discretization scheme, such that levels in the low-energy window [−D * , +D * ] are equally spaced by δ (blue), with a smooth transition to logarithmic spacing ∼ Λ n at large energies (red). (c) The targeted occupation of each individual lead mode, originally encoded in the couplings to the Lindblad baths in (b), can be transferred onto the lead-impurity couplings, such that the resulting two leads now represent particles or holes, and are driven to be completely full or empty respectively. (c') This rotation of the local Liouvillian basis is performed separately for each original (physical) lead level. A general such level is coupled to the impurity with coupling constant v, and to two Lindblad baths, one filling it and the other emptying it, at rates proportional to γ and chosen such that they drive the level towards its equilibrium occupation (determined by chemical potential and temperature). From the Liouvillian description, an auxiliary level is introduced at the same energy, and linear combinations of the two levels are chosen such that one is driven to be completely full (particle) and the other to be completely empty (hole), both at rate γ. (d) The particle and hole leads can be exactly mapped onto nearestneighbor Wilson chains via tridiagonalization, with local dissipators filling one chain and emptying the other. The hopping amplitudes along the chains away from the impurity initially exhibit exponential decay due to the logarithmic discretization at large energies (red), until they cross over into more uniform hopping amplitudes of order δ in the linear discretization regime (blue). (e) The sites in the logarithmic sector, together with the impurity, are numerically integrated out in standard NRG spirit. This provides an effective subspace for the low-energy description in terms of an effective renormalized impurity (RI) with multiple dissipators. (e') This is achieved by collecting all the logarithmic sector sites into a single Wilson chain (via a re-tridiagonalization) for the sake of numerical stability of the subsequent iterative diagonalization by NRG. The fixed number of states coming out of the last NRG iteration constitute the RI low-energy subspace.
continuous system, as shown in Fig. 1(b). Formally, this can be done by dividing the full band [−D, +D] of each lead into consecutive distinct intervals I n . By the bilinear structure of the coupling in Eq. (5), the impurity couples to a particular state in each interval which itself then is coupled to the remainder of the states in that interval. The latter can be integrated out, leaving a single representative level for each interval. Explicitly performing this integration (under the Born and Markov approximations) yields the structure of the Lindblad baths and their couplings to the system. However, we will allow ourselves some freedom in choosing the exact values of the couplings to the Lindblad baths so as to simplify the subsequent simulation of the driven system, while enforcing that the correct steady-state is obtained.
The choice of the intervals I n and corresponding energy levels relates to coarse-graining that depends on a discretization scheme. A common discretization scheme used for treating quantum impurity models is the loga-rithmic discretization scheme, introduced by Wilson as part of the NRG [3,4]. This scheme produces discrete semi-infinite leads with level spacing shrinking exponentially as the lead Fermi energy is approached. It is designed to generate energy scale separation, and subsequently justifies integrating out of high-energy modes via iterative exact diagonalization as part of a logarithmic RG flow. This enables us to accurately and reliably resolve exponentially small energy scales which frequently arise in impurity models due to Kondo-like physics. However, for an open system, e.g., via coupling to a thermal reservoir or the presence of finite voltage bias, energy scale separation ceases to exist below the corresponding energy scale, and the logarithmic RG flow will be cut off. In the nonequilibrium case this gives rise to a dynamical low-energy window described by a reduced bandwidth D * , which is of order of the bias voltage or temperature (see below). For a least-biased numerical approach then, the discretization scheme within this regime should be uniform. Therefore, RL-NESS employs a mixed discretization scheme [29,36]. This consists of a logarithmically discretized region extending from the band edge down to just above the lead bias voltage or temperature, that smoothly crosses over into a linearly discretized region (with uniform level spacing) in the bias window [−D * , +D * ]. Such a scheme allows one to make use of NRG to integrate out high-energy modes (relative to V or T ), in order to obtain an effective low-energy nonequilibrium system, to be simulated in a controlled manner by a DMRG-like approach. This scheme has also been discussed for the setup of two leads with a voltage or temperature bias in a previous work [29], but without the Lindblad driving (previously suggested in Ref. [31]). It is therefore briefly outlined here for completeness.
We define D * , the characteristic energy scale of the leads, as the energy at which the Fermi-Dirac distribution of the lead drops below some pre-selected threshold. For zero temperature this implies D * = V 2 , while for finite temperature the specific value of D * depends on the chosen threshold. The intervals I n , as shown in Fig. 1 , are chosen such that in the range [−D * , +D * ] they are of equal size δ, referred to as the linear discretization parameter, while away from this range they scale exponentially as ∼ Λ n , where Λ > 1 is referred to as the logarithmic discretization parameter. In the intermediate region the interval widths cross over smoothly between being constant and growing exponentially. In each interval a representative energy level ε n is selected with corresponding coupling v αnλ to the impurity λth level. For details regarding the choice the intervals and corresponding energies and couplings see Appendix A. The same intervals are chosen for both leads such that by construction ε n are lead independent, while the coupling constants v αnλ can differ between the leads. The resulting leads and coupling Hamiltonians are: Following through with this procedure, the dissipators can formally be derived. If such a path is pursued, the continuum of states of a specific interval will serve as the environment only of its representative level, thus resulting in a local dissipator for each level: where f α (ε) ≡ f FD (ε; µ α , T α ) is the Fermi-Dirac distribution for lead α (depending on the lead specific chemical potential and temperature), and {γ n } are referred to as Lindblad driving rates. This structure implies that when the leads are decoupled from the impurity, i.e., v αnλ = 0, they are driven to their equilibrium occupation, as expected. The total Lindblad equation is then given by: where H (disc) coupling is the total Hamiltonian, now with discrete leads, and hence effectively of a finite system that becomes an open system by means of the Lindblad driving. As shown in Ref. [31], a wide range of driving rates reproduce the same continuum limit observables. Therefore one is free to choose them, in this range, so as to best suit the numerics. With this in mind, and for reasons to be explained in the following section, the rates will all be chosen to be energy independent, i.e., γ n = γ, and of order of the linear level spacing δ. Let us note that driving of energy modes (exponentially) larger than D * will have negligible effect on the results, since, importantly, these modes start and practically remain in equilibrium throughout the dynamics. Thus, the corresponding couplings can further be tuned, or even completely turned off, in order to enhance numerical stability, as we further discuss in Sec. III and Appendix F. At this point, the Lindblad equation to be solved is fully defined. As a consistency check, note that properly taking the limits of this equation converges back to the continuum limit: In the limit Λ → 1 the discretization scheme collapses to a linear (equal spacing) discretization, which in turn converges to the continuous system in the γ = δ → 0 limit [31].

C. Local Form
The Lindblad driven discretized levels (LDDL) scheme [31] is a set of exact manipulations, applied to the Lindblad equation (11) with the goal of bringing it to a form more favorable for treatment in the framework of tensornetworks. The system Hamiltonian obtained after discretization is formulated in the so-called star geometry, involving diagonal leads, as in Eq. (8), with all levels directly coupled to the impurity, as in Eq. (9). This geometry is non-local (in the sense that all lead levels couple to the impurity), and therefore less convenient in the framework of tensor-networks. The dissipators, on the other hand, are already local in this geometry, with each lead level coupled to its own Lindblad bath, which is a property we would like to retain. A standard procedure, employed for example in NRG, is to perform an exact mapping in terms of a basis transformation from the star geometry to a chain geometry [4]. The bilinear structure of the coupling term in Eq. (9) directly defines the only bath level c † α0 that the impurity couples to. This level constitutes the first site of a nearest-neighbor tightbinding chain, which can be obtained by tridiagonalizing the single-particle basis of the remainder of the lead levels, e.g., by construction of a full Krylov space: (12) Such a basis transformation, however, will result in nonlocal dissipators due to the n-dependent prefactors in Eq. (10), which include the Fermi factors. The LDDL scheme circumvents this problem and yields a Lindblad equation which is local in both the dissipators and the Hamiltonian in the chain geometry. For completeness it will be described here briefly. The idea behind this scheme is to shift the Fermi-Dirac information from the dissipators in Eq. (10) to the lead-impurity couplings in an effective Hamiltonian, still in the star geometry. With an appropriate choice of Lindblad driving rates, the system can then be tridiagonalized into the chain geometry, without loosing the locality of the dissipators. Consider a single discrete lead level with creation operator c † αn , referred to as a physical level, in lead α at energy ε n and coupling constants v αnλ to the impurity levels. We temporarily drop the subscripts αn for readability in what follows. Its dissipator is given according to Eq. (10), meaning it is constantly depopulated and re-populated at a constant Lindblad driving rate γ, weighted by 1 − f (ε) and f (ε), respectively. In the LDDL scheme this single physical level is mapped onto two artificial lead levels with creation operators c † h and c † p , referred to as hole and particle levels, thus effectively doubling the number of levels. The former is constantly depopulated at rate γ and the latter constantly re-populated at rate γ: These two levels both have the same onsite energy ε, yet are now coupled to the impurity with amplitudes that depend on temperature and chemical potentials: Formally this mapping is obtained by introducing an auxiliary level at energy ε which is decoupled both from the impurity and the Lindblad baths, and performing a unitary rotation between the physical and auxiliary levels, thus shifting the Fermi-Dirac information from the dissipators to the lead-impurity couplings, as shown in Fig. 1(c'). For more details, as well as a discussion of the resemblance of this procedure to purification of the level, or the thermofield approach, see Ref. [29]. The described procedure is repeated for each lead level. This replaces each physical lead with a corresponding hole lead and particle lead, as in Fig. 1(c), thus doubling the total number of lead levels. By selecting Lindblad driving rates γ to be energy independent, one obtains dissipators for each of the hole or particle leads which do not depend on the energy index n (i.e., are proportional to to the identity matrix w.r.t. this index). Each such lead can therefore be tridiagonalized separately into a nearest-neighbor chain while the dissipators remain unaltered, resulting in a Lindblad equation which is local both in the dissipators and the Hamiltonian, as desired [see Fig. 1 Two remarks are in order regarding the doubling of lead levels, before the tridiagonalization is performed. The first is that for physical levels lying far from the chemical potential in units of temperature, where f (ε) is 0 (1), the particle (hole) level decouples from the impurity, and can thus be disregarded in subsequent calculations. For zero temperature this holds for all physical levels, and so the described mapping is reduced to relabeling physical levels above (below) the lead chemical potential as holes (particles), with no doubling actually occurring.
The second remark relates to exploiting a left-right symmetry in the lead spectrum. In equilibrium calculations, when both leads have the same energy levels and for each lead level the left and right coupling constants to all impurity levels are proportional, only a specific linear combination of left and right levels couples to the impurity, precisely as defined by the coupling Hamiltonian. The complementary orthogonal combination of left and right levels decouples from the impurity and hence becomes irrelevant for the impurity dynamics. This simplifies the model from a two-lead model to an effective single-lead model. In the nonequilibrium case, the different potentials applied to the left and right leads break this symmetry, and prevent its exploitation. However, once the physical leads are separated into hole and particle leads, the symmetry is reinstated (for holes and particles separately), and can therefore be exploited. Thus for models in which this symmetry exists, the final number of artificial lead levels is actually smaller than the original number of physical levels.

D. Renormalized Impurity -NRG
The LDDL scheme is indifferent to the specific discretization scheme employed, as long as the Lindblad driving rates are kept energy independent. Observe now the implications of the linear-logarithmic scheme on the resulting Lindblad equation. The obtained on-site energies {ε αk } and nearest-neighbor hopping amplitudes {t αk } in the vicinity of the impurity are of the largest magnitude and decay exponentially as the distance from the impurity grows, all the way down to D * . The corresponding sites will therefore be referred to as the logarithmic sector. Below D * , the on-site energies and hopping amplitudes remain of order of the linear level spacing δ and D * , respectively, and will be referred to as the linear sector. Due to the smooth transition in the discretization, the exact boundaries between these two sectors are fuzzy, and in practice are chosen with some fine tuning in order to enhance convergence.
In the chain geometry, the logarithmic sector, including the impurity, can be considered as a mesoscopic system, coupled to the linear sector leads. By construction, the vast majority of the (many-body) energy levels of this mesoscopic system are at energies larger than D * , and so are expected to be indifferent to the voltage bias applied, thus largely remaining in the equilibrium state. They are therefore expected to contribute to the nonequilibrium dynamics only through renormalization effects on the low-energy modes in the linear sector, which in turn actively participate in the dynamics. As argued in Ref. [29], it is therefore sufficient to approximate the mesoscopic system by a renormalized impurity (RI) residing in an effective significantly reduced low-energy basis. This is a controlled approximation, as one can monitor the weight on all states in the RI while time-evolving towards ρ NESS . Note that the chemical potential of this RI is set midway between the chemical potentials of the leads, so that the effective low-energy subspace consists of states with RI particle number which is close to its average occupation in the NESS.
The RI is obtained by the following procedure, as described in Fig. 1(e'): An additional subsequent tridiagonalization is applied to merge the two (particle and hole) chains in the logarithmic sector. This brings them into a single-lead Wilson-chain structure, which is important for NRG, since it keeps correlations at a given energy scale local. An NRG sweep is then applied to the chain -starting from the impurity, at each step a site is added to the chain, the Hamiltonian is diagonalized, and highenergy modes are discarded. At the end of the sweep through the logarithmic sector, the R lowest-lying manybody states are taken as the effective basis of the RI. All operators acting on the impurity, or on sites in the logarithmic sector, are then projected to this effective reduced basis.
The leads in the linear sector, together with the RI, now form the dynamical system under consideration, as shown in Fig. 1(e). The Lindblad equation for this system still consists of a nearest-neighbor Hamiltonian, however with a more complicated local term acting on the RI site. The dissipators are also local in this setup, and again the local terms acting on the RI are more complicated, corresponding to the multiple dissipators acting on the logarithmic sector. Note that although the dissipators on different sites of the chain originally commute, the logarithmic sector dissipators, after being projected to the RI basis, no longer do. Another concern regarding the logarithmic sector dissipators is that because they were not taken into account during the RG flow, they might drive the RI out of the effective low-energy basis. Figure 2. (a) MPDO description of the system density operator -each site is described by a rank-4 tensor (i.e., 4 legs) with two physical indices, of local dimension d for chain sites or R for the renormalized impurity (RI), and two virtual indices, of bond dimension χ, connecting it to its neighboring tensors. The initial MPDO is chosen as a product state, where all particle sites are full, all hole sites are empty, and the RI is in its ground state. In practice this issue can be handled, as discussed below in Sec. III and Appendix F

E. MPDO Solution -tDMRG
The obtained Lindblad equation is solved for the steady-state by real-time evolution, implemented in the tensor-network formalism. The system (mixed) state is represented as a matrix-product density operator (MPDO) [37]. Analogously to the matrix-product state (MPS) representation of wavefunctions, which has a single (physical) index for each chain site, the MPDO has two (physical) indices for each chain site, as shown in Fig. 2(a). For the chain at hand, each of these physical indices is of dimension d corresponding to a single fermionic Hilbert space, except at the RI, where it is of dimension R corresponding to the effective low-energy subspace. It is common practice to combine the two physical indices of each MPDO site into a single effective index of dimension d 2 and simply treat it as an MPS. However, in the case of a large physical index dimension, e.g., for the RI, keeping the indices separate enables more efficient contractions and reduces the computation cost.
The derivation of the Lindblad equation respects the same local symmetries as the original continuous system. In our case the full continuous system conserves U (1) charge (particle number). The discretization procedure does not break any of these symmetries, so that the resulting Hamiltonian still conserves the same charges for the discrete finite system alone. On the other hand, the derived Lindblad operators do not necessarily conserve the charge in the finite system alone, but do respect the symmetry and conserve the charge for the full system, including the baths. Hence, although the Lindblad dynamics does not conserve this charge for the finite system, it does conserve a related quantity [38], which can be exploited in order to decompose the MPDO into symmetry sectors, further reducing the computational cost. We demonstrate this for particle number conservation. Define the super-operators N ± = N ⊗ I ± I ⊗ N , where N is the particle number operator and I is the identity. The Liouvillian super-operator L commutes with N − but not with N + . The dynamics therefore does not conserve particle number N ⊗ I = N++N− 2 . However the conservation of N − suffices in order to decompose the MPDO into particle-number symmetry sectors. It also implies that the parity of N + and thus of N ⊗ I is conserved, which suffices in order to account locally for fermionic signs. This example holds for any abelian symmetry, while a more involved argument is required in the case of nonabelian symmetries.
Following Ref. [29], the system is set in an initial state |ψ (0) which is a product state between the ground state of the decoupled RI, and the steady-state of the decoupled linear sector leads. The latter is defined as the pure product state where all lead particle (hole) sites are full (empty). This initial state can be written either as an MPS, or as an MPDO for ρ (0) ≡ |ψ 0 ψ 0 |, both with bond dimension χ = 1. This starting point is assumed to be sufficiently close to the desired final steady-state solution, so that when the coupling to the RI is turned on, the full system will quickly converge to its steady-state (as our results confirm). One could also initialize the RI in its decoupled steady-state, but in practice this does not improve convergence. Note that the initial setup, together with its transient dynamics, are regarded only as a means to obtain the desired steady-state, so that the specific choice of initial state can be fully based on numerical considerations.
The coupling between the RI and the leads is then turned on, and the system is evolved in time by a variant of tDMRG, formulated to accommodate for Lindblad dynamics. Note that in this work the RI-lead coupling is turned on in an immediate quench, and slow ramping up of the coupling, as employed in Ref. [29], was not necessary. In the spirit of tDMRG, this time-evolution is based on a second-order Trotter-Suzuki decomposition with a sufficiently small time-step τ . Then the propagator can be written as a product of short-time propagators e Lt = Nt i=1 e Lτ , with N t steps required in order to arrive at a time t = τ N t . Each short-time propagator is Trotter decomposed into local and nearest-neighbor gates based on the short-rangedness of the Liouvillian introduced above. The total Hamiltonian can be written as a sum of local two-site operators H = N −1 k=1 H k,k+1 where non-adjacent terms commute. Defining the Hamiltonian two-site super-operators as H k,k+1 ρ ≡ −i [H k,k+1 , ρ], the Liouvillian L can then be written as the sum of these two-site Hamiltonian terms and single-site dissipative (hole/particle) terms defined in Eq. (13): For an exact representation of the super-operators, the Hamiltonian terms commute with all non-adjacent Hamiltonian and dissipative terms, and the dissipative terms all commute with each other. However, inside the RI the fermionic anti-commutation relations of the original fermionic operators are compromised by the NRG truncation, which results in non-commuting terms in its vicinity. The second-order Trotter decomposition adopted here and depicted in Fig. 2(b), is similar to the one discussed in Ref. [39]: .
(16) The two-site Hamiltonian gates are given by: and the dissipative single-site gates translate into Kraus operators [40,41]. In the spinless case they are respectively given for particles or holes by: For spinfull fermions there will be 4 Kraus operators for each site, replacing η → (η, σ), with σ ∈ {↑, ↓}. For terms which are bilinear in the fermionic operators, such as the Hamiltonian or K 1η , fermionic signs arising from the anti-commutation relations can be accounted for locally. The operators K 2η in the dissipative terms, however, act simultaneously on both sides of the density matrix. Hence, they give rise to a global Jordan-Wigner string. In the present MPDO setup, it can be efficiently 'pulled' in locally [42]. This requires that charge parity is fully tracked on all tensors, which is the case here when decomposing the MPDO into U (1) charge symmetry sectors, in the sense discussed above. In the local configuration, as shown in Fig. 2(b.1), the crossing of the red line with the bond index implies that the charge parity operator Z ≡ (−1) q , with charge q, must be simultaneously applied to the bond state space when acting with K 2η .
A quick overall complexity analysis of the method can be performed assuming a fixed bond dimension χ on all MPDO sites. Since the treatment of the RI is clearly the most expensive step, the following considers operations involving the RI. The analysis is completely analogous for all other sites where one just replaces R with the regular local dimension d of a physical site. The cost of the Trotter gate contraction is O d 2 R 2 χ 3 + d 3 R 3 χ 2 , where the two terms correspond to merging the RI tensor with its neighboring tensor and to applying the nearestneighbor Trotter gate, respectively. The SVD back into local tensors then costs O d 4 R 2 χ 3 . Finally, the cost of the Kraus gates contractions is O kR 3 χ 2 , where k is the number of Kraus gates acting on the RI, which is proportional to the number of sites in the logarithmic sector. Ignoring the cost of all other sites in the linear sector, the total cost of the method can be approximated as O N t d 3 R + kR + d 4 χ R 2 χ 2 , where N t is the number of sweeps.
The most important property of the MPDO ansatz is that it can efficiently represent the steady-state, using a relatively small number of parameters. Another important constraint on the represented state is that it must be a physical state, i.e., a positive semi-definite Hermitian operator with finite trace. While the MPDO ansatz does not enforce these constraints, the Lindblad evolution (also after Trotterization) is a completely-positive trace-preserving (CPT) map [41], and thus guarantees that starting from a physical state will always result in a physical state. The only loss of positivity (and trace) can come from truncation of singular values during the tDMRG sweep. This relates to a drawback of the MPDO ansatz -the singular values obtained after a Schmidt decomposition no longer correspond to the singular values of the reduced density matrix (as for an MPS). However, if the singular values drop quickly enough, as is the case for the models analyzed here, only small singular values are truncated, resulting in negligible loss of positivity.

F. Observables
At any point throughout the evolution, single-time correlations can be extracted from the MPDO, with the long-time limit representing the steady-state value. As correlations in the chain geometry are easily obtained, in practice it is convenient to map all observables of interest to this geometry, as shown in Ref. [29]. In this work we focus on the particle current flowing from one lead to the other. The time derivative of the impurity occupation can be separated into contributions I α corresponding to the current flowing from lead α into the impurity: Thus I α is given in the continuum limit (taking e = 1, = 1), and approximated after discretization by: In the steady-state, d dt n d = 0, the current flowing from the left lead into the impurity is equal to the current flowing from the impurity into the right lead I ≡ I L = −I R . From a numerical perspective, the average combination I = (I L − I R ) /2 converges more rapidly, and is less prone to noise. Running simulations for different voltages, a full I-V curve can be obtained, and numerically differentiated in order to produce the differential conductance G (V ) = ∆I/∆V . Note that the numerical derivative is very sensitive to noise, so that the I-V curve must be obtained with high accuracy.

III. ERROR ANALYSIS AND THE RLM
In order to estimate the accuracy of the presented method, we apply it to the exactly solvable noninteracting resonant level model (RLM). This model, with dot Hamiltonian describes a single spinless impurity level with energy ε 0 (e.g., controlled by a gate voltage), coupled to two spinless leads described by Eq. (5) via the coupling Hamiltonian (4). An end-to-end comparison of the steady-state current and the differential conductance, between the exact result of the RLM in the continuum limit and the RL-NESS result, is shown in Fig. 3(a). It displays a good agreement over a wide range of bias voltages and impurity level energies ε 0 , with parameter values given in the caption. The RL-NESS real-time evolution of the current, for a typical case of V = Γ, ε 0 = 0 with finite Λ = 6 and linearly extrapolated to γ → 0, is plotted in Fig 3(b) (blue). It demonstrates several key aspects of the method. After an initial rapid rise in the current over a period of order V −1 , the oscillations (discretization artifacts related to the logarithmic sector) decay exponentially at a rate which is proportional to γ, finally stabilizing on a steady-state value. For further discussion regarding the evolution time scales see Appendix D. The convergence to the steady-state can also be observed in the lower panel, where the truncation error saturates. Throughout the entire evolution, our method displays excellent agreement with the corresponding exact result (shaded blue) of the same driving protocol. After linearly extrapolating also to the Λ → 1 limit, the RL-NESS current (black), once the NESS is reached, displays excellent agreement with the continuum limit exact result (shaded gray).
If the dissipation is initially set to γ = 0, while keeping a finite level spacing, RL-NESS reduces to the NRG-tDMRG scheme [29], in which the state of the system is represented by an MPS (instead of an MPDO), and the real-time evolution is unitary. In what follows this will simply be referred to as MPS evolution. As in the case of finite dissipation, the γ = 0 evolution of the current can be calculated (for the same Λ = 6) either explicitly as an MPS evolution (solid red) or exactly in the single particle basis (shaded red). These results agree with the RL-NESS current in the early transient oscillatory regime, but later residual oscillations persist. Thus only a quasi-steady-state is obtained, whose mean is nevertheless consistent with the γ → 0 limit. Eventually, the current drifts away, due to truncation errors that, without dissipation, do not saturate. Later on, even for the exact solution, this quasi-steady-state will be lost due to reflection off the edges of this closed system at a time t = L/v F = 2π δ , dictated by the finite linear level spacing. In stark contrast, in the case of RL-NESS, for strong enough damping γ δ, the discrete levels become sufficiently blurred out, such that the dynamics truly represents an open system, where reflection off the edges and the accompanying drop in the current no longer occur.
As part of the analysis, steady-state observables of the RLM are calculated exactly by means of Keldysh formalism, both for the continuous system and in an arbitrary discretization (Appendix B). The exact time evolution of single-time observables (in a given discretization) is also calculated by solving a differential continuous Lyapunov equation for the single-particle correlation matrix (Appendix C).
The remainder of this section is dedicated to an analysis of the two major error contributions in the method: the lead discretization error -how well do the discrete system observables represent the continuous system, and the simulation error -how accurately does the tensornetwork method solve for the discrete system steadystate. Generally there is a trade-off between the two contributions, as a finer discretization better reproduces the continuum limit, but is also harder to solve for numerically.
The lead discretization error depends both on the fineness of the discretization grid, controlled by the logarithmic Λ and linear δ discretization parameters, and on the broadening of the discrete levels, controlled by the Lindblad driving rates γ. We introduce the relative measure for the discretization error, disc ≡ max V |1 − I DL /I CL |, as the maximal relative distance over a range of bias voltages V ∈ [0.01, 100] Γ, between the exact discrete leads current I DL (for a specific choice of Λ and ratios δ/V, γ/V ) and the exact continuous leads current I CL . This measure can be explicitly evaluated for the RLM and is plotted in Fig. 3(c) as a function of γ/V and for several values of Λ. The specific choice of δ has only a minor effect, as long as δ ≤ γ, which is required in order to negate finite size effects. Fixing δ/V to a small value results in a smooth curve (solid), while taking δ = γ results in a slightly more noisy curve (shaded) with the same trend. Note that disc is approximately linear in γ/V , down to a lower bound on the error, dictated by Λ. This lower bound in turn is linear in Λ [see inset to Fig. 3(c)]. These two observation justify a linear extrapolation in these two parameters to the continuum limit The simulation error has multiple contributions, listed in ascending order of significance. First consider the Trotter error, arising from the discretization of the Liouvillian real-time evolution. In practice, exploiting secondorder Trotter decomposition, it is numerically feasible to choose sufficiently small time-steps such that this error is negligible compared to the other ones. A second source of simulation error is introduced by the NRG procedure, and controlled by the number of kept states in each NRG iteration. As in equilibrium, the number of required kept states can be reduced by taking a coarser logarithmic discretization, i.e., larger Λ. In practice, only the number of kept states R in the last NRG iteration, dictating the size of the restricted low-energy subspace of the RI, pose a computational bottleneck. The numerical cost in setting up the RI by previous NRG iterations is entirely negligible. Therefore earlier NRG iterations are less harshly truncated, but a larger Λ is still required in any case in order to keep R sufficiently small. The third, and most significant, source of simulation error is the truncation of the MPDO to a fixed bond dimension χ after each time-step, by discarding small singular values. Empirically the singular values decay faster than polynomially with the singular value index, at a rate which decreases with decreasing γ (see Appendix E). This implies that the required bond dimension (for a fixed truncation error) scales exponentially with γ. This exponential scaling can naturally be understood in the γ → 0 limit, in which the entanglement entropy grows linearly in time, thus leading to an exponential blowup in the required bond dimension. Choosing a finite γ sets a time scale 1/γ at which the entanglement entropy stops growing. For any finite γ the steady-state can therefore be represented with a finite (possibly large) bond dimension, which in the small γ limit must grow exponentially with 1/γ in order to match the expected exponential blowup. It is important to stress that even though there is an exponential bound on simulating small γ, this represents the thermodynamic limit, which can be approached by working with finite γ and then linearly extrapolating to small γ.
Finally, let us discuss the issue of whether or not to apply the Lindblad terms coupled to the RI. Physically, since the RI represents the degrees of freedom far above the voltage and temperature bias scales, it is reasonable to expect that they are barely affected by the nonequilib-rium conditions. Thus whether or not the Lindblad terms acting on the RI are applied, we expect to obtain similar results. We demonstrate this for the RLM in Appendix F. Numerically, however, the effort involved in the two approaches (for the same accuracy) is different. The effect on the numerical results becomes more pronounced in the interacting case, considered in the next section. There, for a large logarithmic discretization parameter Λ, the RI spectrum contains nearly degenerate levels, which can be coupled even by weak Lindblad driving at the sites composing the RI. In practice this can drive, and hence affect, high-energy modes in the RI (which in principle should remain in equilibrium), resulting in artifacts which are enhanced in the differential conductance. Taking small values of Λ ∼ 2 could resolve this problem. However, this necessitates increasing the number R of states kept in the RI, and therefore is often impractical. Taking the manageable intermediate value Λ = 3 for the interacting case (instead of extrapolating to Λ → 1), at the cost of a larger R = 64, suppresses these artifacts, but still does not completely eliminate them. For these reasons, in the interacting case it becomes preferable to entirely turn off the Lindblad terms coupled to the RI. This does not adversely affect the physics. To the contrary, it leads to a stable numerical solution without artifacts, with reasonable computational costs.

IV. INTERACTING SYSTEM
We now wish to demonstrate the method on an interacting system, which has no known solution for the NESS current. For this we choose an interacting twolevel model (I2LM), consisting of two interacting dot levels ε 1 , ε 2 with onsite interaction energy U , coupled to non-interacting leads. The dot levels are taken with level spacing ∆ ≡ ε 2 −ε 1 , and can be shifted by changing ε 0 by a gate voltage (taken relative to particle-hole symmetry), such that the dot Hamiltonian is given by: Both levels are coupled symmetrically to the left and right leads, so that the lead and coupling Hamiltonians are given by Eqs. (4) and (5), with equal hybridization Γ αλ = Γ. For simplicity, we take all the dot-lead couplings to be real and with the same sign. Since our main goal is demonstrating the method rather than studying the model, we do not explore the full impurity parameter space, but concentrate on restricted yet representative sets of parameter values. The level spacing and the interaction are fixed to U 5 = Γ = 5∆ such that ∆ < Γ < U , thus having a separation of energy scales in a strongly correlated regime, and the bias and gate voltages are varied. (b) NESS current and (c) its derived differential conductance at finite bias, as calculated by RL-NESS for two gate voltages, corresponding to the valley at ε0 = 0 (blue) and the peak at ε0 = U 2 (red). The low-bias behavior of the current exhibits a linear dependence for the peak, but a cubic dependence for the valley, where the linear response conductance thus vanishes quadratically in V . For comparison the equilibrium spectral function of the I2LM is plugged into the Meir-Wingreen formula for the current and conductance (shaded). Note that it quantitatively captures the small-and large-bias features, but qualitatively misses various physical features in the intermediate bias regime.
First we explore the zero-bias linear conductance for small bias voltage. Due to the Fermi liquid nature of the low-energy fixed point, the T = 0 linear response conductance is determined by the total phase shift, which in turn is set by the Friedel sum rule [43,44], leading to the relation: with G 0 = e 2 h the conductance quantum, and n d the total dot occupation in equilibrium, which can be calculated by equilibrium NRG. We show our results for V = 0.01Γ ∆, Γ, U in Fig. 4(a) (circles) vs. NRG (shaded). At gate voltage ε 0 = 0 the system is particle-hole symmetric, and the impurity is occupied exactly by one electron, such that the linear conductance vanishes. At ε 0 ≈ ∓ U 2 the dot population is close to 1 ± 1 2 respectively, hence the linear conductance features Coulomb blockade peaks with height G 0 (red dot).
Next the NESS current at finite bias is shown in Fig. 4(b) for two gate voltages ε 0 = 0 and U 2 , corresponding to the valley and the peak (in the zero-bias conductance). For ε 0 = U 2 the low-bias behavior exhibits a linear dependence (shaded red), as expected. For ε 0 = 0 however, the linear response term vanishes by symmetry, and as the current is an odd function of the bias voltage, the next term is expected to be cubic in the bias voltage, as is indeed observed (shaded blue). In the limiting regime of very large bias V , exceeding all other energy scales (except for bandwidth D), the current saturates for both cases to the maximal value of 2πΓ, directly corresponding to the two conduction channels at coupling strength Γ each. The differential conductance is shown in Fig. 4(c), with peaks corresponding to conductance channels opening up. For ε 0 = 0 we get very clear peaks, with the first conductance channel opening at ∆ with sequential tunneling and thus fluctuations between the two dot levels, and the second conductance channel opening at U , corresponding to full charge fluctuations in the dot occupation. For ε 0 = U 2 , a single-particle level is midway between the two chemical potentials, so there is already a single channel fully open at zero-bias, resulting in a differential conductance of G 0 . The conductance starts dropping close to V = Γ to about half its value. In the vicinity of V = U there is a shoulder, beyond which the conductance drops to zero since the current saturates.
As an interesting comparison consider an approximate form of the Meir-Wingreen formula for the-steady state current [45]. The exact version of this formula reads for the I2LM where G R,A (ω) are the exact retarded and advanced impurity nonequilibrium Green's functions (2 × 2 matrices for the 2 impurity modes), and Γ (ω) = Γ 1 1 1 1 , corresponding to symmetric and equal hybridization of both modes to the two leads (for details see Appendix B). There is of course no known expression for the nonequilibrium Green's functions. However, one could calculate the equilibrium (V = 0) spectral function, e.g., by fdm-NRG [46], and plug it into Eq. (24). This approximation is valid in the linear response regime, and is expected to produce quantitatively good results for large bias, but in the intermediate regime is uncontrolled. Fig. 4 therefore also shows the steady-state current (b) and differential conductance (c) obtained in this manner (shaded). We see the the equilibrium spectral function results agree quantitatively with our nonequilibrium results in the low-and large-bias limits. Interestingly, for ε 0 = 0 they also capture the charge-fluctuation peak at V = U in the intermediate region, but only hint at the level fluctuation peak at V = ∆. On the other hand, for ε 0 = U 2 , they completely miss the shoulder in the drop of the conductance. Thus we conclude that the RL-NESS method successfully reproduces the current and conductance in the known limits, but also gives numerically convergent results in the intermediate regime.

V. DISCUSSION
To conclude, in this work we have introduced RL-NESS, a new numerically exact algorithm for finding the steady-state of general impurities far from equilibrium. It builds on the power of equilibrium NRG in addressing equilibrium quantum impurities with widelyseparated bare and emergent energy scales, and brings it into the nonequilibrium realm. The method is based on coherently coupling the impurity to appropriately loglinearly discretized leads, which in turn are subject to weak Lindblad driving representing incoherent reservoirs. This model setup corresponds to the physical picture of, e.g., a quantum dot coherently coupled to quantum wires, which are in turn coupled to a classical voltage bias source. The resulting system is numerically simulated by a combination of NRG reduction of the high-energy degrees of freedom, followed by tDMRG-based MPDO Lindblad evolution. We benchmark our approach by presenting results for both noninteracting and interacting models. The accuracy of these demonstrate the power of our method, accompanied with a detailed analysis of all error sources and their treatment.
One can envision different ways to try to improve the algorithm. Having shown that an efficient representation of the steady-state as a tensor-network exists, it would be useful to search for more compact representations. One candidate for such a representation is the locally purified tensor-network (LPTN) ansatz [37,39], which enforces physical constraints on the density operator such as positivity, and as such resides in a smaller manifold. However it is not guaranteed that such an ansatz will efficiently capture the entanglement structure of the steady-state [47], as preliminary investigation seems to suggest for the case at hand. So-called disentanglement schemes for the ancilla index [48] might improve the situation, but require further investigation. Recent works [49,50] claim that the entanglement structure of the chain geometry is not optimal, suggesting that applying time evolution in the star geometry might result in a slower growth of entanglement entropy, thus requiring a smaller bond dimension. Testing this idea together with RL-NESS is left for future work. Another direction which might be worth investigating is directly solving the Lindblad equation Lρ = 0 for the steady-state [51,52], instead of obtaining it by real-time evolution.
It would be interesting to apply RL-NESS to more complicated models, such as the single impurity Anderson model [2], the interacting resonant level model [53], and the I2LM with non-symmetric coupling [54,55], all of which are expected to demonstrate Kondo-like physics. RL-NESS already incorporates a temperature for each lead, and so can immediately be employed for finite temperature calculations, as well as calculating thermal conductance, by assigning a different temperature to each lead. We also plan to leverage the success of RL-NESS in obtaining a stable steady-state solution, in order to extract dynamical properties, i.e., time correlators and spectral functions. In the longer run, we envision the treatment of far-from-equilibrium higher-dimensional correlated quantum systems, using, e.g., the dynamical mean field approach [56][57][58].
We thank M. C. Bañuls, J. I. Cirac, J. Eisert, F. Schwarz and A. Werner for helpful discussions. ML thanks F. Schwarz for making her code, developed for Ref. [29], available to him at the initial stages of this project. This joint work was supported by the German Israeli Foundation As specified in Sec. II.B, the intervals I n are chosen such that in the range [−D * , +D * ] they are of size δ and far from this range they scale exponentially as ∼ Λ n . This choice of intervals is achieved by defining a function f (x) for positive x, which is linear for x < D * δ and has a smooth transition to exponential ∼ Λ x for large |x|: with n running on all integers such that the full band is covered up to the cutoff D. The edge of the last interval is then manually fixed to be D. The parameter z ∈ [0, 1) is referred to as the z-shift parameter (in the NRG jargon), and can be used to shift the lead energy levels.
Since different z-shifts result in different yet equivalent discretizations, it is common practice to average simulations using different z-shifts in order to reduce numerical artifacts due to the discretization [59], especially when calculating spectral functions. However in this work it was sufficient to use z = 0. The intervals for negative energies are taken as a mirror image of the positive intervals. This guarantees particle-hole symmetry for any z, at the cost of the interval closest to 0 not necessarily being of size δ. In each interval a representative energy level is selected, with its energy ε n chosen as the arithmetic mean of the interval boundaries in the linear sector (below D * ) and the geometric mean in the logarithmic sector (above D * ). The coupling Hamiltonian is then integrated over each interval in order to derive the appropriate coupling v αnλ of the new lead level with the impurity λ level:

Appendix B: Exact Solution of the Continuous Noninteracting Case
The exact solution for a quadratic continuous system can be calculated in the Keldysh formalism. For noninteracting leads, all the effects of the couplings to the leads on the impurity are encoded in the hybridization function, defined between the λth and νth impurity levels (λ, ν ∈ {1, . . . , m}) for lead α ∈ {L, R} as: where v αnλ are the coupling constants between the λth impurity level and the nth energy level of lead α (in the star geometry). In the case of continuous leads, the sum over dense levels ε n is understood as an integral over the energies. The total hybridization function is then defined as a sum on the hybridization functions of all leads: For a quadratic dot Hamiltonian H, the retarded and advanced Green's functions of the dressed impurity are then given by: where H and Γ are understood here to be m × m matrices. The NESS current can be obtained by the Meir-Wingreen formula [45], which for equal hybridization functions Γ L (ω) = Γ R (ω) = Γ (ω) /2, simplifies to: where f α (ω) is the lead specific Fermi-Dirac distribution. The Keldysh Green's function equals: (B.5) and the impurity single-particle density matrix can then be obtained by integrating over it: Specifying a box hybridization function for the λth level with half bandwidth D: and taking coupling constants v αnλ = v which are lead, n and λ independent (so that all the elements of Γ (ω) are equal), as is indeed the case for both models under investigation in the continuum limit, Eq. (B.4) simplifies to: which can then be evaluated for any desired bias voltage. The Keldysh Green's function in Eq. (B.5) also simplifies to: resulting in a simple integral for the single-particle density matrix.
impurity current. Furthermore, for quadratic systems, this matrix encodes all information about the state of the system, so that finding P (t) amounts to fully solving the system. In the case of a quadratic Lindblad equation (both in the Hamiltonian and the dissipative terms), the exact evolution, as well as the steady-state solution, can be reduced to a continuous Lyapunov equation for P .
The key parts of this reduction are derived in this appendix, following Ref. [31]. We start from the most general Lindblad equation for fermionic Lindblad operators {c q }: where Λ (1,2) encode the Lindblad driving rates. The time dependence of a general single-time observable A (t) ≡ tr (Aρ (t)) is then given by: Assuming a quadratic Hamiltonian H = mn H mn c † m c n , and substituting A = P rs into Eq. (C.2), results in a differential continuous Lyapunov equation for P : The general solution of this equation, for some initial condition P 0 , is: In this section we analyze the time scales of the current evolution for the RLM with ε 0 = 0, after discretization in the linear-logarithmic scheme, and with energy independent Lindblad driving γ, as discussed in Sec. II.B. The time scale of the initial rise in the current can be characterized by t 1/2 , the time at which the current first reaches half of its final value. This time scale appears to be inversely proportional to the bias voltage V , as can be seen in Fig. 5(a) where t 1/2 · V is of order unity over the full range of explored bias.
The time scale of the decay towards the steady-state can be extracted for a quadratic model from the matrix A, defined in Eq. (C.3). The (negative) real parts of the eigenvalues of this matrix dictate the decay rate of each mode. The ones with the smallest magnitude set a bound on the total decay rate of the system. For sufficiently small Lindblad driving, the imaginary part of the eigenvalues depends mainly on the Hamiltonian, while the real part will depend on the driving rates. Thus for the RLM in the discussed discretization scheme, the real part of the eigenvalues naturally scales with γ, as can be seen in Fig. 5(b) for several choices of γ, and the decay rate is proportional to γ, with the proportionality constant of order 1.

Appendix E: MPDO Singular Value Spectrum
In this section we discuss the dependence of the longtime limit steady-state singular value spectrum of the MPDO on the Lindblad driving rate γ. As an example we plot in Fig. 5(c) the singular value spectrum, taken at the bond connecting the RI to one of the linear sector leads (as indicated in the cartoon), for the RLM with parameters as given in the caption. First note that while the normalization of the wavefunction constraints the squared singular values of an MPS to sum up to 1, the density operator normalization condition does not impose any constraint on the MPDO singular values. Thus the global prefactor is arbitrary, and for clarity the singular values are rescaled such that the largest singular value for each γ is 1. As can be seen in the figure, the singular values decay at a faster than power-law rate, implying an efficient representation of the steady-state as an MPDO with finite bond dimension χ. Moreover, we observe that the decay rate grows monotonically with increasing γ, implying that larger γ requires a smaller bond dimension in order to efficiently represent the state of the system.
A full characterization of the exact functional dependence of the singular values λ j on the index j and the system parameters requires a more detailed analysis than carried out in this work. We do note however, that we can fit it to a log-Gaussian behavior λ j ∝ e −(a log j+b) 2 , with a and b the fitting parameters. We suspect that RI Figure 5. (a) The time t 1/2 (multiplied by V ) at which the current reaches half its final value, plotted for the RLM with ε0 = 0 over a wide range of bias voltages. (b) The distribution in the complex plane of the eigenvalues of the matrix A = −iH − Λ (1) − Λ (2) , defined in Eq. (C.3), for the RLM with several choices of driving rates γ. The imaginary parts, mostly corresponding to the Hamiltonian, are plotted in units of Γ, while the real parts, which are related to the Lindblad driving, are rescaled by γ. The closely bunched points near Re (s) = −γ correspond approximately to the single-particle energies of the Hamiltonian arising from the linear sector. (c) Example of the long-time limit steady-state singular value spectrum of the MPDO bond connecting the RI to one of the linear sector leads (as indicated by the red line in the cartoon). The spectrum was obtained for the RLM with V = 1, δ = 0.025V, χ = 512, R = 32 and several values of γ. The singular values were rescaled such that the largest singular value for each γ is 1, and fitted to a log-Gaussian (solid line). (d) NESS current and differential conductance of the RLM as a function of bias voltage V and at different gate voltages ε0, as calculated with (pluses) and without (circles) Lindblad driving of the sites enclosed in the RI, and compared with the continuum limit (shaded). The current with and without the driving at the RI is calculated exactly for γ/δ = 2, 1 (with δ = 0.05V ) and Λ = 8, 6, and is then linearly extrapolated to γ → 0, Λ → 1.
this specific behavior for an MPDO steady-state is not coincidental, since a similar behavior has been argued to occur for an MPS ground-state [60]. We further observe that the fitting parameter a, which dictates the decay rate, is monotonic in γ and goes to zero in the γ → 0 limit. Thus in this limit the required bond dimension χ diverges. This is to be expected, as the steady-state corresponds to evolution to infinite time without dissipation, and we get the well known exponential growth in entanglement entropy for unitary evolution.
Appendix F: Driving RI Sites Fig. 5(d) demonstrates that Lindblad driving of the RI itself has a negligible effect on the resulting current and differential conductance, with respect to an exact solution (which is attainable for the RLM). As argued in Sec. III, this is because the RI represents energy levels far from the voltage or temperature bias scales. These levels are not expected to be affected by the nonequilibrium conditions and thus only set the (renormalized) stage for the low-energy dynamics. Moreover, the exact solution of the modified Lindblad equation (without driving the RI) is still a valid approximation for the continuous system in the limits Λ → 1, γ = δ → 0. This justifies turning off the driving for the interacting case, thus suppressing numerical artifacts arising due to the interplay between NRG and the dissipative dynamics.