Steady-state spectra, current and stability diagram of a quantum dot: a non-equilibrium Variational Cluster Approach

We calculate steady-state properties of a strongly correlated quantum dot under voltage bias by means of non-equilibrium Cluster Perturbation Theory and the non-equilibrium Variational Cluster Approach, respectively. Results for the steady-state current are benchmarked against data from accurate Matrix Product State based time evolution. We show that for low to medium interaction strength, non-equilibrium Cluster Perturbation Theory already yields good results, while for higher interaction strength the self-consistent feedback of the non-equilibrium Variational Cluster Approach significantly enhances the accuracy. We report the current-voltage characteristics for different interaction strengths. Furthermore we investigate the non-equilibrium local density of states of the quantum dot and illustrate that within the variational approach a linear splitting and broadening of the Kondo resonance is predicted which depends on interaction strength. Calculations with applied gate voltage, away from particle hole symmetry, reveal that the maximum current is reached at the crossover from the Kondo regime to the doubly-occupied or empty quantum dot. Obtained stability diagrams compare very well to recent experimental data [Phys. Rev. B, 84, 245316 (2011)].


I. INTRODUCTION
The understanding of non-equilibrium phenomena in strongly correlated many-body systems may reveal previously unknown fundamental aspects of physics as well as prove crucial for the development of technical applications in the fields of nano or molecular electronics. Currently, combined insight of experiments on nanodevices 1-3 and results from artificial quantum simulators (see e.g. refs. [4][5][6][7][8] are capable of providing coherent insight into the non-equilibrium behavior of the quantum world. These experiments provide both a challenge for theoretical concepts as well as an accurate check for theoretical results. Both nano-devices and condensed matter simulators are often described remarkably well by interacting model Hamiltonians which are in general not exactly solvable. Here we focus on a model of a single quantum dot, the single-impurity Anderson model (SIAM) 9 . This model, incorporating spin and charge fluctuations as well as Kondo correlations, has been studied as an idealized realization of an interacting system by a wide array of techniques in equilibrium (for an overview see e.g. Hewson ref. 10). The model we use here is applicable to generic single quantum dot systems including effects of a finite lead bandwidth. However, the evaluation of dynamic quantities of strongly correlated quantum many-body systems out of equilibrium poses a notoriously difficult problem. A particular challenge for the SIAM is that it is expected to remain in a strong coupling regime, even under the influence of a bias voltage 11 . Additional challenges are posed by the need for particle and/or energy dissipation mechanisms. Techniques which allow the calculation of physical quantities beyond the linear response regime are quite restrictive up to now and no full understanding of the non-equilibrium dynamics or the steady-state under bias are available for the SIAM, apart from some special cases [13][14][15] . Nevertheless, controlled results have been obtained with the analytical Bethe Ansatz for some physical quantities in the interacting resonant level model 16 and for the SIAM 17 . Perturbative calculations 18,19 as well as non-crossing approximation studies 20,21 extend early fundamental work on the non-equilibrium problem [22][23][24][25] and the impurity out of equilibrium [26][27][28] . Depending on the setup of parameters and model details, insight may be gained by semiclassical methods 29 or master equation approaches 30 . Recently moreover, several techniques, which have proven very successful in the equilibrium theory, have been extended to the non-equilibrium case. Among them are many-body cluster methods 31,32 , renormalization group (RG) approaches [33][34][35][36][37][38][39] , flow equation methods 40,41 , real time path integral calculations 42 , out of equilibrium noncrossing approximation (NCA) 94 , generalized slave-boson methods 12,96 , diagrammatic quantum Monte Carlo (QMC) [43][44][45] or QMC methods based on a complex chemical potential [46][47][48][49] . The Gutzwiller approximation has been generalized to the time-dependent case 50 and so has numerical renormalization group (NRG) [51][52][53][54] where however some issues with the use of Wilson chains in non-equilibrium systems have been pointed out by Rosch 55 . Dual fermion approaches 56 have been proposed as well as super operator techniques 57,58 . Some recent work attempts to compare several of these theories [59][60][61] and shed light on the critical issue of time scales involved 62 . Finally, some results for the SIAM are available 63 from numerically exact time evolution by a combination of Density Matrix Renormalization Group (DMRG) 64 and successive time evolution via time-dependent DMRG (tDMRG) 63,[65][66][67] . They are currently limited to small bias voltages and moderate interaction strengths. In the present paper we explore the capabilities of non-equilibrium Cluster Perturbation Theory (nCPT) 31 and benchmark the non-equilibrium Variational Cluster Approach (nVCA) 32 on the SIAM. We obtain the full current-voltage characteristics which we compare to results from a very accurate time evolution by means of Time Evolving Block Decimation (TEBD) 68 and QMC 43 . We also comment on the comparison to other results obtained by tDMRG 69 , perturbation theory 18 and FRG 34 . Because of good agreement with these data, we then proceed to evaluate the single-particle spectrum of the quantum dot in the steady-state. Since this is a dynamic quantity, it is even harder to obtain for most numerical methods. We find a linear splitting of the Kondo resonance 85,86 which depends on the interaction strength. Detailed results for the particle-hole symmetric model where Kondo correlations dominate are presented and supplemented by data in the whole parameter space including an applied gate voltage. We highlight the crucial edge which nVCA gains over nCPT via its self-consistent feedback. Both many-body cluster techniques, nCPT and nVCA are based on the well established equilibrium counterparts Cluster Perturbation Theory (CPT) 70,71 and the Variational Cluster Approach (VCA) 72,73 . They make use of the non-equilibrium Keldysh-Schwinger Green's function technique [74][75][76] . The present paper extends and benchmarks the ideas for adapting VCA to non-equilibrium situations introduced in ref. 32. We attempt to give a comprehensive overview of the current capabilities and shortcomings of nCPT and nVCA for the application to steady-state problems of strongly correlated systems. The SIAM provides an excellent probing ground for our purposes as a model where the effects of correlations are crucial. It is however still relatively simple which permits systematic study and some results are available which allow for comparison. Reasonable results for the SIAM in equilibrium have been obtained previously by CPT as well as VCA 77 . Both cluster methods are approximate but they yield fairly reliable results and are therefore interesting due to their great flexibility and versatility. They are computationally not very demanding and allow in principle to treat a wide range of fermionic / bosonic lattice Hamiltonians out of equilibrium. Possible extensions include electronic multi-band or multi-orbital systems in one, two or three dimensions also including phonons. This paper is organized as follows: In sec. II we sketch the SIAM and describe the setup used. We proceed by outlining nCPT (sec. III A) and nVCA (sec. III B). In sec. IV A, a comparison of steady-state currents as obtained by nCPT, nVCA, DMRG and successive TEBD 78 and QMC 43 is presented. The non-equilibrium local density of states (nLDOS) is examined in sec. IV B. Finally effects of an applied gate voltage on the steadystate current are studied in sec. IV C.

II. MODEL OF A SINGLE QUANTUM DOT
We model the setup, consisting of a single correlated quantum dot in-between two metallic leads, by a singlesite Hubbard model embedded in a one dimensional infinite tight-binding chain. The Hamiltonian of this singleimpurity Anderson model (SIAM) 9 reads (see fig. 1) The electronic annihilation (creation) operators c iασ , f σ (c † iασ , f † σ ) obey the usual anti-commutation relations with spin σ = {↑, ↓} and annihilate (create) electrons in the left or right lead α = {L, R} or the quantum dot, respectively. The particle number operator of the quantum dot is defined byn f σ = f † σ f σ . U represents the on-site Hubbard repulsion. The single-particle energy of the quantum dot ǫ f = − U 2 + V G is comprised of a gate voltage V G and a shift by − U 2 (eq. (1b)), such that the particle-hole symmetric point is reached when the gate voltage vanishes and the bias is applied in an antisymmetric fashion as described below. Non-interacting left and right leads are described by tight binding chains with a nearest-neighbor hopping t (eq. (1c)). The lead on-site energies are denoted by ǫ L/R . A bias voltage V B may be applied to the system in an anti-symmetric fashion by setting ǫ L = µ L = −ǫ R = −µ R = VB 2 , where µ L/R denote the chemical potentials of the decoupled leads (t ′ = 0). Finally, the symmetric tunneling between the non-interacting leads and the quantum dot is denoted by t ′ (eq. (1d)). In this work we adopt units in which , e, as well as the inter-lead hopping t are equal to 1. The lead-dot tunneling is fixed to t ′ = 0.3162 throughout this paper, which implies an Anderson width of ∆ ≡ π t ′2 ρ lead (0) = t ′2 t ≈ 0.1. This choice of parameters furthermore implies an effective lead bandwidth of D = 40 ∆. For higher bias voltages we therefore include additional effects due to finite bandwidth and band-shape as compared to the often used wide band limit. The prefactor for the current is then given by j 0 = 2|e|t . In the following, we always consider the zero temperature case. In order to employ the cluster approaches, the system is divided into three pieces, a left part, a right part and a central region of length L, which includes the quantum dot as well as parts of left and right leads.

III. METHODS
All information needed to evaluate expectation values in the steady-state of an interacting fermionic lattice Hamiltonian is contained in a single object, the singleparticle Green's function in Keldysh space G (eq. (A1)) (for details see app. A). Evaluating this object is in general impossible to do exactly but a handle on G is given by nCPT and nVCA. These approximations are to be discussed in the following.

A. Non-equilibrium cluster perturbation theory
The idea of nCPT is to split an (infinitely) large system H for times τ < τ 0 into a set of small decoupled clusters (described by the Hamiltonianĥ), for which the singleparticle Green's functions can be determined exactly, by numerical meansĤ At time τ 0 the coupling of the individual subsystemŝ T =Ĥ −ĥ is switched on and the solution for the singleparticle Green's function of the total system can be obtained by the CPT equation 70 in Keldysh-, site-and spinspace where G/ g denote the single-particle Green's function in Keldysh space of the total system/split systems. The unit matrix in Keldysh space is denoted by11. Eq. (2) is a strong coupling perturbation theory result and holds up to first order in the inter-cluster hoppingT, as far as the self-energy is concerned. As a consequence, the self-energy ΣĤ of the total system is approximated by the self-energy of the initial system Σĥ. This implies that nCPT/nVCA become exact in the non-interacting limit, since in this case, the self-energy of the total system agrees with that of the initial system, which are both zero. The results of nCPT and nVCA converge towards the exact results with increasing cluster size.
In the present paper we illustrate this approach for the SIAM out of equilibrium. We start out by splitting the infinite chain (eq. (1a)) into three parts: i) an interacting central region of length L consisting of the interacting site (the quantum dot) as well as an equal amount of sites of the left and the right leads, ii) a semi-infinite, noninteracting left region consisting of the remaining part of the left lead and iii) a semi-infinite, non-interacting right region consisting of the remaining part of the right lead (see fig. 1). To proceed, the single-particle Green's functions in Keldysh space have to be determined exactly for those three systems. Results for the left and right part are available analytically by the retarded Green's function of a semi-infinite tight binding chain 82 where υ i,j is the retarded Green's function of the infinite tight binding chain. The single-particle Green's function of the central interacting region can be determined in the Q-matrix formalism by the Band-Lanczos algorithm (see ref. 77 and ref. 83). The advanced component can always be obtained by the relation G A = (G R ) † . Before coupling the three subsystems at time τ = τ 0 , each of them is in equilibrium at different chemical potentials µ L/R/C . Therefore we may evaluate the corresponding Keldysh components by the relation 84 where p FD denotes the Fermi-Dirac distribution function at inverse temperature β: p FD = 1 1+e β(ω−µ) . In the zero temperature limit, (1 − 2 p FD (ω, µ)) may be re-expressed as sign(ω − µ). This is the only expression in which the chemical potential µ enters. It is crucial that this relation does not hold in a non-equilibrium situation any longer. The hermitian part of G K (ω, µ) is zero in equilibrium. The imaginary part consists of contributions due to delta peaks for finite size systems like the central regions. The operatorT just contains the two hopping terms from the left to the central and from the central to the right region. At time τ 0 , the hopping processes between these three regions are switched on and the steady-state single-particle Green's function G is determined using nCPT, eq. (2). As in the equilibrium case, the CPT results can be improved by the variational cluster approach, in which the single-particle part of the initial system is suitably modified. In the following a variational extension (nVCA) of the scheme described above will be presented following ref. 32.

