Tomographically reconstructed master equations for any open quantum dynamics

Memory effects in open quantum dynamics are often incorporated in the equation of motion through a superoperator known as the memory kernel, which encodes how past states affect future dynamics. However, the usual prescription for determining the memory kernel requires information about the underlying system-environment dynamics. Here, by deriving the transfer tensor method from first principles, we show how a memory kernel master equation, for any quantum process, can be entirely expressed in terms of a family of completely positive dynamical maps. These can be reconstructed through quantum process tomography on the system alone, either experimentally or numerically, and the resulting equation of motion is equivalent to a generalised Nakajima-Zwanzig equation. For experimental settings, we give a full prescription for the reconstruction procedure, rendering the memory kernel operational. When simulation of an open system is the goal, we show how our procedure yields a considerable advantage for numerically calculating dynamics, even when the system is arbitrarily periodically (or transiently) driven or initially correlated with its environment. Namely, we show that the long time dynamics can be efficiently obtained from a set of reconstructed maps over a much shorter time.


I. INTRODUCTION
A large fraction of active research in physics and chemistry, both theoretical and experimental, involves characterising or modelling the dynamics of open quantum systems. Such dynamics is, in general, significantly more complicated than that dictated by Schrödinger's equation. Nevertheless, many techniques have been developed to predict how open systems will evolve in time. These include schemes for tomographically determining an unknown open process [1][2][3][4][5], perturbative approaches to theoretically modelling dynamics [6][7][8], and computationally expensive numerical techniques which can exactly simulate an open system under certain circumstances [9][10][11][12][13]. However, the resources required for the latter often scale exponentially with the evolution time, and numerical shortcuts usually rely on assumptions about the dynamics, such as no time-dependent driving [14] or limited memory effects [15][16][17], meaning long-time dynamical simulations are intractable for the most general open processes. Here, among other results, we develop a method for extracting the longtime dynamics for general driven open quantum systems.
Conceptually, approaches to describing open quantum dynamics tend to take one of two perspectives, depicted in Fig. 1: First is the "dilated" or "underlying" picture, where the system in question (S) is represented as evolving with its environment (E) from an initial state ρ SE t0 , with a time evolution superoperator U SE t:t0 [18] that propagates solutions to an equation of the formρ All the physics of the system is encoded in the superoperator L SE t , which is usually taken to evolve SE according to the von Neumann equation [19], but could also incorporate dissipative and decohering dynamics through (possibly timedependent) Lindblad terms [20,21].
The second picture is the operational one, where the entire dynamics is described in terms of quantities that would * felix.pollock@monash.edu † kavan.modi@monash.edu be, in principle, directly measurable by an experimenter with access to the system alone. Such a description contains only that information necessary to predict the evolution of the system, nothing more, and is usually more compact than a corresponding description in the dilated picture. In this picture, the dynamics may often be described in terms of a dynamical map Λ t:t0 , which is completely positive and trace preserving, and evolves the initial system state ρ t0 to its counterpart at time t: ρ t = Λ t:t0 ρ t0 [22]. This description is valid when the system and environment are initially uncorrelated in the other picture. More generally, the dynamics can be described in terms of a superchannel [2,23] M t:t0 , which maps an initial preparation superoperator A-any transformation that could be applied to the system at time t 0 -to the state at later times: ρ t = ρ t,A = M t:t0 [A]. The superchannel guarantees the complete positivity of the dynamics even in the presence of initial correlations with an environment. The time-evolved system state that one would obtain without performing any preparation procedure (the 'freely-evolved' state) is simply given by M t:t0 [I], where I is the identity superoperator. The dynamics, Λ t:t0 or M t:t0 , can be reconstructed operationally by means of quantum process tomography [1,3,24].
While these two pictures are equivalent, it is not always clear how to convert from the operational to the dilated picture (the reverse is straightforwardly achieved by tracing out over the environment, ρ t = tr E ρ SE t ). That is, there is no generally applicable technique which allows the underlying SE quantities, ρ SE t0 and L SE t , to be recovered from operationally reconstructible maps, such as M t:t0 . Where this is possible, even partially, the system can be used as a probe of its environment, extracting information about its Hamiltonian directly [25,26], or otherwise determining the parameters of a phenomenological master equation [27][28][29][30].
Partial knowledge of the global dynamics, inferred from the evolution of the system, can also be used to more efficiently predict future dynamics. This is embodied in the transfer tensor method [14], where short time system dynamics can be used to construct a discretized memory kernel, in the form of the eponymous transfer tensors. These can then be used to propagate the system to later times, simulating the long-time dynamics with exponentially fewer resources than t 0 with its environment, is prepared at time t0 with operation A. The two then evolve together according toρ SE t = L SE t ρ SE t , leading to a final state ρt for the system that depends on the initial preparation. (b) The system state at time t can also be thought of as a function of the preparation operation at time t0. All the uncontrollable parts of the process (enclosed within the dashed lines in the first picture), including the initial system-environment state, are encapsulated in the linear map Mt:t 0 . When ρ SE t 0 is uncorrelated, the dynamics can be equivalently thought of as a completely-positive trace-preserving map from the initial state of the system. many other methods [31,32]. Memory kernels are often used to formulate exact, continuous time open dynamics, most famously in the form of the Nakajima-Zwanzig master equation [6]. However, unless the underlying dynamics is homogeneous in time, and the system is initially uncorrelated with its environment, there is no clear way to construct the memory kernel or transfer tensors operationally, either in experiment or in numerical simulations.
In this Article, we start by solving this problem definitively: First, in Sec. II we derive transfer tensors from first principles for any open quantum dynamics, in particular those with initial correlations and time-dependent generators. We then go on, in Sec. III to demonstrate the power of the transfer tensor technique for simulating long-time dynamics, explicitly recovering the dynamical steady state of a driven, dissipative example system. In Sec. IV, we construct a master equation from the transfer tensors, demonstrating a direct correspondence between tomographically reconstructed dynamics and a generalised Nakajima-Zwanzig master equation. Finally, in Sec. V, we show how this correspondence could be used in an experimental setting to relate directly measurable quantities to properties of the underlying dynamics.
Our results render a large class of difficult numerical problems efficiently solvable and build a new link between underlying physics and experimentally accessible quantities. In addition, they open up the possibility of deriving approximate master equations for the system, through operationally meaningful coarse-graining of the reduced dynamics. Before presenting our main results, we first introduce memory kernels and their relation to system's equation of motion.

