Density-operator evolution: Complete positivity and the Keldysh real-time expansion

We study the reduced time-evolution of open quantum systems by combining quantum-information and statistical field theory. Inspired by prior work [EPL 102, 60001 (2013) and Phys. Rev. Lett. 111, 050402 (2013)] we establish the explicit structure guaranteeing the complete positivity (CP) and trace-preservation (TP) of the real-time evolution expansion in terms of the microscopic system-environment coupling. This reveals a fundamental two-stage structure of the coupling expansion: Whereas the first stage defines the dissipative timescales of the system --before having integrated out the environment completely-- the second stage sums up elementary physical processes described by CP superoperators. This allows us to establish the nontrivial relation between the (Nakajima-Zwanzig) memory-kernel superoperator for the density operator and novel memory-kernel operators that generate the Kraus operators of an operator-sum. Importantly, this operational approach can be implemented in the existing Keldysh real-time technique and allows approximations for general time-nonlocal quantum master equations to be systematically compared and developed while keeping the CP and TP structure explicit. Our considerations build on the result that a Kraus operator for a physical measurement process on the environment can be obtained by 'cutting' a group of Keldysh real-time diagrams 'in half'. This naturally leads to Kraus operators lifted to the system plus environment which have a diagrammatic expansion in terms of time-nonlocal memory-kernel operators. These lifted Kraus operators obey coupled time-evolution equations which constitute an unraveling of the original Schr\"odinger equation for system plus environment. Whereas both equations lead to the same reduced dynamics, only the former explicitly encodes the operator-sum structure of the coupling expansion.


Introduction
The experimental progress in nanoscale and mesoscopic devices has continued to drive the development of theoretical methods to tackle models of nonequilibrium quantum systems that interact with their environment [1]. Open-system approaches based on the reduced density operator are particularly suitable for dealing with the strong local interaction effects that are often central to device operation, for example, in read-out circuits of qubits and many quantum transport devices. These interactions together with the continuum of environment energies often complicate the analysis when the system-environment coupling becomes strong.
In some limits, the evolution of the density operator is known to exhibit a special Markovian semi-group property. Gorini, Kossakowski, Sudarshan [2] and Lindblad [3] (GKSL) have established a form of the underlying quantum master equations that is both necessary and sufficient for the dynamics to be completely 1 positive (CP) and tracepreserving (TP), cf. also Ref. [5]. Thus, in these specific limits only quantum master equations of this form guarantee that the evolved reduced density operator still makes statistical sense which is one of the reasons why GKSL master equations have become popular and successful. More recently, such restrictions on the form of quantum master equations have been extended, see, e.g., Refs. [6][7][8][9].
However, the experimental importance of effects due to strong coupling and non-Markovianity-which are tied together [10]-has spurred progress in a variety of other approaches: Inclusion of parts of the environment into the system [11], time-convolutionless master equations [12,13], stochastic descriptions such as quantum trajectories [14], pathintegral [15], QMC [16], and hierarchical methods [17,18], multilayer multiconfiguration time-dependent Hartree method [19], perturbative expansions [20][21][22][23][24][25][26], resummation [27][28][29] and related techniques [30][31][32], projection techniques [33,34] and real-time renormalization-group methods [35][36][37][38][39][40][41][42]. These are applicable to more general reduced dynamics derived from unitary evolution U of initially uncorrelated states ρ(0) of the system (S) and ρ E of the environment (E): We use unitsh = k B = |e| = 1 and -unless stated otherwise-work in the interaction picture with respect to the system-environment coupling V (t) which generates the unitary evolution U (t). The dynamics (1) is always CP and TP [43] but need not have the GKSL form. Due to the above mentioned complications, it is unavoidable to make approximations at some stage in such approaches. A key challenge is to ensure that these do not uproot the CP and TP properties of the reduced dynamics while simultaneously maintaining a clear view on microscopic contributions to physical processes.

Complementary approaches to reduced dynamics
In this paper we make a further [44,45] step towards clarifying this challenging issue in a general setting by combining complementary approaches from quantum information theory and statistical physics. On the one hand, we are guided by the operator-sum form 1 Complete positivity ensures that even when the dynamics acts nontrivially only on a subsystem it should still yield a valid physical state for any initially entangled joint system. It is thus specific to quantum evolution [4] and does not simply require that the output is positive for every input but something stronger. of the general reduced density operator dynamics, ρ(t) = e K e (t)ρ(0)K e † (t), (2) developed by Sudarshan, Mathews, and Rau [46] and Kraus [47], cf. Ref. [48]. Dynamics of the type (1) can always be written in this form. The quadratic form manifests the CP property whereas TP is enforced by the sum rule e K e † (t)K e (t) = 1 S . The sum (2) averages the evolution over all outcomes e of possible measurements on the inaccessible environment. Although it is well-known how the Kraus operators K e (t) -acting only on the system-can be expressed in the joint unitary U and the initial environment state ρ E , the direct evaluation of these quantities becomes intractable when dealing with microscopic models with a continuous environment and strong local interactions.
Approaches to approximately solve such models have been developed in statistical physics where one instead considers general kinetic equations of the form The time-nonlocal nature of the self-energy or memory kernel superoperator Σ(t, t ) and the time convolution * in Eq. (3) underscore the non-Markovian nature of the general dynamics (1). The kinetic equation (3) is derived equivalently [12,24] either by Nakajima-Zwanzig projection superoperators or by the Keldysh real-time expansion 2 on which we will focus here. In either way what makes the practical evaluation feasible is the exploitation of Wick factorization of environment correlations in the cases of practical interest where the environment is noninteracting and initially in a mixed thermal state. Noting that Wick's theorem can be generalized (cf. Ref. [49] and App. A), we stress that the kinetic method is in general on equal footing with the operator-sum approach. It must thus be possible to express the dynamical map rho(t) = Π(t)ρ(0) which solves Eq. (3) in the form of an operator sum, Π(t) = e K e (t) • K e † (t), letting • indicate the position of the argument. How to actually do this on the level of the microscopic coupling V (t) has remained unclear until recently. In this paper we will find the explicit nontrivial relation between the two fundamental equations (2) and (3), in particular between Σ(t, t ) and K e (t), at every step linking the two complementary approaches. Following up work in this direction [44,45] we develop a transparent set of diagrammatic rules by which all of this can be achieved. This allows for a microscopic understanding of the CP-TP structure of the reduced evolution Π(t) and its memory kernel Σ(t, t ) in terms of quantities that are relevant in practical calculations for open systems with continuous environments. This is a crucial prerequisite for an improved understanding of approximations.

Key ideas and prior work
The key idea of our work is to purify the state of the environment as well as the unitary evolution map U to reveal its internal structure related to Wick's theorem. It is remarkable that the technical advantages of purification and entanglement have been long recognized in many-body approaches such as DMRG [50][51][52] and tensor networks [53][54][55] but seem to have been hardly explored for statistical field-theoretical approaches to the present problem. Thus, despite a great deal of intuition behind Keldysh time-evolution diagrams, their precise operational meaning in terms of physical measurements has remained unclear. Our approach fills in this gap: we identify which sums of Keldysh diagrams are physically meaningful by explicitly relating them to a quantum circuit describing evolution conditioned on partial measurements of the open system's environment. This leads to a diagrammatic expansion of (the operators underlying) the Kraus operators in terms of the microscopic system-environment coupling which allows one to exploit field-theoretical methods in the setting of quantum information theory. Although our analysis focuses on the reduced density-operator method, it relates to similar developments in non-equilibrium Green's functions [56] and quantum-field theory [57] (Cutkosky 'cutting rules').
As mentioned, our work was stimulated in particular by two prior publications [44,45] and aims to complement these, cf. also Refs. [58][59][60][61]. Our formalism extends and generalizes the scope of the important Refs. [44,60,61] whose formulation is most easily connected to the standard density-operator Keldysh formalism [27]. In particular, we allow for multilinear coupling V to both fermionic and bosonic environments which enables switching back and forth between the Kraus description (2) and the quantum kinetic description (3). Ref. [45], on the other hand, ingeniously circumvents positivity problems by using a projection technique to derive the evolution of ρ(t) and squaring the result to obtain a positive density operator. While both approaches at first seem quite distinct, they turn out to be closely related. This is nontrivial to see directly, but becomes clear once they are tied to the basic operational formulation of quantum information theory. Whereas the present paper achieves this for the Keldysh diagrammatic real-time approach related to Refs. [44,60,61], the projection-superoperator formulation of Ref. [45] will be addressed in subsequent work [62]. This underscores the extended scope 3 of the ideas developed here. Throughout the paper we will point out the extensions and advantages of our our approach relative to these works.

