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 small 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 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 successfully converges to the steady state at finite bond dimension due to the introduced dissipation, which bounds the growth of entanglement. We thoroughly test the method against the exact solution of the noninteracting resonant level model. We then demonstrate its power using an interacting two-level model, for which it correctly reproduces the known limits, and gives the full $I\text{−}V$ curve between them.

Figure 1

Overview of the RL-NESS method (for more details see Figs. 2 and 3). (a) Two leads with continuous spectra and at a temperature or voltage bias are replaced by a discrete set of energy levels coupled to Lindblad baths. This discrete set is chosen by logarithmic discretization (red) above the bias scale and linear discretization below it (blue). (b) The system is mapped to a tight-binding chain with a local Lindblad bath coupled to each site. (c) High-energy modes are integrated out by equilibrium NRG on the logarithmic sector, resulting in a renormalized impurity (RI). (d) The system is evolved in real time toward the long-time limit steady state, using a dissipative variant of tDMRG for a matrix-product density operator.

Figure 2

Schematic description of the RL-NESS method leading up to the point of solving the Lindblad equation. (a) 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 the modes in that interval. The discrete set of energy levels together with the impurity thus form a finite system. The remaining lead modes are integrated out, resulting in a Lindblad bath coupled to each discrete energy level. (${\mathrm{b}}^{\prime }$) The width of the intervals is chosen according to a logarithmic-linear discretization scheme, such that levels in the low-energy window $\left[-D^{*},+D^{*}\right]$ are equally spaced by $\delta$ (blue), with a smooth transition to logarithmic spacing $\sim {\mathrm{\Lambda }}^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. (${\mathrm{c}}^{\prime }$) 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 $\gamma$ and chosen such that they drive the level toward 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 $\gamma$. (d) The particle and hole leads can be exactly mapped onto nearest-neighbor 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 $\delta$ 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. (${\mathrm{e}}^{\prime }$) This is achieved by collecting all the logarithmic sector sites into a single Wilson chain (via a retridiagonalization) 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.

Figure 3

