Master equation based steady-state cluster perturbation theory

A simple and efficient approximation scheme to study electronic transport characteristics of strongly correlated nano devices, molecular junctions or heterostructures out of equilibrium is provided by steady-state cluster perturbation theory. In this work, we improve the starting point of this perturbative, nonequilibrium Green's function based method. Specifically, we employ an improved unperturbed (so-called reference) state $\hat{\rho}^S$, constructed as the steady-state of a quantum master equation within the Born-Markov approximation. This resulting hybrid method inherits beneficial aspects of both, the quantum master equation as well as the nonequilibrium Green's function technique. We benchmark the new scheme on two experimentally relevant systems in the single-electron transistor regime: An electron-electron interaction based quantum diode and a triple quantum dot ring junction, which both feature negative differential conductance. The results of the new method improve significantly with respect to the plain quantum maste equation treatment at modest additional computational cost.

A well-established method for treating such open quantum systems is by means of quantum master equations (Qme). [75][76][77][78][79][80] Herein, the environment-degrees of freedom are integrated out and usually incorporated in a perturbative manner. The Qme approach allows a detailed investigation of transport phenomena 45,46 and recent selfconsistent extensions attempt to cure some of its longstanding limitations. 81 In the framework of nonequilibrium Green's functions (NEGF) various schemes exist to approximately calculate the electronic self-energy of the correlated region, see e.g. Ref. 68,82-87. In cluster approaches, such as cluster perturbation theory (CPT) and its improvement, the variational cluster approach (VCA), 88 the whole system is partitioned into parts which can be treated exactly and determine the self-energy. Originally devised for strongly correlated systems in equilibrium, 89,90 both approaches have recently been extended to nonequilibrium situations in the time dependent case 91,92 as well as in the steady-state. 92,93 In previous work we applied the steady-state CPT (stsCPT) to obtain transport characteristics of heterostructures, 93 quantum dots [94][95][96] and molecular junctions 97,98 and obtained good results even in the challenging Kondo regime. 49,94,95 A key issue in the CPT approach is to identify an appropriate many-body state for the disconnected correlated cluster in the central region, as a starting point of perturbation theory, the so-called reference state. Up to now, a common choice in stsCPT is to use an equilibrium state at some temperature T S (often T S = 0) and chemical potential µ S in-between the values of the leads. Such an ad-hoc choice is clearly unsatisfactory. Furthermore, it fails to describe certain quantum interference effects in transport phenomena as for example so-called current blocking effects. 45,46,97 The purpose of the present work is to improve on stsCPT by constructing a consistent and conceptually more appropriate reference state, given by the steadystate reduced many-body density matrixρ S obtained from a Qme in the Born-Markov approximation. Within this quantum master equation based stsCPT (meCPT), the ambiguity in defining µ S and T S for the central region is resolved. The equilibrium case, in which µ S and T S coincide with those of the environment, is automatically included. In contrast to standard Qme approaches, lead induced level-broadening effects are accounted for and the noninteracting limit is reproduced exactly, as in the original stsCPT. In addition, meCPT is able to capture the previously mentioned current blocking effects, as shown below.
Other NEGF/Qme hybrid methods exist in the literature. 67,99,100 For instance, in a recent work 101,102 we have proposed a so-called auxiliary master equation (AME) approach, whereby a Lindblad equation is introduced which models the leads by a small number of bath sites plus Markovian environments. The AME is suited to address steady-state properties of single impurity problems as encountered in the framework of nonequilibrium dynamical mean field theory. 72,87,101,[103][104][105] In contrast, the meCPT presented in this work is more appropriate to treat non-local self-energy effects which cannot be captured by single-site DMFT.
This paper is organized as follows: After defining the model Hamiltonian in Sec. II, the meCPT is introduced in detail in Sec. III. We present results obtained with the improved method for two experimentally realizable devices: i) In Sec. IV A, an electron-electron interaction based quantum diode, ii) and in Sec. IV B, a triple quantum dot ring junction which both feature negative differential conductance (NDC).
For ring systems, extensive Qme results and an explanation of the NDC in terms of quantum interference mediated blocking are available in Ref. 45,46.

