Steady-state and quench dependent relaxation of a quantum dot coupled to one-dimensional leads

We study the time evolution and steady state of the charge current in a single-impurity Anderson model, using matrix product states techniques. A nonequilibrium situation is imposed by applying a bias voltage across one-dimensional tight-binding leads. Focusing on particle-hole symmetry, we extract current-voltage characteristics from universal low-bias up to high-bias regimes, where band effects start to play a dominant role. We discuss three quenches, which after strongly quench-dependent transients yield the same steady-state current. Among these quenches we identify those favorable for extracting steady-state observables. The period of short-time oscillations is shown to compare well to real-time renormalization group results for a simpler model of spinless fermions. We find indications that many-body effects play an important role at high-bias voltage and finite bandwidth of the metallic leads. The growth of entanglement entropy after a certain time scale (proportional to the inverse of Delta) is the major limiting factor for calculating the time evolution. We show that the magnitude of the steady-state current positively correlates with entanglement entropy. The role of high-energy states for the steady-state current is explored by considering a damping term in the time evolution.


I. INTRODUCTION
Over the past decade, experimental control over quantum systems has increased considerably. Possible realizations reach from model Hamiltonians 1,2 using ultracold atoms in optical lattices to experimental setups of nanoscopic devices such as molecular junctions, quantum wires, or quantum dots. 3,4 Many of these systems show remarkable properties, often due to reduced effective dimensionality and many-body interactions. A prominent example is the Kondo effect, 5 which plays an essential role in transport across quantum dots. A theoretical understanding of transport in out-of-equilibrium conditions is highly interesting for applications in nanoelectronics and molecular electronics and even in biological systems.
Electron-electron interactions render the theoretical description of nonequilibrium dynamics one of the most challenging problems in today's condensed matter physics. 6 However, with the advent of efficient numerical techniques to simulate one-dimensional (1D) quantum systems, 7-12 many physical problems are well within grasp of theoretical physicists. Even nonequilibrium setups in regimes where the potential bias is large with respect to the energy scales of the unperturbed systems are now feasible to study. [13][14][15] In this work, we obtain the steady-state charge current of a single interacting quantum dot under voltage bias, modeled by a single-impurity Anderson model (SIAM). 16 This model is commonly discussed in the wide-band limit 17 approximation, tailored towards a universal, low-bias transport description. Here, we extend the discussion to the case of a finite (semicircular) conduction band in the leads, which has not been explored specifically. A particular realization could consist of two one-dimensional leads such as nanowires [18][19][20][21] and a junction between them comprised of a magnetic impurity, i.e., the quantum dot. We use generic one-dimensional tightbinding leads with finite electronic bandwidth which mimic the electronic properties of, for example, carbon nanotubes. 22 In such a device, the electronic density of states (DOS) of the leads would have a bandwidth on the order of 15 eV (Refs. 22 and 23) and effects arising from their specific structure are to be expected when using corresponding bias voltages which are larger than those typically applied in current experiments with nanoscopic devices.
The steady state is obtained by combining density matrix renormalization group 7,11 (DMRG) and time evolving block decimation 9,11 (TEBD) techniques to perform real-time evolution of the system after several different quenches. This technique is known to yield reliable results for a wide parameter range of one-dimensional models [12][13][14][24][25][26][27][28][29][30][31][32][33][34] and to agree with analytical data. 13 We focus on the particle-hole-symmetric point which shows the most pronounced many-body effects. 35 The bias voltage for most of our data is much larger than the equilibrium Kondo temperature (see Sec. V A), so that Kondo correlations should not influence the steady-state current. We show that the same steady-state current is reached independent of the type of quench used and identify quenches which are superior to others when it comes to extracting steady-state data. We investigate quench-induced oscillations in the transients and compare to real-time renormalization group results. We have performed a careful convergence study in all auxiliary numerical and system parameters and found the major limitation to be the truncation of the many-body state space in each iteration. The method is well suited for reaching relevant time scales to study the steady-state current. We find that our approach is capable of yielding unbiased results valid in the thermodynamic limit. Data presented in this work reproduce analytical results in the noninteracting system. In the low-bias region, our results for the current-voltage characteristics agree with previous data (Heidrich-Meisner et al. 14 ). We are able to extend earlier results 14,[36][37][38][39][40][41] to a wider parameter regime and discuss the interplay of finite lead bandwidth and electronic correlations. We find evidence for pronounced many-body effects at high-bias voltages in interplay with finite electronic bandwidths of the leads. 42 Finally, we discuss the role of high-energy states for the steady-state current in low-and high-bias-voltage regimes.
The text is organized as follows: In Sec. II we introduce our model and describe in detail the different quenches to be performed. We present data for the transient response in Sec. III. Results for the steady-state current are presented in Sec. IV where we also outline how to extract steady-state data from time-evolved quantities. We analyze time scales of individual parameter regimes in Sec. V. The role of highenergy states in different bias regimes is discussed in Sec. VI. A detailed convergence analysis is presented in Appendix A.