(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 $\chi$, 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. (b) A single second-order Trotterized time step consists of a dissipative half-time-step sweep (b.1), a Hamiltonian full time-step sweep (b.2: sweep forward and backward at $\tau /2$), and then another dissipative half-time-step sweep in the opposite direction [cf. Eq. (16)]. (b.1) Dissipative evolution half-time-step sweep: each red square corresponds to the set of Kraus gates applied to a specific site, and the sum on all gates at that site is implied by the red contraction line (see text). At the RI, multiple gates are sequentially applied, corresponding to the multiple driven sites incorporated into it. (b.2) Hamiltonian evolution time step: each red bar corresponds to a half-time-step Trotter gate $e^{\pm i{H}_{k,k+1}\tau /2}$ [cf. Eq. (17)] applied to sites $k,k+1$. The gates are applied from left to right and then back from right to left for a full time step.

Figure 4

(a.1) NESS current and (a.2) differential conductance of the RLM (in units of the conductance quantum $G_0=e^2/h$) as a function of bias voltage $V$ and at different gate voltages ${\varepsilon }_0$, as calculated exactly (solid lines) and by RL-NESS (circles). Simulation results are obtained for $\delta =0.1D^{*},D^{*}=\frac{V}{2}$ and after linear extrapolation in $\gamma /\delta =2,1\to 0$ and $\mathrm{\Lambda }=8,6\to 1$, with bond dimension $\chi =256$ and $R=32$ states kept in the RI. (b.1) Comparison of the RL-NESS MPDO evolution ($\chi =256$), after extrapolation to $\gamma \to 0$ (blue), and the NRG-tDMRG MPS evolution ($\chi =1024$) for $\gamma =0$ (red), both with the same discretization, at bias voltage $V=\mathrm{\Gamma }$ and ${\varepsilon }_0=0$. Exact results for both cases are also plotted (shaded). The MPDO result is then further linearly extrapolated from $\mathrm{\Lambda }=8,6$ to $\mathrm{\Lambda }=1$ (solid black) and compared to the exact continuous leads (CL) result (shaded gray). (b.2) Truncation error (sum of singular values squared accumulated over several time steps) plotted for the RL-NESS MPDO with $\gamma =\delta$ (blue) and the NRG-tDMRG MPS (red) evolution. (c) Discretization error ${\epsilon }_{\mathrm{disc}}$ for the RLM at ${\varepsilon }_0=0$ as a function of $\gamma$ for different values of $\mathrm{\Lambda }$, with $\delta =0.0005V$ (solid) and $\delta =\gamma$ (shaded). The inset displays the $\mathrm{\Lambda }$ dependence of the lower bound on the discretization error ${\epsilon }_{\mathrm{disc}}$ (corresponding to the circles in the main figure).

Figure 5

RL-NESS simulation results for the interacting two-level model (I2LM) with level spacing $\mathrm{\Delta }=\frac{\mathrm{\Gamma }}{5}$ and interaction $U=5\mathrm{\Gamma }$, marked by dashed lines. The simulation is run with discretization parameters $\mathrm{\Lambda }=3,\delta =0.1D^{*},D^{*}=\frac{V}{2}$, linearly extrapolated to $\gamma \to 0$ from $\gamma /\delta =4,2,1$, and with simulation parameters $\chi =256,R=64$. (a) Zero-bias linear conductance as a function of gate voltage ${\varepsilon }_0$, calculated by RL-NESS at $V=0.01\mathrm{\Gamma }$ (circles) and compared to the numerically exact equilibrium NRG result (shaded). (b.1) Comparison of the RL-NESS MPDO evolution ($\chi =256$), after extrapolation to $\gamma \to 0$ (solid), and the NRG-tDMRG MPS evolution ($\chi =1024$) for $\gamma =0$ (shaded), both with the same discretization, at ${\varepsilon }_0=0$ (blue) and $\frac{U}{2}$ (red), and bias voltage $V=\mathrm{\Gamma }$. (b.2) Truncation error (sum of singular values squared accumulated over several time steps) plotted for the RL-NESS MPDO with $\gamma =\delta$ (solid) and the NRG-tDMRG MPS (shaded) evolution. (c) NESS current and (d) its derived differential conductance at finite bias, as calculated by RL-NESS for two gate voltages, corresponding to the valley at ${\varepsilon }_0=0$ (blue) and the peak at ${\varepsilon }_0=\frac{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 differential conductance (shaded). Note that it quantitatively captures the small- and large-bias features, but qualitatively misses various physical features in the intermediate bias regime.

Figure 6

(a) The time $t_{1/2}$ (multiplied by $V$) at which the current reaches half its final value, plotted for the RLM with ${\varepsilon }_0=0$ over a wide range of bias voltages. (b) The distribution in the complex plane of the eigenvalues of the matrix $A=-iH-{\mathrm{\Lambda }}^{\left(1\right)}-{\mathrm{\Lambda }}^{\left(2\right)}$, defined in Eq. (C3), for the RLM with several choices of driving rates $\gamma$. The imaginary parts, mostly corresponding to the Hamiltonian, are plotted in units of $\mathrm{\Gamma }$, while the real parts, which are related to the Lindblad driving, are rescaled by $\gamma$. The closely bunched points near $\mathrm{Re}\left(\mathrm{s}\right)=-\gamma$ 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=\mathrm{\Gamma },\delta =0.025V,\chi =512,R=32$ and several values of $\gamma$. The singular values were rescaled such that the largest singular value for each $\gamma$ 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 ${\varepsilon }_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 $\gamma /\delta =2,1$ (with $\delta =0.05V$) and $\mathrm{\Lambda }=8,6$, and is then linearly extrapolated to $\gamma \to 0,\mathrm{\Lambda }\to 1$.