II. MODEL
We consider a model of spin-1 2 fermions, having in mind the electronic degrees of freedom of a contacted nano structure, heterostructure or a molecular junction. The Hamiltonian consists of three parts: i) The "system"Ĥ S represents the interacting central region i.e. the nano device or molecule consisting of single-particle as well as interaction many-body terms. It is described by electronic annihilation/creation operators f iσ /f † iσ at site i = [1, . . . , N S ] where N S is typically small and spin σ = {↑, ↓}. 106 We will specify the particular form ofĤ S in the respective results section. ii) The "environment" HamiltonianĤ E describes the two noninteracting electronic leadŝ where c λkσ /c † λkσ denote the fermion operators of the infinite size lead λ with energies ǫ λkσ and electronic density of states (DOS) ρ λσ (ω) = 1 N λ → ∞ are the number of levels in the leads. The disconnected leads are held at constant temperatures T λ and chemical potentials µ λ so that the particles obey the Fermi-Dirac distribution p FD λ (ω, T λ , µ λ ). 106,107 iii) Finally the system and the environment are coupled by the single-particle hoppinĝ

III. MASTER EQUATION BASED CLUSTER PERTURBATION THEORY
Our goal is to obtain the steady-state transport characteristics of the HamiltonianĤ, Eq. (1a) in a nonequilibrium situation induced by environment parameters like a bias voltage V B or temperature gradient ∆T . The important step consists in evaluating the steady-state single-particle Green's function in Keldysh space G in the well established Keldysh-Schwinger nonequilibrium Green's function formalism. [108][109][110] In generalĤ is both interacting and of infinite spatial extent. Therefore explicit evaluation of G is prohibitive in all but the most simple cases which motivates the introduction of approximate schemes.
One such scheme is CPT, 89,90 in which one performs an expansion in a 'small' single-particle perturbation, for example the system-environment couplingĤ SE of Eq. (1c). The unperturbed HamiltonianĤ S +Ĥ E can be solved exactly. While in the noninteracting case CPT becomes exact, results obtained in the presence of interaction are approximate and depend on the reference state for the unperturbed system. A common practice within stsCPT 92,94-98 is to use a pure state given by the equilibrium ground state |Ψ 0 S of the disconnected interacting system HamiltonianĤ S . In a nonequilibrium situation, this is still ambiguous, as it depends on an arbitrary choice of the chemical potential µ S and/or temperature T S for the interacting finite system.
The goal of this work is to provide an unambiguous and conceptually more rigorous criterion for the choice of the reference state for the interacting central region. Ideally, the reference state is selected such that it resembles best the situation of the coupled system, i.e. for the full Hamiltonian, Eq. (1a) in the steady-state. An appropriate choice in equilibrium is to use the grand-canonical density operator 106ρS gc as reference state. In this case, T S and µ S are uniquely determined by the equilibrium situation. Equivalently,ρ S gc is given by the steady-state solution of a Qme in the Born-Markov approximation (see Sec. III B), when coupling the system to one thermal environment. From this viewpoint a natural extension to the nonequilibrium situation is to make use of a Qme as well in order to obtain a consistent reference state, given then by the steady-state reduced density operator of the systemρ S . In this work, a second order Born-Markov Qme is employed, which yields the correct zeroth order reduced density operatorρ S (adjusted toĤ SE ). 111,112 Subsequently,Ĥ SE is included within the CPT approximation, 89,90 in order to obtain improved results for the Green's function and in turn for the transport observables.
In summary, the meCPT method consists of the following three main steps, analogous to a standard CPT treatment: 1. Decompose the whole system into a small interacting central region (system) and noninteracting leads of infinite size (environment), seeĤ S andĤ E in Eq. (1a).
2. The new step introduced in this work is to solve a Qme for the system in order to obtain the reduced density operatorρ S , which serves as a reference state to calculate the cluster (retarded) Green's function 113 3. Reintroduce the system-environment couplingĤ SE perturbatively, see Sec. III A and Eq. (4), to determine the Green's function of the coupled system.

A. Steady-state cluster perturbation theory
Here we briefly recall the main, well-established CPT concepts and equations, as this is the starting point for the formalism presented in this work. For an in depth discussion of CPT 114 and its nonequilibrium extension we refer the reader to the literature. 91,93,95,97 The central element of stsCPT is the steady-state single-particle Green's function in Keldysh space 115 where R denotes the retarded, A the advanced, and K the Keldysh component. In the present formalism, G R/A/K become matrices in the space of cluster sites and depend on one energy variable ω since time translational invariance applies in the steady-state. As explained above, in order to compute G(ω) within stsCPT one partitionsĤ, Eq. (1a) in real space, into individually exactly solvable parts, in this case, the systemĤ S and the environmentĤ E , which leaves the coupling Hamiltonian H SE as a perturbation. The single-particle Green's function of the disconnected Hamiltonian is denoted by g(ω), which obviously does not mix the disconnected regions. For the noninteracting environment, the respective block entries of g(ω) are available analytically. 94,116 For the interacting part the respective entries of g(ω) are calculated via the Lehmann representation with respect to the reference state. This can be computed e.g. based on the Band Lanczos method. [117][118][119] The full steady-state Green's function in the CPT approximation is found by reintroducing the inter-cluster coupling perturbatively where we denote by the matrix M the single-particle Wannier representation ofĤ SE . CPT is equivalent to using the self-energy Σ of the disconnected Hamiltonian as an approximation to the full self-energy. Therefore, the quality of the approximation can in principle be systematically improved by adding more and more sites of the leads to the central cluster. However, in doing so the complexity for the exact solution of the central cluster grows exponentially. Independent of the reference state, this scheme becomes exact in the noninteracting limit.

