Fermionic duality: General symmetry of open systems with strong dissipation and memory

We consider the exact time-evolution of a broad class of fermionic open quantum systems with both strong interactions and strong coupling to wide-band reservoirs. We present a nontrivial fermionic duality relation between the evolution of states (Schr\"odinger) and of observables (Heisenberg). We show how this highly nonintuitive relation can be understood and exploited in analytical calculations within all canonical approaches to quantum dynamics, covering Kraus measurement operators, the Choi-Jamio{\l}kowski state, time-convolution and convolutionless quantum master equations and generalized Lindblad jump operators. We discuss the insights this offers into the divisibility and causal structure of the dynamics and the application to nonperturbative Markov approximations and their initial-slip corrections. Our results underscore that predictions for fermionic models are already fixed by fundamental principles to a much greater extent than previously thought.


Introduction
The dynamics of open quantum systems is a problem of interest in a range of research fields. Their higher complexity as compared to closed systems evolving unitarily continues to motivate the development of new frameworks and approximation schemes to make further progress. Complementary to this, it has become more important to maximally reduce this complexity within existing well-developed approaches using basic symmetries and other general structures, see, e.g., Ref. [1] and references therein. For closed quantum systems, simplification by exploiting symmetries for some fixed set of system parameters is a highly developed subject and builds on the unitarity of transformations and the corresponding Hermicity of its generators. When turning to dynamics of open systems one runs into interesting new problems because the latter properties are lost in a reduced description.
In this paper we instead consider a different kind of simplification offered by a duality mapping in which the dynamics of a fermionic open system of interest is associated in a simple way to the dynamics of a similar system governed by different parameters.
What is fermionic duality? The idea of the duality mapping is particularly easy to describe for a closed quantum system evolving unitarily with a time-constant Hamilto-nian H. In this case the mapping explicitly constructs the adjoint Heisenberg evolution (superscript H) from the Schrödinger one by a substitution of physical parameters: U H (t) := e iH H t = U (t) † = e iHt = e −iHt =:Ū (t). (1) We will denote such a parameter mapping by an overbar. In the present simple example of a duality, the required relation between the Hamiltonian evolution generators, is achieved by inverting the signs of all energies H → −H: all local energies, all hopping amplitudes and all many-body interactions. To motivate this duality mapping consider the computation of the evolution of an arbitrary state, |ψ(t) = i |u i u i (t) u i |ψ(0) which requires both the right eigenvectors {|u i } of U (t) and its left eigenvectors { u i |}. Equivalently, one needs the right eigenvectors of U (t) and of U H (t) =Ū (t) where we consider the Heisenberg evolution "as" a Schrödinger evolution at different parameter values. The duality mapping (1) makes explicit that these two sets of eigenvectors are related in a simple way through their parameter dependence allowing unnecessary algebra to be bypassed.
Having outlined the key idea, we immediately observe that for closed systems with time-constant H this trick is completely pointless because there is an obvious shortcut: the eigenvectors of U (t) are time-constant and coincide with the eigenvectors of H = H † which are related by taking the Hermitian adjoint, Since [H, U (t)] = 0, time-dependence arises only through the eigenvalues u i (t) = e −ih i t where h i are the constant energy eigenvalues. Also, one may hesitate to work with the Hamiltonian H since it is clearly unphysical : inverting energies destabilizes any physical system which does not have an upper bound on its energy spectrum. Notably, for a fermionic system with a finite number of modes this latter objection is not really an issue since its spectrum is bounded by the Pauli exclusion principle.
However, when considering an open system with evolving density operator ρ(t) = Π(t)ρ(0) the above mentioned shortcut completely breaks down. Although it turns out that non-unitary open-system evolutions can still be generated time-locally [2][3][4][5] as ∂ t ρ(t) = −iG(t)ρ(t), new problems arise because the generator is a time-dependent, non-Hermitian superoperator G(t) = [G(t)] † even though the total system evolution is generated by a time-constant, Hermitian Hamiltonian operator. Physically, these new complications derive from memory (retardation) and dissipation effects, hallmarks of open-system dynamics. They cause the left and right eigenvectors of the generator G(t) to be distinct, timedependent and different from the eigenvectors of the evolution propagator Π(t) which is ultimately of interest: [G(t), Π(t)] = 0. This implies that for an open system the transformation between the Schrödinger and Heisenberg generators is highly nontrivial (Ref. [4], p. 125) unlike the relation (2) for the underlying closed total system. An alternative simple mapping between the Schrödinger and Heisenberg picture evolution would dramatically simplify time-evolution calculations by providing a link between the left and right vectors.
Transposing the simple duality mapping for fermions from a closed to an open system does not seem to be possible at first: it is unclear how to evaluate the average of the simple relation (1) or (2) over the reservoir degrees of freedom (partial trace), even when making specific microscopic model assumptions. This is in contrast to other closed-system duality mappings [6][7][8] which are distinct from the one considered here [9]. It is therefore remarkable that for a very large class of fermionic open systems there does exist a nontrivial and useful extension of the duality which applies to the reduced time-evolution superoperator Π(t). Anticipating its later detailed discussion [Eq. (23)], it provides an elegant formula for the adjoint Heisenberg evolution analogous to Eq. (1) [10]: Here P is a linear transformation involving the fermion parity. Its presence hints at fermion parity superselection-forbidding quantum superpositions of states with even and odd fermion parity-as one fundamental principle on which the duality (3) is based [10]. Γ is the lump sum of microscopic tunneling constants-known by inspecting the underlying model Hamiltonian-and the overbar again denotes a parameter mapping. This generalization of Eq. (1) is truly dissipative. For example, for a resonant level coupled to a reservoir the parameter mapping inverts not only the sign of the level energy ε and the electrochemical potential µ, but also the dissipative tunnel-rate constant Γ. The relation (3) was first derived in Ref. [10] without making weak coupling and / or "Markovian" assumptions, requiring the techniques of Refs. [11][12][13][14] to explicitly consider all orders of the expansion of Π(t) in the system-environment coupling. In the following we will denote this direct consideration of the exact propagator Π(t) as approach (i) to duality. Applied to weakly coupled but locally interacting open systems, the fermionic duality (3) has already provided several interesting insights and predictions [10,[15][16][17]. For example, the time-dependent response of a "kicked" quantum dot with repulsive Coulomb interaction was shown to exhibit effects of electron-attraction. This surprising effect can be nicely understood from the duality mapping which involves the inversion of the local interaction parameter as in Eq. (2). This explains pronounced effects in the measurable time-dependent heat current which is sensitive to interactions. The same formulas are very difficult to understand directly in terms of the real repulsive interaction, but are easily rationalized by electron-pairing induced by the attraction in the fictitious dual system defined by the duality mapping. More generally, the thermoelectric response of a quantum dot-although studied long ago-entails several features that turned out to have a very simple explanation in terms of an effective attractive model that is dual to the repulsive system of interest [15,16]. These conclusions hold even beyond linear response to electro-thermal biases where, e.g., Onsager relations no longer apply, and the effects can be understood by extending the weak-coupling fermionic duality beyond the wide-band limit [16]. In all these cases, the original system is analyzed by a dual system, an effective system with at worst unconventional properties. Conversely, it was also shown that the response of a physically attractive dual system can be understood better by exploiting its repulsive original system [17]. Thus in the weak-coupling limit the duality can also be used in reverse.
Extension to other approaches. So far, these applications were in fact all based on a different formulation of the duality which we will denote by approach (ii) in the following. It differs from Eq. (3) by relying on the time-nonlocal quantum master equation (QME) also called Nakajima-Zwanzig (NZ) [18,19] or time-convolution type QME. By introducing a memory kernel it anticipates the time-convolution structure of the higher-order systemreservoir coupling terms encountered in the microscopic derivation of the duality [10]. Whereas in the weak-coupling limit approach (ii) recovers various other types of quantum master equations, for the interesting regime of stronger coupling it differs in essential points. So far it has remained unclear what the implications of the fermionic duality are in general for other approaches to open-system dynamics. This problem is solved in the present paper: besides extending the propagator approach (i) and the memory kernel approach (ii), we establish the fermionic duality for three additional approaches which are fundamentally different and complementary as we now outline.
(iii) The Sudarshan-Kraus or measurement-operator approach [20][21][22]-ubiquitous in quantum-information theory-also directly addresses the time-evolution superoperator Π(t). However, it is an operational approach which decomposes the evolution into independent physical processes conditioned on possible outcomes of measurements performed on the system's environment. Theoretically, this has the distinct advantage that approximations formulated in terms of these operational building blocks automatically preserve the positivity of quantum states, also in the presence of initial entanglement with a reference system (complete positivity). From these measurement operators acting only on the system one can furthermore compute the evolution of its effective environment and quantify the exchange of information as illustrated in Ref. [23]. Barring special limits where simpler Lindblad equations [24,25] apply (Markovian semigroups), for most systems of interest the insights offered by this approach seem practically impossible to gain in other formulations. The same applies to the so-called Choi-Jamio lkowski state, which is closely related to the measurement-operator sum by a well-known isomorphism [26][27][28]. The microscopic calculation of Kraus operators has received recent attention [29][30][31][32][33] but remains very difficult, motivating our search for analytic simplifications.
(iv) The time-convolutionless (TCL) or time-local quantum master equation approach [2,3,34] has the advantage that it allows the Markovianity of the evolution to be scrutinized more conveniently through the time-local generator G(t) mentioned earlier. It is thus closely tied to the question of the divisibility of the dynamics [35,36]. In practice, time-local QMEs also arise naturally from the time-nonlocal QMEs of approach (ii) when consistently accounting for frequency dependence of the memory kernel in decay problems [37,38]recently generalized in Ref. [5]-or in adiabatic expansions for situations with external driving [9,[39][40][41]. The microscopic calculation of G(t) is, however, very challenging making additional analytic insights valuable [34,42,43].
(v) The closely related jump-operator approach decomposes the time-local generator G(t) into "quantum jumps" with intermittent renormalized Hamiltonian evolution occurring at infinitesimal time steps. Although this is similar to the measurement-operator approach (iii) it provides distinct insights by primarily making the conditions for divisibilityrather than complete positivity-explicit. Importantly, this approach is also at the basis of the successful stochastic simulation method for open-system dynamics [44][45][46][47][48][49][50][51] and includes the familiar Markovian Lindblad QMEs as a special case. In the present work the jump-operator approach is particularly interesting because it most explicitly generalizes the closed-system duality (2) discussed above.
In none of the approaches (iii)-(v) the implications of fermionic duality have been explored. Doing so is of particular interest since these methods are pivotal for the continued fruitful application of ideas from quantum information theory to open-system dynamics [4,50,[52][53][54][55][56]. One should note that although approaches (i)-(v) are exactly equivalent, they define completely different starting points for approximations and formal considerations. Thus, having a formulation of fermionic duality in hand for each case will enable attaining independent insights. This holds true even when applied to the simplest, explicitly solvable transport model of a strongly coupled resonant level as was recently highlighted in Ref. [23] and we will draw on this reference for illustration. Even though this model has been studied for decades [57] the fermionic duality relations presented here went unnoticed. Importantly, our results continue to hold for a large class of much more complicated models whose detailed discussion is however beyond the present scope.
Fermionic duality: Useful but unphysical? Before proceeding it is important to neutralize two potential points of confusion. An immediate worry is that the fermionic duality for open systems maps some physical parameters to unphysical values as noted above. In fact, for open systems the unphysical destabilization of the system by inverting the signs of all local energies discussed after Eq. (1) becomes more prominent. For one, the duality mapping even makes the system-reservoir coupling Hamiltonian anti -Hermitian, H T → iH T . Although no real physical parameters become imaginary, this does invert the sign of all dissipative decay rates. As mentioned, for a resonant level this means that the decay constant is inverted Γ → −Γ. This does not lead to divergent quantities since in the duality relation (3) the negative decay rates are explicitly compensated by an exponentially decaying prefactor e −Γt . Moreover, in the weak coupling limit close inspection reveals [10] that one can use fermionic duality to set up a relation between two dual physical systems which both have nonnegative decay constants but otherwise different physical parameters. This simplification facilitated the applications in the weak coupling limit cited earlier. In the present paper we will, however, focus on the general case of strong coupling where this simplification fails 1 and this non-physicality must be confronted.
We will show that the anti-Hermitian coupling Hamiltonian causes the reduced dynam-icsΠ(t) to violate complete positivity, giving a clear operational meaning to the vague notion of an "unphysical" system. This is important since it will allow us to identify which contributions to the evolution of the dual system are unavoidably unphysical, a question that cannot be answered directly using the original derivation of the duality in Ref. [10]. Instead, by leaving aside the derivation and only considering the duality relation (3) as such, this paper shows that the loss of complete positivity is associated with the fermionic parity transformation P, a key ingredient of the duality. Hence, unlike in the weak coupling limit, the dual evolution has no statistical meaning anymore and one can no longer refer to an effective, physical dual system which simulates the original system. Although this may sound disastrous at first, it will become clear that in none of the discussed approaches these unintuitive features of fermionic duality limit its practical usefulness. Since in each of these approaches the fundamental property of complete positivity is expressedif at all-by different constraints, a careful discussion what is unphysical about the dual equations will be a recurring side-theme. It will emerge that the general unphysicality of fermionic duality is instead of an artifact a key feature unveiling its unconventional insights as compared to ordinary symmetries, see Sec. 6.
To avoid confusion about the domain of applicability we note that the fermionic duality (3) is primarily important for analytical calculations where one obtains some quantities of interest as functions of physical parameters. By a simple substitution of parameters it allows one to bypass very complicated and nonintuitive algebra. The ultimate importance of fermionic duality lies therein that this simplification allows the analysis of physical effects [10,[15][16][17] to be pushed much further. Unlike doing algebra, solution by parameter substitution has the advantage that it preserves the compact form of an expression that has already been calculated, simplified and physically well-understood. For example, it makes explicit which quantities have a similar functional dependence: if one knows that some contribution is an exponential function of time then generically the dual contribution obtained by a parameter substitution is exponential as well. Loosely speaking, one can thus distinguish individual nontrivial contributions to the dynamics (exponential time dependence) from trivial (exponential) ones. The duality mapping also implies the concept of self-dual quantities: Despite being physical, such quantities are mapped onto themselves by the generally unphysical duality relation, and are thereby constrained in ways that are impossible to see with common physical dualities, symmetries or intuition.
Outline. The outline of the paper is as follows. In Sec. 2 we first consider the simplest, exactly solvable open system that exhibits fermionic duality beyond the weak-coupling Table 1: Duality relations for the approaches numbered (i)-(v) mentioned in the introduction. By ! = we indicate equalities that are valid only in the special case where the evolution commutes with its generator, [G(t), Π(t)] = 0, which includes the weakcoupling limit where Π(t) = e −iGt with constant G. Hat denotes the Laplace transform f (ω) = ∞ 0 dte iωt f (t). Bar denotes the duality mapping of parameters which effects H → −H, H T → iH T and µ r → −µ r . I is the identity superoperator.

Finite evolution approaches
Infinitesimal evolution approaches Superoperator approaches limit [16], the resonant level. This provides the simplest yet nontrivial illustration of the general results derived in the subsequent sections. In Sec. 3 we consider the fermionic duality in its most basic form (3) obtained in Ref. [10] as a mapping between the finite-time Schrödinger and Heisenberg superoperators. From this we derive a fermionic duality for the set of Kraus measurement operators using the Choi-Jamio lkowski state associated with the dynamics. In Sec. 4 we consider fermionic duality for the infinitesimal-time generators of the evolution, either via a time-local or time-nonlocal quantum master equation. Whereas the time-nonlocal formulation allows for a solution in the Laplace-frequency domain, the time-local formulation allows for a further decomposition into jump operators. This leads to some unexpected insights into the divisibility of the dynamics and its causal structure. Finally, in Sec. 5 we combine these approaches to gain deeper insight into a generally applicable nonperturbative semigroup approximation [5,23] and its correction by an "initial slip". Here we combine the fermionic duality with a recently found exact functional relation between the time-local generator G and the time-nonlocal memory kernel K [5]. In Sec. 6 we conclude and outline directions for follow-up work. Table 1 provides a guide to the paper by summarizing the duality relations for all discussed approaches. Throughout the paper we set = k B = 1.

Simple example: Fermionic resonant level
The general fermionic duality is best illustrated for the simple model of a resonant level with arbitrary tunnel coupling Γ to a single reservoir at temperature T and electrochemical potential µ, where d and b ω are fermionic annihilation operators. This gives only a very basic description of the time-dependent (dis)charging of quantum dot systems realized in a range of heterostructures including molecular junctions [58][59][60] and atomic impurities [61]. Because it ignores Coulomb interaction effects this model is exactly solvable for strong coupling Γ whose effects are primarily of interest here. Although the solution is known since Ref. [57], several interesting aspects of the dynamics of the density operator ρ(t) were overlooked until recently [23], including measurable effects of the breakdown of two different notions of Markovianity. In the present paper we complement the study [23] by pointing out further interesting properties of this model which are implied by fermionic duality and which, importantly, turn out to be more generally valid. For example, the spectral properties of the various time-evolution quantities of approaches (i), (ii) and (iv) were noted in Ref. [23] to display striking regularities which could not be rationalized by any symmetry of the resonant level model. The measurement-and jump-operators of approach (iii) and (v) were likewise found to exhibit a pronounced pattern and to obey an unexpected sum rule whose form depends only on one microscopic parameter and time. These formal regularities indicate that there is still more to be said about the physical properties of this model, for instance, its degree of (non-)Markovianity as introduced below. As we will see in Sec. 5 this ties in with the practical task of constructing Markovian semigroup approximations to the solution and their corrections. All these points-and more-will be explored step by step in this paper. Physical properties of the dynamics. To appreciate these implications of fermionic duality, we, however, first need to give a summary of the various analytical forms of the resonant level model dynamics for later reference. Further details on the latter are given in Ref. [23]. Explicitly, all nontrivial dependence on the level detuning ε − µ, temperature T and tunnel coupling Γ is captured by three related functions of time: These functions contain pronounced damped oscillatory contributions at low T . Fortunately their complicated precise form given in Ref. [23] is not needed here. We first review how the functions k, g and p encode basic physical properties of the dynamics which will be important later on. For fixed time t, the state-evolution map ρ(0) → ρ(t) = Π(t)ρ(0) has two fundamental properties, namely trace-preservation (TP), Tr Π(t)ρ(0) = Tr ρ(0) for all initial states ρ(0), and complete positivity (CP)-reviewed in App. A-which implies Π(t)ρ(0) ≥ 0 for all ρ(0) ≥ 0. Although TP imposes no restrictions on k, g and p, the CP property is fully encoded in the range of values that the function p(t) can take at any time: Π(t) is CP if and only if |p(t)| ≤ 1. This nontrivial constraint is indeed [23] satisfied by formula (7) for all times t and all physical values of the parameters ε, µ, T, Γ of the model, as it should. Depending on these parameters, the dynamics may be Markovian in two different ways with different observable consequences as we now explain.
(1) The dynamics may be divisible by itself, Π(t) = Π(t−t )Π(t ), for all t ≥ t , which is the case if and only if the function g(t) = constant for t > 0. Despite the simplicity of the model, this semigroup-divisibility always breaks down except at resonance, |ε − µ|/T → 0, where p(t) = 0 (which always holds in the limit of a hot reservoir, T → ∞), or, when the level is completely off-resonance, |ε − µ|/T → ∞, where p(t) = ±Θ(t) is a step function. Thus, the dynamics is virtually never Markovian in the semigroup sense and cannot be Whenever |g(t)| > 1 for some time t the dynamics is not CP-divisible. The black curve marks where the maximal value equals 1. Note that max t≥0 |g(t)| depends only on the ratios (ε − µ)/Γ and T /Γ, whereas g(t) like k(t) depends on all parameters separately.
described by a Lindblad quantum master equation. The breakdown of this property is witnessed by an anomalous transient enhanced level occupation when decaying to a more depleted stationary state [23].
(2) In a less strict sense, the dynamics may still be divisible as Π(t) = Π(t, t )Π(t ) for all t ≥ t by another physical evolution Π(t, t ), a CP-TP map [35,36]. This CPdivisibility turns out to occur if and only if |g(t)| ≤ 1 for all t. This nontrivial condition is mapped out in Fig. 1 as function of level position and temperature relative to the coupling energy. One sees that the dynamics fails to be Markovian in the sense of CP-divisibility whenever the level is off-resonant by more than the tunneling and thermal broadening, |ε − µ| max{Γ, T }. For the resonant level this distinct property can be observed in transport by checking whether there is no reversal of the measured current as function of time for any initial level occupation [23].
Evolution. Having outlined some of the physics of this model, we now describe how the functions (5)-(7) explicitly determine the structure of the exact dynamics of ρ(t). The evolution can be written in three different ways and features the function p(t) [Eq. (7)]: By [H, •] we denote the commutator of the system Hamiltonian 2 H = εd † d with argument •, and the dissipator D η will be defined in Eq. (11). Also, we write A • := Tr A † • and 2 Throughout the paper we consider the action of the map Π(t) on arbitrary initial states ρ(0) since this enables the techniques of Refs. [11,13,14] which lead to the neat exponential form (8), see Ref. [23]. Also, this allows us later on [Eq. (61)] in the jump-operator approach (v) to generalize Eq. (2) of the introduction. If one restricts the action of the map Π(t) to operators ρ(0) which are fermionic states commuting with the parity, [ρ(0), (−1) N ] = 0 (superselection), then the contribution of the system Hamiltonian in Eq. (8) is not relevant in this model. For this restricted map one can also find a simpler set of measurement operators. B := B for operators A, B. Throughout the paper we will label each right eigenvector by its eigenvalue and the corresponding left eigenvector is distinguished by an additional prime as in Eq. (9).
The exponential form (8) is particular to this simple model. Due to the nontrivial dependence of p(t) on time [Eq. (6)] it is not the exponential solution of some Lindblad equation, despite the appearance of the familiar dissipator superoperators where we defined d + := d † , d − := d. The more general spectral decomposition (9) is natural to approach (i). The eigenvalues and their distinct left and right eigenvectors for this model are listed in Table 2. Finally, the Kraus operator sum (10) of approach (iii) always exists and in the present case the measurement operators read with nonnegative coefficients and the shorthand Quantum master equations. The above dynamics is the solution of the exact timenonlocal QME of approach (ii), whose memory kernel features the function k(t) [Eq.
Note that we included the system Hamiltonian H into K and used the normalization t 0 ds δ(s) = 1. Finally, the dynamics is also the solution of an exact time-local QME that defines approach (iv), with a generator that features the function The eigenvalues and eigenvectors in Eq. (18b) are listed in Table 2. Although Eq. (16) and Eq. (18a) look similar to the exponent of Eq. (8), they involve the three very different functions (5)-(7). Table 2: Time-dependent eigenvectors and eigenvalues of the evolution Π(t) and its time-local generator G(t) derived in Ref. [23]. Observe that the eigenvalues are trivial exponentially decaying functions and constants, respectively, fixed by the T = ∞ semigroup limit of the model. All nontrivial time-dependence is due to finite-T effects [23] and enters the dynamics through the functions p(t) and g(t) which appear only in the eigenvectors of Π(t) and G(t), respectively. Note that here the normalization of the eigenvectors of G(t) differs from the normalization fixed in Eqs. (50b) and (50c). The duality relation for the coherences can be seen from Having summarized the exact equations for this model to be discussed, we note that in the weak coupling limit it is not difficult to reveal a simple structure. For example, in the time-local QME (17) one can simplify 1 − ηg(t) ≈ f (η(ε − µ)/T ) where f is the Fermi function and one then directly derives a fermionic duality relation using the formal replacement f ((ε − µ)/T ) → 1 − f ((ε − µ)/T ) as explained in Ref. [16]. However, this is not possible in the case of arbitrarily strong coupling Γ considered here. Thus, despite the simplicity of this solvable model none of the above representations of its exact dynamics seems to exhibit an obvious general structure. In the following we will derive such a structure and illustrate it for each of the above expressions.

Fermionic duality for exact time-evolution
We now extend the scope to the much broader class of models of the form H tot = H +H R + H T where only the following assumptions are made: (I) The multiple fermionic reservoirs described by H R are noninteracting with structureless, infinitely wide bands, each one being separately in equilibrium at the initial time. (II) The coupling to the fermions in the system (indexed by l) is bilinear in the field operators, H T = rl dω t rl d † l c rω + h.c., and independent of the energy ω of the fermionic modes in the reservoirs (indexed further by r). (III) The system Hamiltonian H obeys parity superselection, [H, (−1) N ] = 0, and as a result so does the total system. The only microscopic quantity that explicitly plays a role in the fermionic duality is the lump sum of tunnel-coupling constants over the system and reservoir indices: Here l (r) includes all relevant quantum numbers (spin, orbital moment, etc.) on the system (reservoir) which need not be conserved by H T , unlike the fermion number. Based on these three assumptions the duality was established in Ref. [10]. However, the detailed derivation given there does not lead to the insights reported in the present paper. The conditions (I)-(III) do not help to understand the results obtained here by starting from the established duality relation Eq. (3). We refer to Refs. [10,16] for further discussion of these assumptions and the derivation and to Refs. [10,[13][14][15][16][17]62] for numerous detailed illustrations of how the duality can be technically applied and physically exploited in the weak coupling limit.
No other assumptions are necessary: in particular, the system, described by H, may consist of any finite number of levels with any type of multi-particle interaction of arbitrary strength, including superconducting pairing terms that break particle conservation but preserve parity. Also, the magnitude of the couplings, temperatures and electrochemical biases can be arbitrary assuming that the employed perturbation series converges. Thus, the following results apply to a very large class of actively studied models which are relevant to nonequilibrium quantum-impurity physics, quantum transport and open-system dynamics. We also note that for weak coupling, the duality relation can be generalized beyond the case of structureless wide bands [16].
Of central interest is the superoperator Π(t) describing the state evolution, i.e., the Schrödinger propagator, It is obtained by tracing out the fermionic reservoirs, assuming that each of these is initially uncorrelated with the system and separately in an equilibrium state. The propagator is thus a function of the parameters specifying the system Hamiltonian H, the coupling H T , and the different electrochemical potentials of the reservoirs, collected in µ = (µ r ). This dependence is important in the following and will be denoted by Π(t, H, µ, H T ) when required. The dependence on the different reservoir temperatures T r need not be indicated. The superoperator Π(t) † describes the time-evolution of system observables A, i.e., the Heisenberg propagator, for expectation values. Here the superadjoint of a superoperator, indicated by bold †, is defined by Tr A † (ΠB) = Tr (Π † A) † B and is of central importance in this paper. It is defined relative to the Hilbert-Schmidt scalar product between operators (A|B) = Tr A † B and therefore distinct from the ordinary adjoint † of an operator relative to the scalar product between vectors ψ|Aφ = A † ψ|φ . For superoperators with the special form • → (L • R) of a left and right multiplication by operators L and R, respectively, the two distinct adjoint operations are related in a simple way: In the following these distinctions will be clear in the context and we will talk about adjoints, eigenvectors, and orthogonality without further specification ("super"). Since generally the evolution Π(t) is a not represented by a normal matrix, [Π(t) † , Π(t)] = 0, its left and right eigenvectors are not simply related by taking the adjoint †. As both sets of vectors are required in the analysis of dynamics, this presents a crucial complicating factor in any (semi-)analytical treatment of open quantum systems. This is what fermionic duality addresses.

Evolution superoperator
The fermionic duality establishes a relation between Π H (t) = Π(t) † and Π(t) evaluated at different parameter values which is denoted byΠ(t). By first explicitly evaluating the wide-band limit, this relation can be derived within a renormalized perturbation expansion of all finite-T corrections of the propagator Π(t) around the T → ∞ limit [10,13,14]. For the considered class of models the propagator obeys the fermionic duality relation order-by-order. Here Γ is the lump sum of couplings (19) and the superoperator denotes the left multiplication with the system parity operator (−1) N := exp(iπN ). By the overbar we denote the following parameter substitution of some function X: For example, for the resonant level model of Sec. 2 this parameter mapping corresponds to (ε, µ, Γ) → (−ε, −µ, −Γ) which transforms the functions encoding all nontrivial parameter dependence as follows 3 :k The fermionic duality (23) expresses an exact restriction on the possible parameter dependence of Π(t) based only on the quite generic physical assumptions (I)-(III) mentioned at the beginning of Sec. 3 and two fundamental physical principles, the Pauli exclusion (anticommutation relations) and fermion-parity superselection applied to the total system. One may think ofΠ(t) as a continuation of Π(t)-considered as function of microscopic parameters-from a physical domain to a larger domain of unphysical values. This is not uncommon in physics, cf. for example, the complexification of angular momentum in scattering theory (Regge theory). In the present case, the system-reservoir coupling Hamiltonian is mapped to anti-Hermitian values, H T → iH T . This corresponds 4 to a Wick-rotation t rlω → it rlω together with inversion of the relative sign between the two tunneling terms in H T . The fermionic duality is the imprint left behind in the reduced description (after tracing out reservoirs) of the mentioned physical assumptions and principles (before tracing). It takes the form of a restriction (23) on the continuation beyond the physical parameter domain. It is not required-or to be expected-that the superop-eratorΠ(t) resulting from the parameter substitution (25) should be a physical evolution. The construction as a continuation guarantees thatΠ(t) is still a TP map, but we will see that it is not CP. Nevertheless, approximations that break fermionic duality are inconsistent with the physical assumptions and principles governing the underlying total system, see Sec. 6. After presenting all our results we will compare with other works in the discussion [Sec. 6]. 3 Relation (27), written as (1 − e −Γt )p(t) = g(t) + e −Γtḡ (t), is obtained by inserting (7) on the left and partially integrating using p(0) = 0. Relation (28) follows by taking the overbar of Eq. (27). 4 The substitution HT → iHT means that we treat the conjugate pair of tunnel constants in HT = rl dω t rl d † l crω + t * rl c † rω d l as independent parameters: t rl → it rl but t * rl → it * rl = (−it rl ) * . This inverts the sign of all spectral densities t rl t * r l → −t rl t * r l determining the decay rates, see Ref. [10].

Cross-relation left and right eigenvectors
We now first explain the usefulness of relation (23), extending the analysis of Ref. [10]. It implies that if π i (t) is a right eigenvector of Π(t) with eigenvalue π i (t) then is also an-in general different-eigenvalue, numbered j = i, with left eigenvector Similarly, right eigenvectors are related to left ones by Thus, although Π(t) is not a unitary matrix [cf. Eq. (1)] its left and right eigenvectors are nevertheless related by conjugation up to parity signs (P) and a parameter substitution (25) (overbar). The duality only ensures proportionality of the vectors in Eqs. (29b)-(29c). The proportionality constants were chosen such that binormalization imposed for pair i is preserved for pair j: (π j (t)|π j (t)) = (π i (t)|π i (t)) * = 1. One is then still free to gauge the right hand side of Eq. (29c) by any nonzero time-dependent complex scalar θ j (t) and correspondingly Eq. (29b) by 1/θ j (t). If an eigenvalue happens to be self-dual, π i (t) = e −Γtπ i (t) * , we have i = j in Eq. (29a). In this case the gauge freedom is fixed by binormalization (π i (t)|π i (t)) = 1: with related factors π i (t) P π i (t) · π i (t) P π i (t) = 1. Table 2 shows that for the resonant level model all eigenvalues are indeed cross-related by the duality relation (29a). The nontrivial, non-exponential time-dependence is located in the eigenvectors. The duality relation (29b) now dictates that if the right eigenvector to eigenvalue π 0 (t) = 1 depends nontrivially on time through p(t), then the same must hold for the left eigenvector to eigenvalue π 3 (t) = e −Γt , see Table 2 and Eq. (28). Analogously, the time-constancy of the left and right eigenvectors i = 1 dictates the time-constancy of the i = 2 eigenvectors. Thus, duality provides a fine-grained insight into the location of nontrivial contributions to the dynamics.
In an analytical calculation of the spectrum of Π(t) one may, for example, determine for each dual pair only one eigenvalue and its left and right eigenvector algebraically, and then obtain the remaining eigenvalues and eigenvectors via a mere parameter substitution [Eq. (25)] and parity transform P [Eq. (24)]. This is much simpler and, moreover, preserves the compactness of analytical expressions already obtained. For models only slightly more complicated than the resonant level this already leads to significant simplifications and some surprising insights as shown in the weak coupling limit [10,[15][16][17].

Constraints on evolution of states and observables
We have seen for the resonant level that the duality (29) dictates that terms with qualitatively similar time-dependence in the spectral decomposition of Π(t) occur pairwise on opposite ends of the real part of the eigenspectrum. In the general dynamics, one pair of contributions is of particular interest. The right eigenvector to eigenvalue π 0 = 1 is a time-dependent fixed point 5 , Π(t) π 0 (t) = π 0 (t) , which is guaranteed to exist by the evolution's TP property, π 0 Π(t) = π 0 writing π 0 = Tr. Often the operator π 0 (t) is unique and can then be scaled to a positive, trace-normalized physical state 6 . For simplicity we assume throughout the paper that the eigenvalue π 0 is nondegenerate. The time-dependence of the fixed-point is important and its significance was recently highlighted [23]: If one initially prepares ρ(0) = π 0 (t r ) where the reentrance time t r > 0 is a parameter, then the nontrivial evolution is guaranteed to exactly recover this state at the preset time t = t r , Π(t r ) ρ(0) = π 0 (t r ) , even though the environment state for t = 0 generally differs from the one at t = t r . For the resonant level model this reentrant behavior signals the breakdown of semigroup-Markovianity [23].
The duality cross-relation (29a) now dictates that the dynamics has another fundamental eigenvalue π n (t) = e −Γt with trivial time-dependence at the opposite end of the spectrum. Here we number i = 0, . . . , n where n := d 2 −1 and d is the system Hilbert space dimension. The right eigenvector π n = (−1) N is the time-constant parity operator, We note that this follows directly from the fact that the dual propagatorΠ(t) is also a TP map 7 , 1 Π (t) = 1 . The corresponding left eigenvector can be expressed via the zeroth right eigenvector, π n (t) = Tr{π 0 (t)(−1) N •}, whereπ 0 (t) denotes the self-adjoint operator specifying π 0 (t) . It determines the amplitude in the expansion of the timedependent state: Thus, the nontrivial time-dependence of the coefficient of the fast Γ-decay is also determined by the time-dependent non-decaying fixed-point component π 0 (t) , namely through its functional dependence on parameters. For semigroup dynamics this coefficient is timeconstant, but in general it is time-dependent, even in the resonant level model, see π 0 (t) in Table 2. The result (33) implies that the expectation value of a system observable A can be decomposed into an instantaneous expectation value in the time-dependent fixed-point state plus corrections: The corrections with the fast Γ-decay appear only for observables which overlap with the fermion-parity, Tr{A(−1) N } = 0. Such operators A depend multiplicatively on the occupations of all fermionic orbitals in the open system, i.e., they probe global correlations within the system. The above insight into the general dynamics extends the weak-coupling results of Ref. [16].

Measurement operator sum
We now turn to an entirely different formulation of the same dynamics which is ubiquitous in quantum information theory. We can apply this approach here since we are assured that Π(t)-being the exact evolution-is a CP map 8 . It can therefore be written in the form of a Sudarshan-Kraus operator sum [20][21][22] Without loss of generality we choose to normalize the measurement operators using the Hilbert-Schmidt scalar product, Tr M α (t) † M α (t) = 1. The coefficients m α (t) are then real and positive by CP, see App. A, and the TP property of Π(t) is equivalent to By taking the trace this implies a scalar sum rule: the coefficients must sum to the Hilbert Each term in the operator sum (35) describes a physical process in which outcome α is obtained by a measurement on the environment R in some basis. For each different choice of a basis, there is a set of measurement operators {M α (t)} and thus a different operator-sum representation. We fix this freedom by considering canonical measurement operators which are orthonormal, Tr M α (t) † M α (t) = δ αα . If the m α (t) are nondegenerate, this fixes the set {M α (t)} uniquely up to trivial changes by phase factors which cancel out term-by-term in the sum (35), see App. A for the case of degeneracy. Importantly, for fermionic systems the operators must have a definite parity denoted by (−1) Nα , i.e., since the operators describe measurements 9 .

Cross-relation of Heisenberg and Schrödinger measurement operators
From the operator sum (35) it is easy to find the measurement operators for the Heisenberg evolution Π H (t) = Π(t) † by using Eq. (22), To see the nontrivial implication of fermionic duality (23), we insert Eq. (35) and Eq. (38) and show that the individual terms in the two operator sums must be equal up to a permutation of the summation index α. This follows most elegantly by the Choi-Jamio lkowski (CJ) correspondence for which the fermionic duality is worked out in App. A. We obtain the key result that pairs of orthonormal measurement operators with the same parity obey and their corresponding coefficients fulfill where α = α is allowed. In Eq. (39a) the only freedom left in the relation between the operators M α and M α is a complex phase factor, which we set to 1. The fermionic duality relation (39) implies that if a coefficient is self-dual, m α = e −Γt (−1) Nαm α , the measurement operator is a strongly constrained function: its adjoint must correspond to dual parameters, M † α =M α . In all other cases, for each pair α, α of dual coefficients one needs to determine only one of the measurement operators, obtaining its dual operator for free. Thus, very similar to the relation (29) between left and right eigenvectors of Π(t), the difficult task of analytically finding the measurement operators and coefficients for nontrivial fermionic open systems is significantly simplified.

Additional fermionic sum rule for measurement operators
Since the dual propagatorΠ(t) = αm αMα (t) •M α (t) † is also a TP map, the dual measurement operators also obey a sum rule: αm αMα (t) †M α (t) = 1. Notably, this is not an obvious consequence of the TP sum rule (36) for Π(t): inserting Eq. (39) we instead find 10 that the original measurement operators of the fermionic systems must obey an additional, independent sum rule: This shares with Eq. (32) the remarkable feature of depending only on a single detail of the microscopic model, the lump sum of couplings Γ, independent of interactions and external controls such as temperature, and chemical potentials. Unlike the familiar sum rule (36), the adjoint appears on the right operator and the difference of even and odd parity terms is constrained to a time-dependent operator. The trace of Eq. (40) implies an extra scalar sum rule where we used that Tr M α (t) † M α (t) = δ αα for canonical measurement operators. Together with Eq. (37) we obtain separate sum rules for the coefficients of the even and odd measurement-operators as functions of time: Whereas for t = 0 the even operators must contain all the weight to produce Π(0) = 1 • 1, the even and odd weights coincide in the stationary limit t → ∞, evenly splitting the standard sum rule (37). In other words, the stationary evolution gives equal weight to parity changing and parity preserving processes.
Correspondingly, the even-parity operators (12) are self-dual and the odd ones are dual partners, The self-duality strongly constraints the coefficients of d † d and dd † inside the operator M † 0η : the substitution (ε, µ, Γ) → (−ε, −µ, −Γ) maps the coefficients to their complex conjugates. We stress that without the unphysical inversion of the decay rates Γ → −Γ one cannot explain this puzzling "symmetry" of this exact result. Furthermore, it is by no means obvious from the explicit solutions (12)-(13) that the additional simple sum rule (40) indeed holds. Also, the scalar sum rule (42) implies for fixed level occupation N = 0 or 1 that any nontrivial time-dependence of the coefficients (13) for η = ± must be the same up to a sign. All these structural features of the measurement operators were left unexplained in Ref. [23].

Unphysicality of the duality mapping
Of all the approaches to be discussed, the measurement-operator formulation (39) most clearly reveals that the dual propagatorΠ(t) is unphysical. It is not CP whenever Π(t) and Π H (t) are CP: the fermion-parity signs in Eq. (39b) imply that for operatorsM α with odd-parity the coefficientsm α are strictly negative. This means that due to the inversion of coupling constants, Γ → −Γ,Π(t) cannot correspond to the evolution of any physical system. In the duality relation (23) this is reflected in the parity transformation P.
We stress that this unphysicality in no way obstructs the derivation of Eq. (23) or its useful application to physical problems. On the contrary, it makes the duality mapping particularly interesting: by continuation of parameters to non-physical domains [Eq. (28) ff.] it points out functional dependencies which are not just physically "unintuitive" but even impossible to motivate by strictly physical parameter mappings.

Fermionic duality for exact quantum master equations
We now consider how fermionic duality constrains equivalent exact quantum master equations which generate the evolution Π(t).

Time-local quantum master equation
The dynamics can be described by a time-local QME [5] d dt Importantly, the time-local generator G(t) is in general time-dependent even though it derives microscopically from a time-constant Hamiltonian generator H tot for system plus reservoirs. The generator of the corresponding Heisenberg evolution acts from the left, d dt in order to generate the evolution of observables (21) as d dt A(t) = iG H (t)A(t). As a consequence, it is not 11 simply equal to the adjoint of the generator: 11 The adjoint equation d dt Π(t) † := iΠ(t) † G(t) † suggests to identify G(t) † with the generator. However, since it acts on the right one verifies that it is not the generator in the equation of motion for an observable This difference implies that for open systems one cannot switch from the equation of motion in the Schrödinger picture, d dt ρ(t) = −iG(t)ρ(t), to the Heisenberg picture, d dt A(t) = iG H (t)A(t), without first solving the dynamics. This is a known complication (Ref. [4], p. 125) of the analysis of open-system evolutions not commuting with their generator [65].
This nontrivial problem-specific to open systems-is solved by fermionic duality and will play a role in the construction of approximations in Sec. 5. Only in the simple cases where [G(t), Π(t)] = 0 do we have G H (t) = G(t) † . This includes the familiar case of Markovian semigroup dynamics where a time-constant G generates Π(t) = e −iGt . However, already for the resonant level model we have G H (t) = G(t) † since time-ordering of the generator matters except for special parameters (T = ∞ or ε = µ, see Sec. 2).
The TP property of the Schrödinger evolution, Tr G(t) = 0, by Eq. (47), corresponds to the Heisenberg evolution being unit-preserving or unital, Physically this means that trivial measurements stay trivial.
Taking the time-derivative of relation (23) one obtains the fermionic duality for the time-local generator: This relation is another key result of the paper which we again stress is exact, in particular, it is not based on any time-local approximation. It solves the nontrivial task of obtaining the Heisenberg generator G H (t) directly from G(t), without computing Π(t). Already for the resonant level model this presents a significant simplification: instead of performing quite some superoperator algebra 12 as required by Eq. (47), we obtain G H (t) by the simple parameter substitution g(t) →ḡ(t) given in Eq. (27).

Cross-relation left and right eigenvectors of the generator
The time-local fermionic duality (49) immediately implies that if g i (t) is a right eigenvector of G(t) with eigenvalue g i (t) then is also a-generally different-eigenvalue g j with left eigenvector Similarly, for right eigenvectors: As before [Eq. (29) ff.], the proportionality factors were chosen to ensure that the binormalization of pair i is passed on to pair j: (g j (t)|g j (t)) = (ḡ i (t|ḡ i (t)) * = 1. For self-dual eigenvalues g i (t) = [iΓ −ḡ i (t)] * biorthonormality (g i (t)|g i (t)) = 1 implies where ḡ i (t) PΠ(t) −1 g i (t) · g i (t) Π(t)P ḡ i (t) = 1. It is expected that fermionic duality takes a more complicated form here since the generator G(t) incorporates a great deal of the complexity of the solution Π(t) into the QME (45) in order to eliminate the memory integral of QME (15). In this respect, the simplicity of the eigenvalue duality (50a) is surprising and presents a definite advantage for analytical calculations. It generalizes the relations of Refs. [10,16] which for weak coupling imply that Γ is always the largest decay rate [16]. As expected the cross-relation (50b)-(50c) of the eigenvectors is more complicated due to the involvement of Π(t). Only the special case [G(t), Π(t)] = 0 leads to a simpler duality relation G(t) † = iΓ I − PḠ(t)P that directly relates left and right eigenvectors of G(t): g j (t) = ḡ i (t) P and g j (t) = P ḡ i (t) . This includes the Markovian-semigroup limit with time-constant G, recovering the weak-coupling results of Ref. [16]. Table 2 shows that for the resonant level model, the eigenvalues of G(t) indeed satisfy the cross-relation (50a). Yet, since [G(t), Π(t)] = 0 for this model, the eigenvector relations (50b)-(50c) remain nontrivial: their verification requires the transformation (27) of the functionḡ(t) and some algebra to verify Eq. (49). We note that G(t) and its eigenvectors also satisfy another, simpler relation which is, however, specific to the model and not related to general principles, see App. B.

Constraints on time derivatives of states and observables
Analogous to the fermionic duality (31) for the propagator, its time-local version (49) provides general insight into where nontrivial (non-exponential) contributions occur in the dynamics. In this case it concerns the time-derivative of the state: The prime indicates that we sum over pairs of dual eigenvalues i and j keeping only one term for self-dual ones. Here there is a catch because the evaluation of Eq. (52c) requires that the normalization of g i (t) is known, which we implicitly fixed in the duality relation Eq. (50). This is not an issue for two important contributions which we now discuss. The first one is the missing contribution: the time-dependent zero-mode of the generator, G(t) g 0 (t) = 0. Such a right eigenvector with eigenvalue g 0 (t) = 0 always exists since by trace preservation g 0 G(t) = 0, writing g 0 = Tr = π 0 . At finite times g 0 (t) is distinct from the time-dependent fixed point π 0 (t) of Π(t) [Eq. (31)], even though asymptotically both converge to the stationary state ρ(∞) = g 0 (∞) = π 0 (∞) whenever it is unique 13 . In fact, for the resonant level g 0 (t) even fails to be a positive operator in time intervals where |g(t)| > 1 [ Table 2], which happens precisely in parameter regimes where the dynamics is not CP-divisible shown Fig. 1. In contrast, π 0 (t) is always positive since |p(t)| ≤ 1 for all t, see Sec. 2.
The fermionic duality (50a) implies that there is another fundamental contribution to Eq. (52c) with eigenvalue g n (t) = −iΓ at the other end of the spectrum. Remarkably, it only depends on Γ and its right eigenvector does not depend on any microscopic detail: This follows from g n (t) = e Γt/2 Π(t) (−1) N = e −Γt/2 (−1) N [Eq. (32)]. Note that this also follows from the fact thatḠ(t) generates a TP map, 1 Ḡ (t) = 0. The corresponding left eigenvector is g n (t) = e −Γt/2 Tr{ḡ 0 (t)(−1) N Π(t) −1 •} [Eq. (29b)] whereḡ 0 (t) denotes the self-adjoint operator specifying ḡ 0 (t) . It inherits the parameter dependence of the zero-mode which in general is complicated, For the time-derivative of the expectation value of an observable A the nontrivial timedependent prefactor of the Γ-decay rate, is completely determined by the parameter dependence of the zero-mode of the generator, ḡ 0 , the missing term in the expansion (52c). This remarkable structure is a generalization of the weak-coupling result of Ref. [16] which introduces a new distinction: whereas the fixed-point of Π(t) determines this fast contribution to ρ(t) and expectation values A (t) [Eq. (34)], it is the distinct zero-mode of G(t) that determines the fast contribution to d dt ρ(t) or currents d dt A (t) [Eq. (55)]. Whereas for [G(t), Π(t)] = 0 (including Markovian semigroups) the fixed point and zero mode coincide at any time, they are in general different, π 0 (t) = ḡ 0 (t) . This illustrates that fermionic duality leads to independent insights when formulated in complementary approaches.
For the simple resonant level model this leads to an interesting insight into the transport current by taking A = N , the level occupation operator. Using Eq. (52c) one can verify that the omitted terms in Eq. (55) are zero because (N |g j (t)) = 0 except for j = 3. The normalization of the eigenvector g 3 (t) can then be calculated 14 by taking the trace of Eq. (50c) and we obtain from Eq. (55) Thus, the observable transport current is automatically decomposed in two contributions of which one is a trivial exponential decay depending on the initial state through 1 2 ((−1) N |ρ(0)) = 1 2 − N (0). All nontrivial time-dependence is captured by the single function g(t) from the generator G(t) but evaluated at dual parameters. Note that this relation does not follow from Eq. (34) unless one laboriously uses identities connecting the nontrivial functions g and p.

Jump operator sum
As mentioned in the introduction, a distinct advantage of the previous time-local QME approach is that it connects to the divisibility properties of the dynamics (Markovianity). These are, however, only revealed when decomposing the generators G and G H (t) [Eq. (47)] appearing in the time-local fermionic duality (49) into jump-operator sums, analogous to the decomposition of Π(t) into a measurement-operator-sum.

Causal and anti-causal divisibility
This approach requires some preliminary discussion. We first note that the Hermicity-and trace-preservation properties of the dynamics alone already imply the following structure of the generator [App. C] due to Lindblad, Gorini, Kossakowski and Sudarshan [24,25] The • contain jump operators J α (t) and are weighted with real coefficients j α (t) which we assume to be nondegenerate, see App. C for the degenerate case. Their structure guarantees that the generated dynamics is TP (Tr D α (t) = 0). The coefficients j α (t) need not be positive, in contrast to the coefficients of measurement operators [Eq. (37)]. The effective Hamiltonian H(t) is Hermitian but differs from the bare one, H, which we indicate by the time argument.
Similar to the measurement operators [Sec. 3.2], we eliminate gauge freedom by working with canonical jump operators, which are orthonormal, both mutually Tr J α (t) † J β (t) = δ αβ and to the identity, Tr J α (t) = 0. Importantly, the canonical jump operators J α (t) have a definite parity (−1) Nα , i.e. (−1) N J α (t)(−1) N = (−1) Nα J α (t) and the canonical effective Hamiltonian has even parity, The dynamics is CP divisible, Π(t) = Π(t, t )Π(t ) for all 0 ≤ t ≤ t ≤ ∞, if and only if the condition 15 j α (t) ≥ 0 holds for all α and t ≥ 0 [35,36]. In this case the jump operators have an operational meaning: J α (t) is a measurement operator for outcome α measured on the environment with infinitesimal probability j α (t)δt ≥ 0 during infinitesimal evolution. Divisibility also has a clear operational meaning in terms of a simulation task. The condition states that the full evolution up to time t can be simulated by stopping the evolution earlier at t -decoupling and discarding the environment-and then applying to the output some postprocessing device described by Π(t, t ). Such a physical device exists if and only if the latter is a CP map, see Sec. III of Ref. [23] for a discussion. If such a simulation is possible for every t and every t, then Π(t) is called CP divisible. This indicates that the input-output correlations of the dynamics are weak. For this purpose we only need to inquire into the possibility of such a simulation, not its implementation.
To derive a fermionic duality for jump coefficients and operators, we need to decompose the Heisenberg generator G H (t) appearing in Eq. (49) in a similar way, with Hermitian H H (t). The different structure of the Heisenberg dissipator, • , now ensures that the Heisenberg evolution is unitpreserving [Eq. (48)]. Moreover, the coefficients j H α (t) are distinct from the j α (t) and related to a different type of divisibility: j H α (t) ≥ 0 is the condition 16 for what can be called anti-causal CP divisibility of the state dynamics, Π(t) = Π(t )Π a (t, t ) for all t ∈ [0, t] by some CP-TP map Π a (t, t ) on the right, in contrast to the usual division of the dynamics 15 If G(t) satisfies this condition of CP-divisibility, it implies that Π(t) = Π(t, 0) is CP. Note that if G(t) does not satisfy this condition, it is not known which sufficient conditions the Jα(t) and jα(t) should satisfy to ensure that Π(t) is CP. 16 If G H (t) satisfies this condition, then it implies that Π(t) = Πa(t, 0) is CP. (b) Continuation of max t≥0 |g(t)| from physical (Γ > 0) to dual parameters (Γ < 0) shows the connection of causal and anti-causal divisibility revealed by fermionic duality. In (a) this connection is hidden in the T /Γ → ∞ limit. At resonance (ε = µ) the evolution is a semigroup which is trivially both causally and anti-causally divisible. As one tunes further away from resonance, the evolution first looses the anti-causal and then the causal divisibility. For T < Γ/(2π) (white regions) max t≥0 |ḡ(t)| = ∞, reflecting that the stationary limit of G H (t) does not exist even though the limit Π H (∞) is well defined. This is a peculiarity of the time-local description. by postprocessing to the left. Whereas semigroup dynamics is both causally and anticausally CP divisible, this does not hold for more general dynamics as studied here. For the resonant level model the parameter regime of anti-causal divisibility is mapped out in Fig. 2(a) and does not coincide with the regimes of causal divisibility.
The operational meaning of anti-causal divisibility becomes clear when viewed as a simulation task: The condition states that the full evolution Π(t) up to time t can be simulated by preprocessing its input by some device described by Π a (t, t ), and then afterwards running the evolution Π(t ) only up to time t . Also here, a physical preprocessing device exists if and only if Π a (t, t ) is a CP map. If such a simulation is possible for every t and every t, the evolution can be called anti-causally CP divisible. This indicates that the input-output correlations of the dynamics are weak and additionally that the causal ordering is weak, i.e., the dynamics is robust against interruption at t and reversal of causal ordering. As for causal divisibility, we only inquire into the possibility of such a simulation, not its implementation.
These two types of divisibility are not related in an obvious way. It is in general possible to express H H (t), J H α (t) and j H α (t) in H(t), J α (t) and j α (t) using the measurement operators M α (t) and m α (t) of the solution of the dynamics Π(t). However, just like the relation between G(t) and G H (t), this relation is highly nontrivial whenever [G(t), Π(t)] = 0. As a result, anti-causal divisibility cannot be easily related to causal divisibility. The fermionic duality for jump operators solves this nontrivial problem by relating the two jump operator sums, as we will explain below.

Fermionic sum rule for jump operators
To derive the duality relation for the jump operator sum, we first note a special implication of Eq. (49), the exact fermion-parity zero mode of G(t), Eq. (53). This is equivalent to a fundamental sum rule for the jump operators: This is remarkable since in general the jump operators are not constrained by any sum rule independent of model details, and here the only such detail is the lump sum Γ. Taking the trace, we find that the time-dependent coefficients j α (t) of the odd-parity jump operators sum to a constant, leaving the even parity jump coefficients unrestricted. Although Eqs. (59) and (60) are clearly analogous to the additional sum rules (40)- (41) and originate from the same fermionic duality relation they are not simple consequences of each other. In the resonant level model there are only odd-parity jump operators [Eqs. (11) and (18a)] and the fermionic sum rule for jump operators (59) , which in this case is a multiple of the scalar sum rule (60), η j η (t) = Γ.

Cross-relations between Heisenberg and Schrödinger jump operators
Inserting Eq. (57) into Eq. (49) and using the fermionic sum rule (59) we obtain 17 G H (t) in the form (58) where the effective Heisenberg Hamiltonian equals minus the Schrödinger one evaluated at dual parameters, (61a) We thus explicitly recover the closed-system fermionic duality (2) extended nontrivially by the inclusion of the time-dependent renormalization by the environment (H(t) = H). In close analogy to the measurement-operator duality (39), the Heisenberg jump operators are related pairwise to Schrödinger jump operators at dual parameters: whereas their corresponding coefficients obey This duality relation implies that the distinct anti-causal divisibility of the dynamics (all j H α (t) > 0) can be decided by the parameter dependence of the coefficients determining the causal divisibility properties (all j α (t) > 0). For the resonant level model, this is achieved by simply replotting Fig. 2(a) in units of temperature while varying the coupling Γ as shown in Fig. 2(b). The continuation of the causal divisibility boundary to negative coupling Γ precisely gives the anti-causal divisibility boundary that was shown in Fig. 2(a). The duality relation (61c) tells us precisely when causal and anti-causal divisibility coincide: j α (t) ≥ 0 for all α must imply (−1) Nαj α ≥ 0 and vice versa. This imposes a very strong constraint on the parameter dependence of the dynamics. This always holds when the evolution commutes with its generator, which includes the case of Markovian semigroup evolutions. We then have G H (t) = G(t) † , implying by Eq. (22) that H H (t) = H(t), J H α (t) = J α (t) † and j H α (t) = j α (t). Moreover, Eq. (61a) becomesH(t) = −H(t): the effective Hamiltonian is constrained to change sign under the duality mapping. In this case equation (58) additionally strengthens to a cross relation between the jump operators of G(t) alone: J α (t) † =J α (t) and their coefficients j α (t) = (−1) N α j α (t). This extends the results of Ref. [16] for the weak-coupling generators G in Lindblad form (57).
Beyond this trivial case the two types of divisibility need not coincide, as evidenced by the resonant level model [ Fig. 2(a)]. Thus, fermionic duality strongly suggests that anti-causal divisibility generically differs from causal divisibility by an explicit strong constraint on model parameters. This is in line with the general intuition that this type of divisibility additionally requires weak causal ordering of the dynamics and is thus a more fragile property. This motivates further investigation, for example in relation to recent work on causal ordering in quantum information theory [66,67].
We stress that although the duality relations (39) and (61) have a common origin, for general dynamics one cannot derive the jump-operator duality by using the trick of "linearizing" the measurement operators in Eq. (39) as it is possible in the Markovian semigroup limit [68]. The analogy between (39) and (61) is best seen in the Choi-Jamio lkowski correspondence to G(t) (instead of Π(t)) as discussed in App. C.

Unphysicality of the duality mapping
To conclude we verify thatḠ(t) is the generator of the dual evolutionΠ(t): using Eq. (23) and (49) What is interesting here is thatḠ(t) generatesΠ(t) in the same causal order as the Schrödinger evolution Π(t). On the other hand the dual propagatorΠ(t) is related to the propagator in the Heisenberg picture Π H (t) = e −Γt PΠ(t)P. This implies thatḠ(t) is related to the generator G H (t) acting from the left in the Heisenberg evolution [Eq. (46)] and not to G(t) † acting from the right. This explains why in Eq. (61c) the jump-coefficients j α (t) are related to the coefficients j H α (t) characterizing the anti-causal divisibility of the evolution [Eq. (58)] and not to the j α (t) describing ordinary causal divisibility. Note that the generator G H builds up the Heisenberg evolution in anti-causal order in contrast to G † . These observations are gratifying since they tie the physical divisibility properties of the dynamics to a key step in the derivation of the duality (23), the formal reversal of the causal ordering [Eq. (S-71) of Ref. [10]], within a completely different formalism.
We also observe that the relation between Π H (t) andΠ(t) is reflected by their generators appearing in the duality (49). Written as jump operator sum similar to Eq. (58), the dual generator reads The TP property ofΠ(t) corresponds to TrḠ(t) = 0 which is ensured by Eq. (53). In Eq. (63) this property is ensured by the causal structure of the dual dissipatorsD Even thoughḠ(t) has the causal structure and the TP property of a Schrödinger picture generator, we know that the dynamicsΠ(t) it generates is never CP [Sec. 3.2.3]. This general conclusion is not readily seen from the jump expansion (57) ofḠ(t), which is tailored to reflect divisibility properties. 18 However, it can be seen in the special case where G is time-constant: then Π(t) = e −iGt is CP-TP if and only if j α ≥ 0. Since in this case we also have j α = (−1) N α j α [Eq. (61) ff.] this impliesj α < 0 for all odd-parity jump-operators and thus the generated mapΠ(t) = e −iḠt is not CP.

Time-nonlocal quantum master equation
We now turn to the expression of fermionic duality in the last approach discussed in this paper, which will be particularly important for the application in Sec. 5. We now exploit that the evolution Π(t) is also the solution of the completely different time-nonlocal QME d dt In contrast to the time-local QME (45), its convolution structure matches the one obtained in the microscopic derivation of the evolution [10-14, 18, 19]: the propagator decomposes into a geometric series of convolutions of memory-kernel blocks −iK(t − t ) of duration t − t , giving a self-consistent Dyson equation: denoting (f * g)(t) =  64) and (65) we obtain the time-nonlocal Heisenberg QME noting that under the convolution one may commute 19 Π(t ) and K(t − t ). Inserting the propagator duality (23) and Eq. (64) on the left-hand-side of Eq. (66) we obtain the fermionic duality for the memory kernel For the resonant level model the memory kernel (16) indeed obeys this relation: this follows from the parameter dependence of the nontrivial functionk(t) = −k(t) [Eq. (26)] and the parity transformation of the dissipator PD † η P = −D −η − I.

Complex-frequency representation of dynamics
An advantage of the time-nonlocal QME (64) is that it allows a particularly simple explicit expression of Π(t) in terms of the Laplace transformK(ω) := ∞ 0 dte iωt K(t) of the memory kernel which facilitates further analysis: 18 For time-dependent generators which commute with the evolution, [G(t), Π(t)] = 0, we still have a direct relation jα(t) = (−1) N α j α (t). If Π(t) is CP-divisible, i.e., jα(t) ≥ 0 for all t ≥ 0 thenjα(t) ≤ 0 for all odd-parity operators. As mentioned in footnote 15 this does not allow to infer whetherΠ(t) is CP or not. If Π(t) is not CP-divisible, the signs of jα(t) are unrestricted and we cannot conclude either. 19 The identity K * Π = Π * K follows from the two ways Laplace transforming relation (23) gives the fermionic duality in frequency-domain reported in Ref. [10]: This relatesΠ(ω) and Π (ω) in complex-frequency regions where either both their Laplace transforms converge, or in regions where both are defined by analytical continuation. The mapping of the complex frequency argument reverses the real energy part of ω, while maintaining the sign of the dissipative imaginary part of ω up to a shift iΓ into the upper half plane. The fermionic duality for the frequency-domain memory kernel has the same structure:

Unphysicality of the duality mapping
While the above discussed operational approaches concern algebraic properties at each instance of time, the analytical structure of the memory kernel makes explicit how physical properties evolve in time. This makes fermionic duality in the frequency domain of independent interest (see below). It also reveals another way in which the dual propagator is unphysical as follows. Since a physical evolution Π(t) in general shows oscillations and decay or a combination thereof, its Laplace transform converges only for complex frequencies in the upper half plane. The obtained functionΠ(ω) has a unique extension to the lower half plane where in general it has both poles and branch points. This is illustrated for the resonant level model in Fig. 3. Integrating along any clockwise oriented contour C enclosing the poles and branch cuts (parallel to the imaginary axis) gives the general solution for the real-time evolution [11,12]: where Res f (ω p ) denotes the residue of f (ω) at ω = ω p . In view of our later application in Sec. 5.2, we note that the first term on the right hand side of (71) sums up contributions from two types of poles: those that arise due to the frequency-dependent eigenvalueŝ π i (ω) = i/[ω −k i (ω)] obeying the pole equationk i (ω p ) = ω p wherek i (ω) is an eigenvalue ofK(ω), and the remaining poles which also involve the eigenvectors. Since the parameter map Π(t) →Π(t) appearing in the duality relation (23) inverts the sign of dissipative decay rates, the Laplace transform ofΠ(t) converges only for frequencies above an imaginary cutoff in the upper half of the complex plane which is at least iΓ, the fundamental parity eigenvalue: If Π(t) converges to some stationary value Π(∞) this implies by Eq. (23) thatΠ(t) diverges at least as fast as e Γt which must be suppressed by e iωt in the Laplace transform. Thus, also its analytical structure proves clearly that the time-dependence of the dual propagatorΠ(t) is not physical, complementing the discussion of the algebraic, operational constraints of CP and TP [Sec. 3.2.3] which are independent of time, see also [33].

Cross-relations frequency-dependent left and right eigenvectors
Laplace transforming Π(t) and K(t), analytically continuing and diagonalizing giveŝ At resonance (ε − µ T ) these poles (branch cuts) cancel exactly leaving just a single pole at −iΓ whereas off-resonance (ε−µ T, Γ) they move to the sides where they become suppressed in amplitude. Only in these two limits four poles remain and the dynamics is a semigroup [23]. For T → ∞ the first case always applies.
Inserted into the memory-kernel duality (70) we obtain for the eigenvalueŝ with the duality between left and right eigenvectors Due to the simple relation (68) in the frequency domain the eigenvectors coincide, π i (ω) = k i (ω) and π i (ω) = k i (ω) , with eigenvaluesπ i (ω) = i/[ω −k i (ω)]: The left and right eigenvectors are indeed cross-related as dictated by Eq. (73b). In particular, the (non)trivial frequency dependence of the left (right) eigenvector for the eigenvalue with pole ω = 0 necessarily implies that the right (left) eigenvector for the eigenvalue with pole ω = −iΓ is (non)trivial as well. Thus, also in the frequency domain duality provides fine-grained insight into the location of nontrivial (non-exponential in time) contributions to the dynamics. In particular, the frequency dependence of the eigenvectors through the Laplace transform of k(t) [Eq. (5)], generates infinitely many additional poles at ω = ±(ε−µ)−iΓ/2−iπT (2n+1), n = 0, 1, . . . due to the digamma function ψ. For T → 0 the poles merge to form two branch cuts as shown in Fig. 3. In our application in the next section this analytic structure turns out to provide crucial insights.

Nonperturbative semigroup approximation and initial slip
Finally, we consider an application of fermionic duality where the insights of several of the discussed approaches come together. We consider analytic approximations to the solution of the time-local QME (45), d dt Π(t) = −iG(t)Π(t), constructed from the generator G(t) which we assume to be exactly known (best case). This equation naturally suggests a nonperturbative semigroup approximation which does not rely on any weak-coupling assumption [5,38]: requiring only that the generator converges to a stationary value G(∞), which is diagonalizable. It is not in general clear how accurate this approximation and corrections to it are. We will show that the quality of these approximations can be deeply understood using its exact relation to the corresponding time-nonlocal QME (64) and its memory kernel K(t) combined with fermionic duality. This effort is motivated by two attractive properties of the approximation (77): (i) For the large class of evolutions which are CP-divisible in the stationary limit, i.e., j α (∞) ≥ 0 in Eq. (57), the approximate evolution (77) is both CP and TP. This is in general very difficult to achieve for nonperturbative approximations [29][30][31][32][33]. This class includes dynamics which is not a trivial semigroup described by a Lindblad equation, which is already the case for the resonant level model (except for T = ∞ or ε = µ, see Sec. 2). It also includes dynamics which is not CP-divisible as long as j α (t) ≤ 0 occurs only for finite times.

Fixed-point relation between generator G and memory kernel K
To address this problem, we will make use of a recent exact result [5] which shows that the stationary generator obeys the self-consistent equation [App. E] Here the superoperator G(∞) takes the role of the complex frequency ω in the Laplace transform of the memory kernel K(t). Inserting the spectral decompositions, this implieŝ , and a careful analysis shows that the eigenvalues of the stationary generator G(∞) are eigenvalue-poles 21 ofΠ(ω), i.e., ω p =k j (ω p ) [Eq. (68)]. Equation (78) thus states that G(∞) "samples" the Laplace transformK(ω) of the memory kernel precisely at complex frequencies given by the eigenvalues of G(∞): Herek j i (g i (∞)) = g i where j i denotes the index of the eigenvaluek j (ω) ofK(ω) which equals g i at frequency ω = g i . By trace-preservation, the frequency sampling always includes zero, g 0 (∞) = 0, and thus G(∞) andK(0) have the same stationary state eigenvector, proving property (ii) mentioned above. However, also nonzero complex frequencies are sampled. Thus, only the set of eigenvectors { k j i (g i (∞)) }, collected from different superoperators [K(ω) at different frequencies ω = g i (∞)] provides the full set of right eigenvectors { g i (∞) } of the single superoperator G(∞). This in turn determines the left eigenvectors { g i (∞) } of G(∞) and thus, remarkably, G(∞) can be directly constructed from the memory kernelK(ω) once one knows the eigenvalues g i (∞). For the resonant level model this surprising construction was verified explicitly in Ref. [5].
Ref. [5] focused on the implications of relation (78) assuming that one knows K and aims to compute G. Here we focus in a sense on the converse question: Using a known G to construct an approximation for Π(t), what general insights into the quality of the approximation does K provide, and how does fermionic duality help in this analysis? We consider the simplest approximation beyond the semigroup approximation (77), a corrected semigroup which aims at improving the description of the evolution at long times. To illustrate the simplest type of such a correction we assume in the following that the poles ofΠ(ω) at the eigenvalues of G(∞) are of first order. 22 We start by noting that by Eq. (79) the eigenvalues of G(∞) must be poles in the eigenvalues ofΠ(ω) [first term of Eq. (71)]. Their contribution to the exact evolution can be expressed as Here the sum over i runs over distinct values of g i (∞). This indeed looks like an initial slip correction to the semigroup approximation. We will now discuss the quality of the two approximations Π (1) and Π (2) and then derive an exact restriction on this procedure imposed by fermionic duality.

Nonperturbative semigroup approximation
The quality of the semigroup approximation (77) can be investigated using the result (80a) of the time-nonlocal QME approach. We observe that in the complex-frequency domain the errorΠ(ω) −Π (1) (ω) has in general the same poles as the exact dynamics except for the pole at zero frequency. Thus, the approximation captures only the stationary state systematically. 21 Because we assume that G(∞) exists and Eq. (78) holds, K(gi(∞)) cannot have poles in its eigenvectors. 22 Our assumption is equivalent to dkj i (ω)/dω| ω=g i (∞) = 1. Higher order poles would require a timedependent slip superoperator S(t) = S1 + S2t + S3t 2 + · · · where Sn is constructed from the residues of all n th order poles ofΠ(ω) at eigenvalues of G(∞).

Nonperturbative initial slip correction
Within the time-local QME approach one may set up a nonperturbative correction of the semigroup approximation by an initial slip [69][70][71][72][73] using some superoperator S: This initial-slip approximation is precisely what is obtained in the time-nonlocal QME approach when selecting the exact poles g i (∞) in the contour integration (71) of the inverse Laplace transform. This relation can be used to investigate how well one can do in the time-local approach: we know from Eq. (80) that such a correction can be achieved by Note that this expansion is not the spectral decomposition of S, since the sets of left and right eigenvectors are not biorthogonal. This is in fact the best one can do for an approximation in the long time limit since in the frequency domain this approximation for the resolvent readŝ showing that now the exact poles g i are completely canceled in the errorΠ(ω) −Π (2) (ω) as illustrated in Fig. 4(b). Here we assumed that the eigenvalues g i (∞) are nondegenerate as is the case in the resonant level model, see App. D for the degenerate case. In Fig. 5(a) we illustrate that for the resonant level model the inclusion of the exact initial slip S may indeed lead to a much faster approach to the exact dynamics than the nonperturbative semigroup. In Eq. (80b) this is achieved by first mapping the initial state using S to a possibly nonpositive density operator-reflected in the plot by the unphysical occupation > 1-and then applying the semigroup evolution (77). Indeed, although S and Π (2) (t) are both TP maps 23 , the map S is not CP. Thus, the initial dynamics is unphysical whenever the slip is nontrivial, since by construction Π (2) (0) = S = I [74][75][76]. For the same reason, the slip-approximated dynamics (81) is not a semigroup: even though the decay is exponential, any nontrivial slip obstructs addition of the exponentials: e −iG(∞)t Se −iG(∞)t S = e −iG(∞)(t+t ) S. Although this is not a problem-the exact dynamics Π(t) is not a semigroup either-it may be overlooked that summing isolated pole contributions [Eq. (80a)] in the time-nonlocal QME approach, one does not obtain a Markovian semigroup approximation.
In Fig. 6 we investigate this correction more closely for the resonant level model by plotting the time at which Π (2) (t) becomes CP, a necessary indicator of quality. Interestingly, near resonance, ε = µ, this time may diverge for specific physical values of the strong coupling. Fig. 5(b) illustrates that in their vicinity the initial slip may give a very large correction even though the semigroup approximation is very close to the exact dynamics. This breakdown is easily overlooked when constructing the slip approximation within the time-local formulation in which the frequency-domain structure is not available. This crucial insight is enabled by the fixed-point relation (78) of Ref. [5].
We note that in the time-local approach, one might expect the exact slip superoperator to be expressible as S = lim t→∞ e iG(∞)t Π(t) but this limit may actually fail to exist. In the resonant level model this indeed happens when the coupling exceeds a sharp threshold, Γ ≥ 2πT . However, one can regularize this expression to recover Eq. (82) by taking the zero-frequency residual of its Laplace transform L: S = −i Res L e iG(∞)t Π(t) (0). We further note that Markovian approximations can also be performed starting from the Heisenberg equation of motion d dt A(t) = iG H (t)A(t) as is often done in quantum optics. However, naively following the same steps in this picture leads to different results with problematic features. This further subtlety and its resolution are discussed in App. D.

Fermionic duality for the initial slip
The restrictions imposed by fermionic duality now provide a crucial insight: already for the resonant level model, divergences of the CP-time quality indicator [ Fig. 6] are unavoidable. This can be understood by taking the stationary limit of the duality relation (49) which, provided S −1 exists, simplifies to Here G(∞) applies the duality mapping after the stationary limit [App. D]. This matters since although the generator G(∞) exists for all parameters of the model this is not true for G H (∞) andḠ(∞), a further illustration of the nontrivial relation between the Schrödinger and Heisenberg time-local generators. Importantly, the slip superoperator that appears here obeys a separate fermionic duality which is another key result of the paper: For the resonant level model, the exact relations (84) and (85), together with the trace preserving property of S, completely fix the relevant part of the slip superoperator constructed from the generator G(∞) [App. D]: wherek(i 1 2 Γ) = (2/π) Im ψ 1 2 + [Γ/2 + i(ε − µ)]/(2πT ) [cf. Eq. (76)]. As function of parameters the expressionk(i 1 2 Γ) dictated by fermionic duality is singular at physical parameter points ε = µ ± 0 + , Γ = (1 + 2n)2πT (87) with n = 0, 1, 2, . . ., causing the slip approximation (80b) to break down as in Fig. 5(b). Moreover, the Heisenberg and dual generator diverge here. This illustrates that fermionic duality can be a useful tool for understanding the often highly nontrivial properties of time-local generators of quantum evolutions which complicate analytic approximations.
In the time-nonlocal approach the breakdown can be clearly understood in the frequency representation: Close to the values (87), three poles in the exact result forΠ(ω) [Eq. (75)] approach the same point (ω = −iΓ). One pole comes from an eigenvalue of Π(ω) (lowest cross in Fig. 3) and two others come from its eigenvectors (top two unmarked poles in Fig. 3). The prefactor of the eigenvalue-pole diverges as function of parameters as 1/(ε − µ) but in the exact result this is canceled by corresponding divergence of the prefactors of the two eigenvector poles. However, in the slip approximation the latter two poles are discarded [ Fig. 4], leaving the spurious divergence of the remaining one.

Discussion
In this paper we have shown that for a large class of fermionic open systems the nontrivial relation between state and observable evolution can be completely bypassed by a simple fermionic duality mapping that exploits and generalizes the functional dependence of the two evolutions on the microscopic physical parameters. We have shown that this works for essentially all canonical approaches used in quantum transport, open-system dynamics and quantum information theory without introducing any assumptions (such as weak coupling, high temperature, various Markovian approximations, etc.) except for wideband coupling to the reservoirs. The obtained fermionic duality relations are summarized in Table 1. In the superoperator-based approaches these imply exact parity eigenvalues and eigenvectors [Eqs. (32), (53)] and nontrivial cross-relations for the entire spectrum [Eqs. (29), (50), (73), (74)]. Correspondingly, in the operational approaches we derived additional fermionic sum rules for measurement and jump operators [Eqs. (40), (59)] and their scalars coefficients [Eqs. (41), (42), (60)], and nontrivial cross-relations for these entire sets of operators [Eqs. (39), (61)]. Combining the latter approaches with fermionic duality naturally led us to consider a new type of divisibility of the dynamics: We noted that the Schrödinger and Heisenberg time-local generators through their quantum-jump coefficients encode both ordinary causal and anti-causal divisibility, respectively. Using fermionic duality we showed how the operational condition of anti-causal divisibility can be inferred from the condition of ordinary, causal divisibility. Dynamics which does not commute with its generator may be causally divisible but fail to be anti-causally divisible, as we demonstrated by an explicit example. This provides definite information about the causal ordering of the dynamics, i.e., whether the ordering of dividing the dynamics matters.
Throughout we emphasized both the usefulness of duality relations such as Π(t) † = e −Γt PΠ(t)P and the unphysicality of mappings such as Π(t) →Π(t). The usefulness was illustrated by identifying various "symmetries" in the dynamics of the resonant level model which went unnoticed so far. Importantly, our results apply much more generally and are therefore relevant for a large class of outstanding dynamical open-system problems in regimes of strong coupling, strong interactions, finite temperature and nonequilibrium. This holds in particular for much of the analysis and applications in Sec. 5, which apply to complex models of high interest for which the calculations are of course very complicated.
Here one should remember that the merit of duality lies in simplifying a calculation given a method of choice, not in providing this method. In this application we exploited fermionic duality to provide deeper insight into a nonperturbative semigroup approximation which respects CP-TP for a broad class of dynamics, using the key result from Ref. [5] that relates time-local and non-local QMEs. We showed that inclusion of the initial-slip correction [Eq. (80b)] in the time-local approach (TCL) corresponds precisely to a selection of poles in the time-nonlocal approach (Nakajima-Zwanzig). We found that even for the resonant level model this correction unexpectedly fails at specific physical parameters and that this failure is unavoidable due to fermionic duality constraints. The failure of the slipcorrection is interesting: in the resonant level model it is caused by the poles that form branch-cuts in the limit T → 0. That branch cuts invalidate exponential approximations at T = 0 is not unexpected, but here we found that precursors of branch-cuts already cause havoc at finite temperature T .
The unphysicality of the duality mapping appeared particularly clearly in the operational formulations tailored to quantum information theory which show that the evolution Π(t) at dual parameters violates complete positivity. We highlighted this point since it implies that when tacitly assuming this valid and important restriction in any of the operational approaches one may easily overlook the powerful constraints imposed by fermionic duality. For example, in the operational approach evolution maps are invariably written as Π = α M α • M † α (which we avoided doing) by anticipating positive coefficients in the operator sum and absorbing them into the norm of the measurement operators, i.e., Tr M † α M α = m α . This automatically eliminates the dual superoperatorΠ(t) with its benefits from further considerations since it has no such expression.
Relation to other works. The fermionic duality is most closely related to the extension of PT-symmetry [77,78] to dissipative systems [79,80]. Unlike ordinary symmetries which relate for example the evolution to itself by conjugation with a symmetry transformation, the evolution is related to its adjoint. In the present paper we have emphasized this relation as a connection between the mathematically and physically distinct evolutions of states and observables. In the Refs. [79,80] this was achieved under the strong assumption of Markovianity in the sense of semigroup divisibility (Lindblad). The fermionic duality which was derived in Ref. [10] is instead based on far less restrictive assumptions and is applicable to strongly non-Markovian dynamics such as in the resonant level model.
The involvement of the adjoint sets fermionic duality apart from standard symmetry consideration transposed to Liouville space-called "weak symmetry" in Ref. [81]-where the time-local generator commutes with the symmetry superoperator. Open systems also allow for a notion of "strong symmetry" where the Hamiltonian and jump operators of this time-local generator commute with a symmetry operator. This stronger notion introduced in Ref. [81] for semigroup-Markovian systems played a role in our analysis albeit in modified form (the jump operators may either commute or anticommute with the fermion parity operator, see Eq. (57) ff. and App. C). However, fermionic duality is distinct from both these notions of symmetry.
We furthermore note that after the original derivation of fermionic duality in Ref. [10] subsequent work appeared [82,83] which exploited similar tricks to simplify the calculation of open system evolutions, such as unphysical, non-Hermitian coupling to reservoirs with structureless wide bands. However, in addition to focussing on bosonic systems, these works relate the environment of a system of interest to a simpler, effective environment to reduce the computational complexity. In contrast, the fermionic duality relates a system and its environment to an equally complex dual system and environment, in order to exploit the functional parameter dependence of the evolution.
Finally, our finding of a new type of divisibility of the dynamics raises an interesting question regarding stochastic simulations as different types of divisibility are at the basis of different simulation methods. Whereas CP-divisibility of the dynamics enables an implementation that directly uses the jump operators to generate quantum jumps [51,84,85], P-divisibility allows only for an indirect implementation of the jumps via Diosi's rateoperator [51,86,87]. If the evolution has no divisibility property, its simulation is more complicated, requiring reverse quantum jumps connecting trajectories [51]. However, all these distinctions are based on causal division of the dynamics. Fermionic duality surprisingly enables conclusions about anti-causal divisibility, an apparently new concept with a clear operational formulation independent of fermionic duality. It presents a refined distinction between different types of dynamics. For the resonant level model we found that anti-causal CP-divisibility is lost while the dynamics remains causally CP-divisible and thus efficiently simulateable. It is an interesting open question how to detect the loss of this property on the level of simulated quantum trajectories.
Outlook. Exploiting the results reported here in applications to interacting models with strong coupling requires that one uses an approach that maintains fermionic duality in approximations. The renormalized perturbation theory [13,14] which was used originally [10,16] to derive the duality relation (3) exhibits this feature: It preserves duality order by order in a renormalized temperature-dependent coupling, see Ref. [88] for a recent implementation. Ordinary perturbation theory in the bare coupling Γ [89][90][91] already breaks fermionic duality in the next to leading order Γ 2 . For the Anderson model this breakage seems specifically related to the electron pair-tunneling contributions [90,92,93].
For calculations involving stronger coupling the renormalization-group approach of Refs. [11][12][13][94][95][96][97][98][99] is well-suited. Its original formulation [11,13] is build on the renormalized perturbation expansion, which explicitly preserves fermionic duality [10]. Although in Ref. [13] the first implications of fermionic duality were discovered and applied in a very advanced context, it remains an interesting open question which truncation schemes for this exact hierarchy of the RG equations maintain the fermionic duality. Moreover, the potential advantages of the duality for the more recent E-flow formulation [95][96][97]99] also remain unexplored. Finally, it is of interest to understand how approximations can be formulated within other nonperturbative approaches in a way that maintains fermionic duality, i.e., by which rules. A key step in this direction would be an elementary microscopic derivation of fermionic-duality within these approaches. Our finding that fermionic duality takes a simple form in each of the canonical approaches to quantum dynamics suggests that this is possible.
A Duality for Choi-Jamio lkowski state of propagator Π(t) CP and CJ operator. A dynamical map ρ(t) = Π(t)ρ(0) is completely positive (CP) if and only if it preserves positivity when evolving the system together a non-evolving auxiliary system: (Π(t) ⊗ I)ρ ext (0) ≥ 0 for any initial state ρ ext (0) of the system plus any auxiliary system. This is equivalent to positivity for the worst case ρ ext (0) = 1 d |1 1| of a maximally entangled state |1 = k |k |k on the tensor product of the system Hilbert space with an auxiliary copy of itself. Thus, the so-called Choi-Jamio lkowski (CJ) operator should be positive: choi[Π(t)] := Π(t) ⊗ I |1 1| The operators M α (t) in the operator sum (35) for a Hermicity-preserving superoperator Π(t) are obtained by diagonalizing the Hermitian choi[Π(t)]: When diagonalized with eigenvectors normalized to 1 the real eigenvalues m α (t) are the coefficients in the operator sum (35). If Π(t) is CP they are positive since then choi Π(t) ≥ 0 as noted in the main text. The bipartite eigenvectors |M α (t) = (M α (t) ⊗ 1)|1 = ij i|M α (t)|j |i |j determine the matrix elements of the canonical measurement operators relative to the chosen basis {|i }.
Parity of measurement operators. Fermion superselection for the total system evolution and initial reservoir state implies that the reduced system evolution Π(t) commutes with the parity superoperator (−1) This shows that the eigenvectors |M α (t) have definite bipartite parity and thus the measurement operators must have definite parity: (−1) N M α (t)(−1) N = (−1) Nα M α (t) with N α = even or odd, as claimed in the main text.
Fermionic duality. Using the property that in the maximally entangled state the action of any operator on the system is perfectly transposed to its copy, A ⊗ 1|1 = 1 ⊗ A T |1 , the relation (23) implies the fermionic duality for the CJ state involving the bipartite swap operator S|i |j = |j |i . Thus, if |M α (t) is a right eigenvector of choi[Π(t)] with eigenvalue m α (t), then [S|M α (t) ] * = |M α (t) † is also a right eigenvector with eigenvalue m α (t) = e −Γt (−1) N α m α (t) proving Eq. (39b) in the main text. Using |M α = (M α ⊗ 1)|1 we also establish the fermionic duality (39a) for measurement operators.
This means that corresponding partial operator sums for the dual eigenvalue pair α and α are equal:

B Relation specific to resonant level model
In addition to the generally valid duality relation (49), the time-local generator of the resonant level model obeys another, simpler relation. This may be understood also from the formal similarity of Eq. (8) and Eq. (18). Despite the fact that Π(t) = T exp − i t 0 dsG(t) and the time-ordering T is nontrivial for this model, it holds true that Π(t) = exp(−itG(t))| g(t)→p(t) . Inserting this into the relation (23) and usingp(t) = −p(t) one obtains by comparing exponents where on the right we replace "by hand" g(t) → −g(t) in analogy to the transformation of p(t) under the parameter substitution. As a result, the left and right eigenvectors of G(t) are formally related by taking the adjoint, multiplying with the fermion parity (−1) N and replacing g(t) → −g(t), as one can verify in Table 2. Although relation (94) is simpler and inferred by inspection, it is not valid for the general class of models for which Eq. (50) holds.
C Duality for Choi-Jamio lkowski operator of generator G(t) Jump expansion for G(t). The canonical jump-expansion for a time-local generator G(t) is obtained by constructing its CJ-operator (88) proceeding analogous to Π(t) as in App. A. However, G(t) is not a CP map (unlike Π(t)) and we can only use that Tr G = 0 (instead of Tr Π = Tr), and that −iG is Hermicity-preserving and thus choi[−iG(t)] is Hermitian. The canonical form is obtained by diagonalizing the projection of choi[−iG(t)] on the orthogonal space of the maximally entangled state |1 := k |k |k . This gives the last term in Eq. (95) and the first two terms account for the remaining matrix elements. The projector to this orthogonal space is denoted as Q := 1 − |1 1|/d. Splitting B = Re B + i Im B one checks that the Hermitian part is fixed to Re B = − 1 2 α j α J † α J α by the condition Tr G(t) = 0. The remaining anti-Hermitian part defines the effective Hamiltonian H(t) = − Im B(t) in Eq. (57).
Parity of jump operators. In our case we can also use that fermion superselection, [(−1) The remaining eigenvectors of the projection with Q have definite bipartite parity of either sign, (−1) N ⊗ (−1) N |J α = (−1) Nα |J α for N α being even or odd. These determine the jump operators through |J α (t) = (J α (t) ⊗ 1)|1 for which (−1) N J α (t)(−1) N = (−1) Nα J α (t) as claimed in the main text. Note that this implies that Re B has even parity consistent with the above.
Degenerate coefficients j α (t). If some coefficient j α is a degenerate eigenvalue of Q choi[−iG(t)]Q in the construction of Eq. (95) then similar remarks apply as for the measurement operators. The partial operator sums for the dual eigenvalue pair α and α are equal, we find the result of the main textH(t) = −H(t), J α (t) † =J α (t) and j α (t) = (−1) N α j α (t) for nondegenerate coefficients.

D Duality for stationary generator and slip superoperator
Duality for the slip superoperator S [Eq. (86)]. Using the t → ∞ limit of Eq. (50a), g j (∞) = iΓ −ḡ i (∞) * , and Eq. (69) in Eq. (82) we obtain writing g i for g i (∞) Solving these equations for the resonant level model gives Eq. (86) up to an additional term η α η (t) d η d η with undetermined function α η (t) * =ᾱ η (t). The latter term is irrelevant for the occupation dynamics, in which the breakdown of the initial slip correction occurs. Thus duality fixes the relevant part. Note that the evolution of the coherences is already exact in the semigroup approximation. Setting α η (t) = 0 we obtain the exact result.
Stationary limit of time-local duality (49) [Eq. (84)]. Noting that the left action of G(∞) on the slip operator gives G(∞)S = −i i g i ResΠ(g i ) by Eq. (82b), we insert the with anti-timeordering T → . In Ref. [5] it was shown that for t → ∞ this gives the stationary fixed-point relation (78) used in the main text.
Analogous to Eq. (102) the Heisenberg generator G H (t) = [Π(t) −1 G(t)Π(t)] † also obeys a functional fixed-point equation with the memory kernel K H (t) = K(t) † : The proof of this relation is analogous to the proof of the functional fixed-point equation in Schrödinger picture in Ref. [5]. In neither of the above general relations fermionic duality was used and it is thus important to check that it is consistent with these relations. To see this, substitute (49) and (67)