B. Non-equilibrium variational cluster approach
The idea of enlarging the size of the central region in nCPT is to find a better approximation for the starting point of perturbation theory. As pointed out before, in CPT the self-energy of the total system is approximated by that of the decoupled system. By increasing the size of the central region, the self-energy of the decoupled system converges gradually towards that of the total system. In nVCA an even more suitable initial system is chosen. This can be achieved by making use of the fact that the decomposition ofĤ into clustersĥ and inter-cluster partŝ T is not unique. We are free to add single-particle opera-tors∆(x), which depend linearly on parameters x, to the initial hamiltonianĥ, provided we subtract them again fromTĤ where∆ l is quadratic in the fermion operators. This parametrized single-particle field introduces additional freedom in the starting guess for perturbation theory and can be used for a self-consistent feedback on the clusters. The system described by ĥ +∆(x) , usually referred to as reference system in the context of VCA, has the same structure as in equilibrium VCA. The condition to fix the single-particle parameters x, however, is different in nVCA and VCA. Equilibrium VCA is based on the Selfenergy Functional Approach (SFA) 87,88 and provides a variational principle for the generalized grand potential functional, which is not well defined in the context of nonequilibrium systems any longer. An alternative criterion for fixing the variational parameters x was introduced in ref. 32 and compared to the traditional VCA criterion in ref. 77. It is based on the idea of starting from a system which is as similar as possible to the total, original system in terms of physically observable quantities. We demand the expectation values of the operators∆ l to coincide in the initial (reference) system and the steady-state of the total system, i.e.
For example, adding variational freedom in the on-site energy of the quantum dot, corresponding to ∆ = x ǫ fn f σ yields the self-consistency condition: n f σ g ! = n f σ G . From a more conceptual point of view, it is interesting that these self-consistency conditions follow from the condition 32 which is closely related to the stationarity condition of the generalized grand potential functional Ω in the equilibrium approach 89 . Here τ 1 = 011 110 is a Pauli matrix in Keldysh space and the subscript zero denotes noninteracting Green's functions. From the numerical point of view, one has to find the roots of an n-dimensional set of non-linear equations (where n is the number of variational parameters x). Like in first order Dirac perturbation theory, the influence of the perturbation increases with time (perturbations introduced due to the coupling with hopping elements t at τ 0 are proportional to t · (τ − τ 0 ) since we are considering first order perturbation theory) and one might argue that nCPT is bound to fail in the long-time, steady-state limit, even for small couplings between the clusters 31 . This argument cannot be true in general, as nCPT yields exact results in the case of non-interacting particles, although the initial decoupled systems are far from the steady-state behavior.
The initial system in nCPT is independent of the nonequilibrium situation and this shortcoming is improved to some extent within nVCA, where the information about the applied bias voltage is self-consistently fed back to the initial reference system. As we will see in this paper, there are circumstances under which nCPT already yields reasonable results. In general, however, we observe that nVCA represents a significant improvement over nCPT. This implies, on the one hand, that under steady-state conditions, the self-energy in the central cluster is significantly modified when going from nCPT to nVCA. On the other hand, it indicates that the steady-state situation can be mimicked in an equilibrium system by the auxiliary one-particle terms, determined self-consistently. A point for improvement of this approach is to modify the self-energy so that it is a genuine non-equilibrium one. Details on the nVCA procedure and the particular choice of variational parameters are provided in app. B.

