Open quantum systems beyond Fermi’s golden rule: Diagrammatic expansion of the steady-state time-convolutionless master equation

Steady-state observables, such as occupation numbers and currents, are crucial experimental signatures in open quantum systems. The time-convolutionless (TCL) master equation, which is both exact and time-local, is an ideal candidate for the perturbative computation of such observables. We develop a diagrammatic approach to evaluate the steady-state TCL generator based on operators rather than superoperators. We obtain the steady-state occupation numbers, extend our formulation to the calculation of currents, and provide a simple physical interpretation of the diagrams. We further benchmark our method on a single non-interacting level coupled to Fermi reservoirs, where we recover the exact expansion to next-to-leading order. The low number of diagrams appearing in our formulation makes the extension to higher orders accessible. Combined, these properties make the steady-state time-convolutionless master equation an effective tool for the calculation of steady-state properties in open quantum systems.


I. INTRODUCTION
Open quantum systems constitute a wide research area that permeates both fundamental and applied physics [1,2]. Specific examples include transport phenomena in semiconductor devices [3,4], quantum simulations based on cold atom experiments [5][6][7], as well as quantum information processing with trapped ions [8,9] or with superconducting circuits [10,11]. Given the recent developments in quantum technologies, such systems promise great advances in computing [12,13], simulation [14], and sensing [15][16][17]. Regardless of the specific realisation, open systems can all be broadly described as containing a small system of central interest that is coupled to a large environment. The presence of the environment can fundamentally change the dynamics of the system [18,19], while leaving it sufficiently coherent for quantum effects to be crucial in explaining its behaviour. In electronic mesoscopic transport several phenomena, such as cotunneling [20][21][22] or the Kondo effect [20][21][22][23], fall into this category.
The standard way to cope with open systems is to construct the effective dynamics of the system by integrating out the environmental degrees of freedom. A common phenomenological approach to do so is based on the Lindblad master equation [1], where the most general and physically-allowed evolution of the system is parametrised and then constrained by physical assumptions and experimental data. Conversely, bottomup methods start from a microscopic description of the entire setup including the system, the environment, and their coupling [20,[24][25][26][27]. This latter approach offers greater predictive power by reducing the number of (or even eliminating the need for) fitting parameters [20]. Furthermore, a microscopic description can be readily extended to include higher-order effects.
In setting up such a microscopic formalism, several assumptions have to be made about the environment's state and its coupling to the system. Commonly, we assume an environment that is equilibrated and decoupled from the system in the far distant past [20]. The coupling between the system and the environment then involves a slow switch-on. Technically, this is done with the introduction of a switch-on rate η/ that defines the time scale over which the system and environment are coupled. Such a slow switch-on appears, explicitly or implicitly, in a range of methods, from the functional renormalisation group [28] to simpler perturbative master equations [1,20,29]. The latter can be broadly categorised into three families: (i) Formally exact time-non-local methods, such as the equivalent real-time-diagrammatic (RT) method [30][31][32][33][34][35][36][37], Nakajima-Zwanzig (NZ) master equation [25,26], and Bloch-Redfield (BR) master equation [38][39][40]. (ii) Formally exact time-local master approaches such as the time-convolutionless master equation (TCL) [41][42][43]. (iii) Approximate methods, that include approximations on top of (i) and (ii), in particular Fermi's golden rule [20] and the T-matrix master equation [44][45][46].
The T-matrix approach is often used to generalise Fermi's golden rule [20], however, it predicts unphysical divergences in the switch-on rate /η due to the time non-local nature of the rates [29]. Deep in the perturbative regime, it has been shown that a physically motivated regularisation scheme for the T-matrix [44][45][46] becomes an acceptable approximation when computing currents, but not occupation probabilities [47,48]. On the other hand, the RT or NZ master equation is a time non-local method, that naturally avoids divergences in η [30,31]. The TCL master equation provides a further formally exact description, naturally free of divergences, and produces a conceptually simpler time-local master equation [41][42][43]49]. Furthermore, the TCL has recently been combined with the slow switch-on approximation [49][50][51] such that it can directly be used to compute a perturbative expansion of the steady-state; a development which we call the steady-state time convolutionless STCL master equation.
In this work, we focus on the STCL master equation approach and demonstrate that it serves as a practical arXiv:2010.09838v1 [quant-ph] 19 Oct 2020 tool to compute the steady-state of open quantum systems. We provide a brief overview of current approaches to open systems dynamics and highlight the merits of using the STCL. We then develop a diagrammatic approach to compute the STCL generator and perform the expansion explicitly up to fourth-order for quadratic environments. For practical applications, we extend the STCL to the calculation of currents and again perform the expansion explicitly to fourth-order for quadratic environments. We demonstrate the implementation of our formalism on a non-interacting setup that serves as a test bed.
The paper is structured as follows: in Section II, we briefly highlight and discuss the main results of this work. In Section III, we review the state of the art in the field. We introduce the T-matrix approach in both the usual operator formalism and in terms of superoperators. The real-time-diagrammatic and steady-state time-convolutionless master equations then are directly formulated in the superoperator language. En route, we show that the STCL, which relies on a switch-on process in the distant past, is suitable to compute the steadystate, order by order, whereas it cannot directly be used to compute dynamics without further assumptions. In Section IV, we develop a diagrammatic formulation of the STCL generator S. We use both the operator and superoperator formalisms to minimise the complexity of the diagrams. In Section V, we apply the STCL master equation to setups with quadratic environments and take advantage of Wick's theorem [52]. We then show how the STCL recovers exact results for the occupation numbers in a non-interacting setup. In Section VI, we extend the STCL to compute currents flowing through the system in steady-state, and again show that we recover exact results for the currents in a non-interacting setup. Finally, in Section VII, we summarise the results of our work and give an outlook on future applications for the STCL master equation.

II. MAIN RESULTS
We start with a short tour through our main results and discuss their implications. The goal of the present work is to develop a practical, though exact at every order, method to calculate steady-state probability distributions and transport currents in driven open quantum systems, see e.g., Fig. 1, where we illustrate the electronic mesoscopic setup serving as our physical motivation. We achieve this goal with the help of the steady-state timeconvolutionless master equation [41,42,50,51] ∂ t ρ p (t) = − i L 0 ρ p (t) + S(t, η)ρ p (t), for the (projected) density matrix ρ p (t), with L 0 the Liouvillian of the uncoupled system-environment setup, S the STCL generator and η → 0 a slow switch on rate, see ...
µ n , T n FIG. 1. Open quantum system, device schematic and model. (a) Sketch of a gate-defined quantum dot device. The semiconductor heterostructure defines a two-dimensional electron gas (2DEG) at the boundary (green) between two semiconductor layers (transparent) [21]. The gold top-gates (centre, biased with a set of gate voltages VG) deplete the 2DEG and form a quantum dot (purple), which is coupled to metallic leads on each side. These are, in turn, contacted with wires (gold, sides) that impose a bias voltage (source VS and drain VD) across the device. The latter drives the measurable current I through the dot. A magnetic field B can be applied to the entire setup making the setup spinless for large fields. (b) Generic setup for an open quantum system out of equilibrium, see Eqs. (2) and (3). A system (purple disc) is coupled (via Vr) to the n different reservoirs (grey ovals) constituting the environment. Each individual reservoir r is at equilibrium with a temperature Tr and a chemical potential µr. The environment is usually assumed to be out-of-equilibrium, e.g., µr = µ r and/or Tr = T r . the discussion around Eq. (41). Specifically, we perform three steps: (i) We construct a set of diagrammatic rules to compute the STCL generator S for arbitrary systemenvironment setups, order by order in the systemenvironment coupling V . The results are found in Eqs. (63)- (65) describing the recursive expansion of the STCL generator S in terms of the propagator G presented in Eq. (69). The latter involves the expansion of the evolution operator in Eq. (14). The corresponding diagrammatic representation is depicted in Figs. 4, 5, and 7. (ii) We apply the diagrammatic expansion to setups with quadratic environments, where Wick's theorem greatly simplifies the calculations, and use the corresponding diagrammatic formulation from Figs. 10 and 11 to obtain explicit fourth-order rates for the STCL. These rates are reported in Eqs. (101), (105)-(107), and (B13)-(B15). In Fig. 18, we highlight the result of computing these rates for a non-interacting setup and then calculating the steady state occupation P 1 of a single-level in a non-interacting quantum-dot at equilibrium. We compare this analysis (up to next-to-leading order) to the expansion of the exact result and find perfect agreement.
(iii) We extend the STCL to compute steady-state currents for setups with quadratic environments. We modify the diagrammatic formulation to include so-called current rates. This provides us with the results in Eq. (127) for the currents as expressed through the current generator (141), and the charge transfer propagator (132), see also Fig. 21. In Fig. 24, we again show up to fourthorder, that these results recover the exact results for a non-interacting single-level setup.
The diagrammatic expansion presented here can be systematically expanded beyond fourth-order, e.g., we predict that only three times more diagrams appear at sixth order as compared with what we have computed at fourth-order. There are several effects, such as measurement backaction [53] and the Kondo effect [20], that manifest at sixth order, motivating further development of our diagrammatic scheme. Lastly, we highlight that the STCL lends itself to resummation schemes, similar in spirit to those that have already been developed for the RT method [30][31][32]36], see Ref. [48] for a brief discussion.