Motivation: Approximations
TP approximations. Our work is motivated by the need to perform approximations for the mentioned systems of interest. As it turns out, most approximations formulated in terms of the memory kernel Σ(t, t ) of the kinetic equation (3) have no problem with guaranteeing the trace-preservation but impose no restrictions at all to ensure complete positivity. For such TP approximations on cannot be sure that complete positivity is not lost, implying that negative density operators may 4 come out depending on the input. This may be fatal for a calculation: actual loss of positivity of the state implies that the results cannot even be taken as a 'rough indication' of what is going on since negative probabilities are meaningless. Such loss of positivity is sometimes related to the quality of the approximation indicating that 'not enough processes were taken into account'. However, we will make explicit that this failure instead indicates that what one is actually summing up in such approximations are not physical processes but rather partial mathematical contributions to such processes.
These issues tie in with the current interest from the side of statistical physics and quantum thermodynamics to formulate approximation schemes which guarantee meaningful results for entropic quantities. These are completely undefined unless the quantum 3 The projection approach [4,34] based on the quantum superchannel formalism is more general by accounting for an initially correlated system-environment state. It is an interesting open question how this can be related to the present approach and that of Refs. [44,45,60,61]. 4 Note carefully the difference between CP and positivity-preserving (PP) superoperators. Approximate evolutions that are not CP but only PP still produce positive outputs for any input. Examples of such superoperators are well known in entanglement theory [63] but no simple general characterization is known. Approximate evolutions that are neither CP nor PP may still produce positive outputs but not for all inputs. Also, a positive stationary state does not guarantee there was no breakdown of CP in an approximation: one can construct families of TP superoperators that are physically meaningless but nevertheless pass this test. state ρ(t) is positive for all times. In this context the Keldysh approach has proven useful by its extension to 'parallel-worlds' [26] to compute Renyi entropies. Our operator-sum formulation of the Keldysh approach identifies a variety of approximations that give meaningful entropic quantities, including also the exchange entropy [64,65] associated with the environment [66] which requires the evolution Π(t) to be completely positive.
CP approximations. Indeed, we will find it is not difficult characterize large families of approximations which strictly guarantee complete positivity independent of the details -and thus of the quality-of the approximation. Although they are harder to evaluate and have not been explored much, they also have the advantage of a bounded quantitative measure for the TP error related to probabilities of physical processes 5 unaccounted for. Importantly, the characterization of CP approximations also provides a guide for less ambitious strategies motivated by practical concerns. It indicates how in TP-approximations one can systematically 'improve towards complete positivity' by identifying which diagrams should be added to give a better partial account of a physical process without insisting on strict CP. It must be stressed that failure of CP approximations to guarantee TP can be fatal, too: Although perhaps less important for short-time dynamics, TP ensures that dissipative dynamics achieves a stationary state at large times.
Such approximations guaranteeing complete positivity connect to the interests of quantum information theory where this property is essential to maintain the connection to several fundamental theorems. For example, evaluating Eq. (1) for the evolution Π(t) using some non-CP approximation breaks the chain of reasoning because fundamental insights in this field rely on this property, such as the quantum data-processing inequality [64]. Our approach is centered around the time-evolution of (operators underlying the) Kraus operators and provides a systematic way for going beyond the limited GKSL approach without giving up CP. Importantly, we show how this can be done by exploiting known statistical physics techniques.
CP-TP duality. Thus, neither route is without challenges. In fact, the two approximation strategies and their problems are fundamentally tied together. It has been noticed [67] that there is a kind of duality between complete positivity and trace-preservation: in practice, rigorously fixing CP in some approximation tends to uproot TP and vice versa. Although this is not entirely unexpected based on the correspondence between quantum states and evolutions (App. C), we will identify the microscopic origin of this duality in terms of elementary diagrammatic contributions to the real-time expansion. We will show that in some approximation schemes based on selection of diagrams CP and TP cannot be achieved simultaneously. Approximations that are both CP and TP require more sophisticated schemes explored in Refs. [44,45,[58][59][60][61] which we will also touch upon.
In the present paper we focus on the prerequisites for a clearer understanding of all the above mentioned approximation schemes. This is of common interest to both research fields: although our encompassing formalism remains geared towards applications -by its formulation in terms of microscopic couplings and Keldysh diagrams of Eq. (3)-it is firmly rooted in the operator-sum (2) and defined operationally in terms of quantum circuits.

Outline and guide
The paper is organized as follows. Sec. 2 first applies methods of quantum-information to the problem. We review the operator-sum representation, purify the environment state -expressing finite temperature in terms of entanglement-and also purify the evolution map based on its internal structure related to Wick normal-ordering. These physically motivated steps reveal clearly how measurements condition the time-evolution and connect the CP structure with Wick's theorem.
In Sec. 3 we turn to statistical field theory and connect the obtained operator-sum representation with Keldysh real-time diagrams. This allows us to identify groups of diagrams associated with a physical process conditioned on a measurement outcome. Each such process also corresponds to a Keldysh operator acting on both system and environment which is obtained by 'cutting' the group of Keldysh diagrams 'in halves' and summing them. These are the central objects of interest, not the corresponding Kraus operators acting only on the system which can be obtained by normal-ordering and taking environment matrix elements.
Sec. 4 builds on this partial 'undoing' of the environment trace which is a convenient way to see the connection between microscopic Keldysh diagrams and Kraus operators. Focusing on the generators of infinitesimal evolution in either approach, we derive two equivalent hierarchies of time-nonlocal evolution equations. The hierarchy for the Keldysh operators presents a useful exact unraveling of the original Schrödinger equation for system plus environment. It explicitly encodes the operator-sum structure of the reduced dynamics before tracing out the environment. The result of this analysis is the explicit functional dependence of the memory kernel superoperator Σ = m Σ m [σ 0 , . . . , σ m ] in Eq. (3) on the self-energy operators σ 0 , . . . , σ m generating the Keldysh operators in Eq. (2). It expresses the operator-sum theorem for general infinitesimal reduced evolutions and reveals at the deepest level that CP and TP appear as incompatible constraints because of the difference between time-ordering of the microscopic couplings on one and on two branches of the Keldysh real-time contour.
Finally, in Sec. 5 we illustrate the generally applicable formalism for the simplest possible model of a single fermionic mode interacting with a hot continuum of fermion modes. This Markovian example only serves to highlight the inner workings of the approach and physical principles involved, delegating the technically much more demanding applications to future work.

Quantum-information approach to reduced dynamics
The two approaches (2) and (3) that we combine both address the calculation of the same dynamical map (1). We are ultimately interested in the challenging case where the environment consists of reservoirs with a continuum of modes in mixed thermal states and where particles strongly interact on the system. Microscopic models of this type lead to interesting quantum many-body and transport phenomena which preclude direct evaluation of Eq. (1) and require advanced methods to obtain good approximations to the density operator. In quantum information one often deals with simpler models where such a direct approach is feasible. We now show that the techniques developed there can still guide calculations that address the more complex setups mentioned above. Fig. 1 shows how the evolution (1) can be depicted as a quantum circuit. Taking the partial trace amounts to trashing all information about the environment, i.e., averaging over all possible measurement outcomes. The coupling V (t) between system and environment generating the unitary evolution U (t) = T exp − i t 0 dτ V (τ ) correlates the system with the environment over time. In the following, the central idea is to explicitly keep track of these correlations by extending the environment and the evolution map. A computational approach which makes this explicit at each step will manifestly guarantee the complete positivity of the evolution as we work out in detail for the the Keldysh diagrammatic approach in Sec. 3.

Purification of the mixed environment state
To achieve this, we first need to distinguish correlations between the system S and the environment E built up over time from the correlations already present in the initial, mixed environment state written in its eigenbasis as ρ E = n ρ E n |n n|. The standard way to achieve this in quantum information theory is by purifying the environment state, i.e., expressing it as denotes the purified environment state. This is a pure, entangled state in a doubled 6 environment space EE written in terms of the tensor product of the eigenbasis {|n } of ρ E . We already note that this purification is important for the analysis but will not need to be computed in the end. We can bring Eq. (1) into the equivalent form where we note that the ancilla E is not affected by unitary evolution. Due to unitarity of U the new map W EE : S → S ⊗ E ⊗ E is an isometry, W † EE W EE = 1 S , which together with the form (5) is necessary and sufficient for Π to be a CP-TP superoperator [43,68]. This ensures it corresponds to some quantum circuit built using only unitary components, pure states and projective measurements. In Fig. 2(a) we show this quantum circuit for the form (5) in which W EE is indicated by the shaded part.
Any change of the purified environment state |0 is now certain to be caused by the microscopic coupling V . This becomes explicit when we expand the trace in Eq. (5) in terms of some orthonormal basis {|e } for EE giving the operator-sum representation (2): When a measurement on this extended environment does not find it in its known original state, |e = |0 , this directly probes the entanglement generated between the system and the original environment E. Accordingly, the evolution conditioned on this outcome is K e (t)ρ(0)K e † (t) as shown in Fig. 2(b). We stress that for the original environment E, due to its mixed nature (e.g. thermal noise), it is not possible to keep track of changes incurred by the evolution only from measurements at the final time t. This train of thoughts is of course well known from quantum information theory. However, tangible means to obtain the Kraus operators (6) from the full unitary U ⊗ 1 E acting on a continuous environment in practice require field-theoretical methods. There, it is easier to evaluate Wick-averages of the unitary U after it has been expanded in the microscopic coupling V , and, furthermore, to consider infinitesimal evolutions using kinetic equations. In the following, we will explicitly provide the connection to this fieldtheoretical approach allowing one to preserve complete positivity in a framework better suited for approximate treatments.

Purification of the evolution operator
To do so, the picture of the measurement-conditioned evolution needs to be further refined: the above introduced states |e on the purified environment do not account for the internal structure of the evolution U relative to the environment state ρ E . This is revealed when expanding the unitary into normal-ordered components as discussed in App. B: :k m (t): .
By construction, the operator k 0 is exceptional in that it acts trivially on the environment. The normal-ordering of the remaining terms m = 0, indicated by : : , is a standard procedure reviewed in App. A. It is generally applicable, i.e., it does not only apply to noninteracting environments that we will consider later on. In view of this later application we will refer to m as the modes of the environment. We already note that the normal-ordering : : need not be evaluated in the end but continues to play a simplifying role: it ensures by construction that tracing out the environment in gives a nonzero system operator only if the modes m and m match. We can exploit the structure (7) from the start by explicitly keeping track of the labels of the environment modes m using a superoperator W EE M : S → S ⊗ E ⊗ E ⊗ M which stores these in a register |m in an auxiliary meter space M. This leaves the total evolution unaltered if we additionally trace out the meter: Due to the normal-ordering property (8) the new map is also an isometry, W † EE M W EE M = 1 S , and can thus be regarded as another purification 7 of the original evolution map (1). 7 The purification of the isometry W EE • (W EE ) † = TrM W EE M • (W EE M ) † corresponds to altering  Expanding both the trace over M and EE we obtain a refined operator-sum where by measuring the meter's register we can pick out the normal-ordered component :k m : of the evolution. Here, the Keldysh 8 operators k m are still operators on system and environment, in contrast to the Kraus operators K e m (t) which act on the system only. Keeping the environment E via the k m is at first a technical trick to avoid writing out matrix elements. Later on it will allow for an interesting unraveling of the Schrödinger equation of system plus environment [Sec. 4.2]. By construction, for m = 0 the Keldysh operator factorizes as k 0 = K 0 0 ⊗ 1 E and the only remaining Kraus operator is the environment average of the unitary U . The circuit in Fig. 3 summarizes the above: purifying the environment to an entangled state |0 (to eliminate thermal noise) and purifying the unitary U to an isometry W MEE (using normal-ordering) results in the refined operator-sum (10) for Π. The trace over the purified environment and the meter correspond to discarding the outcomes m (mode labels) and e (state of purified environment) by summing over them. Importantly, we can now choose the level of refinement: if one only resolves the trace over the meter space M one may write the original time-evolution problem as This is the key relation of the paper: in Sec. 3 we show that this exactly corresponds to calculations employing Keldysh real-time diagrams and Wick's theorem. There the Keldysh operators k m (acting on the system and environment) are uniquely represented by groups of Keldysh diagrams for a conditional propagator (acting on the system only), depending on the measurement outcome m. This corresponds precisely to the modified quantum circuit as shown in Fig. 4 and explicitly ensures that Π m is a CP superoperator. It is these objects that one should calculate while formally keeping track of m in order to obtain a CP result for Π = m Π m . We remark that the above purification approach can be completely avoided by directly working in Hilbert-Schmidt or Liouville-space [62] which is useful for, e.g., deriving projection-based approaches that guarantee CP evolution [45]. The above quantuminformation approach, however, emphasizes the clear operational test: if a complicated mathematical expression can be expressed as a physically implementable quantum circuit without initial system-environment correlations, then it must be a CP superoperator.

Trace preservation (I): State-evolution correspondence
So far we have focused on the CP property which is almost trivialized in the quantuminformation approach: when Π is written in the operator-sum form it can be verified term-by-term. However, for the TP property of the evolution this implies the opposite: it requires a nontrivial sum of quadratic operator expressions to match the system identity operator at all times t. For the Keldysh operators this condition translates to This implies that the conditional propagators Π m are trace-nonincreasing (TNI) superoperators, Tr Π m ≤ Tr. Keeping in mind applications with continuous environments these are highly nontrivial 9 conditions to check or guarantee. In this respect, the relation of k m to the Keldysh real-time expansion announced at Eq. (12) is intriguing: in the latter approach TP is a trivial condition to satisfy whereas CP is nontrivial to guarantee. There is thus a fundamental incompatibility in the practical requirements for ensuring CP and TP, no matter in which form one writes the propagators Π and Π m . This CP-TP duality was noted in Ref. [67] and is in fact an expression of the fundamental correspondence between states and evolutions [69] based on the isomorphism of de Pillis [70], Jamiolkowski [71] and Choi [72]. This general point of view is expanded in App. C but we will instead focus on the way this CP-TP duality is expressed on the microscopic level of diagrammatic contributions that one encounters in practice when we return to this issue in Sec. 3.4.

Statistical field-theory approach to reduced dynamics
The considerations of the previous section hold generally -noting the key assumption of an initially uncorrelated system and environment. In the following, we will focus on the Keldysh real-time diagrammatic technique [20] which is a powerful tool for calculating density-operator evolutions beyond the weak-coupling limit [22][23][24]27] when the environment is a composite of noninteracting, continuous reservoirs of fermions or bosons, each in a different thermal state. We note, however, that most of our conclusion can be generalized [62]. Our aim is to reorganize this well-established diagrammatic expansion such that one can easily identify the Keldysh-and Kraus operators in terms of practical Wickaverages. This makes our approach of immediate relevance for various formulations of approximations [27][28][29]33] in terms of such diagrams.