IV. RESULTS
In the following we compare nCPT and nVCA data for the steady-state current with TEBD 68 and QMC 43 results in the particle-hole symmetric model. We elaborate on the non-equilibrium density of states and finally discuss results for the steady-state current away from particle-hole symmetry by applying a gate voltage. Earlier preliminary results obtained for the SIAM out of equilibrium by nCPT and nVCA are available in ref. 92 and ref. 93.

A. Steady-state current
Here we investigate the steady-state current-voltage characteristics of the particle-hole symmetric SIAM. In this parameter region, Kondo correlations are especially important. All nCPT and nVCA results are for an infinite system using a self-energy based on an L ≤ 11 site interacting reference system. Although the length of Kondo correlations scales exponentially in interaction strength, it has been shown in ref. 77 that the VCA approximation, similarly to other approximate methods such as FRG or Gutzwiller wave functions, is capable of retaining qualitatively most features of the ongoing Kondo physics, although, possibly, with renormalized scales. The steady-state current can be evaluated on any link either within the central region or between the central and the neighboring regions, since we find that the continuity equation is fulfilled within at least 10 −6 relative to the steady-state current amplitude. We note again that in the non-interacting case nCPT as well as nVCA become exact so we do not show these data explicitly.
Results obtained by eq. (A2) for the nCPT case are shown in the left column of fig. 2 for various values of the interaction strength U . The linear response result j lin = 2 G 0 V B , where G 0 is the conductance quantum, is fulfilled within nCPT (and nVCA) in the very low bias region. This is because these methods fulfill the Friedel sum rule in equilibrium 77 . Actually, nCPT tends to stay in the linear response regime longer than nVCA and therefore overestimates the current for low (but not very low) bias voltages. We model the leads by one dimensional tight binding chains, yielding a semi-circular density of states (eq. (4)) of bandwidth D = 40 ∆. This implies a vanishing current in the high bias limit at V B = 40 ∆, due to non-overlapping lead density of states. This limit again is trivially fulfilled within nCPT and nVCA. Note that in a wide band limit the current curves would approach a constant roughly at their respective maxima in the data shown. For comparison, numerically exact results 68 as obtained by a real-time evolution with TEBD are also indicated in fig. 2. Note that the TEBD data are obtained from the steady-state plateau in the time dependence of the current. At small and medium bias, the current converges well and the TEBD results provide an accurate benchmark. At higher bias, convergence is low, but the TEBD data still provide a reasonably reliable upper bound for the steady-state current. In the intermediate bias region it is interesting to investigate the behavior of nCPT with increasing size of the central region L. As can be seen from the plots, for any interaction strength increasing L yields mono-tonically improving results. While for low interaction strength U = 4 ∆ the nCPT results almost coincide with the TEBD data, greater deviations arise at higher interaction strengths. For very large interaction strength (see for example U = 20 ∆), the lengths of the central region L considered here are not sufficient. At high bias voltage some spurious finite size effects of the reference system are visible in the form of peaks in the steady-state current.
Next we would like to discuss the performance of nVCA. An illustration is shown in the mid column of fig. 2. Here we compare nCPT and nVCA T for L = 3, 7 and 11 again for increasing interaction strengths U/∆ = 4, 8, 12 and 20. We show a zoom to the low bias region where benchmarking data from other techniques, like TEBD data 68 , tDMRG 69 , FRG 34 , and QMC data obtained by Werner et al. 43 in the wide-band limit are available. One aspect to note immediately is that nVCA, like TEBD, departs sooner than nCPT from the linear response data, which is due to the better reproduction of the Kondo resonance within nVCA. The thinning of this narrow resonance at zero bias with increasing interaction strength is responsible for the departure from the linear response result 57 . Since the Kondo resonance is better accounted for in nVCA, the curves leave the linear response data sooner than the respective nCPT results. The TEBD data 68 and QMC data obtained by Werner et al. 43 in the wide-band limit are also plotted and serve as a benchmark up to V B ≈ 5 ∆ before effects of a different lead density of states become important. Curves obtained by tDMRG 69 / FRG 34 are not shown but lie basically on top of the TEBD / QMC data. Early data from perturbation theory 18 are known to have some additional spurious bumps in the low bias region. The improvement due to the variational feedback in nVCA becomes clear by comparing to the nCPT results for the same L, especially at higher interaction strengths. The right column of fig. 2 shows a comparison of the performance of different variational parameters considered within nVCA for L = 7. For small interaction strength all nVCA parameter sets work equally well. With increasing interaction strength, however, the different variational parameters predict a different behavior of the steady-state current. It should be noted here that some parameters like the ones used in nVCA t b cause almost no deviation from the nCPT result while others improve it appreciably, like e.g. nVCA T . It predicts a two peak structure in the high bias voltage regime, which however is not observed directly in the TEBD benchmark data. No final conclusion can be drawn about the behavior of the current in this bias regime because TEBD can only predict upper bounds for the steady-state current there. The position of the dip in the steady-state current is at V B = U for all interaction strengths. We note that this is where the bare level position of the quantum dot ǫ f = − U 2 is about to stop overlapping with one of the lead density of states (while still overlapping with the other). The calculation for two independent variational parameters (nVCA t,t ′ ) yields similar results as nVCA T in the lower bias region but goes to zero quickly for high bias voltages and does not show a dip.
On the whole, one may say that the self-consistent feedback implemented within nVCA provides a significant improvement over the nCPT results. This has been already observed in ref. 77. However, there this result was deduced from the convergence of results with increasing cluster size and not with a benchmark comparison with alternative numerical methods. For low interaction strengths U 4 ∆, bare nCPT already performs very well (independent of L) while for larger interaction strengths the variational improvement of nVCA becomes important. Motivated by the success of nVCA for the steady-state current, we proceed by evaluating the nonequilibrium local density of states of the quantum dot.