Memory kernel master equations
The most general form for the generator of a memoryless, or Markov, quantum evolution was famously derived by Lindblad [20], as well as Gorini, Kossakowski and Sudarshan [21]. In this case, the derivative of the system state at time t is a homogeneous linear function of the state at that same time, as in Eq. (1). For more general, non-Markovian processes, one natural way to extend the formalism is to consider master equations in which the time derivative depends on the history of the system's evolution: where L t is a time-local generator, K t,s is a superoperator known as the memory kernel, and J t,t0 is an inhomogeneous term, that can depend on the initial condition. Put simply, the memory kernel K t,s quantifies how the system state at time s influences, through interaction with the environment, the evolution at later time t. Knowledge of the functions on the right hand side of Eq. (2) allows for the self-consistent solution of the system's dynamics.
In recent years, several phenomenological master equations with memory kernels have been proposed [33][34][35][36][37]. It is always possible, given a microscopic description, to cast any open dynamics in the form of Eq. (2) using the Nakajima-Zwanzig projection superoperator technique [6] (see Sec. IV and Appendix D), though the solution of the resulting equation is far from trivial and calculating the memory kernel for different models remains an ongoing research problem [38][39][40]. However, there is currently no prescription for determining the elements of Eq. (2) when a microscopic model is unavailable. In Sec. IV we show how this can be done, in a simulation or an experiment, in terms of completely positive dynamical maps, but first we show how a discrete time memory kernel structure arises, for general open dynamics, in the form of transfer tensors.

II. OPERATIONAL DERIVATION OF TRANSFER TENSORS
The transfer tensor method was first proposed as an ansatz in Ref. [14] under certain restrictive conditions: (i) the underlying dynamics is time homogeneous, i.e., Our key result is a first principles derivation of the transfer tensors for a general open evolution, without making any of the above assumptions.
We begin by introducing a set of objects {P s , Q s } which act as such that (P s + Q s )ρ t = ρ t . Subsequent action of these objects further breaks up the evolution, i.e., P s Λ t:s = Λ t:s Λ s :s . Here, Λ t:s is a completely positive dynamical map from time s to time t, Λ t: is a unit trace time-dependent density operator not necessarily related to the actual environment dynamics; we will elaborate on the interpretation of these maps in Sec. V.
Next, iterating the identity (P s + Q s ) = I, we can write the evolution up to time t = t N in terms of the evolution to Here, the T (n) t,s are exactly the transfer tensors defined in Ref. [14], except that we are able to explicitly derive them in terms of, in principle, operationally accessible maps, without appealing to time-homogeneity of the memory kernel. They are given by Chasing this definition leads to the explicit form Finally, the remaining inhomogeneous term in Eq. (5) encodes all memory effects arising from the initially correlated SE state, and is defined as the difference ρ tj ; this can be operationally determined for any initial preparation by reconstructing the superchannel M t:t0 . For long enough times, this term can be neglected [41]; however, in Appendix A we show how this term can be included without direct calculation at the expense of reconstructing a larger set of dynamical maps. There, we account for the initial correlations operationally by replacing the map Λ t:t0 with the completely positive superchannel M t:t0 [2]. We will now show explicitly how the transfer tensors derived in this section can be used to simulate driven open systems.

