Tomographically reconstructed master equations for any open quantum dynamics

TL;DR

This paper introduces a method to reconstruct memory kernel master equations for open quantum systems using quantum process tomography, enabling efficient simulation of complex dynamics without detailed system-environment information.

Contribution

The authors derive a transfer tensor method from first principles, allowing memory kernels to be expressed solely in terms of dynamical maps reconstructed via process tomography.

Findings

Reconstruction of memory kernels from dynamical maps is experimentally feasible.

The method simplifies long-time dynamics simulation for driven or initially correlated systems.

It offers computational advantages for simulating open quantum systems over traditional approaches.

Abstract

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.

Figures

Equations

$${\dot{\rho }}_t^{SE}={\mathcal{L}}_t^{SE}{\rho }_t^{SE}.$$

$${\dot{\rho }}_t={\mathcal{L}}_t{\rho }_t+{\int }_{t_0}^t\mathrm{d}s{\mathcal{K}}_{t,s}{\rho }_s+{\mathcal{J}}_{t,t_0},$$

$${\mathbf{P}}_s{\rho }_t={\Lambda }_{t:s}{\rho }_s\text{mbox}and{\mathbf{Q}}_s{\rho }_t={\rho }_t-{\Lambda }_{t:s}{\rho }_s,$$

$${\rho }_t=$$

$$=$$

$$T_{t:t_{N-1}}^{(1)}={\Lambda }_{t:t_{N-1}}\text{mbox}and$$

$$T_{t:t_j}^{(N-j)}={\Lambda }_{t:t_j}-\sum _{k=j+1}^N{T}_{t:t_k}^{(N-k)}{\Lambda }_{t_k:t_j}.$$

$$T_{t:t_j}^{(N-j)}=$$

$$+\sum _{k=j+2}^{N-1}\sum _{l=j+1}^{k-1}{\Lambda }_{t:t_k}{\Lambda }_{t_k:t_l}{\Lambda }_{t_l:t_j}-\dots .$$

$$T_{t_k:t_{k-l}}^{(l)}=T_{t_{k+c}:t_{k-l+c}}^{(l)},\forall l<k\text{mbox}and\forall t_{k-l}$$

$${\Vert {\rho }_{t_k}-{\tilde{\rho }}_{t_k}^{(m)}\Vert }_1$$

$$\lesssim \sum _{l=1}^m\Vert T_{t_{k-2m}\mathrm{mod}T+2t_m:t_{k-2m}\mathrm{mod}T+l\delta t}^{(2m-l)}\Vert ,$$

$${\mathcal{L}}_t^{SE}{\rho }_t^{SE}=$$

$$+\Gamma \left(L{\rho }_t^{SE}{L}^{†}-\frac{1}{2}\{L^{†}L,{\rho }_t^{SE}\}\right).$$

$$\frac{{\rho }_t-{\rho }_{t_{N-1}}}{\delta t}=\frac{({\Lambda }_{t:t_{N-1}}-\mathcal{I})}{\delta t}{\rho }_{t_{N-1}}+\sum _{j=0}^{N-2}\delta t\frac{T_{t:t_j}^{(N-j)}}{\delta t^2}{\rho }_{t_j}+\frac{{\Xi }_{t:t_0}}{\delta t},$$

$${\mathcal{K}}_{t,t_j}=\underset{N\to \infty }{\lim }\frac{T_{t:t_j}^{(N-j)}}{\delta t^2},{\mathcal{J}}_{t,t_0}=\underset{N\to \infty }{\lim }\frac{{\Xi }_{t:t_0}}{\delta t},$$

$${\mathcal{L}}_t=\underset{N\to \infty }{\lim }\frac{{\Lambda }_{t:t_{N-1}}-\mathcal{I}}{\delta t}.$$

$${\mathbf{P}}_s{\rho }_t={\mathrm{tr}}_E\{{\mathcal{U}}_{t:s}^{SE}{\mathcal{P}}_s^{SE}{\mathcal{U}}_{t:t_0}^{SE}\mathcal{A}{\rho }_{t_0}^{SE}\}\text{mbox}and$$

$${\mathbf{Q}}_s{\rho }_t={\mathrm{tr}}_E\{{\mathcal{U}}_{t:s}^{SE}{\mathcal{Q}}_s^{SE}{\mathcal{U}}_{t:t_0}^{SE}\mathcal{A}{\rho }_{t_0}^{SE}\},$$

$${\mathcal{L}}_t{\rho }_t={\mathrm{tr}}_E\left\{{\mathcal{P}}_t^{SE}{\mathcal{L}}_t^{SE}{\mathcal{P}}_t^{SE}{\rho }_t\otimes x^E\right\}$$

$${\mathcal{K}}_{t,s}{\rho }_s={\mathrm{tr}}_E\{{\mathcal{P}}_t^{SE}{\mathcal{L}}_t^{SE}{\mathcal{G}}_{t,s}$$

$$\times \left({\mathcal{Q}}_s^{SE}{\mathcal{L}}_s^{SE}{\mathcal{P}}_s^{SE}-{\dot{\mathcal{P}}}_s^{SE}\right){\rho }_s\otimes x^E\},$$