B. Non-equilibrium local density of states
The nLDOS in the quantum dot ρ f (ω) as obtained by nCPT and nVCA T is depicted in fig. 3 for three sizes of the central part of the reference system: L = 3, 7 and 11 and a large interaction strength of U/∆ = 12 (corresponding to the steady-state current in fig. 2 (third row)). We plot the nLDOS in a density plot as a function of energy ω (horizontal) and applied bias voltage (vertical). The equilibrium result, consisting of a thin Kondo resonance at the Fermi energy (ǫ F = 0 here) and two broad incoherent peaks located at ≈ ± U 2 with a width of ≈ 2 ∆, can be inferred from a horizontal cut at V B = 0. For finite bias voltage, a splitting of the Kondo resonance is observed in both nCPT (top row) and nVCA T (bottom row). It is well known that the noncrossing approximation (NCA) predicts a splitting of the Kondo resonance into two under voltage bias and that within second order perturbation theory the resonance is not split but suppressed only. 18 A linear splitting and slight broadening of the Kondo resonance with increasing bias voltage is proposed e.g. in ref. 94. Intuitively it is expected that the split Kondo resonances pins at the chemical potentials of the leads. Several other methods yield a splitting with different features: real time diagrammatic 95 and scaling methods 97 as well as the equation of motion technique 20,21,98-101 . Within fourth order perturbation theory the Kondo resonance splits into two, which are located near the chemical potentials of the two leads 18 .
In experiments on nano-devices a linear splitting of the Kondo resonance has been observed 85,86 . Such a linear splitting is also predicted by nVCA (see fig. 4), while nCPT yields a splitting which shows a roughly quadratic dependence on V B . In addition within nVCA one observes an interaction dependent splitting (see fig. 5). In contrast to the prediction of ref. 94 we do however not observe a simple pinning of the Kondo resonance at the chemical potentials of the leads. Our results indicate that the position ω K of the split Kondo resonance depends on the interaction strength U : fig. 5). Here V m is the voltage where the Kondo resonance merges with the high energy part of the spectrum located at ω H ≈ ± U 2 : V m = 15, 17, 21 and 28 ∆ for U = 4, 8, 12 and 20 ∆. The U dependent values of V m have been determined from the respective linear fit to the data. For high voltages, these split peaks merge with the Hubbard bands and saturate, which has also been observed in fourth order perturbation theory calculations 18 . In this bias region a new low energy excitation is observed for L > 3 within nVCA T . This additional peak in the nLDOS has a dominant contribution to the total weight and is responsible for the two-peak structure observed in the nVCA T steady-state current (see fig. 2).