B. Born-Markov equation for the reference state
In the following we outline how to obtain the reference stateρ S by using a Born-Markov-secular (BMsme), or more generally a Born-Markov master equation (BMme). [75][76][77][78][79][80] Although this approach is standard, for completeness we present here the main aspects and notation. We loosely follow the treatment of Ref. 40,78,79. The real time τ evolution of the full many-body density matrixρ is given by the von-Neumann equatioṅ ρ = −i Ĥ ,ρ − . 75 Typically the large size of the Hilbert space ofĤ prohibits the full solution in the interacting case. One thus considers the weak coupling limit |Ĥ SE | ≪ |Ĥ E | and performs a perturbation theory in terms of |Ĥ SE |. 80,120 In the usual way one obtains an equation for the reduced many-body density matrix of the systemρ S (τ ) = tr E {ρ} by working in the interaction pictureρ I (τ ) = e +i(Ĥ S +Ĥ E )τρ (0)e −i(Ĥ S +Ĥ E )τ with respect to the coupling Hamiltonian, Eq. (1c). One then performs three standard approximations: i) Within the Born approximation, valid to lowest order in |Ĥ SE |, the density matrix is factorizedρ I (τ ) ≈ρ S I (τ ) ⊗ρ E I . Furthermore, the environmentρ E I is assumed to be so large that it is not affected by |Ĥ SE | and thus independent of time. ii) The Markov approximation implies a memory-less environment, that is, the system density matrix varies much slower in time than the decay time of the environment correlation functions C αβ (τ ). Upon transforming back to the Schrödinger picture this yields the BMme, which is time-local, preserves trace and hermiticity, and depends on constant coefficients. iii) To obtain an equation of Lindblad form which also preserves positivity one typically employs the secular approximation, which averages over fast oscillating terms, yielding the BMsme. 78,121,122 The system-environment coupling can be quite generally written in the formĤ SE This hermitian form is convenient for further treatment.The tensor product form can be achieved even for fermions by a Jordan-Wigner transformation, 79 see App. B. For our coupling Hamiltonian, Eq. (1c) and particle number conserving systems, the coupling operators take the form In the energy eigenbasis of the system Hamiltonian H S |a = ω a |a , the BMme in the Schrödinger representation reads 113 with where ω ba = ω b − ω a . The Lamb-shift HamiltonianĤ LS and the environment functions ξ αβ (ω 1 , ω 2 ) are defined in App. A. When employing the secular approximation, the terms in the BMsme simplify and in Eq.
Due to the secular approximation the BMsme can only lead to interference between degenerate states. The more general BMme also couples non-degenerate states at the cost of loosing the Lindblad structure of the Qme, see Sec. IV B and Ref. 40.

Single-particle Green's function
As discussed above, for meCPT, the Green's function g(ω) of the isolated system is evaluated from the reference stateρ S . The retarded component Eq. (2) takes the explicit form where i, j denote indices of system sites. The advanced component follows from g A = g R † and the Keldysh component g K of the finite, unperturbed system is not relevant for the CPT equation, Eq. (4). Once g is obtained, the full Green's function is again approximately obtained within CPT by Eq. (4). Notice that for U = 0, G is independent of the reference state, which is why stsCPT, stsVCA as well as meCPT coincide (and become exact) in the noninteracting case.

C. Numerical implementation
From a numerical point of view, the two main steps are to first obtain the reference stateρ S by solving the Qme and then to evaluate the Green's functions using Eq. (8) and Eq. (4). For the solution of the BMme, Eq. (6) one needs to carry out the following: i) Full diagonalization of the interacting system Hamiltonian which is done in LAPACK, making use of the block structure inN andŜ z . ii) Evaluation of the coefficients of the BMme in Eq. (6), which involves coupling matrix elements a|Ŝ α |b and numerical integration of the bath correlations functions, see App. A, C, for which an adaptive Gauss-Kronrod scheme is employed. iii) The steady-stateρ S is finally obtained from the unique eigenvector with eigenvalue zero of Eq. (6), which we determine by a sparse Arnoldi diagonalization. Again, a block structure is related toN and S z . The numerical effort for the exact diagonalization scales with the size of the Hilbert space, and therefore exponentially with the system size N S . In the second major step, the Green's function of the disconnected system is calculated by Eq. (8). Finally, the meCPT Green's function G(ω) is found using Eq. (4). We outline how to evaluate observables within meCPT and the Qme in App. D.