III. EFFICIENTLY SIMULATING LONG-TIME DYNAMICS
The principle advantage of the transfer tensor approach over other simulation methods is that it allows for the determination of long-time dynamics, with an error that does not grow with the evolution time [42] but instead depends on the time beyond which memory effects are neglected. This is in contrast to many numerically exact techniques, where the computational effort to achieve a given numerical precision grows exponentially with the evolution time [13]. In what follows, we will outline how a transfer tensor simulation could be performed and how the error can be bounded.
Fixing δt = t j+1 − t j , which determines the time resolution of the resulting dynamics, we will assume that the dynamical maps Λ t:s are periodic (or transient), such that Λ t+T :s+T = Λ t:s for large enough s, and where the period T = c δt with integer c. This is equivalent to requiring that the underlying generator L SE adjacent time steps and iterating with Eq. (6)). Using Eq. (5), the state at any arbitrary time t k can be found by identifying T (l) t k−l mod T +lδt:t k−l mod T , for l ≤ m, and with 0, for l > m (similarly, Ξ t k ,t0 is set to 0 for k > m), before propagating the initial state. We reiterate that the complexity of this latter stage is independent of the underlying SE model, meaning it could be equally well applied to any open system whose short time dynamics can be determined.
In Appendix B, we show how the maximum distance between the numerically propagated stateρ t k and the exact state ρ t k can be bounded as where X := max ρ,σ tr{σXρ} is the operator norm of the matrix representation of superoperator X that acts on vectorised density operators. In other words, the error grows with the cumulative size of the transfer tensors beyond the first memory time. What this means is that one can guarantee a desired accuracy in the propagated state by varying t m and constructing transfer tensors to time T + 2t m until the right hand side of the above equation is small enough. The upper bound is approximate, insofar as it only considers contributions to the error from the second memory time period (as long as the memory decays faster than arithmetically, then contributions from subsequent memory times will not diverge). However, we find that this bound tends to grossly overestimate the propagation error, which is more closely bounded by the largest norm of the 'longest' transfer tensors: max x≤c T (m) xδt+tm:xδt . On the other hand, if the memory cutoff approximation is too severe, then the propagated dynamics may become unphysical, leading to divergent error in some cases (as in other approximate methods, where throwing away significant terms results in unphysical behaviour [6,7]).
To illustrate the effectiveness of the transfer tensor we now consider a simple example of a driven, dissipative open system. The dynamics of the total system-environment state is governed by a memory-less master equation: While the SE dynamics is Markovian, the dynamics of S alone will be non-Markovian and will have the form of Eq. (2). For concreteness, we take the system-environment to be a pair of qubits. The two evolve according to the Hamiltonian  Example propagated dynamics (a) elements of the system density operator at short and long times for model described in the text. Transfer tensors are reconstructed, using a time independent reference state, between the first nine time steps (red and blue crosses), before being used to propagate the state to later times. The propagated population of the system's excited state (blue dots) is shown alongside exactly calculated dynamics (continuous lines). Real and imaginary parts of the system's coherence are plotted in the inset (the coherence is zero at long times). Also shown (yellow squares) is the trace distance between the exact density operator and its transfer tensor propagated counterpart. (b) Long time propagation error ρt −ρt 1 as a function of memory time tm and time step size δt; also plotted, for reference, is the norm of the memory kernel at tm. Points with arrows correspond to cases where the distance is greater than two (sometimes by many orders of magnitude), indicating non-physical dynamics. In all cases, the bound in Eq. (9) is larger than the error.
In the figure, short time exact dynamics, calculated between the first nine time steps (using a time-independent reference state τ E = tr S ρ SE 0 ), is used to propagate the statẽ ρ t k to much longer times with no noticeable increase in error, as indicated by the yellow dots, which remain close to zero for all times. This can also be observed by comparing the exact dynamics (solid curves) with the propagated points. The dynamics are accurate despite choosing a relatively short memory cutoff time (ωt m = 5). As shown in panel (b), the memory kernel still has a significant norm at this time, which would usually lead one to conclude that we are discarding important memory contributions when making the cutoff.
One of the reasons for the remarkable effectiveness of the transfer tensor approach, in this case, is the coarse sampling (long δt) of the initial dynamics. In general, as seen in panel (b) of the figure, where we calculate the long-time propagation error for a variety of cutoff times and time step lengths, we find that finer grained dynamics require longer memory times for the same accuracy [44]. In every case, we find that the bound in Eq. (9) is extremely loose, meaning that the transfer tensor propagation is more robust than indicated by the conservative analysis presented in Appendix B. In the limit of continuous time dynamics, this approach leads us to a generalised version of the Nakajima-Zwanzig equation, of the form of Eq. (2), as we will now show.