$${\mathcal{J}}_{t,t_0}={\mathrm{tr}}_E\left\{{\mathcal{P}}_t^{SE}{\mathcal{L}}_t^{SE}{\mathcal{G}}_{t,t_0}{\mathcal{Q}}_{t_0}^{SE}\mathcal{A}{\rho }_{t_0}^{SE}\right\},$$

$${\Vert {\rho }_{t_k}-{\tilde{\rho }}_{t_k}^{(m)}\Vert }_1=$$

$$\le$$

$${\Vert {\rho }_{t_k}-{\tilde{\rho }}_{t_k}^{(m)}\Vert }_1\lesssim$$

$${\Vert {\rho }_{t_k}-{\tilde{\rho }}_{t_k}^{(m)}\Vert }_1\le$$

$$+\mathcal{O}\left(\sum _{l=1}^m\Vert T_{(t_k-3t_m)\mathrm{mod}T+3t_m:(t_k-3t_m)\mathrm{mod}T+l\delta t}^{(3m-l)}\Vert \right)$$

$$\simeq$$

$$+\mathcal{O}\left({\int }_0^{t_m}\mathrm{d}s\Vert {\mathcal{K}}_{(t_k-3t_m)\mathrm{mod}T+3t_m,(t_k-3t_m)\mathrm{mod}T+s}\Vert \right)\},$$

$$T_{t:t_j}^{(N-j)}\rho =$$

$$+\sum _{k=j+2}^{N-1}\sum _{l=j+1}^{k-1}{\mathcal{U}}_{t:t_k}{\mathrm{tr}}_E\left\{{\mathcal{U}}_{t_k:t_l}{\mathrm{tr}}_E\left\{{\mathcal{U}}_{t_l:t_j}\rho \otimes {\tau }_{t_j}\right\}\otimes {\tau }_{t_l}\right\}\otimes {\tau }_{t_k}-\dots \}.$$

$$T_{t:t_j}^{(N-j)}\rho =$$

$$=$$

$${\Xi }_{t:t_0}={\mathrm{tr}}_E\left\{{\mathcal{P}}_t{\mathcal{U}}_{t:t_{N-1}}{\mathcal{Q}}_{t_{N-1}}{\mathcal{U}}_{t_{N-1}:t_{N-2}}{\mathcal{Q}}_{t_{N-2}}\dots {\mathcal{Q}}_{t_1}{\mathcal{U}}_{t_1:t_0}{\mathcal{Q}}_{t_0}\mathcal{A}{\rho }_{t_0}^{SE}\right\}.$$

$${\mathcal{U}}_{t_{j+1}:t_j}=$$

$$\simeq$$

$$T_{t:t_j}^{(N-j)}\rho \simeq$$

$$=$$

Full text

1 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, usually in terms of the time evolution of the system’s density operator. 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] (even going beyond density operator evolution to the full multi-time statistics of observables [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 long-time 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 $\rho^{\scriptscriptstyle{SE}}_{t_{0}}$, with a time evolution superoperator $\mathcal{U}^{\scriptscriptstyle{SE}}_{t:t_{0}}$111Superoperator $\mathcal{U}^{\scriptscriptstyle{SE}}_{t:t_{0}}$ need not be unitary in general; however, it must be divisible. that propagates solutions to an equation of the form

[TABLE]

All the physics of the system is encoded in the superoperator $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}$, which is usually taken to evolve $SE$ according to the von Neumann equation222That is, $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}=-\frac{i}{\hbar}[H^{\scriptscriptstyle{SE}}_{t},\rho^{\scriptscriptstyle{SE}}_{t}]$, where $H^{\scriptscriptstyle{SE}}_{t}$ is the $SE$ Hamiltonian., but could also incorporate dissipative and decohering dynamics through (possibly time-dependent) Lindblad terms [18, 19].