III. BACKGROUND
We consider a setup composed of a system and an environment as shown in Fig. 1, and compute its steady-state properties. To this end, we focus on the time evolution of the system after a slow switch-on using a master equation. Expanding the rates that govern this evolution order-by-order in the system-environment coupling V , we obtain a power series for the steady-state properties of interest. In this section, we provide an introduction to the T-matrix [20] and steady-state time-convolutionless [41][42][43]51] master equations. Concurrently, we comment on the relationship between these methods and the realtime diagrammatic approach [1,25,26,[30][31][32][33][34].

A. General properties
The system is assumed to be small, i.e., it is described by a finite Hamiltonian H sys which we can solve exactly, e.g., via numerical techniques, providing eigenstates |i sys and eigenenergies χ i . The environment may be large/infinite, is described by the Hamiltonian H env , and is commonly assumed to be non-interacting such that it is exactly solvable. It has eigenstates |u env with eigenenergies ξ u . The combined unperturbed Hamiltonian describing the decoupled system and environment reads H 0 = 1 env ⊗ H sys + H env ⊗ 1 sys , (2) with eigenstates |u, i = |u env ⊗ |i sys and eigenenergies ξ u + χ i . Throughout this work, we omit system and environment subscripts at times when it improves readability (we use indices i, j, f, g, n, m for system states whereas u, v are reserved for environment states). The environment is composed of multiple reservoirs r that are individually at equilibrium (with corresponding temperatures T r and chemical potentials µ r ) but mutually out-of-equilibrium. Each reservoir may be composed of fermions or bosons (or particles with exotic statistics). They are coupled to the system by individual perturbations V r , the sum of which makes up the total perturbing Hamiltonian The total Hamiltonian H = H 0 + V which governs the physics of the full setup is the sum of the unperturbed part and the system-environment coupling.

Unitary time evolution
The state of the full setup is described by the density matrix ρ, which evolves according to the von Neumann equation [15] Here, we allow the Hamiltonian H to be time dependent to account for the switch-on of the perturbation V . Making use of the unitary time-evolution operator with the time-ordering operator T , the differential equation (4) can be formally integrated, With the density matrix composed of elements ∝ |u, i v, j|, the forward-evolution U acts on the density matrix ρ(t 0 ) from the left and propagates a state |u, i from time t 0 to t; it is commonly referred to as the forward Keldysh branch of the time evolution. Concurrently, the backward evolution U † acts on ρ(t 0 ) from the right and propagates a conjugate state v, j| forward in time, commonly referred to as the backward Keldysh branch. These two branches, shown pictorially in Fig. 2, form the basis for our diagrammatic representation in Section IV. Combining the von Neumann equation (4) with an initial condition ρ(t 0 ) at a time t 0 fully specifies the density matrix of the system-environment setup at a later time t.
The Keldysh contour; graphical representation of Eq. (6). The density matrix ρ(t0) of the entire setup at time t0 is evolved to a later (leftward) time t by unitary time-evolution operators (triple line) for the coupled systemenvironment setup on each of the Keldysh branches. Note that, on the lower branch, the evolution is given by U (t0, t) = U † (t, t0), corresponding to the different directions of propagation on the two branches.

Distant past
We assume that in the distant past t 0 → −∞, the system and environment were decoupled. Such an initial condition implies that the density matrix at time t 0 can be decomposed into the product of the system's density matrix ρ sys (t 0 ) and a locally-equilibrated distributionof the environment Here, the probability distribution P env u describes the thermal-equilibrium configuration of each reservoir in the environment.
The large size of the environment makes the direct calculation of the density matrix ρ(t) from the initial condition (7) and the von Neumann equation (4) intractable. To tackle this problem, two steps are commonly applied: (i) expand the time evolution operator (5) as a perturbation series in V , and (ii) integrate out the environment dynamics to be left with only the system behaviour. As we will see in Eqs. (14) and (15), the expansion is incompatible with the t 0 → −∞ limit for a time-independent perturbation V . To circumvent this problem, we introduce the slow switch-on assumption where η/ is an infinitesimal switch-on rate. The system and environment are in contact for an effective timescale /η → ∞ before the present time t = 0, and thus if a steady-state exists it will have been reached.

Expanding U
The expansion of the time evolution operator (5) can be written in the form where the terms U ν of order V ν can be found using the recurrence relation with the first term U 0 specified by the free evolution of the unperturbed Hamiltonian H 0 in Eq. (2), We perform the integral in the recurrence relation (11), keeping a finite switch-on rate /η and assuming t 0 → −∞ such that This relation is assumed to remain true when taking the second limit η → 0 later on. We work in the eigenbasis |u, i of the unperturbed Hamiltonian H 0 and repeatedly insert identity operators to obtain where U ν now depends on t and η (but not on t 0 → −∞) due to the slow switch-on. There is no clear notion of intermediate times in Eq. (14), however a specific operator ordering for the coupling events V is inherited from the time-ordering in Eq. (11). Note that each time the recurrence (14) is applied, an additional factor of exp(ηt/ ) appears, which is the reason for the ν-fold enhancement of η in the denominator. Furthermore, the rightmost operator in the recurrence is always an unperturbed timeevolution U 0 . The recurrence relation (14) is undefined at η = 0 as the real part of the denominator of the free propagator vanishes for states |u, i that fulfil χ i + ξ u = ω. The resulting divergences motivate the introduction of the slow switch-on (9) with a finite η > 0. Making use of the perturbative expansion (10), we can now integrate out the environment. The T-matrix, usually familiar from scattering theory [20], provides a standard approach to achieve this goal: It is commonly used to construct a rate equation for the probabilities P i of finding the system in a state |i , i.e., the diagonal elements of the density matrix. Such a rate equation is also known as a Pauli master equation, and can be used to study relaxation, i.e., the time scales on which the system decoheres due to transitions between different states [1,20]. In contrast, a full master equation can be used further to study dephasing, where the system state is preserved but the off-diagonal elements of the density matrix decay [1]. Next, we sketch the derivation of the Pauli T-matrix rate equation as it will allow us to effectively highlight its pitfalls.

B. Fermi's golden rule and the T-matrix
Let us consider the probability P sys f |i of finding the system in a state |f at time t, given that in a far distant past the system was in an initial state |i . The environment is in equilibrium at t 0 and its state at time t is irrelevant. To compute the probability P sys f |i , we evolve the state |u, i from t 0 to t, square its overlap with a final state |v, f , and sum over the environment states u, v, weighted by the distribution P env u of initial environment states. This procedure propagates the system probabilities according to and only depends on the original environment distribution, implying that we have integrated out the effect of the environment at all later times. In order to obtain the T-matrix rate equation, we insert the expansion (14) into the expression (16) for the propagated system probability. We differentiate the result with respect to time t which we set to t = 0 without loss of generality. Taking the limit η → 0, we obtain the transition rate where we have introduced δχ if ≡ χ i − χ f and δξ uv = ξ u − ξ v for compactness. The delta function in Eq. (17) ensures energy conservation and arises from two occurrences of the leftmost denominator in Eq. (14) due to the two occurrences of U in Eq. (16). Furthermore, we have introduced the T-matrix T M that is defined by the self-consistency relation with 0 + a positive infinitesimal. Note that in the Tmatrix rate equation the infinitesimals 0 + in each denominator are usually all taken equal, which is justified if Γ T is convergent in the η → 0 limit. However, this convergence is not guaranteed and thus the prefactor ν of η in Eq. (14) is a crucial ingredient that must be tracked in an exact method. At lowest order T M ≈ V and the rates (17) are identical to Fermi's golden gule.
At higher orders in V , the rates (17) are divergent in the limit η → 0.
The rate matrix Γ T and its elements Γ if T relate the change in the probability P sys f (0) (to find the system in state |f at time t = 0) to the corresponding probabilities P sys where Γ ii such that the total probability is conserved. From here onwards we shall drop the "sys" label on the system probabilities for brevity, while the "env" label will remain to avoid ambiguity. Equation (19) with its (time non-local) rates Γ if T does not generate a time-local rate equation, as pointed out in Ref. [29]. The equivalence at lowest order between the T-matrix rates and Fermi's golden rule means that it is tempting to use the T-matrix rate equation (19) as a higher-order generalisation of Fermi's golden rule [20]. To do so, Eq. (19) is turned into an approximate time-local master equation which at lowest non-vanishing order in V is identical to Fermi's golden rule. This path is fraught with difficulties though, even for the seemingly simple task of calculating the steady-state distribution of the system. The standard procedure [20,29,[44][45][46] to do this is to use (20) and assume that the system has reached the steady state at t = 0 such that n Γ nm T P n (0) = 0.
It is then claimed that solving for P n (0) in Eq. (21) gives us the steady-state system-distribution. However, we have in fact solved for the distribution of system probabilities at t 0 = −∞ that leads to the steady-state at time t = 0, cf. Eqs. (19) and (21) or see Ref. [29]. Furthermore, in the limit η → 0, every distribution at t 0 = −∞ leads to the steady-state at t = 0, as an infinite amount of time has elapsed, and the constraint (21) becomes ill-defined. This manifests as divergences in the T-matrix rates Γ T , which signal the breakdown in the approximation (20). Consequently, physically motivated regularisation schemes for Γ T have been developed [44][45][46] that remove the divergences, but produce results which differ from the exact expansion [54]. The STCL, on the other hand, takes the step from Eq. (19) to (20) in a rigorous manner, which is why we choose to use this approach in our present work.

C. Superoperator formulation
Next, we introduce the mathematical tools, specifically the superoperators, which facilitate the calculations of the various master equations [1]. We then briefly rederive the T-matrix in this more formal language and then provide a simplified derivation of the STCL inspired by Refs. [49][50][51]. Additionally, we discuss the relationships between the T-matrix, the STCL, and the RT approaches.

Setup evolution
We start with the von Neumann equation (4) in superoperator notation with the Liouvillian superoperator Lρ ≡ [H, ρ]. Every Hamiltonian (H 0 , V ) is associated with a Liouvillian (L 0 , L V ) (we denote superoperators by caligraphic upper case letters [55]). In analogy to (6), we formally integrate (22) to obtain the time evolution with the time evolution superoperator which can also be written in terms of the unitary evolution Uρ = U ρU † . We expand the superoperator U = ∞ α=0 U α in the Liouvillian L V (equivalent to an expansion in V ) and perform the integrals, in the limit t 0 → −∞, as for the unitary evolution, see Eqs. (5) and (11)- (14). The resulting recurrence relation takes the form which is similar to the one for the unitary time evolution (14). Each time the recurrence is applied an additional factor of exp(ηt/ ) appears, the rightmost superoperator is an unperturbed evolution U 0 , and the superoperator ordering is inherited from the time-ordering in Eq. (24). By expanding the time evolution of the density matrix in the operator-(6) and superoperator- (23) representations and comparing order by order, we obtain the identity Henceforth, a sum over µ+ν = α implies that each of the indices ν, µ runs from 0 to α and in opposite direction for the partner. Note that, the single term U α in Eq. (25) contains α Liouvillians L V , and hence α commutators with V . When written explicitly, these commutators lead to 2 α terms, significantly more than the α + 1 terms on the right hand side of Eq. (26). We thus conclude that Eq. (26) significantly reduces the complexity of computing series expansions of the evolution, and will use this feature in Sec. IV.

Projected space
Having specified the evolution in the superoperator language, we proceed by integrating out the environment part of the evolution. To this end, we define the projector P through its action on a density matrix where Tr env is the partial trace over the environment states |u env . In effect, Eq. (27) projects any density matrix to the space of valid initial conditions, such that The projector P obeys the usual condition P 2 = P, and it is straightforward to verify that [L 0 , P] = 0 as well as where Tr is the trace over all setup (system and environment) states |u, i and O is any setup operator. We define the projected density matrix which carries only the system degrees-of-freedom but resides in the Hilbert space of the full setup, i.e., we can write At this point, we can write down the projected timeevolution superoperator U P = PU P, which directly propagates the projected density matrix from a decoupled initial time t 0 up to time t. This projected evolution will be central to the derivation of master equations in the rest of this work.
To obtain the T-matrix master equation in the superoperator formalism, we differentiate (32) with respect to time t using the identity which follows from the definition of U in Eq. (24). We substitute the result into the projected evolution (32), make use of the commutator [L 0 , P] = 0 and obtain the T-matrix master equation (describing both relaxation and dephasing) in the form see Fig. 3. Here, we have taken the limit t 0 → −∞ and introduced the T-matrix generator in its superoperator form with R = α R α . Its explicit perturbation expansion can be obtained by inserting the expansion for the superoperator U, either from Eq. (25), or the one from Eq. (26) if the operator representation is preferable. At higher orders this leads to divergences in η as for the usual T-matrix formulation, see Section III A.

Pauli projected space
The master equation (34) tracks both relaxation and dephasing, i.e., on-and off-diagonal elements of the system density matrix. On the other hand, Pauli master equations describe only relaxation, and are thus sufficient to calculate occupation probabilities. They require the introduction of the Pauli projectorP, and Pauli projected space ρp, through where Diag sets all non-diagonal elements of a matrix to zero. The projectorP thus limits us to track only the system probabilities which can in turn be used to reconstruct the Pauli projected density matrix Assuming that the initial system condition is completely dephased we can make use of the Pauli projected evolution UP = PUP to construct Pauli master equations.
As an example, the Pauli T-matrix master equation takes on the form which is obtained by replacing each occurrence of the projector P in Eq. (34) by a Pauli projectorP, see Eq. (36). Furthermore, we have used the property L 0P = 0, which can easily be verified, and introduced the Pauli T-matrix generatorR, which is obtained by replacing P →P in Eq. (35). This superoperator formulation of the Tmatrix Pauli master equation is identical to the one in Section III A, cf. Eqs. (19) and (40) with t = 0.

D. Steady-state time-convolutionless approach
The steady-state time-convolutionless [1,41,42,[49][50][51]56] master equation is of the form (see below for the derivation) with the superoperator S(t, η) depending on t and η but not on t 0 , that has been sent to the distant past, t 0 → −∞. Its time-local property naturally sidesteps the issues that appear in the T-matrix approach. As shown in Refs. [50,51], each term S α in the STCL series expansion S = ∞ α=1 S α is convergent for the setups we consider here and in the limit η → 0. The latter limit, in combination with t 0 → −∞ implies that the system and environment have been in contact for an infinite amount of time, such that the asymptotic state has been reached. Specifically, the projected density matrix may display an oscillating behaviour which is basis dependent.
The Pauli projected density matrix, however, reaches a basis dependent but constant steady-state [50]. As for the T-matrix scheme, the STCL master equation has a Pauli equivalent which is of central importance to the present work. Simultaneously, the Pauli STCL superoperator is uniquely defined by its matrix elementsS if , which are the timelocal transition rates between system states |i and |f . Written in terms of (diagonal) matrix elements, the Pauli STCL master equation (42) becomes which we will use for practical calculations in sections V and VI. As the Pauli STCL (42) is time local, a steadystate constraint exists, and is obtained by setting the derivative in Eq. (43) to zero nS nm (0, η → 0)P n = 0, where P n is the steady-state probability distribution of the system [57] that also satisfies n P n = 1. To avoid notational clutter, we denote the steady-state system probability by P n , simply dropping the time-dependence. By expanding Eq. (44) and given an explicit representation for the perturbation series ofS, we obtain a formally exact power series for P n . In the next subsection, inspired by Refs. [49][50][51], we present a simplified derivation of the superoperator S.

Derivation of S
In order to find an expression for the STCL superoperator S, we invert the projected time evolution (32) and insert it back into the superoperator formulation of the T-matrix (34) approach; the result then provides the identification [49] Different master equations in a nutshell. (a) Tmatrix master equation functionality. The full density matrix at time t0, where ρp = ρ, is propagated forward and differentiated with respect to time and then brought into the projected space using the time derivative of the projected evolution ∂tUP . (b) Steady-state time-convolutionless master equation functionality. The projected density matrix at time t and the inverse projected evolution U −1 P are used to construct the projected density matrix at time t0, where ρp = ρ. The result is then fed through the same procedure as in the Tmatrix master equation. see Fig. 3. Starting directly from Eq. (46) provides a significant technical advantage in deriving an expression for the generator S as compared to previous derivations [50,51]. As pointed out in Ref. [58], the operator U P is analytic in time and finite dimensional as it evolves only the projected space density matrix. It is thus invertible at all but isolated points in time, see Ref [48]. This mathematical property guarantees the existence of the STCL for any non-zero finite value of η or t 0 . In the limit |t 0 | /η → ∞, the system and environment have been in contact for an infinite time and the steadystate has been reached. The STCL we use in this work is therefore limited to steady-state calculations and cannot describe dynamical properties.
To obtain a series expansion of the inverse U −1 P , it is convenient to introduce the time-local propagator where I is the identity superoperator, and its Pauli equiv-alentG obtained by replacing all occurrences of P byP. The time-local propagator G encodes the non-trivial evolution in the projected space and its series expansion The terms G α are divergent in the limit η → 0 and must therefore be computed as a Laurent series in η, i.e., a power series that includes terms of negative powers of η.
In the projected space, we can recast the projected time evolution (32) in terms of the time-local propagator Note that this representation of U P is only appropriate when acting on projected density matrices as we have omitted the projector P in the unperturbed evolution. It is then straightforward to write down the inverse of U P as a power series We insert Eq. (50) and the T-matrix generator (35) into Eq. (45) and compare the result, order by order, with the STCL master equation (41). As a result, we obtain the recurrence relation for the expansion of the STCL generator. Both the Tmatrix and STCL generators, R and S, conserve probability (but not necessarily positivity) at each order, i.e., which can be seen by first taking the trace over Eqs. (41) or (34) (or their Pauli equivalents) and then performing an expansion on both sides. The fact that probability must be conserved (52) at each order, puts constraints on each term S α , which we will use to reduce the number of diagrams we compute later on in section V. From equations (51) and (35), we immediately infer that at lowest non-vanishing order the Pauli (PU 0 =P) STCL and T-matrix superoperators coincide [50,51]. Hence, at the lowest order the Pauli STCL is identical to Fermi's golden rule. At higher orders, the T-matrix includes divergent terms, which in the STCL are correctly subtracted by the recursion formula that generates a time-local version A loc (t) = α A loc α from the expansion A α of a propagating projectedspace superoperator A(t, t 0 ). The recursion (53) appears in (51) and was derived in Ref. [51]. In the specific case of the T-matrix and STCL master equations, it can be thought of as a regularisation procedure, which subtracts reducible contributions, and guarantees that each term S α in the expansion of the STCL generator is convergent in the limit η → 0 [50,51]. We stress that the regularisation (53) leads to exact results for the power series of steady-state system properties that are different from those obtained by the more common, physically motivated, regularisation scheme [29,[44][45][46][47]50].
They encode the evolution with and are embodied by the kernel K. This approach is conceptually more complex than the T-matrix or TCL master equations as it involves an integral over all previous times; it is often combined with physical assumptions to generate a Lindblad master equation [1] or other more complex often non-Markovian master equations [59]. Furthermore, a set of powerful tools, including resummation schemes [32][33][34][35][36] and renormalisation group methods [60][61][62][63][64][65] have been developed for these master equations. As for the T-matrix and STCL approaches, Eq. (54) has a Pauli equivalent whereK is the Pauli Kernel. The derivation of either the full K or PauliK kernels, can be found in a breadth of pedagogical texts on the topic, see for example Refs. [1,47,48]. The RT method is commonly used to compute the asymptotic-state of the projected-space density matrix. Here, we focus on the Pauli case where the probabilities reach a steady-state and refer the interested reader to Ref. [48] for a treatment of the off-diagonal terms. We combine Eq. (55) with a steady-state assumption and move the limit t 0 → −∞ of the integral to the far distant past. Making use of the steady-state assumption (56) (also known as the Markov approximation), we add the qualifier steady-state to the RT method (SRT), as we already did for the TCL. We substitute Eq. (56) into the (Pauli) real-time master equation (54) and obtain the conditioñ which when solved for, order-by-order in V , leads to the same steady-state as the STCL condition (44). In Eq. (57), we introduce the Pauli SRT generatorZ, which is obtained by performing the integral over time [30]. The series expansionZ = αZ α of the Pauli SRT generator is closely linked to the expansions of the T-matrix and STCL generators, see Refs. [1,30,47,48] for a derivation.
Here,Ũ α is the Pauli version (u = v and i = j) of the expansion of the full evolution minus secular contributions (I → I −P), such that with the initial conditionŨ 0 = U 0 , cf. Eq. (25) with t = 0. The Pauli SRT generatorZ is thus composed of an unperturbed backward propagation U −1 0 and a full forward propagation minus the secular contributions. The same way the regularising recurrence (53) guarantees that S is finite at each order, the subtraction of secular contributions (through I → I −P) guarantees thatZ is finite at each order. The full SRT generator Z is not as easy to obtain as in the STCL case due to the oscillating off-diagonal elements, see for example Ref. [48].
Both the SRT and STCL methods lead to identical power series for steady-state observables, through related but different generators. This difference between the STCL and SRT methods occurs because the STCL generator S carries a remnant of the dynamics within it. The SRT master equation, instead, includes an explicit steady-state (or Markov) assumption (56) and therefore cannot be straightforwardly extended to an exact timelocal master equation. We choose to work with the STCL over the SRT for three main reasons: (i) its conceptually simpler form, (ii) its direct formulation in terms of U P , which we will use to considerably simplify the computation of rates, and (iii) its potential to be extended to the TCL, a time-local master-equation that can describe the dynamics after a quench. Next, we proceed with the derivation of a diagrammatic expansion of the STCL generator S.

IV. DIAGRAMMATIC EXPANSION
Computing the rates in the superoperator S of the STCL master equation (41) requires evaluating the expressions PL V UU −1 0 P and G order by order in V , see Eq. (51). The derivation in III D does not rely on the limits t 0 → −∞ or η → 0, and Eq. (51) is valid for any initial time as long as t 0 is adapted correspondingly in G and U.
We will now show that, in the t 0 → −∞ limit, the combination PL V UU −1 0 P can be recast in terms of G and thus finding the latter is sufficient for computing S. We then construct a novel diagrammatic approach to compute G, inspired by T-matrix calculations [20,[44][45][46]66] and the simplifications in Ref. [47]. Specifically, we will perform the expansion of G using operators, as it leads to a significant reduction in the number of terms compared with the superoperator formulation, see Appendix A. As G contains terms that diverge in 1/η and higher powers thereof, it is necessary to perform the computations for finite η and then express the result as a Laurent series (including negative powers of η). Our diagrammatic representation of G can then be used to construct a diagrammatic expansion for S. The recurrence (53) then guarantees that negative powers of η are subtracted from the expressions for S.
Our first goal is to simplify the expression (51) for S α , and more specifically the term PL V U α−1 U −1 0 P, by recasting Eq. (51) in terms of G α and L 0 only. To do so, we compare two alternative representations of the derivative ∂ t U of the full time evolution. To obtain the first expression, we expand both sides of the identity (33) order by order in V to arrive at the equality The second form, valid only in the t 0 → −∞ limit, is found by differentiating (25) with respect to time and expanding to obtain Combining the two expressions (60) and (61), we find the relation which can easily be verified using (25). Substitution into Eq. (51) provides us with the result that only involves the expansion of G and and the Liouvillian L 0 . The recursion in (63) is of the form (53) that provides us with a regularised expression and guarantees that each order S α is convergent in the η → 0 limit. The first two terms in Eq. (63) are generated by the time derivative; the first one arises from the derivative of the explicit time dependence of the perturbation V ∝ exp(ηt/ ). The second terṁ is similar to the von Neumann equation (4) but formulated in the superoperator space; it can be thought of as the evolution of the propagator G under the influence of the unperturbed evolution L 0 . The above derivation tells us that the expansion (48) of the generator G is the key to computing the STCL generator order by order. We use the Hamiltonian representation (26) of U to reduce the number of terms in our computation by defining where U νµ ρ = U ν ρU † µ is the contribution to the time evolution of ρ with ν (µ) occurrences of the perturbation on the upper (lower) Keldysh branches. We apply this same procedure to the expansion (48) of the generator G, such that where G νµ = PU νµ U −1 0 P and ν + µ ≥ 1. Similarly, the terms A α in the expansion of any superoperator A generated from G, can be decomposed into terms A νµ . This decomposition will reduce the number of terms that we have to compute, from 2 α to (α + 1), in complete analogy to the diagrammatic grouping in Ref. [47], see also Appendix A.
Next, we construct the terms S νµ by generalising the regularising recursion (63) to account for the two Keldysh branches Here, each occurrence of the subscript α in the definitions (64) and (65) of G andĠ is replaced by α → νµ.

B. Diagrammatic rules
We develop our diagrammatic technique for the elements of the generator S. While these are in fact elements of a four dimensional tensor, we will refer to them as matrix elements to keep the nomenclature consistent between the full and Pauli generators. They are indexed by four system indices i, j, f, g, and are the rates that can then be used to compute the steady-state. Furthermore, they can be used to reconstruct the superoperator through the identity The matrix elements for other projected-space superoperators, such as G, R, and L 0 , are defined in the same way (69) as those for S. Simultaneously, these matrix elements can be used to reconstruct their corresponding superoperators with Eq. (70). The generator S is constructed from the propagator G, see Eq. (63), and we have to consider the expansion of the latter first. We use the decomposition (66) of U α in terms of unitary evolutions (14) in the definition (67) of the expansion of G and insert the result into Eq. (69) to obtain  48) and (67). Note that the lowest order terms contain one occurrence of the coupling V (black dot).
Here, we immediately notice the strength of the operator formulation, where G α contains α + 1 terms, each with α occurrences of V , see Eqs. (67) and (71). On the other hand, the superoperator formulation (25) contains a single term with α occurrences of the Liouvillian L V , which leads to 2 α terms once the commutators in the Liouvillians are written explicitly, see also Appendix A. The Common diagrammatic operations. (a) Internal wiggly lines (or cut) appearing in the multiplication of two projected-space operators, see Eq. (72). The cut signals that there are actually two distinct diagrams which must be multiplied. Alternatively, it is possible to formulate a set of rules for them and retain a single diagram, see main text. (b) The commutator of G with L0, see Eq. (74), is represented by the diagram for G with diamonds at each of the corners, marking the four system energies χ. The ± signs mark positive/negative contributions of the energies χ and are dropped henceforth.
hybrid system-environment state |u, i backwards in time with an unperturbed evolution U † 0 , before propagating it forward in time with ν system-environment scattering events. We use the recurrence relation (14) for U ν to define the diagrammatic rules (in energy space) for the matrix elements v, f | G ν |u, i , see Fig. 4(a): 1. Act with a perturbation V , which mixes the system and environment.
3. Repeat the steps 1 and 2 a total of ν times.
The operator G ν (G † µ ) appears on the upper (lower) Keldysh branches, see Eq. (71) and Figs. 2 and 4(b). The initial environment state |u is the same on each of the Keldysh branches and is summed over, weighted by the probability P env u (encoded in ρ 0 env ) of finding the unperturbed environment in the state |u , see Fig. 4