IV. EQUIVALENCE WITH THE NAKAJIMA-ZWANZIG EQUATION IN THE CONTINUUM LIMIT
While Eq. (2) governs the exact dynamics of the system alone, conventional derivations usually rely implicitly on an underlying model for the SE dynamics. In this sense, the evolution it predicts is only as good as the initial model, which may only approximate the system's true behaviour. We will show how operationally determined quantities, encoded in transfer tensors, can be used to construct discretetime memory kernel equations for general open processes, where the SE model may not be well known, and show that these tend to a generalised Nakajima-Zwanzig equation in the continuum limit.
To do this, we will construct a difference equation for the state at time t, before taking a limit and recovering a derivative at that time. We first divide the interval [t 0 , t] into N parts, as above, and choose the N + 1 times {t k } to each be δt apart: t j+1 − t j = δt := (t − t 0 )/N [45]. By substituting Eq. (5) into the difference between ρ t = ρ t N and ρ t N −1 we get: which resembles, term by term, a discrete version of the memory kernel master equation given in Eq. (2). In the limit N → ∞ we can equate the terms as: In order to relate these quantities directly to the underlying dynamics, in terms of L SE t and ρ SE 0 , we can substitute Λ t:s ρ = tr E {U SE t:s (ρ ⊗ τ E s )} into the left hand sides of these expressions (writing T (N −j) t:tj and Ξ t:t0 in terms of dynamical maps). In this way, we can recast the objects in the decomposition in Eq. (4) in terms of SE projection superoperators as where P SE s ρ SE s = tr E {ρ SE s } ⊗ τ E s and Q SE s = I SE − P SE s , with τ E s the environment state appearing in the underlying representation of the dynamical map Λ t,s ; A is, as before, the operation used to prepare the system at the initial time. We can then take the limit by expanding U SE for small δt: where the propagator G t,s = T ← exp t s ds Q SE s L SE s is a time ordered exponential (as indicated by T ← ) and the value of the unit trace operator x E is unimportant, as it is acted on directly by a projector.
Equation (14) represents exactly the terms appearing in a generalised version of the Nakajima-Zwanzig master equation. In the standard derivation of the Nakajima-Zwanzig equation, a time-independent superoperator P SE and its complement Q SE = I SE − P SE are defined on the joint SE space [6]. In Appendix D, we show that one arrives at Eq. (14) following the standard derivation, but starting with the timedependent projection superoperators defined above. Unlike in earlier approaches utilising time-dependent projectors (see e.g. Refs. [46,47]) we do not place any restrictions on the operator τ E t beyond that it should be a positive, unit trace density operator, and that it and its first derivativeτ E t should be continuous functions of time. The significance of this result is that, even for time-inhomogeneous SE dynamics and initial correlations, the memory kernel can be operationally determined by reconstructing the short time-difference dynamics starting from a series of intermediate time points, though the convergence of the limit in Eqs. (12) may be slow, and will depend on the time difference t − t j in general.
Importantly, there is a freedom in choosing the environment operator τ E t that goes into the definition of the dynamical map Λ t,s and, hence, the projector P SE t . Different choices of projector P SE t lead to different forms for the superoperators in Eqs. (14) but the same solution for the dynamics ρ t . This is exemplified, using the model presented in Sec. III, in Fig. 3, where different (time-dependent) projectors are seen to lead to memory kernels with different norms; in this case, a particular time-dependent choice leads to a memory kernel whose norm decays several times faster than the most obvious time-independent choice (where τ E t = tr S ρ SE 0 ). For simulations, this means that a judicious choice of projector can lead to a much shorter memory cutoff time for the same degree of approximation in the dynamics. While τ E t can be freely chosen in a simulation (exact numerical method permitting), the dynamical maps reconstructed in an experimental setting have an environment operator dictated by the SE dynamics. However, as we will now discuss, there is still an analogous freedom in how the tomography is performed.