IV. RESULTS
In this section we present results obtained from the meCPT approach. In all calculations, except the ones in Sec. IV B, the secular approximation is applied for the reference stateρ S . The main improvements of meCPT with respect to bare BMsme are i) the inclusion of lead induced broadening effects, ii) the correct U = 0 limit and iii) a correction for effects missed by an improper treatment of quasi degenerate states in the BMsme (see below). In comparison to the previous "standard" stsCPT, meCPT also captures current blocking effects, which are discussed in detail in Ref. 39 The leads are in the wide band limit and at the same temperature T .

A. Quantum dot diode
We first discuss a quite simple model system: a quantum diode based on electron-electron interaction effects. Fig. 1 depicts this junction consisting of a single interacting orbital described by a Hubbard interaction and an on-site term to allow for a gate voltage V G : 123 The environment Eq. (1b), consists of two spin dependent, conducting leads. We model both, the left (L) and the right (R) lead by a flat DOS with local retarded single-particle Green's function 116 g R L/R (ω) = − 1 2D ln ω+i0 + −D ω+i0 + +D , with a half-bandwidth D much larger than all other energy scales in the model, mimicking a wide band limit. We keep both leads at the same temperature T L = T R = T and at chemical potentials µ L = −µ R = VB 2 corresponding to a symmetrically applied bias voltage V B . The right lead is fully spin polarized, i.e. tunnelling of one spin species (↓) into the right lead is prohibited while both spin species can tunnel to the left lead. The system is coupled to the two leads via a single-particle hopping amplitude t ′ inĤ SE , Eq. (1c) which results in a lead broadening parameter of , and Γ ↓ R ≡ 0. We use Γ without an argument for Γ(ω = 0) as defined in Eq. (C1). For meCPT we use H SE , see Eq. (1c), as perturbation.
Such a system could be realized in: i) A "metal -artificial atom -half-metallic ferromagnet" 124 nano structure where spin-↑ DOS is present at the Fermi energy while the respective spin-↓ DOS is zero. ii) A graphene nano structure 31,32 with ferromagnetic cobalt electrodes. 33 iii) A one dimensional optical lattice of ultra cold fermions in a quantum simulator 125 where the hopping of spin-↓ particles into the right reservoir is suppressed. For all three systems spin-↓ particles cannot reach the right lead, in the first two due to a vanishing DOS, in the third one due to a vanishing tunnelling amplitude.
We consider parameters such that the junction is operated in a single electron transistor (SET) regime, 37 i.e. temperatures above the Kondo temperature. 49 In this regime we expect an interaction induced -magnetization mediated blocking due to the fact that the system fills up with spin-↓ particles. On the one hand they cannot escape, yielding a vanishing spin-↓ current, and on the other hand they suppress the spin-↑ occupation, at finite repulsive interaction U , resulting also in a vanishing spin-↑ current. 46 Fig . 2 (A) shows the meCPT stability diagram of the interacting system in the V B − V G plane. When applying a particle-hole transformation for all particles, leads and system, along with t ′ → −t ′ we easily find the symmetry properties From the continuity equation it is clear that only spin-↑ steady-state current can flow which limits the maximum current to Γ 2 . The energies ω N of the isolated quantum dot can be labelled by the total particle number N and are for V G = U given by ω 0 = 0, ω 1 = 1 2 U and ω 2 = 2 U . This gate voltage corresponds to the dashed line, marked by (X) in Fig. 2 (A). The corresponding energy differences ∆ 01 = 0.5 U between the single-occupied and the empty dot and ∆ 12 = 1.5 U between double-occupied and single-occupied dot are associated with a further transport channel opening as soon as the bias V B reaches twice their values. The meCPT result for the current exhibits the well known Coulomb diamond 37 close to V B = 0 and V G = 0, where current is hindered because all system energies are far outside the transport window ± VB 2 , see Eq. (D2). At V G = 0 a current sets in at VB 2 = ±| U 2 |, i.e. when transport across the system's single particle level becomes allowed. The point, at which the current sets in, shifts with V G linearly to higher bias voltages. This transition is broadened ∝ max(Γ = 0.1 U, T = 0.025 U ). However, not only the transport window and possible excitations in the system energies determine the currentvoltage characteristics. The particular occupation of the system states may lead to more complicated effects, such as current blocking.
Our first main result is that in contrast to stsCPT the blocking is correctly reproduced in meCPT. The current blocking is visible in Fig. 2 (A) in region (Y), see also the detailed data in subplot (C1). It is asymmetric in V B and therefore responsible for the rectifying behaviour for |V G | > | U 2 |. This feature is easily understood from the plots of the spin resolved densities in Fig. 2 (C2). In the region of interest, for positive V B , n ↓ = 1 which hinders spin-↑ particles from the left lead to enter the system, due to the repulsive interaction U and suppresses the current. For negative V B , the situation is reversed. A direct computation of the current in the framework of BMsme, see App. D 2, also predicts the blocking, which is however not the case if we use stsCPT based on the zero temperature ground state |Ψ 0 S . The blocking is evident in Fig. 2 (B), where we observe that in the blocking regime, the reduced density is ρ S = |↓ ↓|. Independent of the value of U > 0, the blocking sets in at the same values of V B in meCPT and BMsme. Fig. 2 (C1) shows that within BMsme this regime is entered after a U independent hump in the current while within meCPT the hump is broader and weakly U dependent. The current blocking disappears at a bias voltage V B ∝ U in both methods. Immediately apparent are the much broader features in meCPT, which leads to a less pronounced effect in contrast to the total blocking predicted by BMsme. In BMsme the broadening parameter Γ enters merely as prefactor of the current, and broadening is solely induced by the temperature. This temperature induced broadening is correctly taken into account in both methods. For T > Γ the lat-ter dominates and the meCPT results are similar to the plain BMsme solution. A comparison of the three methods is given in Tab. I. In this simple model the blocking can be captured even by a straight forward steady-state mean-field theory in the Keldysh Green's function with self-consistently determined spin densities or in stsVCA. This is not the case for the more elaborate system studied in the next section. In this section we discuss a more elaborate model system: a triple quantum dot ring junction which features negative differential conductance (NDC) based on electron-electron interaction effects mediated by quantum interference due to degenerate states as outlined in detail in Ref. 45,46. Fig. 3 (A) depicts the triple quantum dot ring junction, described by the following Hubbard Hamiltonian 126 In addition to the model parameters described in Sec. IV A, a nearest-neighbour ij hopping t is present. The environment, Eq. (1b) and coupling, Eq. (1c) are now both symmetric in spin. Moreover, we use µ L = −µ R = VB 2 , T = T L = T R and Γ L = Γ R = Γ 2 = π|t ′ | 2 1 2D . Such a junction can be experimentally realized: i) Via local anodic oxidation (LAO) on a GaAs/AlGaAs heterostructure 29 which enables tunable few electron control. 30 ii) In a graphene nano structure. 31,32 Experimentally the stability diagram has been explored 34 alongside characterisation and transport measurements. 29,35,36 The negative differential conductance has been observed in a device aimed as a quantum rectifier. 127 Theoretically the study of the nonequilibrium behaviour of such a device has become an active field recently. 45,46,[128][129][130][131][132][133][134] We investigate transport properties for values of the parameters such that the junction is in a single electron transistor (SET) regime, 37 i.e. temperatures above the Kondo temperature. 49 In this regime we expect an interaction induced -quantum interference mediated blocking as discussed in Ref. 45,46. The rotational symmetry ensures degenerate eigenstates labelled by a quantum number of angular momentum. In situations where these degenerate states participate in the transport they provide two equivalent pathways through the system and lead to quantum interference. 45 The blocking sets in at a bias voltage, where the degenerate states start to participate in the transport. It then becomes possible that a superposition is selected which forms one state with a node at the right lead. In the long time limit this state will be fully occupied while the other one will be empty due to Coulomb repulsion, for reasons very similar to the ones discussed in the previous section. 39,40 The steady-state charge distribution and currentvoltage characteristics of the interacting triple quantum dot are presented in Fig. 3 (B, C) in a wide bias voltage window. The current, depicted in panel (C), in general increases in a stepwise manner and is fully antisymmetric with respect to the bias voltage direction. A blocking effect occurs at V B ≈ 1.5 |t| as can be observed in the BMsme and meCPT data. The previous version of stsCPT based on the pure zero temperature ground state |Ψ 0 S misses this region of NDC. In contrast to the simpler model presented in the previous section, a self-consistent mean-field solution does not capture the blocking effects correctly in this more elaborate system. The BMsme solution shows many more steps in the current than the stsCPT one, which is due to transitions in the reference stateρ S of the central region. The meCPT results in general follow these finer steps, correcting their width to incorporate also lead induced broadening effects in addition to the pure temperature broadening. As can be seen in panel (B1), meCPT predicts a large charge increase at the site connected to the high bias lead. Note that the charge density at site 2, which is connected to the right lead is simply: n 2 (V B ) = n 1 (−V B ). The charge density at site 3 is symmetric with respect to the bias voltage origin.
Next we study the impact of a gate voltage on the blocking. Results obtained by meCPT are depicted as stability diagram in Fig. 4. Upon increasing |V G |, the onset of the blocking shifts linearly to higher V B (Y). We find a Coulomb diamond for 2V G V B − |t| (D). Upon increasing the bias voltage out of the Coulomb diamond, see e.g. line (X), a current sets in but is promptly hindered by the blocking so that the current diminishes after a hump of width ∝ max(T, Γ). Interestingly this device could be operated as a transistor in two fundamentally different modes. In mode (T1), at a source-drain voltage of ≈ |t| the current is on for a gate voltage of V G = 0 and off for V G ≈ 0.5 |t| due to the Coulomb blockade. In mode (T2), at a source-drain voltage of ≈ 1.5 |t| the current is off for a gate voltage of V G = 0 due to quantum interference mediated blocking and on for V G = 0.25 |t|.
Next we discuss the current characteristics in the vicinity of the blocking in more detail, as well as the impact of the interaction strength U . The first row of Fig. 5 shows the total current through the device for different values of U . The blocking region shifts to lower bias voltages with increasing U . As discussed earlier, structures in the BMsme results are only broadened by temperature effects in the steady-state density (compare e.g. the width of the structures in the local density in the second row of Fig. 5), while meCPT additionally takes into account the finite life time of the quasi particles due to the coupling to the leads, given by 1/Γ. This can be seen by solving Eq. (4) for the local Green's function at device sites. Especially for higher lead broadening Γ this gives rise to significant differences in the meCPT results compared to the BMsme data. From the bottom row of Fig. 5 we see that, before the blocking regime is entered, the steady-state changes from a pure N = 2 state to a mixed N = 2 / N = 3 state at the hump in the current. Obviously, blocking arises because the system reaches a pure N = 3 state for U = 2 |t| and U = 3 |t| at V B ≈ 1.4 |t|. For U = |t| the current is only partially blocked, because the contribution of the N = 2 state is not fully suppressed. For all U -values, however, we find NDC. As far as the meCPT current is concerned, the complete blocking at higher interaction strengths, pre- dicted by BMsme, is reduced to a partial blocking due to the lead induced broadening effects in meCPT. Although ρ S ab changes significantly twice in the blocking region (for U = 2 and U = 3), the charge density n i just increases once from n 1 ≈ 0.75 to n 1 ≈ 1.
Details of the steady-state dynamics are provided in Fig. 6. Before the blocking region is entered (V B = 0.4 |t|) the system is in a pure state with N = 2, which corresponds to the zero temperature ground state |Ψ 0 S in the N = 2 sector. Here the transmission function T (ω), Eq. (D3), of meCPT agrees with the one of stsCPT. A small current is obtained due to the N = 2 → 3 excitation at ω ≈ 0.55 |t|. Increasing the bias voltage has no influence on the reference state in stsCPT, which therefore remains in the N = 2 particle sector. Consequently, the transmission function in stsCPT does not change. Only the transport window increases linearly with increasing V B . For V B = 1.4 |t| it includes the peak at ≈ 0.7|t| and results in a significant increase in the current obtained in stsCPT (see stsCPT result in Fig. 3). This is in stark contrast to the BMsme current, depicted in Fig. 5, which exhibits perfect blocking for V B = 1.4 |t|. The reason for the current-blocking is that only two states, both in the N = 3 sector and doubly degenerate, have significant weight in ρ S ab . The meCPT solution is based on the modified density matrix and therefore the current is diminished, since the next possible excitation is at ω ≈ 0.9 |t| (N = 2 → 3), which is outside the transport window W (ω) ≈ (−0.7|t|, 0.7|t|), Eq. (D2). Due to the lead induced broadening of T (ω) and the temperature induced broadening of the transport window, the current is however only partially blocked. For V B = 2.4 |t| this excitation falls into the transport window and the current is no longer blocked. In this case, the state ρ S ab is a mixture of N = 2, 3, 4. The dominant excitation responsible for this current is again the ground state excitation at ω ≈ 0.55 |t| from N = 2 → 3. This is why in this regime the stsCPT current, based on the pure two particle state is again similar to the meCPT current.