Standard Keldysh real-time expansion
The Keldysh real-time diagrammatic approach is based on the formal expansion in the coupling V of the unitary U (t) = T exp − i Wick's theorem expresses the environment trace as the sum of products of paircontractions of environment field operators appearing in the different V 's. These are drawn as connecting lines in Fig. 5 for the example of a bilinear coupling. We assume that the operator V has been normal-ordered with respect to the initial environment state, cf. App. A. Thus, any partial contraction of fields appearing in the same vertex V gives zero 10 . Any nonzero environment average of the coupling is thereby already absorbed into the definition of the system Hamiltonian H(t) + V E (t) → H(t). Rules for translating diagrams to expressions are summarized in App. D but the following considerations do not require these details, the diagrammatic representations suffice. More model-specific diagrammatic rules [22-24, 27, 73] and their equivalent formulations [25,35,38,42] make this a very efficient technique for higher-order calculations [20, 22-24, 26, 28, 29, 73].
From the resulting diagrammatic series one infers the self-consistent Dyson equation, using standard field-theoretical considerations: the self-energy kernel Σ is defined by the sum of two-branch-irreducible diagrams [27]. For this distinction one needs to draw Keldysh diagrams in the standard form shown in Fig. 5(b): such a diagram is (ir)reducible when it can(not) be split up by a vertical cut without hitting a contraction line. Taking the time-derivative of the Dyson equation (14) gives the time-nonlocal kinetic equation (3) noting that π = I in the interaction picture. A crucial problem of the Keldysh diagrammatic expansion in the above standard form is that it does not exhibit explicitly the structure of an operator-sum (10) which it must have by the CP property of the original problem (1). This is especially perilous when considering that approximations are often made by omitting and resumming diagrams: by doing so one may neglect contributions that are essential for the operator-sum structure to . We note that it is common to denote the adjoint U (t) † = U (−t) the 'backward' evolution, imagining to traverse the contour first forward along the upper branch, then across the dashed connector 'backward' on the lower branch. Here, this is misleading: time always runs in the forward direction, also on the lower contour as indicated in (b).
guarantee CP evolution. Moreover, such approximations are frequently made on the level of the self-energy Σ which inherits an even more complicated structure [Sec. 4.3] from the operator-sum form of the propagator Π considered so far.

Reorganized Keldysh real-time expansion:
Cutting and pasting rules We therefore reorganize the diagrammatic series guided by the quantum-information considerations of Sec. 2 and identify which groups of diagrams rigorously correspond to a physical measurement processes. We write the normal-ordering expansion (7) as which is constructed in App. B in detail. Here, k mq...m 1 is the evolution U partially contracted such that only the modes labeled by m = m q . . . m 1 remain. We ignore the ordering of the modes, i.e., k mq...m 1 is invariant under mode permutations and the summation m q . . . m 1 is restricted so as not to double count terms. The special case k 0 denotes the part of U in which no environment operators are left, i.e., the environment average Our central formula (12) explicitly ensures that the CP-TP evolution of interest, which are CP-TNI because they are conditional evolutions depending on measurement outcomes m q . . . m 1 indicating which environment modes the system has interacted with. In Fig. 6 we explain how each Π m 1 ...mq corresponds to a precise group of Keldysh diagrams. In (a)-(b) we first sketch how the standard Keldysh expansion is obtained using the now more pertinent left-right form of diagrams introduced in Fig Thus, the elementary physical processes contained in the standard real-time expansion for the propagator Π when computed in terms of Wick contractions are the conditional propagators Π mq...m 1 . However, once the mode-sums for these contractions have been performed the CP property is hidden and approximations may inadvertently break up the CP structure. We now discuss the two computational stages that the above implies.
Stage 1, Eq. (17b) First, the conditional propagator Π m needs to be computed by eliminating E. This can be done in two equivalent ways: (A) Compute Π mq...m 1 directly. An advantage of our approach as compared to Refs. [44,45,60,61] is that it does not necessarily require a fundamentally new formulation: the standard rules [see App. D] of the existing Keldysh approach merely have to be complemented with a diagrammatic CP-rule which instructs one to keep track of m q . . . m 1 : Rule for the conditional propagator Π mq...m 1 : Sum all diagrams that have a fixed number q of contractions between the different branches of the Keldysh contour with a fixed set of mode indices m q . . . m 1 assigned to these lines in every possible order. This includes summing all possible intra-branch contractions separately and summing over the mode indices of these contractions. Perform all two-branch ordered time integrations.
Note that this rule requires both reducible and irreducible diagrams in the two-branch sense to be summed. In Fig. 7(a) we illustrate why this is technically important and make explicit the intricate difficulty of relating time-ordering on one and two branches of the Keldysh contour. To obtain Π m (t) we sum up terms on the right which are single functions of the final time t obtained by two-branch time-ordered integrations. If the CP-rule given a selection of diagrams on the right hand side, one can verify it has an operator-sum form by first cutting horizontally and summing up the halves. Pasting: subsequently pasting the halves back together should recover the selection. Leaving out, for example, the center two diagrams on the right hand side, this test of 'completing the square' will fail. When directly selecting half -diagrams on the left hand side, every possible approximation (select only green, only blue or select both) by construction produces one of the three groups indicated by dashed lines and shading on the right, each of which is a CP superoperator.
is obeyed, this sum is able to produce a single term that factorizes 11 into two separate functions of time, . Here, the time-integrations are performed independently inside k m on each branch separately. Note that this 'completing of the square' of ordered time-integrations is not possible when summing two-branch irreducible diagrams only, and the latter can therefore never correspond to physical processes [Eq. (26)].
(B) Compute Π mq...m 1 from k mq...m 1 . Alternatively, one can evaluate the right hand side of Eq. (17b). The required Keldysh operators can be obtained from the value of a Keldysh diagram that is cut in half, i.e., from a single time-branch.
Cutting rule for k mq...m 1 : Sum all half-diagrams that have q of external fields acting on a fixed set of environment modes m q . . . m 1 over which we symmetrize. This includes summing over all possible intra-branch contractions and their mode indices. Perform all one-branch time-ordered integrations, including the times of the external vertices.
The full diagrammatic rules for calculating the Keldysh are summarized in App. E and are conveniently similar to the standard Keldysh rules [App. D]. The required value for Π mq...m 1 is then obtained by pasting two halves together as follows: Pasting rule for Π mq...m 1 : Fully contract the external vertices appearing on both sides of in all possible ways. We note that the normal-ordering in Eq. (17b) is only an instruction to simplify the evaluation by discarding all internal Wick contractions of k m and k † m . This is in contrast to the Kraus operators 12 K e m = e|(:k m : ⊗1)|0 for which the normal-ordering has to be explicitly evaluated. It should also be noted that when using Wick contractions, the most elementary physical process Π m is never described by a single or even a finite number of Keldysh real-time diagrams. Typically, guaranteeing CP in approximate dynamics calls for methods that are nonperturbative in the system-environment coupling.
Stage 2, Eq. (17a) It remains to sum the computed Π m over the modes 13 m q . . . m 1 and their number q = 0, 1, . . . in order to eliminate the meter M. Our example in Sec. 5 will illustrate that truncating q, the number of modes, in stage 2 correlates with the time to reach the stationary state relative to the time-scale that is set by the time-dependence of Π mq...m 1 (t) obtained in stage 1 [Eq. (50)]. Because the expansion index m q . . . m 1 is physically motivated, the truncation of stage 2 never jeopardizes CP.

CP approximations
The above transparent rules translate the operator-sum theorem to the framework of the standard real-time expansions on one and two branches of the Keldysh time-contour. Above all, this provides explicit means of translating problems and solutions between the very different formulations employed in quantum information theory and statistical physics on the microscopic level of the system-environment coupling V .
As expected, evaluating the objects obtained by these simple rules is by no means easy. The importance of our general framework is that at least the large variety of existing approximations can be formulated, compared and systematically improved. In Sec. 4 we make a further step towards practical applications by deriving useful self-consistent and infinitesimal evolution equations. The reorganized expansion [Eq. (17), Fig. 6] already shows how CP approximations can be constructed generically. Fig. 7(b) illustrates this for approximations based on selecting complete Keldysh diagrams 14 . In order to guarantee CP, such a selection must pass the test of horizontally cutting the Keldysh diagrams using their standard representation [ Fig. 5(b)] and summing up their halves to obtain k m [left hand side in Fig. 7(b)]. By the pasting rule, one can check that k m reproduces the groups of two-branch diagrams originally selected [right hand side in Fig. 7(b)]. This amounts to checking that one can diagrammatically 'complete the square' to the operator sum (17b).
Clearly, directly selecting groups of two-branch diagrams, each group corresponding to one conditional propagator Π m , passes this test. This amounts to selecting physical processes, i.e., possible measurements performed on the meter M and therefore preserves CP. This is illustrated by the example later on [Eq. (52)]. The approximate propagator is a TNI 15 superoperator [Eq. (13)] expressing that some probability was lost by omitting physical processes.
A more advanced scheme is to select complete one-branch diagrams for the Keldysh operators k m directly [left hand side in Fig. 7(b)]. One can then infer which distinct twobranch diagrams they generate and equivalently select only these in the standard Keldysh expansion. This amounts to a modification of the propagator Π m by selecting a subset of diagrams from that group. Because this is done without breaking the CP operatorsum structure the modified contribution still corresponds to a physical process occurring with a positive probability, cf. Sec. 3.5, provided the superoperators are TNI which is not automatic 15 .  (18). (b) Reorganized expansion: brackets indicate groups of diagrams that belong to the same physical process generated by one Keldysh operator k m . Applying the simple rule of (a) reveals a nontrivial network of cancellations between diagrams in different brackets corresponding to different Keldysh operators. This shows that CP and TP cannot be ensured simultaneously by selection of complete Π m -groups of diagrams.