V. EXPERIMENTAL INTERPRETATION OF PROJECTOR CHOICE
One of the applications of the operational derivation of the Nakajima-Zwanzig equation, presented in the previous subsection, is to relate experimentally accessible quantities to properties of the underlying SE Hamiltonian and initial state (through Eqs. (12) & (14)). In this case, it is important to identify exactly which choices of environment state τ E t and hence dynamical map Λ t,s (or equivalently projector P SE t ) correspond to which experimental reconstruction procedure, if any. (a) After decorrelating the system from its environment (e.g., by projectively measuring and forgetting the outcome) at some intermediate time tj, the subsequent dynamics is described by a CPTP map Λt:t j . This can be operationally determined through normal quantum process tomography, i.e., by repreparing the system in various states ρ and determining the corresponding ρ at later times. In general, operations E = {Et k } could be performed on the system at a series of earlier times {t k < tj} (in addition to the initial preparation operation A), leading to different Λt:t j . (b) One choice of E which leads to a self-consistently solvable reconstructed master equation is Et k = Bσ k , where Bσ k ρ = σ k . (c) The map Λt:t j corresponds to joint evolution with an environment initially in the history dependent state ρ E t j , E .
In the simplest scenario, the joint SE evolution is left unperturbed before tomographically reconstructing Λ t:s (by discarding the state at time s and preparing a fresh one, before measuring the system at time t). The state of the environment at time s, in this case, would be the freely evolved reduced state: ρ E s = tr S {U SE s:t0 ρ SE t0 }. While this is perfectly reasonable from an operational standpoint, it is somewhat limiting from the perspective of extracting information about L SE t from the memory kernel. τ E t must be known to separate out properties of the generator from the memory kernel, and the time-dependent state of the environment is (usually) at least as difficult to determine as ρ t .
In order to get around this problem, one could choose to act on the system with a series of superoperators E = {E t k } (any physically allowed transformation of the system) between times t 0 and s, leading to different environment states ρ E s, E and, correspondingly, different families of intermediate dynamical maps Λ t:s (where we have left the E dependence implicit for notational convenience). This procedure is depicted in Fig. 4 and represents the operational analogue of the freedom to choose P SE t in Eqs. (14). A particularly convenient choice (from the perspective of calculating ρ E tj , E ) is the one where entanglement-breaking operations E tj ρ = B σt j ρ = σ tj are applied at times t j < s, where σ tj is some fixed state for each time step. In this case, the evolution of ρ E tj , E effectively decouples from the system in the limit that these operations are applied infinitely quickly. That is, the reduced state of the environment evolves according to an average generator: In an experiment, this could be achieved by, for instance, dynamically decoupling the system [48] or using the quantum Zeno effect to freeze the dynamics through frequent strong measurement [49]. In many cases, there will be choices of σ t such that ρ E t0 is a stationary state with respect to L E t , leading to time-independent projectors P SE and Q SE .

VI. DISCUSSION
In this Article, we have presented a derivation of the transfer tensor method, leading to a scheme for relating an operationally meaningful description of any open quantum process, in terms of completely positive dynamical maps, to a Nakajima-Zwanzig equation that depends on the underlying system-environment dynamics. We have also generalised the latter from the usual case to include time-dependent projectors. In addition to providing a fundamental connection between two different pictures of open dynamics, our result opens up the possibility for efficient simulation of the longtime evolution of driven systems or those with initial correlations.
Instead of directly solving the dynamics of a system to some long time t, dynamical maps can be reconstructed to a point where the memory kernel has decayed (from the end of the first driving period). These can then be used to calculate transfer tensors which will propagate the system to long times with an error that depends only on the accuracy of the initial reconstruction and the memory cutoff approximation used. This could be used to, for instance, determine dynamical steady states of driven systems [50] or probe signatures of many-body localisation [51].
Using Eq. (11), an experimentally or numerically reconstructed set of transfer tensors could also be used to create an approximate memory kernel through interpolation or form fitting. The resulting master equation would be guaranteed to give physically sensible solutions. Since its accuracy would depend on the smoothness of the memory kernel rather than, e.g., the strength of the SE coupling, it would provide an alternative to the usual perturbative approaches to approximate dynamics.
Finally, since the memory kernel contains information about the SE generator, its operational reconstruction will allow for the extraction of information about the underlying dynamics. That is, the scheme presented here could be used to probe an unknown environment by observing the dynamics of the system alone (cf. Refs. [25][26][27]). In fact, the full set of intermediate completely positive dynamical maps constitutes the maximum possible amount of dynamical information imprinted on the system without considering higher order multitime correlations (see Refs. [3,4]). The latter we will consider in a similar context in subsequent work.