The second picture is the operational one, where the entire dynamics is described in terms of quantities that would 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 $\Lambda_{t:t_{0}}$, which is completely positive and trace preserving, and evolves the initial system state $\rho_{t_{0}}$ to its counterpart at time $t$: $\rho_{t}=\Lambda_{t:t_{0}}\rho_{t_{0}}$333Throughout this Article we denote operators and superoperators acting on the system-environment space with superscript $SE$. We do not put a superscript on operators and superoperators acting on the system alone.. 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 [20, 2] $\mathcal{M}_{t:t_{0}}$, which maps an initial preparation superoperator $\mathcal{A}$—any transformation that could be applied to the system at time $t_{0}$—to the state at later times: $\rho_{t}=\rho_{t,\mathcal{A}}=\mathcal{M}_{t:t_{0}}[\mathcal{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 $\mathcal{M}_{t:t_{0}}[\mathcal{I}]$, where $\mathcal{I}$ is the identity superoperator. The dynamics, $\Lambda_{t:t_{0}}$ or $\mathcal{M}_{t:t_{0}}$, can be reconstructed operationally by means of quantum process tomography [1, 21, 3].

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, $\rho_{t}={\rm tr}_{E}\rho^{\scriptscriptstyle{SE}}_{t}$). That is, recovering the underlying $SE$ quantities, $\rho^{\scriptscriptstyle{SE}}_{t_{0}}$ and $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}$, from operationally reconstructible time-dependent maps, such as $\mathcal{M}_{t:t_{0}}$, is neither uniquely constrained nor easy to achieve in practice; though, in the case of smooth dynamical maps, the underlying objects have been shown to be continuous and smooth themselves [22, 23]. Where such a recovery is possible, even partially, the system can be used as a probe of its environment, extracting information about its Hamiltonian directly [24, 25], or otherwise determining the parameters of a phenomenological master equation [26, 27, 28, 29].

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 other methods when the exact propagation from an underlying model is computationally complex [30, 31]. 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. 2 we derive transfer tensors from first principles for any open quantum dynamics, in particular those with initial correlations and generators that have arbitrary periodic or transient time dependence. We then go on, in Sec. 3 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. 4, 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. 5, 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 [18], as well as Gorini, Kossakowski and Sudarshan [19]. 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:

[TABLE]

where $\mathcal{L}_{t}$ is a time-local generator, $\mathcal{K}_{t,s}$ is a superoperator known as the memory kernel, and $\mathcal{J}_{t,t_{0}}$ is an inhomogeneous term, that can depend on the initial condition. Put simply, the memory kernel $\mathcal{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 [32, 33, 34, 35, 36]. 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. 4 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 [37, 38, 39]. However, there is currently no prescription for determining the elements of Eq. (2) when a microscopic model is unavailable. In Sec. 4 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.

2 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., $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}=\mathcal{L}^{\scriptscriptstyle{SE}}$; (ii) the initial $SE$ state is uncorrelated: $\rho^{\scriptscriptstyle{SE}}_{t_{0}}=\rho_{t_{0}}\otimes\rho^{\scriptscriptstyle{E}}_{t_{0}}$; (iii) $\rho^{\scriptscriptstyle{E}}_{t_{0}}$ is a stationary state. 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 $\{\mathbf{P}_{s},\mathbf{Q}_{s}\}$ which act as

[TABLE]

such that $(\mathbf{P}_{s}+\mathbf{Q}_{s})\rho_{t}=\rho_{t}$. Subsequent action of these objects further breaks up the evolution, i.e., $\mathbf{P}_{s^{\prime}}\Lambda_{t:s}=\Lambda_{t:s^{\prime}}\Lambda_{s^{\prime}:s}$. Here, $\Lambda_{t:s}$ is a completely positive dynamical map from time $s$ to time $t$, $\Lambda_{t:s}\rho={\rm tr}_{E}\{\mathcal{U}^{\scriptscriptstyle{SE}}_{t:s}(\rho\otimes\tau^{\scriptscriptstyle{E}}_{s})\}$, where, as in the introduction, $\mathcal{U}^{\scriptscriptstyle{SE}}_{t:s}$ is the time evolution operator of the underlying dynamics and $\tau^{\scriptscriptstyle{E}}_{s}$ is a unit trace time-dependent density operator not necessarily related to the actual environment dynamics; we will elaborate on the interpretation of $\Lambda_{t:s}$ and $\tau^{\scriptscriptstyle{E}}_{s}$ in numerical and experimental contexts in Secs. 3 and 5 respectively; for now, $\Lambda_{t:s}$ can be thought of simply as a two parameter family of dynamical maps.

Notably, $\mathbf{P}_{s}$ and $\mathbf{Q}_{s}$ are not superoperators, rather they can be seen as similar kinds of objects to the time ordering operator; that is, as consistently defined prescriptions for rewriting the subsequent expression based on its time indices. $\mathbf{P}_{s}$ replaces an object at any $t>s$ (either a dynamical map from $s^{\prime}<s$ to $t$ or the system state at time $t$) with that at time $s$ acted on by $\Lambda_{t,s}$; $\mathbf{Q}_{s}$ gives the complement. These can be seen simply as a shorthand for the decomposition $\rho_{t}=\Lambda_{t:s}\rho_{s}+(\rho_{t}-\Lambda_{t:s}\rho_{s})$ (and a corresponding decomposition for maps), which we will now iterate to further break up the dynamics. Specifically, by repeatedly using the identity $(\mathbf{P}_{s}+\mathbf{Q}_{s})=\mathcal{I}$, we can write the evolution up to time $t=t_{N}$ in terms of the evolution to $N-1$ earlier times $\{t_{j}\}$

[TABLE]

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

[TABLE]

Chasing this definition leads to the explicit form