Trace preservation (II): Microscopic state-evolution correspondence
We now address the constraints imposed by TP [Sec. 2.3] and discuss the apparent duality with the constraints imposed by CP on the microscopic level of the coupling V .
In Fig. 8(a) we first show how in the standard Keldysh real-time diagrammatic approach TP is ensured by pairs of diagrams with the latest (leftmost) coupling vertex V appearing on opposite branches of the Keldysh contour. When taking the system trace, the cyclicity of the total trace allows us to move these latest vertices to either side, where • denotes the initial state as well as all earlier vertices from the n-th and n -th order terms of U and U † , respectively. The resulting terms are already contained in different orders (n±1, n ∓1) of the double expansion in powers of V up to an opposite sign. Thus, all terms cancel pairwise except the zeroth order n = n = 0 term [first diagram in Fig. 8(a)] which guarantees that Tr Π = Tr. Notably, since the total order of V is n + n for the canceling terms, TP can be guaranteed easily order-by-order in the coupling-expansion.
On the other hand, CP can neither be inferred by inspecting just a few terms in the standard Keldysh expansion nor be guaranteed order-by-order in V .
In the reorganized Keldysh expansion the difficulty of guaranteeing CP and TP is reversed. Yet, we can make use of the diagrammatic rule (18) to reveal the nontrivial interrelation of the Keldysh (and Kraus) operators that is implied by the nonlinear constraint (13b) involving all Keldysh operators: In Fig. 8(b) we sketch a group of two-branch diagrams for some process Π mq...m 1 (indicated by the middle bracket). Consider the two possible situations for the latest (leftmost) vertex indicated in red. This vertex is either already contracted within the same branch (green); then, the TP rule (18) demands that we also take into account a corresponding diagram whose latest vertex lies on the opposite branch and thus contributes to a process Π m q+1 ...m 1 with one additional external vertex (right bracket). Otherwise, the vertex is contracted with the opposite branch (blue); then, the TP rule (18) requires a term in the process Π m q−1 ...m 1 with one less external vertex (left bracket). This implies that the TP condition for the Keldysh operators making up each process is guaranteed by the pairwise contribution of coupling vertices to 'adjacent' operators, schematically Thus, the TP condition on the Keldysh and Kraus operators that appears nontrivial [Eq. (13)] is microscopically clear. These pair cancellations do not form a simple chain but rather a network correlating 16 all terms in the operator-sum Π = m Tr E :k m :(•⊗ρ E ) :k † m :. Importantly, the CP structure requires one to include all possible two-branch timeorderings of the individual vertices [cf. Fig. 7(a)]. In particular, we must include both terms in the middle bracket in Fig. 8(b) which then rigidly connects two independent TP constraints (18). Because this applies to every pair of diagrams in the middle bracket, all groups of diagrams connected to it by a TP constraint must be included as well.
Continuing this argument to these other groups one finds that eventually all diagrams have to be included. This is the ultimate microscopic ground for the CP-TP duality as it appears in practical approximations, complementing the general considerations of App. C. We stress that this does not rule out the possibility of CP-TP approximations by diagram selection, see below.

CP vs. TP approximations
Having understood the conflicting demands of CP and TP microscopically, consider again the first CP approximation scheme discussed in Sec. 3.3. Let the approximated propagator Π qmax = qmax q=0 mq...m 1 Π mq...m 1 be obtained by discarding terms with more than q max inter-branch contractions. The TP error for such a truncation is quantified by the trace distance of the operator-constraint (13b), computed from the adjoint of the approximated evolution Π † qmax . This error will be monotonically reduced by increasing q max since the added terms are CP-TNI operators. Physically, this means we account for more modes that have interacted with the system up to the considered time t. One expects that more modes are required when the system interacts longer with the environment. We stress that this truncation is not a naive short-time expansion 17 of Π(t): the time-dependence of Π m (t), obtained in stage 1, defines a cutoff on the number of modes to be kept in the stage 2 expansion depending on the evolution time t. This will be illustrated for a simple example in Sec. 5.
This should be compared with the situation encountered in the frequently used TP approximations. Under the TP constraint, discarding diagrams without further precautions will uproot the operator-sum structure of the expansion, and thus, violate CP. If one has computed the evolution map Π(t) (not the evolution of a specific initial state or only the stationary state) one can check this by calculating the associated Choi operator which is obtained by acting with the superoperator Π on the maximally entangled state |1 = k |k ⊗ |k on the tensor product of the system-space with itself, see App. C. Violation of CP is then detected 18 by finding a negative eigenvalue 19 of the operator choi(Π), but this is hardly ever discussed for TP approximations.
Also, for TP approximations there seems to be no way of systematically 'improving towards CP' -see however Eq. (36) below-other than numerically approximating the exact result better. It is now clear that this breakdown of CP in TP approximations is not due to 'missing physics' but instead indicates that one has uprooted the statistical structure: one is not summing up physical processes but mathematical contributions.
The above indicates that approximations which simultaneously guarantee TP and CP require schemes more sophisticated than the diagram selections considered so far. The difficulty of ensuring TP in CP-approximations was noted in Refs. [45,67]. Building on Refs. [44,58,59] very recent work constructed systematic CP-TP approximations for a broad class of systems [60,61]. There it is noted that the CP constraint is trivially satisfied, but the TP constraint requires an ingenious truncation of environment correlation functions and the detailed proof of TP is nontrivial. It can be shown that this truncation on the level of the Kraus operators used in Refs. [60,61] corresponds to a selection of Keldysh real-time diagrams for Π which allow TP to be verified more easily using the rule (18), illustrating its usefulness. The novel feature of this diagram selection is that it entails a modification of the pasting rule (17b) which we kept fixed in our above considerations, see Ref. [62] for details.
Finally, besides avoiding the dual fatalities of CP-and TP-loss mentioned in Sec. 1.3, approximated dynamics which is both CP and TP has another important feature known as a Stinespring dilation [43]: it can in principle always be expressed in terms of some effective joint unitary evolution U where some effective, initially uncorrelated environment ρ E is traced out as in Eq. (1). Direct construction of the quantities U and ρ E from U and ρ E , respectively, might be an interesting alternative route towards CP-TP approximations.

Hierarchies of self-consistent and kinetic equations
With so much room for approximations, it is of practical interest to have equivalent selfconsistent and infinitesimal formulations of the measurement-conditioned evolution.
In Sec. 4.1 we show how to compute Π m directly from memory kernels Σ m . This approach is closely related to well-developed techniques in statistical physics and is therefore discussed first. It avoids all use of the Keldysh operators k m (and Kraus operators) and instead implements stage 1 through the task of computing the self-energy super operators Σ m .
In Sec. 4.2 we address the calculation of the Keldysh operators k m and by analogy introduce their memory kernels σ m . This is possible since the k m are closely related to the Π m by the cutting rule. This instead implements stage 1 through the task of computing the self-energies σ m -operators on system and environment-and determines the Kraus operators of interest in quantum information theory.
The two approaches are graphically summarized in Fig. 9 and are equivalent, being related by Π m = Tr E :k m :(• ⊗ ρ E ) :k † m :. In both approaches it remains to perform stage 2, the sum over modes m, collecting all CP contributions Π = m Π m . The equations also share the same forward-feeding hierarchical structure which makes them of practical importance in view of the success of numerical methods [74] based on similar hierarchical 18 Checking the evolved state ρ(t) for all possible inputs ρ(0) only establishes positivity-preservation of the evolution but not the stronger condition of complete positivity, see also footnotes 1 and 4. 19 Physically, this implies that the entropy production in the purified environment [65] becomes ill defined. formulations. However, the approaches are dual in that the relative difficulty of guaranteeing CP and TP are opposite. These two hierarchies allows us in Sec. 4.3 to derive an infinitesimal form of the operator-sum theorem by explicitly determining the functional relation between the families of memory kernels Σ m and σ m .

Hierarchy for conditional evolutions Π m
We start from the unraveling (17a) of the reduced evolution into conditional evolutions Π m , each of which is a CP-TNI superoperator. Since Π m is the sum of all diagrams with q fixed contractions involving modes m = m q . . . m 1 between opposite branches, one introduces corresponding conditional self-energy superoperators Σ m by restricting all contractions to be two-contour irreducible [ Fig. 5(b)]. The full selfenergy Σ in the standard Dyson equation (14) -the sum over all two-branch irreducible diagrams-is obtained by summing over all modes m: In the following we use simplified notation for the labels of the environment modes: As sketched in Fig. 9, the expansion (24) unravels the Dyson equation (14) into a hierarchy of equations (23) that is physically motivated: each tier is labeled by the measurement outcome m, listing the environment modes that interacted with the system in the process Π m . In the lowest tier of the hierarchy, q = 0, no modes interacted with system. We obtain an independent self-consistent equation where π denotes the free system propagator which reduces to π = I in the interaction picture. Importantly, this the only self-consistent step in the hierarchy. Its solution Π 0 = π + π ∞ l=1 ( * Σ 0 * π) l amounts to a renormalization of the free system propagator π, giving a CP-TNI superoperator, see Eq. (31) below 20 . Although Σ 0 involves all orders of the microscopic coupling V , it does not contain all diagrams of the standard Σ in Eq. (14).
In the hierarchy (26) it is implicitly understood that one symmetrizes 21 the mode indices appearing on the right hand side. Also note that on the right hand side individual contributions -or partial sums of them-do not correspond to a physical process: only their full sum makes up Π m which is a CP-TNI superoperator due to its operator-sum form (17b). Moreover, these terms consist of sequences of irreducible blocks Σ m and Π 0 and are thus reducible. This goes against the widespread intuition that contributions to irreducible kernels such as Σ (and its components Σ m ) correspond to physical processes.
To further highlight how the physical processes Π m influence each other in time, we cast the self-consistent hierarchy (26) into an equivalent hierarchy of time-nonlocal kinetic equations for infinitesimal conditional evolutions: in the interaction picture with the initial conditions Π 0 (0) = I and Π m (0) = 0 for m = 0. Each column of the matrix of self-energy superoperators sums to the standard self-energy Σ = m Σ m , and by summing the column-entries of Eq. (27) to obtain Π = m Π m , we reproduce the standard kinetic equation (3), d dt Π = Σ * Π, in the interaction picture. Thus, the hierarchy presents an unraveling of the time-nonlocal kinetic equation (3). The lower-triangular form of the matrix ensures that the coupled kinetic equations feed forward and can be iteratively solved for Π 0 , Π 1 , . . ., the formal solution being given by Eq. (26).
An advantage of the hierarchies over the standard (kinetic) equation Eq. (14) [Eq. (3)] is that both TP and CP have become explicit, even though the duality of their constraints persists [Sec. 2.3, 3.4]. Although the individual conditional self-energies do not preserve trace, i.e., Tr Σ m = 0 because the Π m are TNI, we can still see that Π is TP: summing the elements in each column of the matrix appearing in Eq. (27) and taking the trace gives Tr Σ = m Tr Σ m = 0 because contributions from subsequent Σ m cancel pairwise by our earlier arguments 22 .
As we know, CP is guaranteed for each Π m separately by Eq. (17a). In both hierarchies this is expressed by the forward-feeding structure which ensures that neglecting complete higher tiers (corresponding to larger number of modes m q . . . m 1 ) does not affect the lower tiers: if the latter are solved using the required finite set of exact Σ m , the solutions are guaranteed to be CP. Thus, although there are no Kraus operators to be seen in Eq. (26)- (27), we can fully exploit the implications of the operator-sum theorem by unraveling the standard approach of statistical physics based on Eq. (3). At the same time, we can take advantage of the diagrammatic technique to compute the required (sub)set of conditional self-energies Σ m using rules which are very similar to the standard ones [App. D].
To set up more refined CP approximations one may consider additionally modifying the Σ m but doing so without uprooting CP is more complicated. In Sec. 3.3 we discussed a second approximation scheme which achieves this more easily by selecting one-branch diagrams for the Keldysh operators k m . In the next section we set up a hierarchy tailored to this purpose.