Quasi-degenerate states
Next we study the reliability of the secular approximation in the case of quasi degeneracy of the isolated energies of the system and benchmark its applicability to create a reference state for meCPT. To this end we apply a second gate voltage that couples only to the third orbital, see Fig. 3 (left), and leads to an additional term V G,3n f 3 in the system Hamiltonian. This lifts the degeneracy of states present at V G,3 = 0 and therefore requires a treatment within the BMme, see Ref. 40.
In the following we discuss the same parameter regime as above. In Fig. 7 we present results obtained using meCPT (solid lines) and Qme results (dashed lines) for the BMsme (A) and for the BMme (B). The meCPT results of each panel are obtained using the respective Qme. In the BMsme data a very small |V G,3 | has a drastic effect on the current-voltage characteristics. The blocking present at V G,3 = 0 is immediately lifted by very small |V G,3 | and the current jumps to a plateau. For larger |V G,3 | the current stays on this plateau until further transport channels open up. This "jump" at small |V G,3 | arises due to the improper treatment of quasi degeneracies in BMsme. MeCPT results based on BMsme show a smooth change of the current-voltage characteristics. BMme on the other hand correctly accounts for the coupling of the quasi-degenerate states and also exhibits a smooth dependence on V G,3 . For meCPT based on BMme we find qualitative similar results to meCPT based on BMsme, which emphasizes the robustness of the meCPT results in general. From this it is apparent that meCPT is capable of repairing the decoupling of quasi-degenerate states in the BMsme to some degree. However, to study blocking effects at quasi degenerate points it is of advantage to make use of the BMme in meCPT.
As discussed below in Sec. IV C, the BMme is not of Lindblad form and does not necessarily result in a positive definite reduced many-body density matrix ρ S ab in general. Using a not proper density matrix in Eq. (8) may result in non-causal Green's functions when the steady-state ρ S ab is obtained from the BMme. This can be avoided by using a modified reference state ρ S ab → ρ S ab Θ(∆ − |ω a − ω b |), with Θ(x) the Heaviside step function and ∆ a small quantity, being e.g. ≈ 10 −6 , in Eq. (8), which renders the Green's functions causal. This is somewhat an ad-hoc approximation and should be seen simply as a way to explore the effects of continuously breaking degeneracy in the problem.