C. Finite gate voltage: steady-state current and stability diagram
It is interesting to investigate the steady-state current away from the particle-hole symmetric point and the region where Kondo correlations are present. In this section the current through the quantum dot under bias is analyzed as a function of an applied gate voltage V G and applied bias voltage V B at fixed interaction strength U . Results are presented for L = 3. In this case no two-peak structure has been found in the nVCA T current in fig. 2, which is corroborated by fig. 3 (bottom left). The dependence of the current on the gateand bias voltage as obtained by nCPT and nVCA T is depicted in fig. 6 for interaction strengths U/∆ = 4 and U/∆ = 20. First we discuss the so called stability diagram (differential conductance as a function of biasand gate-voltage) which is shown in the third column of fig. 6. The dashed squares mark the region which is typically accessed experimentally. In the nVCA T result, the Kondo region (vertical line in the center) is reproduced extremely well, as are the Coulomb blockade regions (∝ V B ) above and below. Our data compare very well to recent experimental data of ref. 2 (figure 5 therein). Effects at higher bias voltage arise from the finite bandwidth used here and are typically not seen in experiment.
The steady-state current as a function of bias and gate voltage is shown in the first and second column of fig. 6. It is interesting to observe that the largest  current is obtained exactly at the crossover points from the Kondo to the empty or doubly occupied quantum dot (these regions are marked by black-dashed lines in middle and right panel of fig. 6). This aspect can be seen more clearly in fig. 7, where the current is plotted as a function of gate voltage. For U = 4 ∆ there is not much difference between the nCPT and nVCA results as can be inferred from fig. 6. For U = 20 ∆, however, we see that the feedback mechanism in VCA has a significant impact and leads to smoother j − V B curves due to suppression of finite size effects originating from the reference system. The sharp jumps at the crossover point in fig. 7 originate from abrupt changes of the particle number in the ground state of the central cluster L = 3. We expect that this step smoothens with increasing L. Concerning the reliability of the methods used, our data suggest that outside the parameter region where Kondo correlations dominate the single occupied quantum dot (in between the horizontal dashed black lines in fig. 6) nCPT and nVCA perform almost equally well. Results are significantly easier to obtain outside the Kondo plateau (which was mentioned before e.g. in ref. 77 within exact diagonalization/CPT/VCA or in ref. 63 within DMRG). Thereby the convergence within system size is greatly enhanced with respect to the Kondo region and very accurate results may be obtained already with small systems. Therefore, we may argue that for the steady-state of the SIAM, nCPT and nVCA perform quite well outside the Kondo region for any interaction strength as well as in the Kondo region for small interaction strengths. On the other hand, in the Kondo region nVCA outperforms nCPT for higher interaction strengths.