(b) and
Eq. (71). As the environment is traced out in (71), the final environment state |v is the same on the upper and lower Keldysh branches as well, see Fig. 4(b). We define a simpler diagrammatic form for G ijf g νµ in Fig. 4(c), where we merge the system and environment lines for simplicity. Note that the decomposition (26), means that scattering events on each Keldysh branch have an operator order inherited from the temporal order, however, there is no inherited ordering relation between the upper and lower branches.
To construct a diagrammatic expansion of the STCL generator S, see Eq. (68), we require two further elements. The first is the composition AB of two projectedspace superoperators A and B, which in matrix element Diagramatic representation of the STCL generator S, see Eqs. (46), (51) and (63) (68) for Sνµ, where the first two terms refer to G νµ andĠνµ respectively. The remaining lines encode the recursive regularisation (68). The hatched parts again denote lower-order versions of S ν µ , which themselves must be evaluated using the recursion, leading to diagrams with an increasing number of cuts.
form is Diagrammatically, we replace this product by a cut, see Fig. 5(a), which can be understood as reset of the free propagators Π ± 0 from Eq. (15). Namely, the counter ν for the prefactor of η in the propagators Π ± 0 (ω, νη) returns to one, and the energy ω becomes that of the system state at the cut χ n or χ m plus the environment energy ξ u of the unperturbed environment state |u (which are summed over, weighted by the unperturbed distribution P env u ). The second new element involves the commutator of a projected-space operator, e.g., G with the unperturbed Liouvillian L 0 . In the projected space, the matrix elements of the unperturbed Liouvillian are which follows from the definitions of L 0 in (22), of the projector P in (27), and of the matrix elements (69). To obtain the matrix elements of the superoperator commutator [G, L 0 ] from Eq. (65), we combine Eqs. (73) and (72), leading to see Fig. 5(b) for a simple diagrammatic representation. The diagrams for S νµ are composed of three parts, see Eq. (63) and Fig. 6. The first two are G νµ andĠ νµ , which are obtained by multiplying G νµ by (ν + µ)η and −i(δχ ij − δχ f g )/ respectively. The last part is the recursion scheme, which for convenience, we rearrange according to the depth of the recursion, see also Ref. [56]. We thus arrive at rules to generate S νµ directly from G νµ : 1. Create every distinct composite diagram by introducing cuts into the diagram.
2. The prefactor of the subdiagrams is multiplied by 1 (−1) for an even (odd) number of cuts.
3. Transform the leftmost term G in each diagram into the sum of its time differentiated versionsĠ and G .
Combining these three rules and the rules for G in Fig. 4, we generate all terms which contribute to S νµ . We exemplify this scheme with a representation of S 33 in Fig. 7.
Its generalisation provides us with the diagrammatic expansion of the full STCL generator S.