Hierarchy for Keldysh operators k m
Because the Keldysh operators k m are constructed from the Π m by the cutting rule [Sec. 3] they obey an analogous hierarchy of equations (29) illustrated in Fig. 9(b). In this case, we start from the unraveling (15) of the unitary evolution into its average k 0 and operators which are explicitly normal-ordered: We now introduce σ m as the sum of all one-branch irreducible Keldysh half-diagrams representing k m , i.e., irreducibility is understood with respect to the time-ordering on the upper branch of the Keldysh time-contour only. The Keldysh operators k m are then generated by the self-energy operators σ m through the following hierarchy 23 that inherits its structure from Eq. (26): Here u is the free evolution of system and environment which reduces to u = 1 S ⊗ 1 E in the interaction picture. As before, we symmetrize over the mode indices on the right hand side [Eq. (26)]. We note that the sum of Eqs. (29) results in the self-consistent equation 22 TP is preserved through cancellation of pairs of diagrams appearing in different conditional evolutions: Eq. (20) implies that terms in Πm q ...m 1 cancel with those in 'adjacent' Πm q±1 ...m 1 . Since the argument leading to Eq. (20) also holds when removing Π0 at the latest time, we find that contributions to the irreducible block Σm q ...m 1 cancel with those in other blocks Σm q±1 ...m 1 . 23 Refs. [44,60,61] derive a hierarchy for Kraus operators acting only on the system. Here we instead have a hierarchy of Keldysh operators which still act on the environment. Moreover, the hierarchy is expressed in terms of self-energy generators σm which provides some advantages.
for the operator k(t) := m k m (t) with σ(t, t ) := m σ m (t, t ). Importantly, this is not the self-consistent form of the Schrödinger equation due to the lack of normal-ordering relative to Eq. (28). Only after explicitly normal-ordering the hierarchy equations, the Schrödinger-Dyson equation, U = u+u * (−iV )U , with a time-local coupling V is recovered, as it should. This nontrivial consistency check is worked out in App. B.2.
The zeroth tier Eq. (29a) is again an independent self-consistent equation. In this case, the self-energy σ 0 sums up all one-branch irreducible contractions without external fields and therefore acts trivially on the environment just as k 0 . The solution k 0 = u + ∞ l=1 u( * σ 0 * u) l of Eq. (29a) can then be related to the solution of Eq. (26a), noting that the zeroth operator [Eq. (11)] factorizes, k 0 = K 0 0 ⊗ 1 E = U E ⊗ 1 E , and that the free system evolution is π = Tr E u(•⊗ρ E )u † . Although this separates the evolution into operators on the left and right, in general K 0 0 (t) is not generated by an effective nonhermitian Hamiltonian because σ 0 (t, t ) is in general a time-nonlocal operator [cf. Eq. (42)].
Higher tiers determine the k m which act nontrivially on several modes m of the environment. When one has solved the hierarchy for these operators one obtains the conditional propagators Π m by the pasting formula (17b). The Kraus operators K e m (t) = e| :k m (t): ⊗1 E |0 are determined by the solutions for k m (t) obtained from Eq. (29).
Taking the time-derivative of (29) we obtain a hierarchy of time-nonlocal kinetic equations for the Keldysh operators in the interaction picture, with the initial conditions k 0 (0) = 1 and k m (0) = 0 for m = 0. As before, these equations enable an iterative solution of k 0 , k 1 , k 2 , etc., given the conditional self-energy operators σ 0 , σ 1 , σ 2 , . . .. Importantly, this allows one to implement the second kind of approximation strategy discussed in Sec. 3.3. We note that summing the column-entries of Eq. (32) gives a single time-nonlocal evolution equation, which is the differential formulation of Eq. (30), but not the Schrödinger equation due to the lack of normal-ordering in the quantity k. Accounting for the normal-ordering, Eq. (32) [Eq. (29)] constitutes a time-nonlocal unraveling of the (Dyson) Schrödinger equation for system plus environment which has the advantage that it can be solved step-by-step. An advantage relative to Eq. (27) [Eq. (26)] is that the difficulty of guaranteeing CP and TP has been reversed. Indeed, in Eq. (29)- (32) we are free to formulate any approximation in terms of the generators σ m without ever compromising CP because they feed into the Keldysh operators k m constituting the CP operator-sum. Importantly, such k m -approximations will not only select conditional evolutions Π m but also modify them [cf. Sec. 3.3].
Such a one-branch diagrammatic selection does not guarantee TP because the network of pairwise cancellations between different Keldysh operators [Eq. (20), Fig. 8] implies constraints that tie all σ m together. This makes TP harder to implement than the twobranch propagator hierarchy (27).

Operator-sum theorem and memory kernel of the quantum master equation
The  (27) [(26)]. This implements the state-evolution correspondence [App. C] while explicitly separating out the normal-ordering procedure. The nontrivial reversal of the difficulty of imposing the constraints of CP and TP -despite the structural similarities of the two types of hierarchies-results from this duality and will now be traced microscopically to the difference in irreducibility and time-ordering on one and two time-branches: Whereas the self-energy superoperators Σ m were constructed from two-branch irreducible diagrams, the self-energy operators σ m are irreducible on a single branch, ignorant of time variables on the opposite branch 24 .
To clarify this, we explicitly construct Σ m from the σ m by working out the nontrivial time-ordering constraints 25 in the Kraus operator-sum theorem (for finite time evolution) in order to transpose it to the level of the Nakajima-Zwanzig memory-kernel (for infinitesimal evolution). Here in particular, our diagrammatic approach shows its advantage: The functional form of the self-energies Σ m [σ 0 . . . σ m ] can be precisely determined starting from the prescription provided by the irreducible version of the pasting -or purificationequation (17b): Inserting k m [σ 0 . . . σ m ] in terms of its generators, it remains to work out the two-branch irreducibility constraints on the one-branch time-integrations hidden inside k m and k † m . We now show that this amounts to a 'dressed' real-time expansion on the Keldysh contour in terms of effective vertices σ 0 , σ 1 , . . . with an intricate internal time-ordering structure.
In Fig. 10(a) we sketch the case for Σ 0 where the two-branch irreducibility can be achieved in only two ways using one-branch irreducible building blocks. Either a generator σ 0 (green) spans over one whole branch (first diagram); then the time-integrations on the other branch remain unrestricted and the full one-branch propagator k 0 appears (bold line). Otherwise, the two-branch irreducibility can only be maintained by an alternating sequence of generators σ 0 on the upper and σ † 0 on the lower branch with temporal overlaps, filled up with k † 0 resp. k 0 propagations in between (further diagrams). We thus see that the basic blocks of the memory kernel on two branches are four-time superoperators 26 describing irreducible time-evolution from time t + to t + on the upper branch and reducible evolution from t − to t − on the lower branch, with t + ≥ t − ≥ t − ≥ t + and conversely 24 This is a fundamental difference between the Liouville-von Neumann and Schrödinger equation which by our considerations ties in with the state-evolution correspondence. 25 In special limits where the evolution is a dynamical semi-group this goes unnoticed: there one can easily apply the operator-sum theorem directly to Π(t + dt) = Π(dt)Π(t), leading to the GKSL form of the quantum master equation [68]. Here we instead consider the more complicated general evolutions Π(t) where this property fails. 26 Similar objects are considered in Refs. [44,60,61].  Fig. 7(a)]. An explicit construction of Σ 2 is provided in App. F. The irreducible fragments Λ m obtained in step (iii) can be complemented by the generalization Λ 0 (t + , t − , t + , t − ) of Eq. (36a) to a four-time object, see rule (iv).
with t − ≥ t + > t + ≥ t − . Two-branch irreducibility is enforced by letting the one-branch irreducible generator σ 0 (σ † 0 ) start earlier and end later than their reducible counterparts k † 0 (k 0 ) on the opposite branch. In terms of these auxiliary objects the functional form of Σ 0 can be written down explicitly: Here, the two-branch convolution accounts for independent one-branch convolutions imposing the two-branch time-ordering constraints of Eq. (35). Once Σ 0 has been constructed in this way, the remaining tiers of the self-energy Σ m can be obtained from a dressed expansion governed by the following diagrammatic rules illustrated in Fig. 10(b). We start from a Keldysh contour where each branch, running from the external times t to t, is dressed to the Keldysh operator k 0 [k † 0 ]. To obtain Σ mq...m 1 (t, t ) for fixed modes m q . . . m 1 perform the following steps: (i) Distribute in all possible ways the q external vertices acting on the fixed set of modes m q . . . m 1 among self-energy blocks σ µ [σ † µ ] with µ = 0 on the upper [lower] branch. (ii) Integrate over all internal time-arguments of these blocks, but restrict them to definite two-branch time-orderings, see the first three contributions in Fig. 10(b). Let the earliest (latest) self-energy block start (end) at t (t). Note that a single block may extend over a whole branch. Sum over the finite number of possibilities of such orderings.
(iii) Contract the external vertices in all possible ways to obtain two-branch irreducible fragments Λ µ (t + , t − , t + , t − ) as shown for the first three contributions in Fig. 10(b). Analogous to Eq. (35), the four time-arguments indicate propagation times on the upper and lower contour respectively .
(iv) Insert the four-time generalization of Eq. (36a), Λ 0 (t + , t − , t + , t − ), at all positions where the two-branch irreducibility is not yet ensured. Such cases only occur for tiers with more than one mode (q ≥ 2), see App. F. This explicitly enforces two-branch irreducibility.
(v) Add all irreducible extensions of the constructed diagrams by attaching Λ 0 on the left or on the right or on both sides as shown in the last contributions of Fig. 10(b).
Due to this generic structure, the two-branch irreducible fragments Λ µ [σ 0 . . . σ m ] form yet another hierarchy for the conditional self-energy superoperators, where the two-branch convolution takes care of the intricate time-ordering required for CP. Here ∆(t + , t − , t + , t − ) = δ(t + − t + ) δ(t − − t − ) is a shorthand for 'no insertion'. Except for the ∆'s, Σ m inherits a structure similar to that of the hierarchy (26) for Π m . Also here symmetrization over the mode labels is understood implicitly.
The above diagrammatic rules allow different terms of the functional to be explicitly written down as illustrated in App. F for a nontrivial example in second order of the dressed expansion. Importantly, it is again only the zeroth tier Σ 0 that requires summation over an infinite set of diagram contributions via the relation (36a). The remaining contributions -dressed up by the so obtained Σ 0 -represent a finite number of terms to evaluate. As of course expected for this general problem, evaluating these terms is challenging. In the limits where the GKSL form applies the arising complications disappear and the Lindblad jump operators automatically emerge from the conditional self-energies σ m (t, t ) of the Keldysh operators with external fields (m = 0) as the example in Sec. 5 illustrates.
The functional form (36) precisely encodes the nontrivial structure of the memorykernel of the Nakajima-Zwanzig quantum master equation imposed by the Kraus operatorsum theorem. This means that whatever (approximated) σ m one uses, a memory kernel Σ m (or Σ = m Σ m ) with the functional form (36) will generate a CP evolution Π m (Π) through the kinetic equation (27) [Eq. (3)]. It is this structure of the kernel -enforced by Eq. (34)-that guarantees that the solution can be expressed as Π m = Tr E :k m : •ρ E :k † m :. The general importance of the nontrivial functional (36) is that it establishes a microscopic framework for approximations that go beyond the restrictive limits of the GKSL approach without giving up CP. This is by no means easy, but we have found non-Markovian examples for which GKSL does not apply but which are nevertheless solvable with our approach [66], see also examples in Refs. [45,60]. Knowing the explicit structure (36) of the self-energy is also a prerequisite for setting up functional renormalization-group (RG) approaches. It may thus provide a starting point for connecting the exact hierarchies of RG equations formulated on one [75] and two [35,42] branches of the Keldysh contour while explicitly preserving the CP structure of the density operator evolution. This would allow one to assess the implications for CP of truncation schemes for such RG hierarchies. Interestingly, some RG approaches [35,38,42] also entail a two-stage approach similar in spirit to ours but based on an RG flow starting from the Markovian T = ∞ limit [25,38] that we study next in Sec. 5.
We stress that the functional (36) is also relevant for less ambitious efforts. Consider, for example, a TP-approximation based on a selection of complete diagrams [27][28][29][30][31][32] such that it breaks CP [Sec. 3.5]. Then the functional (36) allows one to systematically identify which groups of diagrams are missing, what their time-ordering structure is, and provides a guide for adding a some -but not all-terms to improve towards 27 a complete operatorsum form. This amounts to partially "completing the square" of the solution Π, while working directly on the level of the self-energy of the kinetic equation.