V. CONCLUSIONS
We have presented results for the steady-state of the single-impurity Anderson model. We have applied non-equilibrium Cluster Perturbation Theory and its variational extension, the non-equilibrium Variational Cluster Approach to this model for a single quantum dot under bias. Both methods make use of the Keldysh Green's function formalism and are capable of working in the thermodynamic limit which is necessary to account for particle and energy dissipation mechanisms. Results for the particle-hole symmetric model, which is dominated by Kondo correlations, have been compared to Time Evolving Block Decimation and quantum Monte Carlo. At low values of interaction strength they show excellent agreement already for non-equilibrium Cluster Perturbation Theory. For higher values of interaction strength, the self-consistency implemented within the non-equilibrium Variational Cluster Approach proves crucial in order to obtain reasonable results. Both methods coincide with the low bias linear response data for the steady-state current. Both methods furthermore become exact in the non-interacting limit. The non-equilibrium local density of states of the quantum dot exhibits a linear and interaction dependent splitting of the bias voltage within the non-equilibrium Variational Cluster Approach which is not visible in the non-equilibrium Cluster Perturbation Theory. At a certain (interaction dependent) bias voltage we find that this split Kondo resonance merges with the high energy incoherent part of the spectrum. When applying a gate voltage and thereby leaving the Kondo regime, calculations become a lot easier and non-equilibrium Cluster Perturbation Theory and the non-equilibrium Variational Cluster Approach appear to perform very well, which can be inferred from the convergence of our data. The highest current amplitude is obtained at the crossover from the Kondo to the empty or doubly-occupied quantum dot. Experimental stability diagrams are reproduced very well within the variational approach. They show a clear Kondo regime and a Coulomb blockade region. We may conclude that the non-equilibrium Variational Cluster Approach is a promising method for the evaluation of steady-state quantities of strongly correlated model systems. Dynamic quantities become available in the whole complex plane and, in principle any fermionic or bosonic lattice model may be treated including multi-band or multi-orbital systems. It is both fast and versatile and in principle all parameter regimes of the model are accessible. Subjects for further investigation are the self-consistency criterion, which is not unique, as well as the use of different variational parameters and an optimized out of equilibrium variational principle.
the reference system are also subjected to these on-site energies. All calculations are done for sizes of the central region of L = 3, 7 and 11-sites, which corresponds to symmetric central regions. We note that L = 5, 9 suffer from a finite-size gap which closes with increasing L. This gap arises due to an even amount of sites to the left as well as to the right of the quantum dot. We do not consider even L since this would correspond to an unequal amount of sites of the left and right lead in the central region. Thereby, the geometric symmetry of the problem would be spoiled and the application of a bias voltage which respects particle-hole symmetry would not be possible. In order to reduce finite size effects of the calculated quantities and to make an extrapolation to infinite central cluster sizes easier we scrutinized the idea of averaging over boundary conditions as outlined in ref. 90 or ref. 91. In CPT/VCA one has the freedom to add single-particle terms toĥ and subtract them again in T. Therefore we may add hopping terms which connect the first and last site of the central cluster and which have an arbitrary phase. The spectral function of the impurity site for a central region with complex periodic boundary conditions is in general no longer particle hole symmetric, which makes it necessary to average G φ (ω) over φ to retain the symmetry. The result, however, was disappointing as the expected faster convergence towards the thermodynamic limit has not been confirmed. In this work we do not consider the long-range part of the Coulomb interaction which leads to additional charging effects in real devices.