C. Pauli STCL
The Pauli STCL master equation (42) does not track the off-diagonal elements of the system density matrix in ρ p . Its generatorS is obtained by replacing all occurrences of the projector P →P in the STCL generator Eq. (68), leading tõ where we have used thatPL 0 = 0. As for the full STCL generator, we decompose the Pauli generatorS according to the number of perturbations ν, µ on each of the Keldysh branches The Pauli STCL generator (76) should be contrasted with the Pauli T-matrix generator which is not convergent in the η → 0 limit. Note that the matrix elementsG if νµ of the Pauli superoperatorG are indexed by only two system indices i and f , i.e., they are obtained by enforcing Diagrammatic construction of the STCL generator, exemplified on S33. By introducing every possible set of cuts in the diagrams, we implement the recursion (68). At first the number of diagrams in S33 appears large, however we note that all contributions that contain a cut are obtained from products of lower order diagrams, such that a single sixth order object (G33) must be computed. The prefactor is negative (positive) for diagrams with an odd (even) number of iterations of the recursion, i.e., the number of introduced cuts.
see Eqs. (36) and (48), as well as Fig. 8(a). In contrast, the matrix elementsS if of the Pauli STCL generator contain multiple projectors and are obtained by combining the Pauli propagatorsG according to Eq. (76) and cannot be immediately obtained from S ijf g . In the diagrammatic representation, we differentiate between G and its Pauli counterpartG by a single star in the starting environment distribution. On the other hand, S andS additionally differ by a star at every cut, see Fig. 8 The additional constraint (78) immediately implies that theĠ contributions vanish for the Pauli master equation, see Fig. 8(c) and Eqs. (74) and (78).
There are several properties of the Pauli STCL generator that can be used to reduce the number of diagrams that need to be computed. First, note that matrix elements related by an inversion of the number of scatterings on the two Keldysh branches are complex conjugates of each otherS if νµ =S if * µν . Next, we show how the conservation of probability can be used to avoid computing the diagrams with a vanishing index ν = 0 or µ = 0. Just as for the full master equation, the conservation of probability (52) implies that Here,S ii has the usual interpretation as the inverse lifetime of the system state |i . We use Eqs. (71) and (76) to conclude thatS if α0 ∝ δ if (and similarly forS if 0α ), and use this conclusion to split (79), such that Note that here, unlike in (79), the sum over the system states f is arbitrary and does not exclude diagonal elements. Hence, the sum of diagrams with a zero index follows from the diagrams without vanishing indices, and at lowest order we obtain the simplificationS 01 +S 10 = 0. The rates of the formS ii α0 orS ii 0α are purely diagonal and, as we will see in section VI, they are not associated with a physical process, such that we denote them as probability conserving rates.

D. Discussion
We have developed a diagrammatic method for the matrix elements of the STCL generator S and its Pauli coun-terpartS, which extends upon results from Refs. [30-34, 37, 47, 50, 51, 56]. Within our scheme, there is no inherited time-order between the two Keldysh branches, i.e., the vertices on each branch are ordered but are free to move left or right with respect to the other branch. This should be contrasted with the more common approach, where scattering events on different Keldysh branches obey a specific (left-right) ordering with respect to each other [30,31]. This simplification drastically reduces the number of diagrams for G α to be evaluated, from 2 α to α + 1. These results for G α can then be used in combination with products of lower order diagrams to generate the expansion S α of the STCL. We achieved this simplification by using Eq. (26), which recasts the evolution U of the full density matrix in terms of pairs of the unitary evolution U , see Appendix A. We thus develop a generalised STCL equivalent of the diagrammatic simplification that was presented in Ref. [47] at fourth-order for the RT method. To illustrate the method, we draw the first-and second-order Pauli diagrams in Fig. 9. In the next section, we apply this formalism to setups with quadratic environments, where Wick's theorem allows for further simplifications, and compute explicit rates up to fourth-order for the Pauli STCL.

V. QUADRATIC ENVIRONMENTS
In the following, we focus on setups with quadratic environments. We restrict ourselves to the case of systemenvironment couplings that involve a single environment . . . particle, though our diagrams straightforwardly generalise to multi-particle couplings (such as in the Kondo model [20]). We use the Pauli STCL as it is sufficient to calculate the steady-state distribution in-and the current through-the system.
The strength of the (Pauli) STCL and our diagrammatic approach is demonstrated at fourth-order in V , with a natural and correct regularisation and only five diagrams that need to be computed. The fourth-order ratesS if 4 are convergent in the limit η → 0, as opposed to the T-matrix ratesR if 4 which diverge, see, e.g., Eq. (101) below. The ratesS if 2 andS if 4 provide the exact results for the first two terms P (0) n and P (2) n in the expansion of the steady-state system probabilities P n . We validate our analysis through a comparison with the exact solution of the non-interacting resonant level.

A. Setup
Our setup consists of an arbitrary system, a quadratic environment, and a system-environment coupling that involves only a single environment operator at each vertex. This last assumption does not affect the diagrams from Sec. IV. It does, however, influence the specifics of how Wick's theorem is applied, as seen later in this subsection [67]. While we work with a fermionic example, our approach is generic. To keep the notation concise, we group all indices of the environment into the index κ = (λ κ , k κ ). Here, the indices of all discrete degrees of freedom, such as particle/hole, reservoir r and spin σ are encoded in the discrete multi-index λ = ±(r, σ, . . . ), while k indexes momentum. The sign of the index λ (sgnκ ≡ sgnλ κ ) indicates whether we are considering par- where c κ (c −κ ) creates (annihilates) a particle with quantum numbers κ and energy κ . Note that the creation and annihilation operators obey c † κ = c −κ , and we de- A many-body environment eigenstate |u env is constructed by repeatedly applying creation operators κ > 0 (in proper order) on the vacuum where |0 env is the environment vacuum, and u is a set of positive indices κ. In a non-interacting environment, the single-particle energies are additive, with the eigenenergy of the many-body state |u given by In contrast to the environment, we keep the system general, i.e., its Hamiltonian is where N sys is the number of eigenstates in the system. Furthermore, for simplicity, we assume that the microscopic coupling between the system and environment involves a single environment particle, allowing us to write . . 10. Diagrammatic rules applied to the specific systemenvironment setup described in (81)-(85). (a) The diagrams remain structurally the same as those in Fig. 4 and we have explicitly included the system states i, j, f, g, n, m for clarity. For each application of the perturbation V , a (anti-) particle with quantum numbers κ is added to the environment. The environment indices κ1 . . . κµ+ν are chosen to run clockwise (dashed arrow) through the diagram. (b) A factor in the perturbative series on the upper branch, that adds a particle with quantum numbers κ to the environment, while changing the system state from |m → |n . The operator cκ acts on the environment, whereas all other factors are numbers. The energy δ − β tracks the total number of particles (on the upper branch) that have been added to the environment. (c) Similar to (b) but for the lower branch of the Keldysh contour. We use the fact that the perturbation is Hermitian κ Vnmκcκ = κ V * mnκ c † κ . Note, that the Hermitian conjugation on the lower branch (see Fig. 4) both flips the order of the environment operators and conjugates them individually c † κ → cκ. As a result, the environment correlator assumes the form shown in Eq. (86).
where the amplitudes V * nmκ = V mn−κ of the perturbation guarantee Hermiticity. Last, we assume that the amplitudes V nmκ do not depend on the continuous (momentum) part of the subscript κ. We can thus interchange V nmκ and V nmλ whenever it is convenient.
We insert the explicit representation of the setup (81)-(85) into the expansion (71) for the propagator G, see Fig. 10. Each system-environment scattering event V now changes the system state and describes the emission or absorption of an environment particle, i.e., adding or removing a particle in the environment, respectively. As usual [30,37,47], we embody this emission (absorp-tion) process by a thin line attached to each vertex V in the diagrams 10(a). The additive property (83) of the single-particle environment energies κ leads to a simplifications in the matrix elements G ijf g νµ , cf. Figs. 10(b)-(c) and 4(b). The denominator of the recurrence relation (14) for the time evolution, which appears in (71), now changes by one single-particle environment energy κ after each application of the perturbation V , see Fig. 10(b). We thus use κ to update the energy denominator due to the environment δ ± β , which tracks the total number of particles added to the environment after β applications of the perturbation V on the upper (−) or lower (+) Keldysh branches. We use a clockwise labelling of the environment indices, see Fig. 10

Wick's Theorem
The evaluation of locally equilibrated environment correlators to describe the steady-state, is a great strength of the STCL methodology. It is a consequence of the combination of the back and forth propagation in Eq. (46). We can therefore apply Wick's theorem [52] when evaluating Eq. (86), which is thus recursively reduced by The sign in Eq. (87) indicates fermionic (−) or bosonic (+) environment modes [68]. Note that the normal-order contribution appearing in Wick's theorem vanishes because the environment is locally-equilibrated in the distant past [52,69]. We thus reduce the α-point correlator (86) into a sum over products of two point correlators where n A is the Fermi-Dirac (A=F) or Bose-Einstein (A=B) distribution, with µ λ the chemical potential (0 for massless bosons) and T λ the temperature of the discrete degree of freedom λ, δ is the Kronecker delta, and we again use the fact that the environment is locally equilibrated at t = −∞. Furthermore, in the far distant past, the environment is in an incoherent superposition P env u of states |u env which each have well defined particle numbers, see Eq. (82). There is thus no coherent superposition of states with different numbers of particles, and expectation values c κ = 0 vanish (as well as other odd correlators). Within a diagrammatic formulation, Eqs. (87) and (88) are implemented by: 2. Crossing fermionic lines produce additional minus signs.
In Fig. 11, we provide a specific example for the calculation of G 31 . Pairs that reside on the same Keldysh branch do not influence the environment state, when considered together. On the other hand, a pair connecting the upper and lower branch is a physical process where the system emits a particle (or anti-particle) into the environment, i.e., these diagrams will contribute to the current flow in the steady-state. Diagrammatically, the (thin) Wick contraction lines act in the same way as the (weighted) trace line with an open circle that appears on the right of every diagram, see also Fig. 4.
Using the diagrammatic rules for the Wick decomposition, we build each superoperatorG νµ (or G νµ ) as the sum over its Wick contributions where we have introduced the Wick index w, see below for an example. The Wick index w shows up in the same way in the calculation of other superoperators such asR νµ (77) orS νµ (76). At second order, there is only one contraction c κ1 c κ2 and the Wick index can be dropped. At fourth (second non-vanishing) order the four-point correlator decomposes as where the ± sign accounts for the bosonic (+) or fermionic (−) nature of the particles, see Fig. 11(b). We denote the three contractions in Eq. (90) by the Wick indices w = a, b, c, respectively. Note that, the linear additivity of energies (83) coupled with Wick's theorem implies which we will use to simplify expressions for G νµ later. At this point it is useful to recall the hierarchy of the indices α, ν, µ, and w as we will use them extensively for explicit calculations of the rates. Furthemore, we note that in the Pauli STCL rates, complex conjugation is equivalent to flipping the upper and lower Keldysh branches. This means thatG if νµ =G if * µν , and that for each contraction w there exists a contraction w such thatG if νµw =G if * µνw . In some cases we have w = w , e.g., at second order where there is a single contraction.