[TABLE]

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 $\Xi_{t:t_{0}}=\rho_{t}-\sum_{j=0}^{N-1}T^{(N-j)}_{t:t_{j}}\rho_{t_{j}}$; this can be operationally determined for any initial preparation by reconstructing the superchannel $\mathcal{M}_{t:t_{0}}$. For long enough times, this term can be neglected [40]; 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 $\Lambda_{t:t_{0}}$ with the completely positive superchannel $\mathcal{M}_{t:t_{0}}$ [2]. We will now show explicitly how the transfer tensors derived in this section can be used to simulate driven open systems.

3 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 444Assuming that the short time exact dynamics is reconstructed accurately. 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] (of course, if memory effects persist beyond the time taken to converge on a steady state, then the numerically exact techniques will do just as well). In what follows, we will outline how a transfer tensor simulation could be performed and how the error can be bounded.

Fixing $\delta t=t_{j+1}-t_{j}$ determines the time resolution of the resulting dynamics, this is an external choice that does not affect the accuracy of the approach, but rather limits the degree to which fine-grained dynamical features can be resolved. To guarantee all such features are captured, one could choose it to be smaller than the time scale set by the largest energy in the system to be simulated. Once the temporal resolution is set, we will assume that the dynamical maps $\Lambda_{t:s}$ are periodic (or transient), such that $\Lambda_{t+T:s+T}=\Lambda_{t:s}$ for large enough $s$, and where the period $T=c\,\delta t$ with integer $c$. This is equivalent to requiring that the underlying generator $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}$ and reference state $\tau^{\scriptscriptstyle{E}}_{t}$ are transient or periodic with periods $a\,\delta t$ and $b\,\delta t$, where $c$ is the greatest common denominator of $1/a$ and $1/b$ ($\delta t$ and $\tau^{\scriptscriptstyle{E}}_{t}$ can always be chosen in a simulation to be compatible with the generator in this way). In this case, more general than that considered in Ref. [14], we have periodicity of the transfer tensors, such that

[TABLE]

sufficiently large that transient effects are negligible. Therefore, only transfer tensors $T^{(l)}_{t_{k}:t_{k-l}}$ with $k-l<c$ need to be calculated explicitly (in the transient case, all transfer tensors starting within the transient period must be computed).

Furthermore, in most real open systems, the memory kernel $\mathcal{K}_{t,s}$ decreases in magnitude (often exponentially) with the time difference $t-s$ [6]. That is, the system’s state at one time is continuously forgotten until it no longer affects the future evolution. Consequently, the transfer tensors $T^{(l)}_{t_{k}:t_{k-l}}$ also become negligibly small for sufficiently large $l$. Similarly, the influence of the initial state encoded in $\mathcal{J}_{t,t_{0}}$ or $\Xi_{t,t_{0}}$ commonly decays with $t-t_{0}$. What this means is that, beyond a certain time $t_{m}=m\,\delta t$, memory effects can be neglected. That is, $T^{(l)}_{t_{k}:t_{k-l}}$ with $l>m$ can be effectively set to zero without substantially affecting the dynamics. Combined with the time translation symmetry described above, this means the dynamics at any time, given the finite resolution set by $\delta t$, can be reconstructed using only a finite set of $cm$ transfer tensors ($m$ transfer tensors for each of the $c$ starting points in a driving period). As such, decreasing $\delta t$ while leaving other time scales fixed increases the computational cost quadratically.