II. SETUP
In this section, we define our notation for the SIAM which we use to model a quantum dot. We are interested in electron transport across the quantum dot after each of the several quenches to be described in detail in the following. We explain how we calculate the ground state using DMRG and the realtime evolution using TEBD.

A. Single-impurity Anderson model
We consider a model for a quantum dot including charge as well as spin fluctuations described by the SIAM, consisting of an interacting site connected to a bath of noninteracting electrons. We choose a setup where the quantum dot is located in the middle of a one-dimensional chain of tight-binding electrons. The dynamics is governed bŷ (see Fig. 1) where U parametrizes the onsite interaction strength on the quantum dot, t α ,α ∈ {L,R} is the coupling strength between the quantum dot and the left and right leads. Lead α is characterized by intralead hopping t and onsite potential α . Particle-hole symmetry is enforced for all chosen parameters. When needed, the onsite energy of the quantum dot will be denoted by f .
We choose t = 1 and symmetric couplings t L = t R = 0.3162 t [Eq. (1d)] for all simulations, yielding a bandwidth of D = 4 t of the leads and an equilibrium Anderson width 5 of ≡ π t 2 α ρ reservoir (0) = t 2 α t ≈ 0.1 t. We choose t = 1 and symmetric couplings t L = t R = 0.3162 t [Eq. (1d) for all simulations. This yields a bandwidth of D = 4 t of the leads and an equilibrium Anderson width 5 where the reservoir DOS at the chemical potential is denoted by ρ reservoir (μ). We will display all energies in units of (h,k B and e = 1). We restrict ourselves to the zero-temperature case. Real time will be denoted by τ . In Appendix A we show that within the simulation time τ accessible, the finiteness of the leads does not affect our results.

B. Quench preparation
We are interested in the steady-state current 43 of Eq. (1a) under a finite-bias voltage V B . 44,45 Our strategy to obtain the steady state is by quenching the Hamiltonian parameters x 0 = {U,t,t α , α } at τ = 0 from some initial to their final valueŝ H(x 0 ) →Ĥ(x) and evolve an initial state | 0 withĤ(x). | 0 is chosen to be the ground state of the initial Hamiltonian H(x 0 ) in the canonical ensemble at half-filling with total spin projection S z = 0.
It has been shown that the steady state is independent of the quench rate. 30,31 We apply all quenches at an instant of time, i.e., without a ramp. It could, however, be interesting to study the entanglement growth as a function of the quench ramp.
We consider three different quench types (see Fig. 1) which will be explained in detail below. Unless stated otherwise, we choose a system of L = 150 sites with the quantum dot located at site 75. To drive the system out of equilibrium, a bias voltage V B is applied by setting the respective onsite energies of the leads in an antisymmetric fashion to L = − R = V B 2 . For all quenches, the final parameters are x = {U,t = 1, The initial setup is quench-type (QT) dependent (see Fig. 1):

QT I: Hybridization quench to both leads t
For τ < 0 we take x 0 = {U,t,t α = 0, α = ±V B /2}, i.e., no quantum dot-to-leads coupling, but the bias voltage is already applied. We prepare the ground state ofĤ(x 0 ) at half-filling in the left and right leads and a single up electron on the quantum dot. At τ = 0 the tunneling t α is quenched to its finite value. Note that due to the splitting into three disconnected parts (t α = 0), S z is not zero on the quantum dot and on the right lead initially.

QT II: Quenching the bias voltage α = 0 → ±V B /2
At τ < 0, x 0 = {U,t,t α = 0.3261 t, α = 0}. The system is prepared in the ground state | 0 at half-filling with overall S z = 0 zero. At τ = 0 the bias voltage is quenched to its desired value. As compared to QT I, this setup has the advantage that no subsystems with finite values of S z exist in the ground state. Furthermore, correlations between the three 045132-2 regions are already present in the ground state. Note that the initial state is much more complicated than for QT I. This type of quench has also been used by the authors of Ref. 25.

QT III: Quenching the hybridization t R = 0 → 0.3162 t to the right lead
The initial parameters are chosen x 0 = {U,t,t L = 0.3261 t, t R = 0, α = ±V B /2}, and the system is again solved for the ground state | 0 at half-filling. At τ = 0, we quench t R = 0 → 0.3162 t and evolve | 0 with the quenched Hamiltonian.

C. Methods
To prepare the system in the ground state of the initial Hamiltonian, we employ the DMRG (Ref. 7) algorithm in its two-and single-site formulations. Our implementation exploits conservation of spin projection (S z ) and charge (N), which is crucial for obtaining high-precision data. Time evolution is done using the TEBD (Ref. 9) algorithm, within a second-order Suzuki-Trotter decomposition of the propagator where N τ = T δτ is the number of time slices, T is the total simulation time, and δτ the length of a single time step. The operatorsĤ e andĤ o act on even and odd bonds of the bipartite lattice, respectively. Unless stated otherwise, we use a TEBD matrix dimension of χ TEBD = 2000 and a Trotter time step of δτ = 0.05 t −1 . For additional details including studies of convergence in system size L and all auxiliary numerical parameters, we refer the reader to Appendix A.
The calculations carried out in this work set very high computational demands (≈one million CPU hours) and were only possible due to a parallel code [46][47][48] which respects quantum-number (N , S z ) conservation.

III. TRANSIENT RESPONSE
In this section, we present results for the transient current response of the three quenches. We discuss individual bias regimes and identify oscillations in the time evolution of the current which are reminiscent of results for an interacting resonant level model of spinless fermions. We show that QT II leads to much larger initial oscillations than the other two QTs.

A. Low-, medium-, and high-bias regimes
In our simulations, we identify three regimes of bias voltage V B with qualitatively different behavior. Within each regime, the general features of the time evolution of the current are only moderately dependent on interaction strength. For that reason, we first discuss results for U/ = 12 only (see Fig. 2).
For low-bias voltages [V B / ∈ (0,18)], a steady-state current plateau 6,49 is reached within τ ≈ −1 . In a region of medium-bias voltages [V B / ∈ (18,28)], we observe a fast increase in current over a time scale of τ ≈ 0.3 −1 followed by a rather slow decay which, for some model parameters, is too slow to reach a steady-state plateau within accessible simulation times (see below). In the high-bias region [V B / ∈ (28,40)], the time evolution of the current shows a sharp peak followed by fast decrease of the current into a steady-state plateau within τ ≈ −1 .
Our data indicate that within a simulation time of τ = 3 −1 , approximately one particle is transferred from the left reservoir to the right one. As discussed in detail in Sec. IV, all three QTs eventually approach the same steady state, although in quite different manner. QT II, for example, leads to the largest transient current spike, which is one reason for the lower accuracy in determining the steady state for this quench. We also find that quenching the hybridization(s) (QT I or III) yields much cleaner steady-state plateaux as compared to quenching the bias voltage (QT II), which leads to more pronounced oscillations in these plateaux.
In Fig. 3, we plot τ C (U,V ) as a function of interaction strength and find remarkable agreement with rtRG results at higherbias voltages. The period was extracted from the TEBD time evolution data in three ways: (i) by fitting a sine function, (ii) by identifying the dominant Fourier amplitudes, and (iii) by identifying local maxima. These data were combined and their individual uncertainty taken into account. Error bars in Fig. 3 (not shown) would be sharply growing below V B = 25 . In the lower-bias regime, our data are not significant for a reliable extraction of the period.

IV. STEADY-STATE CURRENT
In this section, we present the current-voltage characteristics of quantum dot. We outline a scheme to extract the steady-state current and investigate the dependence on the type of quench used. The current-voltage characteristics in the low-bias region is compared to existing data obtained with other methods. Furthermore, we present a detailed comparison between an interacting and a noninteracting quantum dot for finite as well as infinite lead bandwidth.

A. Extracting the steady-state current
We identify the steady-state current as the mean value of the time-dependent current taken over a suitable time domain [τ S ,τ E ] over which the current shows an almost constant behavior (apart from small oscillations). τ S typically depends on the model parameters and was chosen by hand, and τ E is taken to be the largest time for which simulations yield reliable results (see Fig. 2). Beyond τ E the current becomes numerically unreliable, resulting in an artificially decaying 045132-3 current (see Appendix A for discussion). We find that in most of the parameter regions, the transients have decayed at τ S ≈ −1 . On the other hand, the end point of the plateau strongly depends on the parameter region under consideration. We define it by two distinct measures. One is the time τ (1) E for which the truncated weight [see Eq. (A3)] reaches a threshold of c = 3 × 10 −6 at any bond (marked by + in Fig. 2). The second definition (τ (2) E marked by • in Fig. 2) is given by the time for which two different definitions of the current, namely the expectation value of the current operator [Eq. (A1)] and the time derivative of the particle number [Eq. (A2)], deviate by more than 7 × 10 −4 , the latter being more susceptible to accumulation of errors. Both times are in good agreement with each other and can be combined into an effective simulation Fig. 2). We choose a value of α = 0.1. Results do not depend on this particular choice. It turns out that this procedure is very robust and does also agree with the point at which the TEBD current starts to deviate from the exact time evolution in the noninteracting system (see Appendix A 5).
The steady-state plateaux obtained in this way usually show oscillations and/or small, parameter-and quench-dependent drifts. We quantify the quality of convergence within the plateau region [τ S ,τ E ] by the slope of a linear fit to the current. A large slope indicates that it is not possible to reach the steady state within the given simulation time τ E , i.e., the physical relaxation time is too long or the reached simulation time is too short. This is further discussed in Sec. V. For these parameter values, we can only provide a likely upper bound for the steady-state current, given by the current at the last reliable simulation time. This is justified because we find the current to always decrease as a function of time (apart from small oscillations). Note that although for some of these parameters the current in some QTs may appear converged but is still 045132-4 considered not converged according to our strict criteria. We consider the current to be converged when the relative slope is below a threshold of ≈5 × 10 −2 . Each curve in addition was inspected by hand for convergence. When we consider the steady-state current converged, we estimate its error as three times the standard deviation taken over the data points in the fitting interval [τ S ,τ E ] (plotted as dashed lines in Figs. 2 and 4). This coincides most of the time with the maximal deviation of the time-dependent current from its mean value.
As an important test, we obtained the current-voltage characteristics for the noninteracting case and compared it to analytical results 51 (see Fig. 4), finding excellent agreement (see also Appendix A 5). Another indication for the reliability of the scheme outlined above is that all three types of quenches investigated yield the same steady-state current within the uncertainty. We note that this is not a priori clear since quench-dependent steady states have been reported in different systems. 52 As noted in Appendix A, the position of the plateau is also stable with respect to variations of technical parameters of the simulation. The quality of the steady-state plateau, however, depends strongly on the values of interaction and bias voltage and may be obscured by initial oscillations or shortened at the end by the truncated weight breakdown.
The behavior of the spin current strongly depends on the quench type and it is even identical to zero for QT II. In this respect, the steady-state charge current does not depend on the properties of the spin current since all three quenches yield the same steady state for the charge current. This turns out to be very advantageous since the time scales in the spin sector are much larger than in the charge sector. 53,54 From our calculations, we find QT I and QT III to yield more reliable data for the extraction of the steady-state current than QT II. Reasons for this behavior are (i) the much more pronounced oscillations in the data of QT II which enlarge the statistical uncertainty of steady-state values and (ii) the much higher transient spike in QT II accompanied by a slightly higher initial entanglement and shorter τ E . Entanglement growth is in general parameter dependent and converges towards the same value for all quench types. 54 In the following, we will present steady-state data extracted from QT I and QT III.

B. Current-voltage characteristics
The current-voltage characteristics of the quantum dot for interaction strengths of U/ = 0,4,8,12, and 20 are shown in Fig. 4. We plot data as obtained from QTs I and III (other QTs would give the same results but with larger error bars, as discussed in Sec. IV A). Results for the noninteracting case agree with analytic results for an infinite system. 51 In some regions, only a likely upper bound for the steady-state current can be provided. This region does not lie on the extreme end of the parameter space. It shows nontrivial dependence on U and V B , which is discussed in detail in Sec. V. The current-voltage characteristics has an approximately semicircular shape, with decreasing maximum as a function of interaction strength U . At small bias V B , the current is linear in V B and agrees with the linear response result j lin = 2G 0 V B (see also Fig. 5). At higher bias, it departs from the linear response result. With increasing U , this departure occurs already at smaller bias V B , which can be attributed to an exponential thinning of the Kondo resonance with increasing U .
In intermediate-bias regions, we observe a flattening in the current-voltage curve. The maximum steady-state current is obtained in a bias regime from V B ≈ 15 to 19 . Increasing the interaction from U = 0 to 12 appears to shift the position of the maximum to higher-bias voltages. For larger values of U our data are not significant to conclude on the behavior of the position of the maximum. We find the maximum current to decrease quadratically with increasing interaction strength: j max 1.675 − 0.003( U ) 2 . Note that these features will likely depend on the actual reservoir DOS.
The decrease of the steady-state current for high-bias voltages can be attributed to the diminishing overlap of the DOS of the two reservoirs. 49 Both have a semicircular DOS with a bandwidth of D = 40 . In the wide-band limit, the curves behave similarly inside the low-bias regime but should 045132-5 saturate as a function of V B for higher-bias voltages (see Fig. 6).
We discuss three simple limits. The TEBD results for the current respect the linear response (j lin ) for very lowbias voltages which gives the conductance quantum G 0 . Furthermore, they respect the high-bias voltage band cutoff where the current has to go to zero (here at V B = 40 ) due to diminishing overlap of the DOS of the reservoirs. The third limit is the noninteracting case (nontrivial for the used numerical method), where we obtain perfect agreement with analytical results for the thermodynamic limit.

C. Comparison to previous results
In the low-bias region, results from other techniques are available for the SIAM out of equilibrium. In the following, we discuss our results for various values of interaction strength U together with data previously obtained (see Fig. 5) by diagrammatic quantum Monte Carlo (QMC), 36 fourth-order Keldysh perturbation theory, 37 time-dependent DMRG, 14 TEBD for temperature T = 0 (this work), nonequilibrium functional renormalization group (fRG), 38 nonequilibrium cluster perturbation theory, 39  summation of real-time path integrals, 41 and the linear response result for the Kondo regime j lin = 2G 0 V B . All methods work at or close to zero temperature. Some of the methods use a wide-band limit and others a semicircular reservoir DOS which (for equal ) become comparable in the shown low-bias region (see Fig. 6 for a comparison). The newly obtained TEBD results agree very well with the unbiased dQMC (Ref. 36) and quasi-exact tDMRG (Ref. 14) data. An earlier comparison including more details but fewer techniques is available in Ref. 56.

D. Comparison to a noninteracting device: Identifying correlation effects from the steady-state charge current
To gain further understanding of the role of correlations, we compare the steady-state current of the interacting quantum dot (U , onsite potential f = −U/2) to the one of a corresponding noninteracting (resonant level) device with U = 0 and onsite potential f = − U 2 (see Fig. 6). Data for the resonant level device are obtained analytically. 51 From the plots in Fig. 6, one can see clear differences in the low-bias region between the noninteracting and interacting device for all interaction strengths, which can be attributed to the presence of the low-energy Kondo resonance in the interacting case. For low bias, the Kondo resonance fixes the linear response current to a U -independent constant and causes a higher current than for a noninteracting quantum dot at the same onsite potential. Furthermore, the curvature of the current-voltage characteristics in the low-bias region is negative in the interacting case as compared to positive in the noninteracting system. For larger values of U = 12 and 20 , the negative curvature turns into a positive one in the low-bias region.
For low values of interaction strength (see data for U = 4 ) we observe deviations in both the low-and high-bias regions. For the latter, this hints at possible many-body effects which may also be important in the high-bias regime. Data in the medium-bias region are almost indistinguishable from the noninteracting case. For high values of interaction strength, the picture changes and many-body effects are present in the whole bias regime.
Summing up, we find that effects of interaction are most pronounced in the low-and also in the high-bias regime, where a larger current is obtained than in the noninteracting device.
Because of the small remaining overlap of the DOS of the leads, this larger current may be due to some low-energy spectral weight in the interacting device, consistent with lowenergy excitations observed in Ref. 39 using a nonequilibrium variational cluster approach calculation.

V. DISCUSSION OF TIME SCALES
In the following, we argue that Kondo correlations do not influence the steady-state charge current in the parameter regime under study (large bias V B compared to Kondo scale). However, our simulations show that depending on bias voltage and interaction strength, the steady-state charge current can not always be reached within the simulation time τ E (see Sec. IV A), due to (i) weak spots of the method (i.e., small τ E ) and/or (ii) long physical relaxation times. To obtain insight into physical mechanisms as well as the parameter dependence of the performance of TEBD, also relevant for future studies, it is desirable to disentangle these two effects. We identify parameter regimes with such long physical time scales to be at U + V B > D (low charge-current regime), where we find our method to perform well, as opposed to parameter regimes with high currents, where only smaller times τ E can be reached, as shown in Appendix B.

A. Finite simulation size/time and Kondo correlations
At the particle-hole-symmetric point of the SIAM, Kondo correlations are especially pronounced. In equilibrium they introduce a characteristic energy scale, the Kondo temperature 5 T K which translates into a length scale of the Kondo singlet ξ K , given by Bethe ansatz 57 Due to the exponential dependence on interaction strength, these spin correlations can not fully develop on a finite-size system, 58,59 already for moderate interaction strength. For the parameters used in this work, the equilibrium Kondo correlations have a spatial extent (screening cloud) of approximately ξ K ≈ 50 sites for U = 4 , ξ K ≈ 200 sites for U = 8 , ξ K ≈ 900 sites for U = 12 , and ξ K ≈ 16 000 sites for U = 20 [see Eq. (2)]. These amount to equilibrium Kondo temperatures of T K ≈ 3 × 10 −1 , 9 × 10 −2 , 2 × 10 −2 , and 1 × 10 −3 , respectively.

045132-7
For very small bias voltages V B T K (and large U ), the Kondo effect introduces a large time scale. In this work, however, we focus on parameters for which V B T K . (An exception is U = 4 and V B < 10 , where the Kondo cloud does fit into our finite-size system. 60 ) For the parameter regime under study, recent numeric 61 and analytic 62-64 studies provide strong indications for suppression of the equilibrium Kondo effect.
It is argued in literature that one expects a splitting of the Kondo resonance, possibly a pinning at the lead potentials 65,66 and/or a suppression 61 of the Kondo effect similar to the effect of temperature 4,61,63 or magnetic field. 67 Renormalization group studies concluded that bias voltage is a relevant energy scale in the problem. [62][63][64] Recent results for the electron dynamics in the steady state indicate a splitting of the Kondo resonance away from zero with bias voltage which further supports our observation that the Kondo induced time scale is not relevant for charge transport at large bias voltages. 39 Note that even in the presence of Kondo correlations, charge relaxation should be orders of magnitudes faster than spin relaxation. 53 From our current simulation we made the observation that an initial system with Kondo correlations (to be precise, their finite-size remnants) as in QT II yields the same steady-state charge current (after a short, and different transient regime) as an initial system without them as in QT I. This indicates that in QT II the Kondo correlations are washed away by bias voltage. We thus conclude that although finite-size systems are not able to capture the full equilibrium Kondo singlet, 58 the steady-state transport in the charge channel is not noticeably affected in the parameter regime under investigation.

B. Time scales in the high-bias regime
We find that relaxation times in the model under discussion are strongly parameter (U , V B ) dependent. These relaxation time scales are estimated by the slope of a linear fit to the plateau region [τ S ,τ E ] (see Sec. IV A). In particular, we identify three regions [see Fig. 7 (top)]: region I is characterized by short physical relaxation times and region II exhibits longer relaxation times. Region II overlaps with the regime in which TEBD restricts us to small final simulation times τ E (high steady-state current regime, see Appendix B for discussion). In region II, we did not obtain a converged steady-state current. In region III, the current is small and the maximum reachable simulation time (see Appendix B) was large enough to determine the steady-state current.
We proceed by providing an intuitive single-particle picture of the transition from region I to II in a Hubbard-I-type description [ Fig. 7 (bottom)]. Then, the leads (assuming infinite reservoirs) are described by semicircular bands of bandwidth D, asymmetrically shifted against each other with increasing bias voltage V B . The quantum dot consists of a single (noninteracting) level, located at the single-particle energy − U 2 . We find that the transition occurs when this single-particle level of the quantum dot leaves the overlap region of both lead DOS (blue line, U trans ≈ D − V B ). We conclude that the existence of an appreciable spectral weight in the overlap region of the lead DOS leads to faster relaxation.

VI. ROLE OF HIGH-ENERGY STATES
To study the role of high-energy states during the time evolution we add a damping term to the propagator which gradually reduces the contribution of high-energy states. In Fig. 8, the effects of damping of high-energy modes on the current is visualized. We show results for very low-bias voltage (V B = 2 ) as well as high-bias voltage (V B = 32 ). The different influence of overdamping (dashed lines) on low-bias setups in contrast to high-bias setups yields insight into the role of high-energy states in the two respective cases. In low-bias settings, strong overdamping (here = 10 ) leads to lower current while in the high-bias case it leads to higher current with respect to the true one. This indicates a qualitatively different role of high-energy states for these two settings.
This result can be made plausible by a simple argument. In the case of small-bias voltage (V B t), the dominant energies should be the kinetic ones and neglecting high-energy states amounts to eliminating those with highest kinetic energy. Such states contribute much to the current and neglecting them leads to a lower total current. On the other hand, for very high-bias voltage (V B t), potential energy is expected to dominate. High-energy states are then those with a lot of particles in the 045132-8 high-bias reservoir. Eliminating them reduces the available state space for hopping of particles back to the side of high potential. Therefore, the current is increased due to less back flow. From a technical point of view, such an approach may reduce entanglement growth (the limiting quantity in real-time evolution using matrix product states), thus reducing the required matrix dimensions of the MPS. Using such an ansatz, however, suffers from two drawbacks. (i) On the one hand, we have just seen that high-energy states can be important for the steady-state current and, on the other hand, estimating a priori a suitable magnitude of the damping is not straightforward since it should in principle be dynamically adjusted during the time evolution taking into account energies and truncated weight. Due to these reasons, we refrain from using such an approach in general. However, we show in Fig. 8 that by choosing a phenomenologically good value for the damping ( = ), one can indeed somewhat prolong the stable time evolution.

VII. CONCLUSIONS
We studied the single-impurity Anderson model out of equilibrium beyond the linear response regime by means of density matrix renormalization group. Real-time evolution was performed making use of the time-evolving block decimation algorithm which allows us to access relevant time scales to reach the steady state. Within this framework, we investigated three different quenches: (i) quenching the hybridization with already applied bias voltage, (ii) quenching the bias voltage, and (iii) quenching the hybridization at one side only.
Calculated current-voltage characteristics agree very well with established results which are available in the low-bias region. We find that the period of characteristic oscillations in the time evolution of the charge current is already very well described by renormalization group results for a different model, the interacting resonant level model of spinless fermions. After an initial transient regime, where on the order of one particle is transferred through the quantum dot, the steady-state current agrees among the three quenches investigated. For the identification of steady-state plateaux in time-dependent quantities, the type of quench is however very important. We show that quenching the lead-dot tunnelings is the most suitable one, contrary to expectations whereas quenching the bias voltage results in large initial oscillations of the current. We furthermore show that limitations of the method such as its inherent finite size do not pose a problem for simulations of the setup discussed here within accessible times. Our findings indicate that the steady-state charge current is not influenced by finite-size effects, hinting that incompletely developed Kondo correlations in the spin channel do not influence charge transport noticeably. We find that a large entanglement entropy correlates positively with a large steady-state current amplitude. By studying a damped time evolution, we find that high-energy states have very different significance in the low-and high-bias regimes, respectively.
Aside from reproducing the universal low-bias physics, we open up new perspectives for devices in which a large-bias voltage is combined with a finite electronic DOS of the leads, such as nanotubes. For such devices, we predict that effects of electron-electron interactions are important even at high-bias voltages.
Interesting extensions within the presented approach may be the application of a gate voltage to study stability diagrams, evaluation of spin correlations which could hint on Kondo correlations, to study effects of asymmetric couplings, the interplay of bias and magnetic fields as well as to investigate correlated leads. 34 On the technical side, it would be interesting to evaluate whether more gently ramped quenches over a finite-time interval further decrease oscillations or even entanglement and further improve the extraction of steady-state data.
from numerical precision. In addition, our setup contains leads of finite size, with two effects in principle. First, this finite size affects the ground state at time zero. We will show in the following that the effect on the current is negligible; it converges already at much smaller lead size than used here.
Second, the finiteness of the reservoirs means that no energy or particle dissipation occurs and eventually the system will show oscillatory behavior. We note in passing that during our simulation time only approximately one particle traverses the quantum dot. The earliest time at which the current can be affected by the finite system size arises from a perturbation which propagates after the quench to the end of a lead and back to the quantum dot. The velocity of this signal is limited by the Lieb-Robinson bound, 68 up to exponentially suppressed parts, and in our case is v ≈ 2t, which can also be clearly seen in the time evolution of local charge expectation values. 54 The perturbation will hit the left and right ends of the chain and return back to the quantum dot after a time of about τ ≈ 2(L/2)/(2t) = L/(2t), i.e., τ/ −1 ≈ L /(2t) ≈ 7.5 for L = 150. This is far beyond the times τ E (see Sec. IV A) up to which we calculate the steady-state current, which is therefore not affected. This conclusion is confirmed by the convergence of the current with respect to system size L, discussed below. The measured current may, however, be affected by other possible errors within our approach: (i) the procedure to measure it, (ii) the Trotter error, and (iii) the limited matrix dimension χ TEBD (i.e., truncated weight).
In the following, we will show that the major uncertainty arises from the limited matrix dimension χ TEBD , while other sources are negligible. The definition of the time intervals from which the steady-state current is evaluated (Sec. IV A) is also relevant. A similar conclusion has been drawn before in the framework of adaptive tDMRG (Ref. 25) and for different systems in the framework of TEBD. 49

Obtaining the current
Within the TEBD time evolution, the steady-state current may be obtained via the expectation value of the current operator at each time step 24,69 where i and j denote adjacent sites and a iσ and a † iσ are annihilation and creation operators for fermions onsite i with spin σ which depend at time τ and t ij is taken to be real. To obtain the current through the quantum dot, a symmetrized version of the inflow and outflow is used: where c L endσ , c †L endσ denote operators on the last site of the left reservoir (number 74 in Fig. 1) and c R 0σ , c R † 0σ denote operators on the first site of the right reservoir (number 76 in Fig. 1).
Another way of computing the current is by calculating the time derivative of the total particle number to the left of the site under consideration: Again, a symmetric combination of the dot's ingoing and outgoing current yields the current under consideration The current through the dot may be evaluated at each TEBD time step using Eq. (A1) or by computing a finite-difference approximation to the differential Eq. (A2) every two successive time steps. Aside from the expected additional source of error by evaluating the time derivative numerically, this method is expected to perform less well due to the influence of all sites in the system on the result for the current, the occupation number of each site having its own limited accuracy. A comparison of the current evaluated by means of Eqs. (A1) and (A2) for various values of interaction strength U and applied bias voltage V B as well as all QTs (I, II, III) shows good agreement in the beginning of the time evolution [see Fig. 9 (left)]. Due to an accumulation of errors in the particle-number expectation values of the individual sites, the results start to deviate at some time τ (2) E . We do not use results beyond τ (2) E (see the discussion in Sec. IV A). Numerical values of all steady-state currents will be obtained using the current operator [Eq. (A1)] which yields a much more stable estimator.

Finite-size effects: L
In this section, we discuss the dependence of the results for the current on the length of the system L. 49,70 We quench both dot-lead tunnelings (i.e., QT I). The qualitative behavior for the other QTs (II and III) is virtually identical. Results for the steady-state current for system sizes of L = 20,40,60,80,100,120, and 150 sites are shown in Fig. 9 (right) in the noninteracting case. We find that the final results for the steady-state current agree with the analytically available results for an infinite system in all cases within the numerical error. This ensures a reliable determination of steady-state properties even on finite-size systems. As mentioned before, the system size limits the maximum simulation time due to signals back-propagating from the borders. In the main part of this work, all calculations are performed for a system size of L = 150 to provide a nice long plateau (maximum simulation time) in the steady-state current. It has been noted in Ref. 71 that in the particle-hole-symmetric half-filled model, the steady-state current is independent of system size. A detailed discussion of finite size and time scales in a model of spinless fermions can be found in Ref. 70.
For completeness, we note that it is possible to extend the available simulation time, when it is limited by the hard boundary conditions of the leads, by applying modified boundary conditions. 34,69,72,73 Exponentially decreasing the matrix elements of the Hamiltonian towards the end of the 045132-10 reservoirs ultimately corresponds to a Wilson chain with logarithmic discretization. 73 In this work, we do not apply any modified boundary conditions because our simulation time is not limited by the size of the chains but the TEBD matrix dimension χ TEBD .

Trotter error: δτ
The Trotter error grows only linearly with simulation time, 25,74 and can be controlled by choosing sufficiently small δτ . Therefore, usually the contribution to the total error arising due to the Trotter approximation is negligible with respect to other approximations. We investigated the influence of the Trotter decomposition on the current. Results for δτ/t −1 = {0.01,0.05,0.1} were found to agree to within 5 × 10 −5 . We do not plot the results because they all lie on top of each other. A good value for the time step was found to be δτ/t −1 = 0.05 which was used in the main section of the paper.

MPS matrix dimension: χ
The quality of the TEBD results is predominantly determined by the maximum matrix size χ TEBD used. A bigger χ TEBD leads to fewer discarded states (i.e., less truncated weight of the reduced density matrices) during the truncation and therefore to a systematically better approximation. 12 The truncated weight is defined by 74 where λ 2 γ denote the eigenvalues of the reduced density matrices. This quantity is zero if no truncation is done. The computational cost of the TEBD algorithm scales essentially like 74 where L is the length of the chain and d = 4 the size of the local fermionic Hilbert space. Therefore, it is essential to keep χ TEBD as low as possible. During the simulations we noticed that at a certain time (long before signals propagating back from the ends of the chain would reach the quantum dot), the truncated weight starts to grow quickly and the results become unstable, 75 causing a decaying current. The effects of enlarging χ TEBD are shown in Fig. 9 (center). As the data indicate, the effect of increasing χ TEBD is to make larger simulation times accessible, before the simulation breaks down due to accumulation of truncated weight. Remarkably, no spurious quasi-steady state is entered when χ TEBD is relatively small. The overall shape of the current appears to be unaffected by enlarging χ TEBD , making reliable predictions for χ TEBD = 2000 possible. We checked our results in all parameter regions for convergence. In the main part of the paper, we always used χ TEBD = 2000 as a good compromise between run time and accuracy.

Comparison to analytical results
In the noninteracting setup U = 0, we compare TEBD data to results from an exact time evolution (see Fig. 10). The exact time evolution was obtained for the same system parameters by time evolving the single-particle density matrix. We find that the TEBD time evolution is reliable up to a systemparameter-dependent time. This time (triangles) again is in accordance with the criterion for the maximum reachable simulation time as defined in Sec. IV A and has a nontrivial dependence on bias voltage and interaction strength.
The noninteracting system is nontrivial for the TEBD method. Our data for the entanglement entropy and the truncated weight at low-, medium-, as well as high-bias voltages for increasing interaction strength 54 indicate that indeed the U = 0 case does not exhibit a particular low entanglement or truncated weight in comparison with higher interaction strength. Since we reproduce the exact analytic steady-state current in the noninteracting case, we conclude that the agreement with exact results is not a peculiarity of the noninteracting system and our way of data extraction can be applied to finite values of U .

Setup
Based on the above considerations, all data in the main part of the paper were obtained for the following parameters: (i) The ground state was obtained by DMRG using a matrix size χ DMRG = 400 and undergoing 10 sweeps of two-site DMRG before switching to 40 runs for single-site DMRG. (ii) The model consists of L = 150 sites. Upon performing one of the three above-described quenches (I, II, or III), we used bias voltages in a range of V B / = (0,42). We always started from an overall half-filled system in the canonical ensemble with total S z = 0 and alternating up and down spins are chosen from left to right. (iii) The time evolution was performed using a TEBD matrix size of χ TEBD = 2000, a trotter step of δτ/ −1 = 0.005 and evolving for 1000 time steps which yielded a final simulation time of T / −1 = 5. Requiring a maximum truncated weight of c = 10 −15 , we dynamically adjusted the size of the TEBD matrices with a maximum matrix size of χ TEBD . We measured observables at each time step.

APPENDIX B: CORRELATION OF ENTROPY AND STEADY-STATE CURRENT
The major limiting factor for time evolution using TEBD is the increase of bipartite entanglement 11 whereρ L/R denotes the reduced density matrix to the left (L) and to the right (R) of a bipartition at bond i. Using a maximum matrix dimension χ TEBD , we stop the simulation (for Fig. 11) whenever the truncated weight at any bond exceeds a threshold value of c = 5 × 10 −5 , which defines our maximum simulation time τ (1) E (see Sec. IV A). In Fig. 11, we plot τ (1) E as a function of U and V B (left) and compare it to the magnitude of the steady-state current for the same parameters (right).
From our data we conclude that reachable simulation times due to accumulation of entanglement (and thus truncated weight) are nonmonotonic in U and V B but can be characterized roughly by the magnitude of current in the system. We find this behavior to be generic to all investigated QTs.