Simple example:
Fermionic mode coupled to a hot fermion bath In this final section we illustrate the developed ideas on a simple example for which all calculations can be explicitly followed, both on the level of the Keldysh operators k m and the conditional propagators Π m . The aim is to illustrate how CP is guaranteed explicitly by the dual hierarchies when systematically constructing the self-energies and then solving each of their tiers. Although other methods also solve this model [see Eq. (52)], these do not show how CP is guaranteed in the general case of ultimate interest. It suffices to consider a fermion mode coupled to a continuous reservoir of noninteracting fermions held at infinite temperature (T = ∞). The total Hamiltonian of this model is given by and we work in Schrödinger picture. The system (field d) is couped bilinearly to the modes in the environment (field b ω ) by tunnel rates Γ that are independent of their energy ω (wideband limit) denoting ω := ∞ −∞ dω. We let the particle-hole index η indicate a creation (η = +) or annihilation (η = −) operator 28 In this example, the modes of the environment are thus labeled by particle-hole indices and continuous energies: (m q . . . m 1 ) = (η q ω q . . . η 1 ω 1 ).
The combined T = ∞ and wideband limit leads to simplifications whenever we contract environment modes and sum the mode energy, This limit therefore causes strictly Markovian semi-group dynamics and renders the model solvable. Here and below, all δ-functions are adapted to working with time-convolutions: they are defined such that τ 0 dτ δ(τ − τ ) = 1 instead of 1/2. We note upfront that the T = ∞ and wideband limit needs to be taken twice: once, when we compute the conditional propagator Π m in stage 1 by summing one-branch-diagrams to obtain k m which is then pasted together with k † m using Eq. (17b); we take the limit a second time in stage 2 where we integrate over the environment energy ω contained in the mode index m = ηω. Importantly, the pasting result (17b) of stage 1 cannot be simplified without this stage 2 energy-integration. As a result, the stage 1-2 distinction is easily blurred when simplifying expressions. We will first illustrate this for k m in Sec. 5.1. This becomes a more prominent issue when computing Π m directly -without the use of Keldysh operatorswhich we illustrate in Sec. 5.2. Figure 11: Hierarchy of the Keldysh-operators (left) and conditional propagators (right) for the model (37). The combined infinite-temperature and wideband limit forces all contractions to be instantaneous in time (drawn as vertical) suppressing all vertex correction diagrams. In combination with the bilinearity of the coupling V this causes all self-energies σ m and Σ m of tier m ≥ 2 to vanish identically, but only after integrating over the environment mode energies ω i , thereby partially performing stage 2.

One-branch hierarchy for k m
Starting with the lowest tier of the hierarchy Eq. (29), we first have to compute the generators for the Keldysh operators. Only the single bubble-diagram depicted in Fig. 11(a) contributes which translates to We next construct the Keldysh operator by iterating the self-consistent equation (29a): In contrast to the general case [Eq. (31) ff.], the zeroth tier evolution is generated by a nonhermitian Hamiltonian for system and environment, H 0 ∞ := H 0 − i Γ 2 1 S ⊗ 1 E , because σ 0 is time-local in this model.
As sketched in Fig. 11(b), the next tier self-energy operators are where we explicitly write the mode index 1 = m 1 = (η 1 ω 1 ). Using this we can then evaluate Eq. (29b) for the Keldysh operator without the need for iteration: For higher tiers q ≥ 2, the self-energies vanish σ mq...m 1 = 0 because the wideband and T = ∞ limit suppresses all one-branch vertex corrections. Thus, only repeated Figure 12: Energy-integrated functional Σ 1 [σ 0 , σ 1 ], cf. Fig. 10(b). Due to the energy integration, the two-branch convolution reduces to a δ-function in time.
convolutions of k 1 remain for the q-th tier equation and we obtain with implicit symmetrization over the mode-indices. The time-ordering on the integrals has been lifted with the prefactor 1/q! accounting for the corresponding double-counting. We see that the decay rate Γ inside the exponential (42) appearing in Eq. (45) is generated initially when solving the lowest hierarchy-tier for k 0 (t). Through the hierarchy it sets the timescale for all higher k m (t) before we have integrated out the meter M by performing the sum over m. We stress that to further simplify the results for Π m in stage 1 we must partially perform stage 2, m Π m , in order to exploit the wideband and T = ∞ limit for a second time: only when integrating Π m = Tr E :k m :(• ⊗ ρ E ) :k † m : over the energies ω i of the modes m i (leaving the η i untouched) can the time-integrations in Eq. (45) be reduced to further δ-functions in time.

Two-branch hierarchy for Π m
First, the zeroth tier conditional propagator is obtained directly by pasting together the Keldysh operators (42): Alternatively, one computes Σ 0 from (36a), where due to time-locality only the first type of diagram in Fig. 10(a) contributes: For the next tier, we can again proceed in two equivalent ways. The direct approach inserts Eq. (44) into the pasting formula Π 1 = Tr E :k 1 : ρ E • :k † 1 : to obtain Eq. (49) below. Alternatively, one can bypass computing k 1 altogether by computing Σ 1 and then generating Π 1 through the two-branch hierarchy. Using our general functional (36b) to compute Σ 1 [σ 0 , σ 1 ], only the single pair of diagrams shown in Fig. 12 contributes. Its construction from k 0 [Eq. (42)] and σ 1 [Eq. (43)] is straightforward with the two-branch convolution: As pointed out above, we are now forced to integrate over the mode energies to exploit the wideband and T = ∞ limit a second time. This causes the general functional (36) to vastly simplify due to the time-locality of the contractions between the opposite time-branches. Equivalent manipulations are also required in the evaluation of the direct approach.
Analogous to Eq. (44), we can then evaluate the first tier equation (26b) without iteration for the conditional propagator: Importantly, the time convolution produced the linear factor t with time-scale Γ ∼ Σ 1 . All higher self-energies identically vanish only when integrated over mode energies ωq...ω 1 Σ mq...m 1 ≡ 0 for q ≥ 2 where the wideband and T = ∞ limits suppress vertex corrections generated by contractions between two branches. For the q-th tier equation integrated over the energies only repeated convolutions of Σ 1 remain, where we used Π 0 * Σ 1 * Π 0 = Π 1 and inserted Eq. (49). This agrees with the direct approach which pastes k m [Eq. (45)] to obtain Π m . For each set of modes integrated over their frequencies, we get a single Kraus operator 29 : The prefactor (Γt) q /q! is generated only in Eq. (50) where we integrate over the modeenergies of the inter-branch contractions, i.e., we partially evaluated stage 2. Doing so results in the δ-functions (second wideband limit), and the time-ordered integral reduces to t q /q!. As we noted at Eq. (45), this cannot be seen in k m or Π m obtained in stage 1. We complete stage 2 by summing over the particle-hole indices η i , and find the GKSL generator Although we could have derived a GKSL master equation ∂ t ρ(t) = Lρ(t) -which is exact for this model-to more easily obtain the result (52), all steps and physical arguments that we discussed can also be applied for general evolutions -for which no GKSL equation can be derived. Clearly, in the latter case the evaluation steps will be technically more difficult.
Discussion. The two distinct occurrences of the coupling constant Γ in the simple results for the time-dependence of the Keldysh operator (45) and the Kraus operator (51) illustrate the two fundamental stages of the real-time expansion. As we anticipated in Sec. 3.3, the physically motivated hierarchy index m = ηω (measurement outcomes) 'counts the powers' in stage 2 of the expansion. This yields the series in powers q of the prefactor Γ in Eq. (52) which is CP order-by-order.
The Γ in the exponential is instead generated in stage 1 and controls how many terms in stage 2 should be kept for a good approximate result. Truncating terms beyond q = q max introduces a bounded TP error [cf. Sec. 3.5], which quantifies the neglected process in which more than q max interactions with the reservoir took place. The contribution of the q-th order process is maximal for times t = qΓ −1 . Due to the complete positivity of this approximation, the error (53) monotonically decreases as q max is increased. The required value of q max for the TP error to be negligible depends on the time scale of interest: for Γt 1 a small q max suffices, while for Γt 1 a large q max is required. This is expected since TP guarantees that the Π(t) has a stationary state. In contrast, CP tends to be more critical at short times, see the discussion in Sec. 1.3.
Finally, we note that the above illustrates how easily limiting procedures -even in an exactly solvable problem-can hide the CP structure (17) of the real-time expansion. We stress that even though the energy integrations -partially accomplishing stage 2-were necessary to analytically solve the problem, they were never needed to reveal CP: in our formulation at each step the conditional propagators have the operator-sum form. Note that Eq. (50) written without energy integrations requires the higher self energies Σ m to recover the explicit operator sum Π m = Tr E :k m : ρ E • :k † m :. Then, even in this simple model, vertex corrections are relevant and the Σ m are not zero for m = 0, 1.