B. Explicit rates
We focus on the Pauli STCL as it is sufficient for the purpose of calculating steady-state occupations and currents. First, we show how to recover Fermi's golden rule as encoded in the leading (or second) order in V Pauli STCL master equation with its associated ratesS if 2 , see Eq. (96) and Fig. 12. These rates involve the exchange of a single particle between the system and an environment reservoir, and are commonly dubbed sequential tunnelling rates in mesoscopic research [20,21], see Figs. 12(b) and (c). We proceed with the fourth-order in V Pauli STCL, where the ratesS if 4 are associated with the exchange of up to two particles (potentially from different reservoirs) between the system and the environment. These rates include co-and pair-tunnelling [20], see Figs. 13(b) and (c), as well as virtually-assisted sequential-tunnelling which are responsible for level broadening and renormalisation [31,47], see Fig. 14(b). The strength of the STCL already manifests at this order: while such fourth-order rates diverge in the T-matrix formulation [29,45], the STCL, by construction, systematically removes these divergences [50]. For illustration, we will demonstrate a perfect agreement of the STCL rates with exact results for a single-particle level, see Fig. 18.

Fermi's golden rule
Fermi's golden rule [70] successfully describes transitions in open quantum systems with weak couplings to large environments; it commonly includes solely the diagonal of the system's density matrix (Pauli) but can be extended to describe the system's coherences as well. We consider the Pauli STCL (42) at lowest non-vanishing (second) order in the perturbation V given in Eq. (85). The first order term has an odd number of environment operators c κ in the correlator (86) and thus vanishes, see Fig. 11(c). At second order, we have three diagrams forG. We apply the rules from Figs. 4, 5, 10, and 11 to find that we display in Fig. 12(a) (recall δχ nm ≡ χ n − χ m ). TheG 11 diagrams are associated with a sequential tunnelling event, see Fig. 12(b), where a particle tunnels in or out of the system. TheG 20 andG 02 diagrams on the other hand are associated with probability conserving events, see Eq. (80), where an excitation tunnels back and forth between the system and the environment along one of the Keldysh branches, with the setup finally returning to the initial state. Note that in our discussion, we will alternate between the superoperatorsS,G and their matrix elementsS if ,G if , according to convenience. We substitute Eq. (93) into Eqs. (76) and (78) to compute explicit values for the second-order Pauli STCL (equivalently, Fermi's golden rule) rates. We recall that κ = (λ κ , k κ ) to convert the sum over momentum k κ in Eq. (93) into an integral over energy while maintaining the sum over discrete degrees of freedom λ. The density of states D λ ∆ λ ( ) associated with λ is composed of two parts. The first D λ is energy independent and is a typical scale of the density of states, while the second ∆ λ ( ) is a dimensionless function that encodes the dependence on energy. Furthermore, we evaluate the expectation value of the (fermionic) environment operators using (88) and introduce the real and dimensionless spectral function We can then write the second-order Pauli STCL (Fermi's golden rule) rates in the form which is the sum of the three diagrams shown in Fig. 12(d). The first term (S 11 ) describes sequential tunnelling and is identical to the rates obtained from Fermi's golden rule. The second term is a manifestation of conservation of probability, see Eq. (80). The explicit rates for the lowest-order full (as opposed to Pauli) STCL for quadratic environments is not required for our discussion but can be found in Ref. [48].

Co-and pair-tunnelling
Fourth-order processes are those that arise from matrix elements with ν +µ = α = 4. In contrast to the T-matrix ratesR 4 ∼ 1/η associated with the same physical processes, the STCL ratesS 4 ∝ η 0 are convergent in the limit η → 0. We start by computing Pauli STCL co-and pair-tunnelling ratesS 22 , before moving to the ratesS 31 Illustration of the co-and pair-tunnelling processes associated with the diagramsG22a andG 22b . A particle or antiparticle (black dot) tunnels into the system from the environment simultaneously with a second particle or anti-particle (white dot) from the same or a different reservoir. If both particles tunnelling to the system are electrons (or holes) this process is known as pair tunnelling. On the other hand, if one particle is a hole and the other is an electron, the process is known as a cotunnelling event. An elastic cotunnelling event is one which leaves the system state unchanged while in an inelastic process the initial and final system states differ. (c) Illustration of processes associated with the diagramG22c. A particle or anti-particle tunnels back and forth between the system and the environment along each of the Keldysh branches. This process leaves the environment unchanged, while it can modify the state of the system. (d) Diagrammatic representation of the three contractions forS22 with all indices dropped.
andS 13 which renormalise the lower-order sequentialtunneling rateS 11 from Fermi's golden rule. We avoid computing the ratesS 40 andS 04 , as they are merely a manifestation of conservation of probability (80). The matrix elements required to computeS 22 , contains two occurrences of the perturbation on each of the Keldysh branches. We combine Eq. (98) with the contractions (90) to obtain the matrix elements of the superoperatorsG 22a ,G 22b , andG 22c , see Fig. 13(a). These matrix elements correspond to new processes beyond sequential tunnelling, where the initial and final state differ by an even (possibly 0) number of particles. Co-and pair-tunnelling arise from the dia-gramsG 22a andG 22b , see Fig. 13(b). In a cotunnelling process, an electron or hole tunnels onto the system while at the same time a second particle of the same kind tunnels out of the system. These two particles may originate from the same or different reservoirs, and a cotunnelling process may change the system state, but not its charge.
If the initial and final system states are identical such a process is termed elastic cotunnelling. A pair-tunnelling process on the other hand, occurs when two electrons (or holes) from the same or different reservoirs simultaneously tunnel into the system, thus changing the state and charge of the system. TheG 22c diagrams gives rise to processes where two particles (electrons or holes), one on each Keldysh branch, tunnel in and out of the system, leaving the environment unchanged, see Fig. 13(c). The ratesS 22c associated with these processes vanish, as we show in Eq. (107). Non-crossing co-and pair-tunnelling. We first consider the type a Wick contraction from Eq. (90). We note that the associated T-matrix co-and pair-tunnelling rates, evaluated by substituting Eq. (77) into Eq. (98) for the a contraction, scale as and thus diverge in the limit η → 0 [50]. On the other hand, combining the recurrence relation (76) for the STCL with the same contraction provides us with the ratesS which contains corrections from two consecutive sequential-tunnelling eventsG im 11G mf 11 , see Fig. 13(d). These corrections cancel the diverging part of (99) and leave the Pauli STCL ratesS if 22a convergent in the limit η → 0. Diagrammatically, corrections due to products of lower order occur whenever a diagram forG can be cut vertically without slicing a contraction line.
We substitute the matrix elementsG if 22a [combining Eqs. (98) and (90)] andG if 11 from (93) into Eq. (100), transform the sums over momenta to integrals over energy and assume constant density of states (94). We then evaluate the expectation values (88) for environment operators and take the η → 0 limit to obtain the expression with dimensionless integrals I and J defined in Eqs. (102) and (103). To arrive at this formula, several lengthy but straightforward steps are required. These are outlined in detail in Appendix B. Eq. (100) is a typical example of the subtraction of reducible contributions (53), where two cancelling diverging terms in 1/η appear (one fromG 22a and one from the productG 11G11 ). In the course of this calculation, we have expanded the expressions in Eq. (101) as Laurent series in η before taking the η → 0 limit; in doing so, we have replaced the derivatives with respect to η arising from this expansion by derivatives with respect to χ.
The first line (101a) in the cotunnelling rateS if 22a is a convergent contribution which is attributed to tunnelling through two different system states n = m on each of the Keldysh branches, see Fig. 13(a). The same contribution appears in the T-matrix rateR if 22a . The second and third lines (101b) correspond to a co-or pair-tunelling process where the intermediate states on both Keldysh branches are the same, i.e., n = m. The second line coincides with the contribution obtained in the T-matrix approach using the phenomenological regularisation scheme developed in Ref. [45], see Ref. [48]. The third line (101c) contains additional corrections which are of the same order as (101b) but are missing in the regularisation of Ref. [45], see also Refs. [29,44,47,66]. For a detailed discussion comparing Eqs. (101) to the corresponding result using the SRT approach, we refer the reader to Ref. [48].
In Eqs. (101), we introduced two dimensionless integrals I and J. To evaluate them, we assume constant densities of states ∆ λ ( ) = 1 for each environment degree of freedom and further assume a constant temperature T across the entire environment in the distant past. Under these assumptions, the first integral I is (see Appendix C) where n B is the Bose-Einstein distribution at temperature T and n ψ (δ) = ψ(1/2 + iδ/2πT ) is a generalised distribution in terms of the digamma function ψ. Note that, the latter is related to both the Fermi-Dirac and Bose-Einstein distributions, see Appendix C. Furthermore, we have dropped the explicit dependence on temperature T in all distribution functions n A , as we assume all temperatures in the setup to be equal. The second integral J is (see Appendix C) where Λ λ is an ultraviolet cutoff for the continuum variable k in the environment reservoir λ. Typically, in electronic systems, this is the bandwidth of the Fermi reservoirs. Due to the derivatives with respect to χ acting on the integrals J in Eq. (101c), the ultraviolet cutoffs do not affect the ratesS 22a in the wide band limit Λ → ∞, see Appendix C.
Crossing co-and pair-tunnelling. Next, we consider the b contraction, see Eq. (90), which we combine with Eqs. (76) and (77) This process includes crossing contraction lines, see No-tunnelling. We now turn to the third contraction c for the rates associated withG 22 . These are associated with processes that do not change the state of the environment. Substituting the contraction c into the expression for the matrix elements ofG 22 , see Eq. (98), we find that the associated T-matrix rates (77) scale as This expression diverges in the η → 0 limit for equal initial and final states and vanishes otherwise. It therefore does not contribute to the rate equation (or current rates see Section VI) and was therefore discarded in the phenomenological regularisation scheme of Ref. [45]. The corresponding STCL rates, see Fig. 13 and vanish in the η → 0 limit, further justifying the omission ofR if 22c in Ref. [45]. The two contributionsS 22a andS 22b to the STCL rates are associated with co-and pair-tunnelling, which are distinct from the sequential-tunnelling process at second order. In the following, we study further contributions at fourth-order, that renormalise the lower-order rates and are interpreted as virtually-assisted sequential tunnelling. These contributions arise from theG 31 andG 13 superoperators (67) and are as important as the cotunnelling rates (101) and (105) in an exact expansion.

Virtually-assisted sequential tunnelling
We now compute the virtually-assisted sequentialtunnelling ratesS 31 andS 13 in the Pauli STCL master equation, see Eqs. (76) and (78). To this end, we use our diagrammatic rules forG from Sec. V A to writẽ We combine (108) with the contractions (90) to obtain the three propagatorsG 31a ,G 31b , andG 31c , see Figs. 14(a) and (b), which can be understood as sequential-tunnelling processes where one of the Keldysh branches takes an indirect path to the final state. We use the recurrence relation (76) 14. Virtually-assisted sequential-tunnelling processes.
(a) Fourth-order diagrams forG31. (b) Illustration of the physical processes associated with the diagrams in (a): A particle or anti-particle (black dot) tunnels into (full arrow) the system on each of the Keldysh branches. Simultaneously, a second (anti-) particle (white dot) enters (dotted arrow) the system before returning to the environment (along the upper Keldysh branch). These processes renormalise the sequential tunnelling (96) and are thus often ignored [46]. (c) Diagrams forS31.

contractions (90) to obtaiñ
see Fig. 14 , c), (b, b) and (c, a). As forS 22b (104), the b contraction (109b) does not involve products of the lower-order terms as this diagram cannot be split. The explicit calculations of the rates (109) involve similar methods as the calculation of the cotunnelling ratesS 22 and can be found in Appendix B.
We have now computed all fourth-order rates that are associated to physical processes. The remaining fourth- Probability conserving diagrams Sα0 at fourthorder α = 4. (a) Fourth-order diagrams forG40. These terms each contain a Kronecker delta, which enforces identical initial (i) and final (f ) system states. (b) Illustration of processes described by the diagrams in (a) resembling those forG22c in Fig. 13(c), except for the fact that now the entire setup, rather than just the environment, initial and final states must be the same. The different diagrams in (a) correspond to different time orderings of the events in (b). (c) Diagrams forS40. order contributionsS 40 andS 04 ensure conservation of probability. Their sum can be obtained from the super-operatorsS 22 ,S 31 , andS 13 or calculated explicitly as for the other fourth-order rates. We take the former route though, for completeness, we now briefly discuss the general structure ofS 40 andS 04 .