APPENDICES Appendix A: Avoiding the inhomogeneous term
Many techniques for simulating open quantum systems, exactly or approximately, rely on a product assumption for the initial state: ρ SE t0 = ρ t0 ⊗ ρ E t0 . When this is not the case, there is always an inhomogeneous term in Eq. (2), and, while techniques have been developed to absorb it into the homogeneous part [52,53], these rely on knowledge of the underlying dynamics. This makes the explicit calculation of ρ t = M t,t0 [A], and hence Ξ t:t0 , problematic for an initially correlated SE, at least when the preparation A does not break those correlations.
However, when a preparation does break the correlations between system and environment, the subsequent dynamics can be simulated in the usual way (with the environment initially in state tr S {Aρ SE t0 }). We now use the property that the superchannel is linear, and that any preparation operation can be written as a linear combination A = α c α A (α) , where {A (α) } is a set of entanglement breaking completely-positive (but not necessarily trace-preserving) maps, and c α ∈ C [3]. This means that the time-evolved state for any preparation can be written as a linear combination of time-evolved states from some finite set of preparations where there are no initial correlations: ρ t,A = α c α ρ t,A (α) .
In the dilated picture, this is equivalent to writing the initial SE state as ρ SE t0 = α c α X (α) ⊗ τ (α) , where {X (α) } is a set of d 2 S linearly independent system operators. By solving for the dynamics with each of the d 2 S uncorrelated initial environment states τ (α) , the need to calculate the inhomogeneous term can be entirely circumvented.
Appendix C: Expanding transfer tensors in terms of system-environment projection superoperators In order to reduce notational clutter, we leave out the explicit SE label on superoperators P, Q, U and L in this section. Since we are assuming underlying dynamics such that Λ t:s ρ = tr E {U t:s ρ ⊗ τ E s }, we can write the transfer tensors in terms of system-environment quantities; the action of the superoperator in Eq. (6) becomes (C1) This can be simplified by introducing the time-dependent projection superoperators P t and Q t = I SE − P t , defined by the action P t X = tr E X ⊗ τ t . In terms of these, the transfer tensor acting on the system is U t:t k P t k U t k :t l P t l U t l :tj P tj ρ SE − . . .
where ρ SE is any state satisfying tr E {ρ SE } = ρ. Noting that M t:t0 [A] = tr E {U t:t0 Aρ SE t0 }, a similar expansion for the inhomogeneous term leads to (C3) In the limit considered in the main text, where t j+1 − t j = δt := (t − t 0 )/N and N is very large, the time evolution operator between two adjacent time points can be expanded in powers of δt as follows: ρ tr E P t (I + δtL t N −1 + . . . )Q t N −1 (I + δtL t N −2 + . . . )Q t N −2 . . . Q tj+1 (I + δtL tj + . . . )P tj ρ SE where we have used the fact that Q t Q s = Q s , P t Q s = 0 and Q t P s = P s − P t . In the limit that N → ∞, δt → 0 and the terms inside the central brackets resum to a time-ordered exponential (this limit is guaranteed to converge uniformly, as long as the generator is bounded and continuous [54]): ds Q s L s Q tj L tj P tj −Ṗ tj ρ SE = δt 2 K t:tj ρ.
The convergence of this equation, for the model presented in the main text, is exemplified in Fig. 5. Again, the same procedure can be performed for the inhomogeneous term in Eq. (C3), leading to t,t 0 for different values of t and N = t/δt, using the same model and parameters as in Fig. 3. Convergence is slower with N for longer evolution times, though, in this case, the same N corresponds to a longer δt. Note that both methods can be used to exactly calculate dynamics, even outside of the continuous limit where they converge.