Summary
Measurement-conditioned evolution. In this paper we studied the dynamics of the reduced density operator ρ(t) = Π(t)ρ(0) by combing the operator-sum representation -ubiquitous in quantum-information theory -with the Keldysh real-time expansion of the statistical field theory. We decomposed the reduced evolution Π into elementary conditional evolutions Π m which are completely positive because they correspond to possible measurement outcomes m in a quantum circuit [ Fig. 3]. This way we recovered the standard real-time expansion in a crucially reorganized two-stage form, revealing it actually consists of two intertwined expansions [Eq. (17)]. Based only on fundamental operational principles of measurements this extends the idea of an expansion in alternating 'life-time broadened evolution' and 'quantum-jumps' to general CP-TP dynamics (1), in particular to dynamics that cannot be described by GKSL quantum master equations or its generalizations.
Cutting diagrams -unraveling the Schrödinger equation. By cutting the Π m superoperator diagrams on the Keldysh double-time contour into halves we obtained what we called Keldysh operators k m on a single time branch. These 'lift' the Kraus operators K e m (t) = e| : k m (t) : ⊗1|0 to act on both system and environment, encoding their timeevolution [Eq. (10)] and the Wick normal-ordering. From the obtained diagrammatic expansion we inferred memory kernels for these quantities and derived a hierarchy of self-consistent evolution equations for the Keldysh operators [Eq. (29)] . This hierarchy generates time-nonlocal, non-Hamiltonian dynamics on one time branch.
The corresponding hierarchy of kinetic equations for the Keldysh operators [Eq. (32)] constitutes a time-nonlocal unraveling of the Schrödinger equation conditioned on the measurements m. The advantage over the original Schrödinger equation is that these onebranch kinetic equations already explicitly encode the operator-sum decomposition that guarantees CP before having integrated out the environment completely.
Pasting diagrams -unraveling the quantum master equation. By pasting the obtained Keldysh / Kraus operator diagrams back together -reconnecting the time branches-we obtained an equivalent self-consistent hierarchy [Eq. (26)]. This expresses the operatorsum theorem on the Keldysh real-time contour of statistical field theory.
The corresponding hierarchy of kinetic equations [Eq. (27)] unravels the general (Nakajima-Zwanzig) quantum master equation. We expressed its two-time-branch memory kernel superoperator explicitly in terms of the one-time-branch memory kernels of the Keldysh operators that also determine the Kraus operators [Eq. (36)]. Irrespective of which approximations to the latter are inserted, a Nakajima-Zwanzig memory-kernel with this functional form generates solutions which are guaranteed to be CP. This puts the celebrated GKSL characterization of time-evolution generators on a more general footing.
The tiers of our hierarchies are labeled by a physical quantity, the modes which have interacted with the system as determined by a possible measurement on the environment, and this guarantees CP order-by-order. We have illustrated how this hierarchy can be truncated based on the dissipative time-scale generated at an earlier stage [Eq. (51)]. It is no surprise that our simple rules in general lead to quantities that are difficult to calculate. However, the merit of the approach is that it indicates new routes to systematically address the open problems of devising and comparing CP approximations, improving positivity in existing TP approximations, and further exploring the most challenging issue of CP-TP approximations [44,45,60,61]. The very recent advances reported in Ref. [60,61] indicate that there is still room for further progress.
The role of time. A recurring theme has been the complementarity between the approaches from quantum-information and statistical field theory that we combined. This was not entirely unexpected given the celebrated correspondence between reduced quantum evolutions and bipartite quantum states [70][71][72]. However, we showed how this entails a duality in approximating reduced time-evolutions -strictly fixing CP tends to uproot TP and vice versa-and our two types of self-consistent and kinetic equations explicitly encode this. By descending to the level of the microscopic coupling V (t) -instead of staying on the formal level of the joint unitary U (t) = T e − t 0 dτ V (τ ) -we exposed the crucial role of time as the common thread parametrizing both evolution operators and super operators. The corresponding difference between time-ordering on one-and two-branches of the Keldysh contour led to two notable mutual insights.
The operational (input-output) approach of quantum information theory ignores time and hides the one-branch time-ordering inherited from the unitary U inside the Kraus operators. This trivializes the CP structure through the operator-sum theorem, but makes TP seem a nontrivial 'global' constraint in practice. In contrast, the statistical physics approach, formulated on the two branches of the Keldysh contour, trivializes the TP structure because the trace connects the two branches, but thereby makes CP seem a nontrivial 'global' constraint. On the microscopic level of the coupling V we found that these two approaches resolve each other's enigmas: the nontrivial TP constraint for Kraus / Keldysh operators living on one time branch is achieved by pairwise cancellations [Eq. (20), Fig. 8] when considering the double time-branch. Likewise, the nontrivial CP constraint on Keldysh diagrams is achieved by grouping their two-branch diagrams according to one-branch time-ordering [Eq. (17), Fig. 7(b)].
Funding information V. R. was supported by the Deutsche Forschungsgemeinschaft (RTG 1995).

A Normal-ordering
The normal-ordering : : of a string of field operators X 1 . . . X q is defined relative to some average . In the main text we mostly consider • = Tr E •ρ E where ρ E is a grand canonical equilibrium state of noninteracting fermions or bosons. In the following we focus on this case as well, but note that partial contractions and normal-ordering can be generalized to interacting environments [49,76].
The construction of normal-ordered expressions follows an iterative reduction into paircontractions X 1 X 2 ≡ X 1 X 2 . This is based on the observation that the total average X 1 . . . X q decomposes into pair contractions by Wick's theorem. With the assumption X i = 0, one recursively defines X 1 = :X 1 : (54a) X 1 X 2 X 3 = :X 1 X 2 X 3 : +X 1 X 2 :X 3 : +X 1 :X 2 : X 3 + :X 1 : X 2 X 3 (54c) . . . and inductively verifies that the average of :X 1 . . . X q : is zero for each q. Moreover, each partial average is zero, for example, :X 1 X 2 X 3 : = 0. For this to work, the partial contractions must include a sign for disentangling in the case of fermions, X 1 :X 2 : X 3 ≡ −X 1 X 3 :X 2 :. Note that for X i = 0 the recursive construction must be extended by substituting X i → X i − X i 1 in Eq. (54) together with :1: = 0.
The sole purpose of normal-ordering is to reorganize the way contractions of X 1 . . . X q are performed: operators within a normal-ordered sequence can only be contracted with operators from a different sequence. Thus, :X 1 . . . X q : :Y 1 . . . Y q : is evaluated by Wick's theorem in the usual way except that the normal-ordering restricts all pair contractions to be of the type X i Y j . The key point used in Eq. (8) is that such an average is nonzero only when the number of operators and the modes on which they act match. Moreover, such an average does not require the calculation of the operators :X 1 . . . X q : and :Y 1 . . . Y q :.
In the main text and in App. B we use this property of normal-ordering to decompose any operator A m acting only on a definite set of environment modes m = m q . . . m 1 into its environment average plus all its normal-ordered partial contractions which systematically account for all possible fluctuations: Here the shorthand notation l ⊂ m denotes a subset of modes on which the remaining uncontracted fields in A l m act and we sum over all possible subsets excluding the empty set l = 0 (not indicated).
In the last line we then collect into U m all strings of repeated couplings which together act on the definite set of 30 environment modes m = m q . . . m 1 . We do not distinguish the order of modes, i.e., U mq...m 1 is symmmetric. Next, we decompose 31 U by expanding each component Here, the Keldysh operator k m is defined for m = 0 as the partially contracted evolution operator U such that only environment modes m ⊂ m remain: The special case denoted by m = 0 is the part of U in which no environment operators are left: it collects the full contraction of U , i.e., its environment average We note that the Keldysh operators are themselves not normal-ordered: k m = :k m :, i.e., k m has nonzero partial contractions for m = 0 while k 0 = k 0 ⊗ 1 E and :k 0 : = 0. Yet, partial contractions of these objects never appear in the theory because the pasting rule Eq. (17b) explicitly normal-orders the Keldysh operators. 30 The discussion applies to any normal-ordered coupling V . Only for the special case of bilinear coupling, does the number of modes count the coupling power, Um 1 = u * Vm 1 u, Um 2 m 1 = u * Vm 2 u * Vm 1 u, and so on. If V is instead bi-quadratic as, for example, in Appelbaum-Schrieffer-Wolf Hamiltonians describing Kondo physics, then Um 2 m 1 may be of first order in the coupling V . 31 We are not normal-ordering U , i.e., altering U → :U :. In Eq. (58) we merely add and subtract canceling partial contractions in the real-time expansion of U to separate out different kinds of fluctuations.