Probability conservation at fourth-order
The probability conserving processes at fourth-order involve the propagator Fig. 15(a). Combined with the contractions (90), the matrix elements (110) correspond to physical processes where the system goes through three intermediate states on one of the branches of the Keldysh contour, before returning to the initial state. This can be thought of as two particles, that tunnel in and out of the environment only to return to the original state, see Fig. 15(b). We combine these fourth-order terms with the pairs of second-order terms according to the recurrence rule (68), see Fig. 15, to obtaiñ where we note thatS if 40w = δ ifS ii 40w for each contraction w. Instead of computing the terms in (111), we make use of the conservation of probability (80) to compute the sum (112)

Steady-state system probability distribution
We now have all expressions needed for the computation of the Pauli STCL rates up to fourth-order. The latter then give us access to the steady-state probability distribution of finding the system in state |n . Here, P (2α) n denote the contributions of order 2α to the probabilities P n (with odd terms vanishing, see Fig. 11). This is the first observable quantity that we can compute and showcases the strength of the STCL. We insert the formal expansions for P n andS nm into the steady-state condition (42) and solve it order by order (all odd orders vanish) to obtain the conditions In order to test our scheme and for illustration, we calculate the STCL rates up to fourth-order for a noninteracting model and compare the results of Eq. (114) with the expansion of the exact solution.

C. Non-interacting system
In non-interacting setups the steady-state probability distribution of the system states and current between the reservoirs through the system can be computed exactly, e.g., using equations of motion for Green's functions or scattering matrices [20,[71][72][73], see also Appendix D. The non-interacting setup, see Fig. 16(a), is thus an ideal platform to test the STCL. 16. Non-interacting lead-dot-lead setup and processes.
(a) The setup involves two reservoirs/leads at temperature T and chemical potentials ±µ/2, that are connected by a hopping amplitude J/ √ 2 to a single non-interacting fermionic mode (the dot-level) with energy 0, see the discussion around Eqs. (115)-(117). (b) Addition spectrum of relevant processes at equilibrium, combining two reservoirs into a single one (grey shaded area). In a sequential tunnelling process, a particle (or hole) enters the system (thin full arrow) and leaves a hole (particle) behind. In the probability conserving trivial processes (S20 andS02), a particle (or hole) visits the environment (dotted arrow) and returns back. A fourth-order, virtually-assisted sequential tunnelling combines the two preceding processes coherently. (c) Cotunnelling process in an out-of-equilibrium setup. A particle or hole exits the system on one side, while another enters on the other (or same) side. This process leaves the system unchanged and its contribution to probabilities is compensated by a probability conserving rate (112).

Setup
The non-interacting resonant level, i.e., a single fermionic level (or quantum dot) coupled to Fermi leads, see Fig. 16, can be realised in gate-defined electronic devices, see Fig. 1(a). Its Hamiltonian is composed of the same three parts as our generic model in Sec. III A with a quadratic environment as described in Sec. V, see Fig. 16(a) for a sketch. The system involves a single fermionic mode (with energy 0 ) and is described by the Hamiltonian where d † 0 (d 0 ) creates (annihilates) an electron in the system. The environment is composed of left (l = L) and right (l = R) leads and is described by where c † lk (c lk ) creates (annihilates) a particle with momentum k and energy k in lead l. The left and right environments are assumed to be at the same temperature T and have chemical potentials µ/2 and −µ/2, respectively, with the total bias given by µ. The system-environment coupling is given by the tunnelling Hamiltonian that couples the system to each reservoir in the environment with hopping amplitude J/ √ 2. We choose this tunnelling amplitude such that in the case µ = 0 the setup is equivalent to a single mode coupled to a single reservoir with amplitude J.

Diagrams
We consider the STCL rates up to fourth-order in the perturbation, including the expressions (96) for Fermi's golden rule, co-and pair-tunnelling (101), (105), (107), virtually assisted sequential tunnelling (B13)-(B15), and the fourth-order probability conserving rates (112), see Fig. 16. As there are only two system states for the non-interacting level, empty |0 or full |1 , there are only a small number of diagrams contributing to the STCL rates, see Fig. 17. The sequential tunnelling ratesS 11 , familiar from Fermi's golden rules change the occupation of the dot from full (|1 ) to empty (|0 ) or viceversa, see Figs. 16(b) and 17(a). During a virtuallyassisted sequential-tunnelling process, an (anti) particle is removed from the environment and changes the system state on both Keldysh branches, while another (anti) particle is removed from the environment before returning on one of the two Keldysh branches, see Figs. 16(b) and 17(b). There are only elastic cotunnelling rates in the non-interacting setup, i.e., rates that leave the system unchanged and therefore do not contribute to the rate equation (42), see Figs. 16(c) and 17(c). These elastic processes do, however, contribute to the steady-state current through the system as we will see in Sec. VI C.

Equilibrium
We start with the equilibrium situation µ = 0. In the steady-state, the exact probability P 1 of the level being occupied is [74] (see also Appendix D) (a) STCL diagram for a sequential tunnelling event, where an electron from the environment enters the system, changing the latter from empty |0 to full |1 . (b) STCL diagram for virtually-assisted sequential tunnelling events where an electron from the environment enters the system. Simultaneously this electron leaves and returns to the environment (potentially in a different location). The b contraction does not contribute for the non-interacting setup. (c) Diagram for an elastic cotunnelling rate, where an electron enters the system before leaving again. These rates do not contribute to the rate equation as they leave the system invariant (elastic process) and are compensated by a probability conserving rate (112). They do however contribute to the current as we will see in Section VI. As for the virtually-assisted processes, the b contraction does not contribute to the non-interacting setup rates, while theS22c contributions always vanish.
where we have introduced the width Γ 0 = 2π|J| 2 D. The exact probability (118) can be expanded in powers of V ∝ J, or equivalently Γ 0 ∝ |J| 2 as odd terms vanish, and thus serves as a benchmark for the STCL results, which we now compute. We insert the model (115)-(117) into our expressions for the sequential rates (96) and obtaiñ where the temperature in the Fermi-Dirac distribution n F is implicit. These rates are familiar from Fermi's golden rule and correspond to an electron entering (leaving) the system on both the upper and lower Keldysh branches, see Fig. 17(a). The second-order diagonal elements of the rate matrixS 00 2 = −S 01 2 andS 11 2 = −S 10 2 are obtained directly from the conservation of probability (79). We solve the lowest-order steady-state constraint (114a) and obtain where P  At this order, the results of the STCL and T-matrix approaches coincide and both give rise to the same result as the one obtained from Fermi's golden rule (120), see Fig. 18(a). Furthermore, the result coincides with the expansion of the exact result (118), see Appendix D. For an infinitely sharp non-interacting level, the result (120) is exact, however, the coupling to the environment broadens the levels [20], which manifests at fourth-order in the rates.
At fourth-order there are three types of rates, see Figs. 13, 14, and 15. The virtually-assisted sequential tunnelling ratesS 31 andS 13 , change the configuration of the system, see Fig. 17(b). We use the formulas (B13)-(B15) in Appendix B to compute the virtually-assisted sequential-tunnelling diagrams, see Fig. 17(b), for the non-interacting setup. We obtain the fourth-order cor-rectionS to the rate for changing the system state from empty to occupied, where ψ is the trigamma function. The rateS 10 4 for the reverse process can be found by replacing 0 → − 0 in Eq. (121). The probability conserving (S 40 +S 04 ) rates at fourth-order are again found by enforcing Eq. (80). We insert (121) and (119) into the constraint (114) for the next-to-leading order steady-state probability distribution correction P (2) and obtain as well as P as required by conservation of probability. The expression (122), is identical to the one obtained from an expansion of the exact result (118) and is plotted in Fig. 18 for the specific case Γ 0 /T = π. The result (122) decays as ∝ Γ 0 / 0 for large onsite energies (such that 0 exp(− 0 /T ) Γ 0 ) and thus embodies the broadening of the level due to its coupling to the environment, see for example Refs. [20,74] for a detailed overview.
We conclude that the STCL approach provides a useful tool to perturbatively compute the steady-state distribution of open systems. It defines a time-local master equation with proper ratesS if that are convergent in the η → 0 limit. In a next step, we take the setup out of equilibrium µ = 0 and investigate the current flowing through the system.