C. Current conservation
Finally we comment on conservation laws in meCPT. Within BMsme and BMme the current conservation (continuity equation) is always maximally violated in a sense that the current within the system is zero. This is due to the zeroth orderρ S as discussed in App. D 2. In BMsme the inflow from the left lead into the system however always equals the outflow from the system to the right lead. Without the secular approximation the quantum master equation (BMme) is not of Lindblad form and the final many-body density matrix is not guaranteed to be positive definite. 121,122 This in turn can lead to slightly negative currents in regions where they are required to be positive by the direction of the bias voltage 78 . Furthermore, the inflow can be slightly different from the outflow. In the noninteracting case, meCPT fully repairs the violation of the continuity equation present in the reference state. For increasing interaction strength, the violation of the continuity equation typically grows also in meCPT. In particular, the overall symmetry of the current stays intact (in our case, inflow equals outflow), while the current on bonds between interacting sites does not exactly match the current between noninteracting sites. This typically small violation of the continuity equation can be attributed to the violation of Ward identities 135,136 in the non-conserving approximation scheme of CPT. 137,138

V. SUMMARY AND CONCLUSIONS
We improved steady-state cluster perturbation theory with an appropriate, consistent reference state. This reference state is obtained by the reduced many-body density matrix in the steady-state obtained from a quantum master equation. The resulting hybrid method inherits beneficial aspects of steady-state cluster perturbation theory as well as from the quantum master equation.
We benchmarked the new method on two experimentally realizable systems: a quantum diode and a triple quantum dot ring, which both feature negative differential conductance and interaction induced current blocking effects. meCPT is able to improve the bare quantum master equation results by a correct inclusion of lead induced level-broadening effects, and the correct noninteracting limit. In contrast to previous realizations of the steady-state cluster perturbation theory, meCPT is able to correctly predict interaction induced current blocking effects. It is well known that the secular approximation (BMsme) is not applicable to quasi degenerate problems, which is corroborated by our results for the steady-state current. However, meCPT based on the BMsme density, is able to repair most of the shortcomings of BMsme. The results are very close to those obtained by meCPT based on the density of BMme, where the quasi-degenerate states are treated consistently.
The computational effort of meCPT beyond that of the bare quantum master equation scales with the number of significant entries in the reference state density matrix but is typically small. In the presented formulation the new method is flexible and fast and therefore well suited to study nano structures, molecular junctions or heterostructures also starting from an ab-inito calculation. 139 Here we provide the detailed expressions for the coefficients in the BMme and BMsme of Eq. (6) and discuss the equations governing the time evolution into the steady-state.
The Lamb-shift Hamiltonian is defined asĤ LS = ab Λ ab |a b|, with Λ ab = 1 2i αβ c λ αβ (ω bc , ω ac ) c|Ŝ β |b c|Ŝ α |a * . The environment functions ξ αβ and λ αβ in Eq. (A1) and Eq. (7) are determined by the time dependent environment correlation functions where the Heisenberg time evolution in the environment operators isÊ α (τ ) = e +iĤ E τÊ α e −iĤ E τ . For the BMme, ξ αβ and λ αβ are given by a sum of complex Laplace transforms whereas for the BMsme (ω 1 = ω 2 ) the expressions simplify to the full even and odd Fourier transforms 78 The coupled equations for the real time evolution of the components of the reduced system many-body density matrix ρ S ab = a|ρ S |b according to the BMsme reaḋ The equations simplify further for system Hamiltoni-ansĤ S with non-degenerate eigenenergies ω a . Then the diagonal components φ a = ρ S aa decouple from the offdiagonals and one recovers the Pauli master equation for classical probabilitieṡ with simplified coefficients Ξ ab := Ξ ab,ab = αβ ξ αβ (ω b − ω a ) a|Ŝ β |b a|Ŝ α |b * .
In this case the dynamics of the decoupled off-diagonal terms (a = b) is given bẏ where the simplified Lamb shift terms are Λ a := Λ aa = 1 2i αβ c λ αβ (ω a − ω c ) c|Ŝ β |a c|Ŝ α |a * .
where ξ i and η j denote local spin-1 2 degrees of freedom at the system and environment sites respectively and the overall ordering of operators is important. L S /L E denote the size of the system / environment. Reordering Eq. (1c) where the minus sign arises due to the fermionic anti-commutator. Plugging in the Jordan-Wigner transformed operators leads tô where in the last line we have defined new operators counts the particles between system site i and environment site j for spin σ depending on the ordering of the environments λ. It is straight forward to show that the bar operators fulfil fermionic anti-commutation rules. Furthermore [f iσ ,c λiσ ] − = 0, which allows us to write the coupling Hamiltonian in a tensor product form. Note that in general [f iσ ,c λ ′ jσ ] − = 0 for i = j which is however not relevant for the coupling Hamiltonian where only the same i couple.
The new operators in hermitian form are given in Eq. (5) by replacing c →c and f →f . Next we show, by examining the BMsme, that in most cases the additional phase operator inc drops out of the calculations and we are even allowed to use the original f and c operators instead of the barred ones. The operatorsc only enter the equations in the environment correlation functions C αβ (τ ) as defined in Eq. (A2). Plugging in the barred operators we obtain for normal systems which preserve particle number where we required that [Ĥ E ,P ij ] − = 0. The dropping out of the phase operators implies that for normal systems where the disconnected environments conserve particle number we can omit the Jordan-Wigner transformation and do all calculations as is with the original environment creation/annihilation operators in hermitian form.