In the first step of the simulation, numerical process tomography is performed by exactly solving the dynamics to a time $t_{m}$ beyond each time $x\delta t$, $x<c$, in the first driving period (or the transient period) for an informationally complete set of initial states. These dynamics can be reconstructed for any choice of environment state $\tau^{\scriptscriptstyle{E}}_{t}$ permitted by the exact technique used to solve for the short time dynamics (though different choices can lead to longer or shorter memory times, as discussed in Sec. 5). The resulting set of dynamical maps is then used to reconstruct the transfer tensors on this interval (starting with those between adjacent time steps and iterating with Eq. (6)). Importantly, despite the fact the transfer tensors are constructed starting from different points of the driving cycle, the exact dynamics only needs to be calculated for an evolution time $t_{m}$ from each point, meaning the efficiency of the technique does not depend on the driving period. Using Eq. (5), the state at any arbitrary time $t_{k}$ can be found by identifying $T^{(l)}_{t_{k}:t_{k-l}}$ with $T^{(l)}_{t_{k-l}\,\mathrm{mod}\,T+l\delta t:t_{k-l}\,\mathrm{mod}\,T}$, for $l\leq m$, and with [math], for $l>m$ (similarly, $\Xi_{t_{k},t_{0}}$ is set to [math] 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 $\tilde{\rho}_{t_{k}}$ and the exact state $\rho_{t_{k}}$ can be bounded as

[TABLE]

where $\|X\|:=\max_{\rho,\sigma}{\rm tr}\{\sigma X\rho\}$ 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: ${\rm max}_{x\leq c}\|T^{(m)}_{x\delta t+t_{m}:x\delta 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:

[TABLE]

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 $H^{\scriptscriptstyle{SE}}_{t}=\frac{\omega}{2}\sigma_{z}^{\scriptscriptstyle{S}}+\frac{\omega^{\prime}}{2}\sigma_{z}^{\scriptscriptstyle{E}}+g[\sigma_{x}^{\scriptscriptstyle{S}}\otimes\sigma_{x}^{\scriptscriptstyle{E}}+\cos(\Omega t)\sigma_{y}^{\scriptscriptstyle{S}}\otimes\sigma_{y}^{\scriptscriptstyle{E}}]$, and the Lindblad operators $L=\mathbbm{1}^{\scriptscriptstyle{S}}\otimes|0\rangle\!\langle 1|^{\scriptscriptstyle{E}}$ act on the environment alone. This generates dissipative dynamics for the environment; incoherently pumping it to the excited state with rate $\Gamma$. This seemingly simple model has many of the features of more complex open systems; in fact, it is approximately equivalent to a class of spin-boson models [41].

The simulated dynamics, with correlated initial state $\rho^{{\scriptscriptstyle{SE}}}_{0}=\frac{3}{4}|0\rangle\!\langle 0|\otimes|+\rangle\!\langle+|+\frac{1}{4}|1\rangle\!\langle 1|\otimes|-\rangle\!\langle-|$, where $\mathinner{|{\pm}\rangle}=(\mathinner{|{0}\rangle}\pm\mathinner{|{1}\rangle})/\sqrt{2}$, is shown in Fig. 2 along with the accompanying error, quantified by the trace distance $\|\rho_{t_{k}}-\tilde{\rho}_{t_{k}}\|_{1}={\rm tr}\{|\rho_{t_{k}}-\tilde{\rho}_{t_{k}}|\}$ (in all figures the parameters are set as $\omega=\Gamma=\omega^{\prime}/16=g/2=\Omega/2$).

In the figure, short time exact dynamics, calculated between the first nine time steps (using a time-independent reference state $\tau^{\scriptscriptstyle{E}}={\rm tr}_{\scriptscriptstyle{S}}\rho_{0}^{\scriptscriptstyle{SE}}$), is used to propagate the state $\tilde{\rho}_{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 ($\omega 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 $\delta 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 555We could have also reduced the error by choosing a time-dependent reference state (since the memory kernel norm decays much faster). However, in this case, a longer memory time would need to be chosen to allow transients in the time-dependent state to decay.. 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. However, beyond a certain point, even this loose bound on the error seems to decrease exponentially with $t_{m}$. In the limit of continuous time dynamics, our approach leads us to a generalised version of the Nakajima-Zwanzig equation, of the form of Eq. (2), as we will now show.

4 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 discrete-time 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 $\delta t$ apart: $t_{j+1}-t_{j}=\delta t:=(t-t_{0})/N$ 666Equidistance between time steps is not necessary for our derivation and is only chosen for convenience. It may be that other choices lead to faster convergence in practice.. By substituting Eq. (5) into the difference between $\rho_{t}=\rho_{t_{N}}$ and $\rho_{t_{N-1}}$ we get:

[TABLE]

which resembles, term by term, a discrete version of the memory kernel master equation given in Eq. (2). In the limit $N\to\infty$ we can equate the terms as:

[TABLE]

In order to relate these quantities directly to the underlying dynamics, in terms of $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}$ and $\rho_{0}^{\scriptscriptstyle{SE}}$, we can substitute $\Lambda_{t:s}\rho={\rm tr}_{E}\{\mathcal{U}^{\scriptscriptstyle{SE}}_{t:s}(\rho\otimes\tau^{\scriptscriptstyle{E}}_{s})\}$ into the left hand sides of these expressions (writing $T^{(N-j)}_{t:t_{j}}$ and $\Xi_{t:t_{0}}$ 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

[TABLE]

where $\mathcal{P}^{\scriptscriptstyle{SE}}_{s}\rho^{\scriptscriptstyle{SE}}_{s}={\rm tr}_{E}\{\rho^{\scriptscriptstyle{SE}}_{s}\}\otimes\tau_{s}^{\scriptscriptstyle{E}}$ and $\mathcal{Q}^{\scriptscriptstyle{SE}}_{s}=\mathcal{I}^{\scriptscriptstyle{SE}}-\mathcal{P}^{\scriptscriptstyle{SE}}_{s}$, with $\tau_{s}^{\scriptscriptstyle{E}}$ the environment state appearing in the underlying representation of the dynamical map $\Lambda_{t,s}$; $\mathcal{A}$ is, as before, the operation used to prepare the system at the initial time. We can then take the limit by expanding $\mathcal{U}^{\scriptscriptstyle{SE}}$ for small $\delta t$: $\mathcal{U}^{\scriptscriptstyle{SE}}_{t+\delta t:t}\simeq\mathcal{I}^{\scriptscriptstyle{SE}}+\delta t\mathcal{L}^{\scriptscriptstyle{SE}}_{t}+\delta t^{2}{\mathcal{L}^{\scriptscriptstyle{SE}}_{t}}^{2}+\dots$. We do this in detail in Appendix C, finding

[TABLE]

where the propagator $\mathcal{G}_{t,s}=T_{\leftarrow}\exp\left[\int_{s}^{t}\rm{d}s^{\prime}\,\mathcal{Q}^{\scriptscriptstyle{SE}}_{s^{\prime}}\mathcal{L}^{\scriptscriptstyle{SE}}_{s^{\prime}}\right]$ is a time ordered exponential (as indicated by $T_{\leftarrow}$) and the value of the unit trace operator $x^{\scriptscriptstyle{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 $\mathcal{P}^{\scriptscriptstyle{SE}}$ and its complement $\mathcal{Q}^{\scriptscriptstyle{SE}}=\mathcal{I}^{\scriptscriptstyle{SE}}-\mathcal{P}^{\scriptscriptstyle{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 time-dependent projection superoperators defined above. Unlike in earlier approaches utilising time-dependent projectors for open systems (see e.g. Refs. [42, 43]) we do not place any restrictions on the operator $\tau_{t}^{\scriptscriptstyle{E}}$ beyond that it should be a positive, unit trace density operator, and that it and its first derivative $\dot{\tau}_{t}^{\scriptscriptstyle{E}}$ should be continuous functions of time. It is worth noting the similarity to recent projection operator approaches to describing non-linear many-body dynamics [44, 45].

The significance of our 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, alleviating the need to explicitly calculate or approximate the time-ordered exponential integrals implicit in Eq. (14). However, 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 $\tau_{t}^{\scriptscriptstyle{E}}$ that goes into the definition of the dynamical map $\Lambda_{t,s}$ and, hence, the projector $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$. Different choices of projector $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$ lead to different forms for the superoperators in Eqs. (14) but the same solution for the dynamics $\rho_{t}$. This is exemplified, using the model presented in Sec. 3, 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 $\tau_{t}^{\scriptscriptstyle{E}}={\rm tr}_{\scriptscriptstyle{S}}\rho_{0}^{\scriptscriptstyle{SE}}$). 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 $\tau_{t}^{\scriptscriptstyle{E}}$ 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.

5 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 $\tau^{\scriptscriptstyle{E}}_{t}$ and hence dynamical map $\Lambda_{t,s}$ (or equivalently projector $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$) correspond to which experimental reconstruction procedure, if any.

In the simplest scenario, the joint $SE$ evolution is left unperturbed before tomographically reconstructing $\Lambda_{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: $\rho^{\scriptscriptstyle{E}}_{s}={\rm tr}_{S}\{\mathcal{U}^{\scriptscriptstyle{SE}}_{s:t_{0}}\rho^{\scriptscriptstyle{SE}}_{t_{0}}\}$. While this is perfectly reasonable from an operational standpoint, it is somewhat limiting from the perspective of extracting information about $\mathcal{L}^{\scriptscriptstyle{SE}}_{t}$ from the memory kernel. $\tau^{\scriptscriptstyle{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 $\rho_{t}$.

In order to get around this problem, one could choose to act on the system with a series of superoperators $\vec{\mathcal{E}}=\{\mathcal{E}_{t_{k}}\}$ (any physically allowed transformation of the system) between times $t_{0}$ and $s$, leading to different environment states $\rho^{\scriptscriptstyle{E}}_{s,\vec{\mathcal{E}}}$ and, correspondingly, different families of intermediate dynamical maps $\Lambda_{t:s}$ (where we have left the $\vec{\mathcal{E}}$ dependence implicit for notational convenience). This procedure is depicted in Fig. 4 and represents the operational analogue of the freedom to choose $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$ in Eqs. (14).

A particularly convenient choice (from the perspective of calculating $\rho^{E}_{t_{j},\vec{\mathcal{E}}}$) is the one where entanglement-breaking operations $\mathcal{E}_{t_{j}}\rho=\mathcal{B}_{\sigma_{t_{j}}}\rho=\sigma_{t_{j}}$ are applied at times $t_{j}<s$, where $\sigma_{t_{j}}$ is some fixed state for each time step. In this case, the evolution of $\rho^{\scriptscriptstyle{E}}_{t_{j},\vec{\mathcal{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: $\dot{\rho}^{\scriptscriptstyle{E}}_{t}={\rm tr}_{S}\left\{\mathcal{L}^{\scriptscriptstyle{SE}}_{t}\sigma_{t}\otimes\rho^{\scriptscriptstyle{E}}_{t}\right\}=\mathcal{L}^{\scriptscriptstyle{E}}_{t}\rho^{\scriptscriptstyle{E}}_{t}$. In an experiment, this could be achieved by, for instance, dynamically decoupling the system [46] or using the quantum Zeno effect to freeze the dynamics through frequent strong measurement [47], though these may be difficult in practice and represent only two possible choices of control that could lead to a known $\rho^{E}_{t_{j},\vec{\mathcal{E}}}$. In many cases, there will be choices of $\sigma_{t}$ such that $\rho^{\scriptscriptstyle{E}}_{t_{0}}$ is a stationary state with respect to $\mathcal{L}^{\scriptscriptstyle{E}}_{t}$, leading to time-independent projectors $\mathcal{P}^{\scriptscriptstyle{SE}}$ and $\mathcal{Q}^{\scriptscriptstyle{SE}}$.

6 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 long-time 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 [48] or probe signatures of many-body localisation [49].

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. [26, 24, 25]). 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.

Acknowledgments—The authors would like to thank Jared Cole, Guy Cohen, and César Rodríguez-Rosario for stimulating discussions regarding this work. KM is supported through ARC FT160100073.[50, 51]

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: $\rho^{\scriptscriptstyle{SE}}_{t_{0}}=\rho_{t_{0}}\otimes\rho^{\scriptscriptstyle{E}}_{t_{0}}$. 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 [50, 51], these rely on knowledge of the underlying dynamics. This makes the explicit calculation of $\rho_{t}=\mathcal{M}_{t,t_{0}}[\mathcal{A}]$, and hence $\Xi_{t:t_{0}}$, problematic for an initially correlated $SE$, at least when the preparation $\mathcal{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 ${\rm tr}_{S}\{\mathcal{A}\rho^{\scriptscriptstyle{SE}}_{t_{0}}\}$). We now use the property that the superchannel is linear, and that any preparation operation can be written as a linear combination $\mathcal{A}=\sum_{\alpha}c_{\alpha}\mathcal{A}^{(\alpha)}$, where $\{\mathcal{A}^{(\alpha)}\}$ is a set of entanglement breaking completely-positive (but not necessarily trace-preserving) maps, and $c_{\alpha}\in\mathbb{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: $\rho_{t,\mathcal{A}}=\sum_{\alpha}c_{\alpha}\rho_{t,\mathcal{A}^{(\alpha)}}$.

In the dilated picture, this is equivalent to writing the initial $SE$ state as $\rho_{t_{0}}^{\scriptscriptstyle{SE}}=\sum_{\alpha}c_{\alpha}X^{(\alpha)}\otimes\tau^{(\alpha)}$, where $\{X^{(\alpha)}\}$ is a set of $d_{\scriptscriptstyle{S}}^{2}$ linearly independent system operators. By solving for the dynamics with each of the $d_{\scriptscriptstyle{S}}^{2}$ uncorrelated initial environment states $\tau^{(\alpha)}$, the need to calculate the inhomogeneous term can be entirely circumvented.

Appendix B Simulation error

Since the expression for the time-evolved density operator in Eq. (5) is exact, errors in a simulation using transfer tensors can arise only if terms in the expansion are neglected, or otherwise from errors in the reconstruction of the maps $\Lambda_{t:s}$. Here, we assume the reconstruction is accurate, and determine the error associated with neglecting memory effects beyond some cutoff time. In what follows, we will quantify the ‘size’ of an operator with the trace norm $\|X\|_{1}={\rm tr}\{|X|\}$, where $|X|$ is the operator formed by taking the magnitudes of the singular values of $X$ and the same singular vectors. For superoperators, we will use $\|\mathcal{X}\|:=\max_{A,B}\{{\rm tr}[A(\mathcal{X}B)]\}$, which is the largest singular value (or operator norm) of the matrix representation of $\mathcal{X}$ which acts on vectorised operators.

The error in the final state $\tilde{\rho}^{(m)}_{t_{k}}$ due to neglecting memory effects beyond a time $t_{m}=m\delta t$ is given by

[TABLE]

where we have used the triangle inequality and the definition of superoperator norm given above. Referring back to Eq. (12), we can see that, for small $\delta t$ (and $t_{m}\gg\delta t$),

[TABLE]

That is, the error decreases with both the length of the time step $\delta t$ and the size of the tail of the memory kernel (beyond the considered memory time $t_{m}$).

Importantly, under the assumption that the decay over one memory time $t_{m}$ (beyond the first) is significant for both $\left\|T^{(l)}_{t_{k}:t_{k-l}}\right\|$ and $\left\|\Xi_{t_{k}:t_{0}}\right\|_{1}$ (or, in the small $\delta t$ limit, $\|\mathcal{K}_{t_{k},{t_{k}-s}}\|$ and $\|\mathcal{J}_{t_{k},t_{0}}\|_{1}$), we have that

[TABLE]

which at leading order does not depend on $t_{k}$ itself, but rather the phase with respect to the driving cycle, $(t_{k}-2t_{m})\,\mathrm{mod}\,T$. In other words, beyond a certain point, the error does not grow with the evolution time.

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 $\mathcal{P}$, $\mathcal{Q}$, $\mathcal{U}$ and $\mathcal{L}$ in this section. Since we are assuming underlying dynamics such that $\Lambda_{t:s}\rho={\rm tr}_{E}\{\mathcal{U}_{t:s}\rho\otimes\tau^{\scriptscriptstyle{E}}_{s}\}$, we can write the transfer tensors in terms of system-environment quantities; the action of the superoperator in Eq. (6) becomes

[TABLE]

This can be simplified by introducing the time-dependent projection superoperators $\mathcal{P}_{t}$ and $\mathcal{Q}_{t}=\mathcal{I}^{{\scriptscriptstyle{SE}}}-\mathcal{P}_{t}$, defined by the action $\mathcal{P}_{t}X={\rm tr}_{E}X\otimes\tau_{t}$. In terms of these, the transfer tensor acting on the system is

[TABLE]

where $\rho^{\scriptscriptstyle{SE}}$ is any state satisfying ${\rm tr}_{E}\{\rho^{\scriptscriptstyle{SE}}\}=\rho$. Noting that $\mathcal{M}_{t:t_{0}}[\mathcal{A}]={\rm tr}_{E}\{\mathcal{U}_{t:t_{0}}\mathcal{A}\rho^{\scriptscriptstyle{SE}}_{t_{0}}\}$, a similar expansion for the inhomogeneous term leads to

[TABLE]

In the limit considered in the main text, where $t_{j+1}-t_{j}=\delta 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 $\delta t$ as follows:

[TABLE]

Substituting this into Eq. (19), we have

[TABLE]

where we have used the fact that $\mathcal{Q}_{t}\mathcal{Q}_{s}=\mathcal{Q}_{s}$, $\mathcal{P}_{t}\mathcal{Q}_{s}=0$ and $\mathcal{Q}_{t}\mathcal{P}_{s}=\mathcal{P}_{s}-\mathcal{P}_{t}$. In the limit that $N\rightarrow\infty$, $\delta t\rightarrow 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 [52]):

[TABLE]

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. (20), leading to

[TABLE]

Appendix D Nakajima-Zwanzig with time-dependent projection operators

Here, we briefly derive the Nakajima Zwanzig equation given by Eqs. (2) & (14) from the underlying equation of motion given in Eq. (1). Starting from the definition of the projection superoperators $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$ and $\mathcal{Q}^{\scriptscriptstyle{SE}}_{t}=\mathcal{I}^{\scriptscriptstyle{SE}}-\mathcal{P}^{\scriptscriptstyle{SE}}_{t}$ (where $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}X^{\scriptscriptstyle{SE}}={\rm tr}_{E}\{X^{\scriptscriptstyle{SE}}\}\otimes\tau_{t}^{\scriptscriptstyle{E}}$), we first note the following identities:

[TABLE]

Next, we consider the equations of motion of the so-called ‘relevant’ and ‘irrelevant’ parts of the $SE$ density operator, $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}$ and $\mathcal{Q}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}$ respectively. Using Eq. (1), we find

[TABLE]

and

[TABLE]

where we have used that $\mathcal{I}^{\scriptscriptstyle{SE}}=\mathcal{P}^{\scriptscriptstyle{SE}}_{t}+\mathcal{Q}^{\scriptscriptstyle{SE}}_{t}$. Formally solving Eq. (30) gives

[TABLE]

As in the main text, we have taken the post-preparation $SE$ state $\mathcal{A}\rho^{\scriptscriptstyle{SE}}_{t_{0}}$ as our initial condition. Substituting this into Eq. (29) and taking the partial trace over $E$ (since $\dot{\rho}_{t}={\rm tr}_{E}\left\{{\rm d}/{\rm d}t\mathcal{P}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}\right\}$) leads to

[TABLE]

In the first line, the derivative-dependent term vanishes, since ${\rm tr}\{\dot{\tau}^{\scriptscriptstyle{E}}_{t}\}=0$, and hence ${\rm tr}_{E}\{\dot{\mathcal{P}}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}\}=0$, when $\tau^{\scriptscriptstyle{E}}_{t}$ is constrained to be a unit trace density operator (as we will always assume). Since $\mathcal{P}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}=\mathcal{P}^{\scriptscriptstyle{SE}}_{t}\rho_{t}\otimes x^{\scriptscriptstyle{E}}$ and $\dot{\mathcal{P}}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}=\dot{\mathcal{P}}^{\scriptscriptstyle{SE}}_{t}\rho_{t}\otimes x^{\scriptscriptstyle{E}}$ for any unit trace operator $x^{\scriptscriptstyle{E}}$, we can identify each line of Eq. (D) with one of the terms in Eqs. (14).

Our equation differs from that in Refs. [42, 43], since we do not assume that $\dot{\mathcal{P}}^{\scriptscriptstyle{SE}}_{t}\rho^{\scriptscriptstyle{SE}}_{t}=0$ $\forall t$. If we were to do so, then the derivative term in the last line of Eq. (D) would vanish, and we would recover the form derived in the aforementioned references.