VI. CURRENTS
The most common transport observable in mesoscopic research is the electrical current [19,21,[75][76][77], where electric charge is flowing in or out of leads attached to the system. Other types of currents can be defined that are associated with the environment reservoirs, e.g., in the case of a spin-full electronic system, each lead l is further split into two reservoirs r = l, σ by the spin degree of freedom σ. It is then possible to consider both electrical and spin currents. Here, we extend the diagrammatic formulation of the STCL rates S from sections IV and V to include current flow out-of-equilibrium.

A. Definition and Derivation
In our discussion, we consider currents to (from) the reservoirs that contain a mean number of particles Tr(ρN λ ), where denotes the reservoir number operator and we use λ = +r, and −λ = −r (with the + describing particles) instead of just r to keep the notation compatible with previous sections. The particle current ∂ t Tr(ρN λ ) can always be thought of as an anti-particle current flowing in the other direction ∂ t Tr(ρN −λ ) = −∂ t Tr(ρN λ ). Furthermore, the particle currents ∂ t Tr(ρN λ ) in and out of the different reservoirs, multiplied by the charge q carried by each particle, give rise to physical and measurable currents I λ . For electrical currents, the charge of each particle (electron) is −e, with the elementary charge e, while for spin currents the 'charge' of each particle is ± /2, depending on the reservoir, i.e., the charge q r may depend on the reservoir (or even the momentum k in the case of a heat current). This leads us to define the charge operator which is the charge operator associated with the reservoir λ and the single particle charge q λ . Note that the number (and charge) operators commute with the unperturbed Hamiltonian [H 0 , N λ ] = 0. In fact, our considerations apply for any environment operator Q env that satisfies [Q env , H env ] = 0. This criterion, allows us to unambiguously promote the operator N λ to a superoperator when desired, see Appendix E. Furthermore, the vanishing commutator [H 0 , N λ ] = 0 implies that no currents flow when the system and environment are decoupled, i.e., when the Hamiltonian V → 0 vanishes. The current I λ (t) is defined as the time derivative of the expectation value of the charge where the density matrix ρ describes the coupled systemenvironment setup. Working in the Schrödinger picture, the time derivative acts solely on ρ(t). An alternative path, particularly common in the RT approach, uses the Heisenberg equation of motion, and thus [H, Q λ ] to encode the time dependence with a static ρ.
In the limit η → 0 and t 0 → −∞ the system is in the steady-state at t ≈ 0 (∂ t ρ p = 0), though out-ofequilibrium with a finite current flow across the system (∂ t ρ = 0). The large environment inhibits a direct solution of Eq. (125), cf. the similar challenge in solving the von Neumann equation (4). Ostensibly, we would like to replace ρ by the projected density matrix ρ p in Eq. (125). However, the latter does not include the evolution of the environment as it has been projected out, i.e., Tr[Q λ ρ p ] = Q 0 λ is constant for all ρ p . Furthermore, the projected space has reached the steady-stateρ p for t ≈ 0 and therefore does not encode any dynamics. Either of the latter two arguments is sufficient to show that a simple substitution ρ → ρ p in Eq. (125) leads to vanishing currents and contains no information The determination of currents thus starts from the full density matrix ρ which includes the non-trivial evolution of the environment. Even though the projected space has reached a steady-state, the full setup has not. Common treatments of currents (as in the T-matrix approach) rely on transition rates, in the form of a rate matrix, acting on steady-state populations [20]. These transition rates change the environment charge and produce finite currents into and out of the reservoirs. In the following, we obtain a similar description for the STCL. Specifically, we find a set of current rates S λ , such that where we add the reservoir index λ to the generator; we will show below that the rates S λ can be obtained directly from the STCL rates S by filtering/weighting different processes.

The current generator S λ
We follow a similar program as in the derivation of the STCL generator S in section III D. First, we construct the full density matrix at time t from the projected one at time t 0 using the evolution U from Eq. (23). We then use the inverse projected evolution (50) to obtain the full density matrix at time t from the projected one at the same time Here, U −1 P = U −1 0 (I +G) −1 propagates the projected density matrix ρ p (t) back in time to t 0 , where ρ p (t 0 ) = ρ(t 0 ), cf. Eqs. (32) and (50). The full density matrix ρ(t) is then forward propagated using the time evolution U. Note that, unlike in Sec. III, we do not immediately apply the projector again once we arrive at time t.
We substitute (128) into Eq. (125) to obtain a timelocal expression for the currents in terms of the projected density matrix where we used TrP = Tr, see Eq. (29), to insert a projector before the trace and thus define the projected space superoperators PQ λ UU −1 P . At this point, it seems that it were sufficient to compute the matrix elements (PQ λ UU −1 P ) ijf g , order by order in V , to reconstruct the current I λ . Unfortunately, these matrix elements diverge in the limit η → 0 (irrespective of the derivative ∂ t ), even though the trace over them, and therefore the current, is finite. We overcome this problem by a subtraction of superoperators whose matrix elements diverge identically to PQ λ UU −1 P but that vanish when traced. We will thus be left with a set of convergent rates S ijf g λ . Furthermore, we will show that the resulting rates S ijf g λ evaluated for quadratic environments are obtained by adding a set of Kronecker deltas to the diagrams used to compute the STCL rates S ijf g . This filter in the calculation of diagrams implies that the currents I λ can be obtained with little additional effort as compared with the determination of the distribution ρ p .
The charge transfer propagator. As the current is the time derivative of the expectation value of the charge, we can subtract the (vanishing) time derivative of the charge at time t 0 where, different from (129), the charge operator Q λ is now applied prior to the forward propagation in time U. Subtracting (130) from (129) we obtain where we have used TrP = Tr to insert the projector P, made use of the definition (50) of U −1 P , used [Q λ , U 0 ] = 0, and introduced the charge-difference propagator Physically, Eq. (131) can be understood as propagating the projected density matrix ρ p backward in time to the distant past, before evolving it forward again. On the way forward, the charge transfer propagator in Eq. (132) tracks the number of particles in reservoir r that have been created minus the ones that have been annihilated. This physical interpretation of the propagator G λ is commonly used in T-matrix calculations of currents [20].
Applying the derivative. We apply the time derivative ∂ t in Eq. (131), and use the simplification (62) (which is valid at t 0 → −∞). The expressions then all become time independent in the limit η → 0 and we thus use the steady-state projected density matrix to obtain Here, we have introduced the time-differentiated counterparts G λ andĠ λ of the charge transfer propagator G λ . These two terms are obtained via the expan- of the projected space superoperator in Eq. (133) are not well defined in the limit η → 0 (even though the current is), just as was the case for (129). The easiest way to verify that (134) is ill defined in the limit η → 0 is to evaluate the matrix elements explicitly at fourth-order, which can be done using the methods developed in Sec. V and Appendix B.
Adding zero. In order to obtain convergent rates, we add a compensating term to Eq. (133), that vanishes when traced over, but contains divergent matrix elements that cancel the diverging part of Eq. (134). In Eq. (133) the projected space superoperators that propagate backwards (I+G) −1 do not contain any occurrences of N λ , whereas those that propagate forward with G λ +Ġ λ do. This asymmetry between the backward and forward propagation is the root causing the diverging matrix elements (134). We therefore add the vanishing term to the current (129), where we have introduced the backward charge transfer propagator It will become clear later, when adapting our diagrammatic approach to the currents, how the term (135) regularises the matrix elements (134). To verify that Eq. (135) is indeed vanishing, we rewrite G andĠ in terms of G using Eqs. (64) and (65), and expand all occurrences of G in powers of V . All of the resulting terms take on the form where A is composed of a product of charge transfer propagators G λβ and regular propagators G β . To verify that Eq. (137) vanishes, we rewrite where we used TrP = Tr, see Eq. (29). We then use the invariance of the trace of an operator under unitary transformation to write We sum Eqs. (133) and (135), expand the result in V , and compare it order by order with (127). We thus arrive at the series expansion for the STCL current generator which can be used to reconstruct the STCL current generator The expansion of the current generator (140) is very similar to the expansion (63) for the STCL generator S but with an additional regularising term in the sum. Constructing the formally exact convergent current rates (140) is one of the main results of this work.

Diagrammatics
We now develop a diagrammatic approach to compute the (divergent) matrix elements of G λνµ and thereof the (convergent) matrix elements of S λνµ . Here, as for G νµ , the νµ indices indicate the number of scattering events on the upper and lower branches respectively, and thus cf. Eq. (92). We insert the expansion (26) of U in terms of the expansion (10) of the unitary evolution U into Eq. (132) and obtain . . Next, we use the cyclic nature of the trace and the fact that the charge operator commutes with the unperturbed environment distribution [N λ , ρ 0 env ] = 0, to obtain which is identical to G ijf g νµ in Eq. (71) up to the commutator with N λ .
Diagrammatically, the commutator in G λνµ is obtained by taking the difference of two terms, where we insert N λ to the far left on one of the Keldysh branches, and the second with N λ inserted to the far right on the same Keldysh branch, cf. Figs. 10(a) and 19(a). On the other hand, mathematically, we could deal with the commutator [N λ , U µ U −1 0 ] by separately evaluating the two terms N λ U µ U −1 0 and U µ U −1 0 N λ before subtracting them. Here, we make use of the structure of U µ to simplify the commutator before performing any explicit calculation. We insert the expansion (10) into the commutator from (144) and obtain an expression of the form where Π 0 are the free propagators, see Eq. (15). The number operator N λ commutes with the unperturbed Hamiltonian H 0 and therefore also with the free propagator [N λ , Π 0 ] = 0, allowing us to rewrite Eq. (145) as a sum of terms each with one occurrence of the perturbation V replaced by the commutator [N λ , V ], see Fig. 19(b). Diagrammatically, Eq. (146) means that we go from G νµ to G λνµ by drawing µ copies of G νµ and replacing one of the systemenvironment scatterings (V ) on the lower Keldysh branch by [N λ , V ], see Fig. 19(c). We could equivalently choose the upper Keldysh branch and create ν copies, leading to the same final result.
Convergent matrix elements of S λνµ . Thanks to the addition of the vanishing terms (135), the relationship between S λνµ and S νµ is identical to the relation between G λνµ and G νµ , cf. Figs. 19(c) and (d). Transforming from the usual superoperator to the current (λ) superoperator, we again draw µ copies of the diagram and replace scattering events V on the lower branch with commutators [N λ , V ] (without the vanishing terms (135), commutators [N λ , V ] would only replace scattering events to the left of the leftmost cut, see Fig. 20). As the dependence on the slow switch-on parameter η lies entirely within the unperturbed propagator Π 0 , the exact form of the scattering events does not influence the convergence in the limit η → 0. Thus if the matrix elements S ijf g νµ are finite, the same will be true for the matrix elements S ijf g λνµ of the current generator.

Quadratic environment and Wick's theorem
We now consider the specific case of quadratic environments with linear coupling, see Sec. V. Here, the commutator [N λ , V ] assumes a particularly simple form which can be thought of as a filter associated with the environment operators c κ , see Figs. 21(a). Each contribution to the filter is a Kronecker delta which constrains the discrete degree of freedom λ κ associated with a scattering event to either λ or −λ. In the diagram for G λ , we have to insert this commutator subsequently for every scattering event on the lower branch of G λνµ , see Fig. 19(c). As a result, we obtain a simple diagrammatic rule that takes us from G to G λ via the substitution in the unperturbed environment correlator (86), with the filtering function  Hermitian conjugation on the lower branch, see Fig. 4 in Sec. IV.
Wick's theorem. The environment correlator (148) for the charge transfer propagator G λνµ is identical to the one for G νµ , Eq. (86), up to the prefactor δ λ composed of a set of Kronecker deltas for the discrete degrees of freedom λ s associated with the scatterings on the lower Keldysh branch. Hence, we can apply Wick's theorem in the same way as before and write G λνµ in terms of a sum over contractions with Wick index w, As a further simplification, we note that the filters associated with two contracted scattering events κ s and κ t on the lower Keldysh branch, such that κ s = −κ t , vanish (δ λ−λs −δ λλs + δ λ−λt −δ λλt ) c κs c κt . . . = 0, (151) see Fig. 21(c). Thus, the only filters that contribute to G λνµ are those associated with contractions that connect the upper and lower branch.
To unify the calculation of both G λνµw and G νµw , we introduce the constrained propagator G νµw (λ p , λ q , ...), see  Fig. 21(d). This constrained propagator is identical to the propagator G νµw except for the fact that the sum over discrete degrees of freedom λ p , λ q , ... associated to contractions that connect upper and lower Keldysh branches are not performed. Once G νµw (λ p , λ q , ...) has been computed one can immediately obtain either the charge difference propagator or the usual propagator see Fig. 22.
The relationships between the rates S νµw (λ p , λ q , ...), S νµw and S λνµw are identical to the relationships between the corresponding propagators G (and similarly for Pauli equivalents). This conclusion follows from equivalence of the rules for rates and propagators shown in Figs. 19(c) and (d). 22. The filtering scheme that produces Gνµw and G λνµw from the constrained propagator Gνµw(λp, λq, ...). Of all the sums over κ1, κ2, ...κν+µ, we make explicit the sums over the discrete part λp, λq, ... of contractions that connect upper and lower Keldysh branches, while all other sums (over discrete and continuous degrees of freedom) are implicit.
The filtering rules (148), (151), and (152) that distinguish the Pauli ratesS if from the current ratesS if λ are the same ones that are commonly employed in Tmatrix calculations [20,78] or the RT approach [32] and provide an intuitive physical description of the processes that contribute to currents: one simply counts the change in the environment's charge during a process to identify its contribution to transport.

Steady-state currents between reservoirs
We proceed with the calculation of the steady-state current I λ flowing out of the reservoir λ using the current ratesS nm λα and the steady-state probability distribution P n as obtained from the ratesS nm α in section V. We remind the reader that we work with the index λ > 0 instead of the reservoir index r in order to keep the notation consistent. The Pauli equivalent to (127), written in matrix element form, reads We expand this expression order by order and. obtain which allows us, along with P (0) n and P (2) n from Eq. (114), to compute the current from a reservoir λ in any setup with a quadratic environment up to fourth-order.

B. Current rates
We now apply the diagrammatic rules for the current generators, see Figs. [19][20][21][22][23], and obtain explicit expressions for the associated current rates; the latter are FIG. 23. Constrained Pauli diagrams, up to fourth-order, that contribute to the current for quadratic environments and linear coupling.S11(λ1): Sequential tunnelling rates involving an environment particle in reservoir λ1. The probability conserving ratesS20 andS02 do not contribute to the current as the two vertices lie on the same Keldysh branch and thus the initial and final environment states are the same, see Fig. 21(c).S31(λ1): The virtually-assisted sequential tunnelling diagrams that contribute to the current. These rates only have one leg connecting the upper and lower branches and therefore only one constrained environment index λ1. S22(λ1, λ2): The co-and pair-tunnelling diagrams that contribute to the current. Note that each leg that connects the upper and lower branches counts once, leading to two contributions per diagram.
closely related to the Pauli STCL rates we computed in Sec. V. Fermi's golden rule. At second order, there are three diagramsS 11 ,S 20 , andS 02 for the STCL rates. The first of these three is Fermi's golden rule, while the latter two contributions ensure conservation of probability (80). Summing up these three terms gives rise to the second-order STCL contributionS 2 . The current ratesS if λ only contain contributions from constrained diagrams that connect the upper and lower Keldysh branches. At second order, there is only one such dia-gramS 2 (λ 1 ) =S 11 (λ 1 ), see Fig 23(a). We compute the constrained diagram for Fermi's golden rule rates by not performing the discrete environment sum in Eq. (96) and obtainS We then use the filtering scheme (152) to obtain the second order current rates which describes sequential tunnelling where a particle of type λ tunnels out of (first term) or into (second term) the environment. Fourth-order. At fourth-order there are three types of constrained contributions,S 31 (λ 1 ),S 13 (λ 1 ), andS 22 (λ 1 , λ 2 ). The first two contain exactly one contraction that connects the upper and lower branch, irrespective of the Wick index w, see Fig. 23. These three diagrams are identical to the ones forS 31w orS 13w , except for the fact that we do not sum over the discrete degrees of freedom that connect upper and lower Keldysh branches, as in Eq. (155). We determine these constrained rates in Appendix B, along with the corresponding STCL rates and can use them to reconstruct the current rates and similarly forS if λ13 . The co-and pair-tunnelling current ratesS if λ22 arise only from the a and b contractions, see Fig. 23. These diagrams have two contractions connecting the upper and lower Keldysh branches, in contrast to the c contraction which has none, see Fig. 13. We constrain the contractions (λ 1 and λ 2 ) that connect the upper and lower branches in Eqs. (101) and (105), for w = a, b. Here the crossed out sum indices λ 1 , λ 2 indicate that we do not sum over any discrete environment degree of freedoms. We then again use the filtering scheme (152) to obtain the co-and pair-tunnelling current rates Summing the contributions fromS 31 ,S 13 , andS 22 , we obtain the total fourth-order current ratẽ We are now in a position to calculate both the steadystate system probability distribution and the environment currents to next-to-leading order. The former is calculated as detailed in Sec. V and then used along with the current rates (156) and (160) in the expression (154) for the currents up to fourth-order.

C. Non-interacting model
As a first application and test of the formalism, we focus on the non-interacting resonant-level (115)-(117) and take the setup out-of-equilibrium (µ = 0), implying that a steady-state current will flow across the device, see Fig. 24. As for the (equilibrium) probability distribution (118), it is possible to compute the resulting (outof-equilibrium) current exactly. This can be done using a scattering matrix or Green's function approach [20,73] and produces a current into the right (R) reservoir with e the unit charge, see Appendix D for a brief sketch of the calculation. The setup is fully described by the level width Γ 0 = 2π|J| 2 D, the single level energy 0 , the chemical potential difference µ between the leads, and the temperature T . The expression (161) can be expanded in powers of V (or equivalently Γ 0 ), leading to at lowest order. This same result can be obtained from the STCL. First we compute the steady-state probability P (0) 1 (P (0) 0 ) of finding the level occupied (empty) as in section V C, but under finite bias conditions. We then use these probabilities and the sequential current rates (156) in the expressions (154a) for the lowest-order current. As for the probabilities in Sec. V C, this lowest-order result (162) is exact in the case of infinitely sharp levels. At higher orders, see Refs. [20,78], the broadening of the level due to its coupling to the environment manifests. A straightforward but lengthy calculation using the fourthorder STCL rates provides the current I

VII. CONCLUSION AND OUTLOOK
In this work, we developed an operator-based diagrammatic approach to the steady-state time-convolutionless master equation. We greatly reduced the number of diagrams to be computed at any given order α when compared to a superoperator formulation, from 2 α to α + 1. We thus kept the complexity of STCL calculations as low as in the T-matrix approach, while remaining formally exact as in the real-time master equation. Going beyond the analysis of a steady-state system distribution, we extended the STCL to perturbatively evaluate the steadystate current through the system. We then applied our diagramatic approach to setups with non-interacting environments and a bilinear system-environment coupling for both steady-state distributions and currents. As an example, we verified our methodology on a noninteracting setup, a single-level coupled to leads, where we demonstrate perfect agreement between our expansion and the expansion of the exact result. These results show that the STCL is a versatile tool which can be used for practical perturbative calculations. The non-interacting resonant-level out-ofequilibrium, for Γ0/T = π and µ/T = 6. The inset provides a sketch of the configuration. The chemical potentials in the two leads are offset by µ. The level moves with 0, and passes through the bias window, when it lies between the two chemical potentials. The main figure shows the current flowing through the level when the system is driven out-ofequilibrium, lowest order (blue dashed), first correction (green dotted), their sum (black dotted-dashed) and the exact result (black full). The corresponding (blue and green) crosses are obtained from an expansion of the exact result and have to be matched by any formally exact perturbative method.
We identify a number of future research topics based on the work presented here, with three key examples: extensions to higher orders, resummation schemes, and dynamics. A sixth order analysis is well within reach as there are only 45 new diagrams to be computed. Specifically, these areG 33w ,G 42w , andG 51w for the fifteen Wick contraction contractions w at sixth order. The sixth order Pauli STCL ratesS 6 , where Kondo signatures are expected to become apparent [20], can then be constructed from these (and lower-order) terms. This is also the order at which certain backaction effects are expected to appear [53]. Furthermore, the STCL provides a formally exact expansion for the steady-state observables, and can therefore be used as a base for resummation schemes. While existing resummation schemes for the SRT method [31][32][33] can be directly applied to the superoperator formulation of the STCL, our operator-based approach enables the resummation of different sets of diagrams. Finally, while this work was exclusively concerned with steady-state properties, our operator-based simplifications can readily be applied to the TCL description of dynamics after a quench [43,48]. We conclude that the time-convolutionless master equation still holds many surprises and a large potential for future applications. The counter for the prefactor of η is shared between the upper and lower branches and runs from 1 to |z|. Here the specific index is z = (↓, ..., ↑, ... ↑, ↓). Note the lack of Hermitian conjugation on the lower branch (cf. Fig. 4), the minus sign on the lower branch scatterings, and the free propagator to the left of the leftmost scattering.
(a) 26. Ordered diagram for Gz, to be contrasted with the diagrams for Gνµ, see Fig. 10. As in Sec. V, the Wick index w is added by creating all pairings of the different environment operators. (a) The full diagram, which is related to Fig. 25 in the same way as Fig. 10(a) is related to Fig. 4(c). The scattering events V (black dots) pick up a contraction line (thin vertical lines) which are indexed clockwise from the bottom right. The system indices along the free propagation (dashed lines on upper and lower branch) are identical to those in Fig. 10(a), while the counter from 1 to α is shared between the upper and lower branch. (b)-(c) Factors that contribute to the diagram in (a) with a scattering event on the upper (b) or lower (c) Keldysh branch. Unlike in the unordered case of Fig. 10 the environment energy counter δ β is shared between the branches. Note the negative sign when the scattering event lies on the lower branch (c).

Remove the Hermitian conjugate operation on the lower branch.
These diagrammatic considerations carry over directly to the specific example of quadratic environments which we considered in section V, see Fig. 26. In contrast to our unordered diagrams, the diagrams in Fig. 26 contain a single environment energy counter δ , which is incre-mented by scattering events both on the upper and lower branches. We can further apply Wick's theorem to the diagrams for G z to obtain the diagrams G zw in exactly the same way as for the G νµ . The STCL generator terms S zw , are generated from the diagrams G zw in exactly the same ways as S νµw are generated from G νµ . Due to the ordering relation between upper and lower branches in the ordered diagrams, any cut must be vertical.
Appendix B: Fourth-order STCL rates In this appendix, we show the steps that lead to the explicit expressions for the fourth-order rates. We will take care to include the η dependence and only take the η → 0 limit for convergent expressions. We start with theS 22a rate, as it includes most of the technical difficulties which arise. We first write down the diagram from Fig. 13(a) as where the second term in the bracket (∝ δ nm ) arises from the diagram with a cut. As a next step, we change the environment sums into sums over the discrete parts of the environment and integrals over the continuous parts. We then perform a change of variables such that both terms in the bracket have the same denominator We then split this expression into two, evaluate the Kronecker delta, and perform some manipulations on the (convergent) n = m sum In the n = m sum, we have performed the 2 integral in the η → 0 limit and taken a partial fraction decomposition of the remaining denominator. The remaining integral is common in higher-order rate equation calculations. In Appendix C, we outline one way of performing it analytically which leads to the expression in Eq. (101a). The second line (B4) with m = n is also convergent, which we will now show, before computing it explicitly. We rewrite Eq. (B4) in the functional form where we have extracted a 1/η factor to guarantee that the integral is convergent, irrespective of the function g. The latter is given by We find the convergence condition on F by multiplying (B5) by η and use L( , η) = πδ( ) + O(η) to obtain which confirms that Eq. (B4) is convergent in the η → 0 limit. As a consequence of L'Hôpital's rule and the fact that F is convergent in the limit η → 0, we can write where we have further made use of the product rule of differentiation and L( , η) = πδ( ) + O(η). We substitute the expression (B6) for g back into this latest expression and obtain (B4) = 2π mλ1λ2 |V imλ1 | 2 |V mf λ2 | 2 ∂ η η C λ1 ( )C λ2 (δχ f i − ) (δχ im + ) 2 + η 2 + C λ1 (δχ mi )C λ2 ( ) (δχ mf + ) 2 + η 2 − C λ1 ( )C λ2 (δχ f m ) (δχ im + ) 2 + η 2 d , (B9) which is valid in the η → 0 limit. The first integrand in the last expression is the same contribution as the regularised T-matrix integral, see Ref [48], whereas the next two are corrections. In Appendix C, we show how these can be cast into the expressions in Eqs. (101b) and (101c).
where in the second line, we performed the 2 integral, a partial fraction decomposition, and relabelled 1 → . The ± sign is + for bosons and − for fermions. The remaining integrals, after a partial fraction decomposition, are again standard and an analytic expression for them can be found in Appendix C. We thus obtain the expression in Eq. (105). The lastS 22 contribution arises from the c contractioñ and vanishes in the η → 0 limit. We can now repeat the same steps as above but for theS 31 andS 13 diagrams. For the first contraction a we obtaiñ where we have again used L'Hôpital's rule and the integrals from Appendix C. The b contraction is again simpler, leading toS where the ± differentiates bosons (+) and fermions (−). The lastS 31 contribution arises from the c contraction. It requires the use of L'Hôpital's rule and the integrals from Appendix C to obtaiñ TheS 13 =S * 31 terms are obtained by complex conjugation. We do not compute theS 40 orS 04 terms as their sum can be obtained easily from the conservation of probability, see Eq. (80).
The last terms in each ofS if 31a andS if 31c contain a J integral which is not differentiated. They thus include ultraviolet divergences Λ, which we will show cancel in physical setups in the wide-band limit. We sum the two contributions and focus on the terms that are linear in Λ to obtain the condition where K is a constant. If we assume that the environment is approximately particle-hole symmetric, the two ultraviolet cutoffs Λ λ ≈ Λ −λ are related. We can then further simplify the constraint (B16) to n |V inλ | 2 + |V niλ | 2 = |V λ | 2 , ∀ i, λ, where V λ are reservoir dependent constants. This is the case in any setup constructed purely from second quantised operators, such as the non-interacting level considered in this work or Anderson's impurity model. It is, however, possible to make the ultraviolet cutoff relevant. One method to do so is to introduce a large (of order unity compared to the bandwidth) particle-hole asymmetry, which breaks the assumption allowing us to go from Eq. (B16) to Eq. (B18). Alternatively, we can discard system states that are considered to be at very high energies. In the Kondo model, for example, high energy intermediate system states are not included, as their presence is purely encoded in an effective coupling. This leads to a bandwidth-dependent Kondo-temperature [20]. On the other hand, a Kondo temperature computed directly from the Anderson impurity, where the empty and doubly-occupied states are included does not depend on the bandwidth [79].

Current rates
Finally we find the current rates for theS 13 contributions. We use the diagrammatic rules from Fig. 21 and thus immediately conclude that we must replace the sums over the reservoir index that connects the two Keldysh contoursS 31 : The out-of-equilibrium( µ = 0) current across a level is an archetypal observable quantity. For a non-interacting system it can be obtained [20,74] using scattering matrices or the Green's function formalism. Alternatively it can be derived from the more general Meir-Wingreen formula [73], which is valid for both interacting and noninteracting systems and takes the Green's function as an input. For the non-interacting level, the current is given by the exact expression [78] see for example equation (9) of Ref. [73]. Physically, the transport occurs in an energy window of width µ, reduced by the temperature T . Within this window the tunnelling probability is proportional to the spectral function A of the non-interacting level, i.e., how much weight the noninteracting level makes available at a given energy. Using the integral Eq. (C4) in the limit Λ λ → ∞ and a partial fraction decomposition on A, we arrive at the exact expression (161) for the current across the non-interacting level.

Appendix E: Combining operators and superoperators
In Sec. VI we obtained expressions which contain both superoperators and the number operator N λ . To maintain an unambiguous order, such expressions require large numbers of brackets. For example the series AABBρ of operators A, B and superoperators A, B acting on a density matrix ρ can be bracketed in a multitude of ways including but not limited to (AA)(BB)ρ, or {A[A(BB)]}ρ, or A{A[B(Bρ)]}, (E1) which in general may produce different results. In this Appendix, we explain the notational shorthand we use to avoid this large number of brackets.
Whenever an operator N λ appears to the right of a projector P (27), we define where A is an arbitrary superoperator and we have used the fact that N λ is an environment operator and the cyclic property of the environment trace in P. Similarly, whenever an operator N λ appears to the left of a projector we define where B is an arbitrary superoperator and we used both the fact that N λ is an environment operator and the fact that [N λ , H 0 ] = 0 (which in turn implies [N λ , ρ 0 env ] = 0).