Appendix C: Bath correlation functions
In the wide band limit, analytical expressions for the bath correlation functions are available in Ref. 39. For arbitrary environment DOS, explicit evaluation of the environment correlation functions becomes convenient for hermitian couplings, Eq. (5) as outlined in App. B. 78 Essentially the environment functions can all be obtained via integrals of the environment DOS ρ(ω). Care has to be taken when going to very low temperatures and solving the integrals with finite precision arithmetic to avoid underflow errors.

Born-Markov master equation
Within the Qme, basic single-particle observables are available in terms of the reduced system many-body density matrixρ S . The single-particle density matrix κ reads where a and b denote eigenstates of the system Hamiltonian. Note that within the BMme/BMsme κ ijσ is purely real and therefore does predict zero current. However, an expression for the current to reservoir λ can be found by making use of the operator of total system chargeQ and total system particle numberN , where q denotes the charge of one carrier λ j λ (τ ) = d dτ Q (τ ) = q tr Nρ S (τ ) .
Takingρ S (τ ) from the Qme we obtain j λ = q abc n c − 1 2 n b − 1 2 n a Ξ λ ca,cb ρ S ab , and for non-degenerate systems, in the Pauli limit we find from the BMsme j λ non-deg = q ab (n a − n b ) Ξ λ ab φ b .