B.2 Consistency with the Schrödinger equation
The main text noted that summing the hierarchy (29) results in a self-consistent equation, with the sum of all self-energies σ ≡ m σ m . Importantly, this is not the self-consistent form of Schrödinger equation because of the lack of explicit normal-ordering, k(t) = U (t) and also σ(t, t ) = −iV (t)δ(t − t ).
To check that the hierarchy (29), when properly normal-ordered, indeed unravels the Schrödinger equation, we must decompose k = k 0 + (k − k 0 ) and separately consider the operator k − k 0 acting nontrivially on the environment. From Eq. (61a) we find that k also obeys the self-consistent equation Subtracting k 0 , we rewrite the right hand side by inserting the zeroth hierarchy equation (29a), k 0 = u + u * σ 0 * k 0 , and again substitute Eq. (62) in the resulting term: which involves the following detailed diagrammatic observations. First, the expansion (55) is applied to Eq. (65a). The first term in Eq. (65b) is the average which sums up all full contractions of fields in V n = :V n : with those in :k n :. The corresponding diagrams can be factorized into two parts. One factor is irreducibly contracted in some time interval [t, t ] such that summing all contributions to it gives the self-energy component σ 0 (t, t ). The remaining factor for the time interval [t , 0] contains all possible contributions without external modes and therefore sums up to k 0 (t ). Since in k n (t) we integrate over all internal times, we also integrate over t and obtain the convolution in the first term in Eq. (65c).
The second term in Eq. (65b) involves 32 (−iV n :k m :) l where a subset l ⊆ nm of modes in V n and :k m : remains uncontracted. Denote by c the -possibly empty 33 -subset c ⊆ n of modes in V n that are contracted with a corresponding subset c ⊆ m of modes in k m . Again the diagrams factorize. One factor contains all diagrams irreducibly connected to V n in some time-interval [t, t ] leaving modes n from V n and modes m from k m uncontracted. Summing over all modes c, we get by definition σ n m (t, t ) representing all such irreducible contributions. The factor at earlier times [t , 0] contains all possible contributions with further uncontracted modes m from k m and therefore sums up to k m (t ).
Thus, we have decomposed the modes in V n as n = n c and in k m as m = cm m such that the uncontracted modes n m and m are attributed to the irreducible and reducible factors of each diagram, respectively. In total, the modes l = n m m thereby remain uncontracted. Independently summing over the modes n m and m -including the empty sets-we obtain σ * k [Eq. (61)] for the second term in Eq. (65c), completing the derivation.

C CP-TP duality: State-evolution correspondence
The duality between CP and TP constraints noted in Ref. [67] can be understood in a general form by combining insights from quantum-information theory and statistical field theory. In the main text we focused on pin-pointing the microscopic origin of this CP-TP duality, i.e., on the level of individual diagrammatic contributions to the Keldysh expansion, rather than using the formal arguments that we now summarize.
In general, the state-evolution connection is provided by the de Pillis-Jamiolkowski-Choi map (22) which turns a superoperator Π into an operator on the double system: Here |1 = k |k ⊗ |k is the maximally entangled state on the double system. Thus evolutions are mapped to states, providing many helpful insights [69]. Kraus form. When the propagator is written in Kraus form, Π = m K m • K † m , the Choi operator reads choi(Π) = m |K m K m | with vectorizations of the Kraus operators |K m := (K m ⊗ 1)|1 . One independently shows that Π is a CP evolution superoperator if and only if choi(Π) is a positive (semidefinite) operator and thus represents some quantum state of the double system. The Kraus form of Π translates CP into a property which can be assessed by inspecting individual terms: physically, it is clear that the statistical mixing of individual pure states |K m K m |, always leads to a positive state choi(Π).
On the other hand, one verifies from Eq. (66) that Π is a TP superoperator if and only if the partial trace of the operator choi(Π) over the first subspace of the tensor product gives the identity on the second: Physically speaking, this requires that the Choi state (which is already a mixed state) has the maximally mixed state as a marginal state, i.e., Eq. (67) constrains the correlations on the bipartite system. This is clearly a nontrivial property to enforce on a decomposition into pure states: indeed, inserting the Kraus form one obtains m Tr 1 |K m K m | = 1 which is equivalent to the condition m K † m K m = 1 mentioned in the main text [Eq. (13a)] and puts a constraint on all terms.
Keldysh Dyson expansion. Now consider the same superoperator Π written in the form of a Dyson series Π = π + ∞ k=1 (π * Σ * ) k π with respect to the physical free evolution π which is CP-TP. The TP property of π guarantees that Π is TP as well, provided that the remaining terms have zero trace, Tr Σ = 0. This condition is guaranteed by simple pairwise cancellation of diagrams discussed in Sec. 3.5.
However, the CP constraint is nontrivial due to the complicated form choi(Π) = choi(π) + ∞ k=1 choi (π * Σ * ) k π of the Choi operator inherited from Π: although choi(π) is positive due to the complete positivity of π, the remaining terms do not make clear that choi(Π) is indeed positive and, consequently, that Π is CP.
We thus see that the state-evolution correspondence (66) changes the difficulty of ensuring CP at the expense of ensuring TP and vice-versa. This explains the noted [67] 'incompatibility' of simultaneously imposing these constraints in practice.
The diagrammatic rules for the contributions to Eq. (69) are the following: (a) Distribute n [n ] vertices on the left (right) in all possible ways over the k = n + n ordered times t k ≥ . . . ≥ t 1 and sum over these diagrams.
(b) Pairwise contract all fields appearing in the vertices V and sum over all possible contractions. The sum over all modes is included in the vertices themselves, V = m =0 V m . For fermions, the sign of each term is given by the parity (−1) nc of the number n c of crossing contraction lines.
(c) Perform the ordered time integrations.
(d) To obtain the conditional propagator Π mq...m 1 , restrict in (b) the q inter-branch contractions to a definite set of modes m q . . . m 1 while summing over all modes of intrabranch contractions. For Σ mq...m 1 restrict this to two-contour irreducible diagrams.
Remarks: (i) The times t i label vertices on either contour: in the Keldysh technique each twobranch time-ordered integral is explicitly represented by a separate diagram. This timeordering allows the transition to a Laplace frequency conjugate to the physical time in ρ(t) which is of practical use in calculations [20,21,23,24,27,29] and crucial for the formulation of RG schemes [35][36][37][40][41][42]. The conversion between single-contour time-ordering (τ i ) used in Eq. (68) to two-contour time-ordering (t i ) in Eq. (69) as is discussed in the main text [Sec. 4.3] has also been useful in other contexts [24].
(ii) In the density-operator technique the Keldysh contour is not drawn closed at latest time [cf. Fig. 5(b)], in contrast to the Keldysh contour for Green's functions [56] where this closure indicates an additional trace over the system to obtain correlation functions.
(iii) For the enviroments of interest here, the modes are labeled by discrete indices (particle type, electron spin, etc.) and a continuous orbital index or equivalent energy [35,42], see Eq. (39). Thus in the above rules, the CP structure requires suspending integrations over the continuous environment energies for the inter-branch contractions but performing them for intra-branch contractions, see Sec. 5.
In practice, to compute the superoperators represented by the diagrams one needs to exploit the known energy eigenbases of H and H E , the system and environment Hamiltonians, respectively. It is then convenient to revert to the Schrödinger picture which only changes the translation rules for the diagrams by inserting the free evolution π(t, t ) = e −iT t t dτ H(τ ) • e +iT t t dτ H(τ ) with formal (anti)time-ordering T (T ) between every two consecutive occurences of vertices at times t ≥ t -irrespective of the branch they occur on.

E Diagrammatic rules for Keldysh operators k m
A key result of the paper [Eq. (17b)] is that the conditional propagator Π mq...m 1 , determined by the rules in App. D, can be 'cut into two halves', representing the Keldysh operator k mq...m 1 and its adjoint, respectively. The Keldysh operators can therefore be obtained directly from a diagrammatic 'cutting and pasting' based on similar rules.

E.1 Cutting rules
The contributions to the Keldysh operator k mq...m 1 (t) are obtained from the formal expansion of the auxiliary operator using the following rules: (a) Pairwise contract all fields appearing in the vertices V and sum over all possible contractions, leaving q fields with a definite set of mode indices m q . . . m 1 uncontracted. For fermions, the number of crossing contraction lines n c,intra introduces a sign (−1) n c,intra , where the crossings with the external contraction lines must be included.
(b) Assign the ordered times τ q ≥ . . . ≥ τ 1 to the uncontracted vertices and integrate over them.
(c) Independently sum over all orderings of the mode indices m q . . . m 1 .
To obtain the self-energies σ mq...m 1 , restrict to one-branch irreducible diagrams. The central objects are the Keldysh operators k m and their self-energies σ m . The Kraus operators can be obtained from these, where |e = ij e ij |i ⊗ |j denotes any desired basis of bipartite vectors for the purified environment EE andê = ij e ij |i j| is a corresponding basis of operators acting on E.

E.2 Pasting rules
The Keldysh real-time diagrams for Π mq...m 1 can be obtained from the above constructed Keldysh operators k mq...m 1 by 'pasting' them together with the vertically flipped diagrams corresponding to k † mq...m 1 . This diagrammatic pasting entails the following rules: (a) Contract in all possible ways the external vertices of k mq...m 1 with those of k † mq...m 1 having the same mode index, remembering that m q . . . m 1 may contain repetitions of the same mode index. For fermions, the number of crossing inter-branch contraction lines n c,inter introduces an additional sign (−1) n c,inter . This produces the correct fermionic sign in App. D since n c = n c,intra + n c,inter .
(b) Assign the ordered times τ 2q ≥ . . . ≥ τ 1 to the 2q consecutive external vertices -irrespective of the branch they occur on-and integrate over them. This produces the two-branch time-ordering in App. D and yields Π mq...m 1 .
(c) Sum over all modes m q . . . m 1 of the inter-branch contractions to obtain Π.
F Functional relation of self-energies Σ m and σ m The general rules (i)-(v) of Sec. 4.3 for constructing the conditional self-energies Σ m from the memory-kernels σ m are illustrated using Fig. 13. As pointed out in the main text, the functional form (36a) for Σ 0 is the only contribution requiring an infinite set of diagrams to be summed up. All following tiers require only a finite set of dressed diagrams constructed from Σ 0 . Here we explicitly discuss the next-to-leading order Σ 2 which covers all difficulties to be encountered when applying the rules to these higher orders.
(i) We first distribute in all possible ways the external vertices among generators σ µ with µ = 0 on each branch. In Fig. 13, the colored boxes group the three possible distributions for Σ 2 : The first group (red) has two σ 1 blocks on each branch. The second group (blue) has two σ 1 blocks on one branch but a single σ 2 on the opposite branch. The final possibility (violet) has a single σ 2 on each branch.
(ii) Within each of these groups, we sum over diagrams with definite two-branch timeordering and integrate over all internal time-arguments. We stress that the number of such orderings is finite, although we refrain from providing all contributions for Σ 2 .
(iii) We contract the individual generators σ µ in all possible ways between the branches, including 'exchange' diagrams that are explicitly shown only in the first (red) group.
(iv) To proceed, we have to distinguish two types of diagrams that we have constructed so far. Two-branch reducible contributions [first three diagrams in Fig. 13] are constructed from multiple two-branch irreducible fragments Λ µ (t + , t − , t + , t − ) 'glued' together to enforce the two-branch irreducibility. This is achieved using the four-time generalization, Λ 0 (t + , t − , t + , t − ), of the self-energy Σ 0 (t, t ) [cf. Eq. (36a)] indicated in gray. Together with the two-branch convolution , we obtain terms of the form for these contributions. The remaining irreducible contributions shown in Fig. 13 are made up of already irreducible fragments and require no such 'fixing': (v) Finally, all constructed diagrams have to be complemented by their possible extensions with Λ 0 blocks on either side or both, as discussed in the main text [Eq. (36)]. These are not shown in Fig. 13.