Completely positive master equation for arbitrary driving and small level spacing

TL;DR

This paper introduces a generalized, completely positive Markovian master equation suitable for arbitrary driving and small level spacing, expanding the applicability of open quantum system modeling.

Contribution

It derives a time-dependent coarse-grained master equation from the Redfield equation, valid for fast driving and small level spacings, with rigorous error bounds.

Findings

The new master equation is completely positive and locally generated.

It applies to arbitrarily fast time-dependent Hamiltonians.

Dynamical decoupling can extend coherence times in a Markovian framework.

Abstract

Markovian master equations are a ubiquitous tool in the study of open quantum systems, but deriving them from first principles involves a series of compromises. On the one hand, the Redfield equation is valid for fast environments (whose correlation function decays much faster than the system relaxation time) regardless of the relative strength of the coupling to the system Hamiltonian, but is notoriously non-completely-positive. On the other hand, the Davies equation preserves complete positivity but is valid only in the ultra-weak coupling limit and for systems with a finite level spacing, which makes it incompatible with arbitrarily fast time-dependent driving. Here we show that a recently derived Markovian coarse-grained master equation (CGME), already known to be completely positive, has a much expanded range of applicability compared to the Davies equation, and moreover, is locally generated and can be generalized to accommodate arbitrarily fast driving. This generalization, which we refer to as the time-dependent CGME, is thus suitable for the analysis of fast operations in gate-model quantum computing, such as quantum error correction and dynamical decoupling. Our derivation proceeds directly from the Redfield equation and allows us to place rigorous error bounds on all three equations: Redfield, Davies, and coarse-grained. Our main result is thus a completely positive Markovian master equation that is a controlled approximation to the true evolution for any time-dependence of the system Hamiltonian, and works for systems with arbitrarily small level spacing. We illustrate this with an analysis showing that dynamical decoupling can extend coherence times even in a strictly Markovian setting.

Equations

$$\text{system-bath coupling strength}\times \text{leavevmode}\sqrt{\text{bath correlation time}/\text{system level-spacing}}.$$

$$\text{system-bath coupling strength}\times \text{bath correlation time},$$

$$H_{\text{tot}}=H\otimes I_{\text{b}}+A\otimes B+I\otimes H_{\text{b}}.$$

$$A_{\omega }=\sum _{mn:E_{mn}=\omega }{A}_{nm}|n⟩⟨m|=\sum _{mn:E_{mn}=\omega }{\Pi }_{n}A{\Pi }_m=A_{-\omega }^{†},$$

$$A(t)=\sum _{\omega }{A}_{\omega }{e}^{-i\omega t}.$$

$$\mathcal{C}(t)$$

$$\gamma (\omega )$$

$$f(\omega )={\int }_0^{\infty }C(t)e^{i\omega t}dt=\frac{1}{2}\gamma (\omega )+iS(\omega ),S(\omega )=\frac{1}{2i}\left[f(\omega )-f^{\ast }(\omega )\right]=S^{\ast }(\omega ).$$

$$S(\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathcal{P}\left(\frac{1}{\omega -{\omega }^{\prime }}\right)[\gamma ({\omega }^{\prime })]d{\omega }^{\prime },$$

$$C(t)=\mathrm{Tr}[{\rho }_{\text{b}}B(0)B(t+i\beta )].$$

$$\gamma (-\omega )=e^{-\beta \omega }\gamma (\omega )\forall \omega \ge 0.$$

$$\frac{1}{{\tau }_{SB}}$$

$${\tau }_B$$

$$\gamma (\omega )=|\gamma (\omega )|\le {\int }_{-\infty }^{\infty }|C(t)|dt=2{\int }_0^{\infty }|C(t)|dt=\frac{2}{{\tau }_{SB}}.$$

$$|{\gamma }^{\prime }(\omega )|\le {\int }_{-\infty }^{\infty }|t||C(t)|dt=2{\int }_0^{\infty }t|C(t)|dt\stackrel{T\to \infty }{=}2\frac{{\tau }_B}{{\tau }_{SB}},$$

$${\gamma }^{\prime }(0)=\frac{1}{2}\beta \gamma (0).$$

$$\beta \gamma (0)\le 4\frac{{\tau }_B}{{\tau }_{SB}}.$$

$${\dot{\rho }}_R(t)=-i[H,{\rho }_R(t)]+(A{\rho }_R{A}_f-{\rho }_R{A}_{f}A+\text{h.c.}),$$

$$A_f={\int }_0^{\infty }\mathcal{C}(-t)A(-t)dt=\sum _{\omega }{A}_{\omega }{\int }_0^{\infty }\mathcal{C}(-t)e^{i\omega t}dt=\sum _{\omega }{f}^{\ast }(-\omega )A_{\omega },$$

$${\dot{\rho }}_D(t)$$

$$H_{\mathrm{LS}}$$

$${\dot{\rho }}_C(t)=-i[H+H_{\mathrm{LS}},{\rho }_C(t)]+{\int }_{-\infty }^{\infty }dϵ\left(A_{ϵ}{\rho }_C(t)A_{ϵ}^{†}-\frac{1}{2}\left\{{\rho }_C(t),A_{ϵ}^{†}{A}_{ϵ}\right\}\right),$$

$$A_{ϵ}$$

$$=\sum _{\omega }f(ϵ,\omega )A_{\omega },f(ϵ,\omega )=\sqrt{\frac{\gamma (ϵ)T_a}{2\pi }}\text{sinc}[T_a(ϵ-\omega )/2],$$

$$H_{\mathrm{LS}}$$

$$F_{\omega {\omega }^{\prime }}$$

$${\dot{\rho }}_C(t)=-i[H(t)+H_{\mathrm{LS}}(t),{\rho }_C]+{\int }_{-\infty }^{\infty }dϵ\left(A_{ϵ}(t){\rho }_C{A}_{ϵ}^{†}(t)-\frac{1}{2}\left\{{\rho }_C,A_{ϵ}^{†}(t)A_{ϵ}(t)\right\}\right),$$

$$A_{ϵ}(t)=\sqrt{\frac{\gamma (ϵ)}{2\pi T_a}}{\int }_{-T_a/2}^{T_a/2}{e}^{iϵt_1}A(t+t_1,t)d{t}_1,$$

$$H_{\mathrm{LS}}(t)=\frac{i}{2T_a}{\int }_{-T_a/2}^{T_a/2}d{t}_1{\int }_{-T_a/2}^{T_a/2}d{t}_2\mathrm{sgn}(t_1-t_2)\mathcal{C}(t_2-t_1)A(t+t_2,t)A(t+t_1,t).$$

$$T_a=\sqrt{{\tau }_{SB}{\tau }_B/5}.$$

$$\Vert {\rho }_{\text{true}}(t)-{\rho }_R(t){\Vert }_1\le O\left(\frac{{\tau }_B}{{\tau }_{SB}}{e}^{12t/{\tau }_{SB}}\right)\text{ln}\left(\frac{{\tau }_{SB}}{{\tau }_B}\right),$$

$$\Vert {\rho }_{\text{true}}(t)-{\rho }_D(t){\Vert }_1\le O\left(\left(\frac{{\tau }_B}{{\tau }_{SB}}+\frac{1}{\sqrt{{\tau }_{SB}\delta E}}\right)e^{12t/{\tau }_{SB}}\right),$$

$$\Vert {\rho }_{\text{true}}(t)-{\rho }_C(t){\Vert }_1\le O\left(\sqrt{\frac{{\tau }_B}{{\tau }_{SB}}}{e}^{6t/{\tau }_{SB}}\right).$$

Full text

1 Introduction

Modeling experiments requires taking into account that physical systems are open, i.e., not ideally isolated from their environments. Usually an environment contains many more degrees of freedom than the system itself, but is not in any interesting phase of matter where the detailed modeling of each of those degrees of freedom is required. This makes the problem of modeling open systems tractable. More precisely, the problem is tractable under a list of assumptions on the parameters of the bath (environment) and its interaction with the system [1, 2, 3].

One natural approach is to consider the case when this interaction is weak compared to all other energy scales in the problem. In this limit, Davies showed [4] that the dynamics of the system are given by a Markovian master equation with what is now called Davies generators. Rigorous bounds on the error between the solution of this equation and the true evolution can be obtained [presented here in Eq. (264)]. The largest contribution to the error comes from making the rotating wave approximation (RWA).111This approximation is stated explicitly in Ref. [3], page 87, around Eq. (3.6.67). After the system evolves for a time set by a relevant relaxation timescale, the error in its density matrix is given by:

[TABLE]

This quantity, in general, becomes exponentially large in the system size due its denominator, so that the coupling strength or the bath correlation time need to be sufficiently small for the use of Davies generators to be justified, up to some upper system size limit. In practice, for convenience and also because it may apply in a range that is larger than the rigorous but conservative bounds imply, the Lindblad master equation [5] with Davies generators is often used outside its formal range of applicability in modeling of experiments (e.g., Ref. [6, 7, 8, 9]). In such cases it should be understood as a semi-empirical model of open system dynamics, that possesses all the physical properties of the dynamics without reproducing them exactly. Unfortunately, since we naturally wish to apply modeling tools to problems for which we do not know a priori what the correct answer is, this approach cannot provide correctness guarantees. However, other methods are even more susceptible to this problem, e.g., a path-integral (non-master equation) approach based on integrating out the bath [10, 11, 12]: this involves integration over a very high-dimensional space, so methods [13] that bring it back into the realm of numerical feasibility usually involve uncontrollable approximations.

One of the important requirements of open system dynamics is complete positivity [14], under which density matrices whose eigenvalues are by definition non-negative are mapped to other such matrices, even when applied to a subsystem of a larger system.222It should be noted that complete positivity can be relaxed without losing physicality; see, e.g., Ref. [15]. This property is equivalent to the Markovian master equation being in Lindblad canonical form [5, 16]. There are numerous other master equations (e.g., Refs. [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]), some in Lindblad form and some not, whose domains of validity and ranges of applicability were not always studied in detail. One promising candidate is a coarse-grained master equation (CGME) [40], which is in Lindblad form and thus automatically completely positive (CP). In this work we show that a suitably generalized variant of this equation, that includes an arbitrarily time-dependent system Hamiltonian (hence called the “time-dependent CGME”), has a much milder error than (1):

[TABLE]

which does not diverge with the system size $n$ for local observables, and has only a polynomial [$O(n^{\alpha})$ with $1/2\leq\alpha\leq 1$] prefactor for nonlocal observables. We note that the full form of the error estimate (2) is not uniform over the evolution time, unlike some of the error bounds developed elsewhere [39, 41]. We also show that this new master equation is an improvement over the Lindblad equation with Davies generators in that it has a local structure, compatible with approximate simulation methods with asymptotically smaller costs [$O(n)$ vs. $O(\exp(n))$], and that the error estimate (2) is valid for an arbitrary time-dependence of the system Hamiltonian. In particular, the time-dependent CGME is compatible with the assumption of arbitrarily fast gates, often made in the circuit model of quantum computing, e.g., in the analysis of fault-tolerant quantum error correction [42]. This assumption was shown to be incompatible with the derivation of the Davies master equation [43].

This paper is structured as follows. We summarize our main results in Sec. 2, where we present the main master equations involved in this work: Redfield and Davies-Lindblad (well known), and coarse-grained, both time-independent and time-dependent (new). We also present bounds that provide the range of applicability of each type of master equation (new). We show that these bounds can be expressed entirely in terms of only two timescales, namely the fastest system decoherence timescale, and the characteristic timescale of the decay of the bath correlation function. The reader interested only in the results and not in the details of derivations can safely read only this section and then skip to the conclusions in Sec. 7. Derivations begin in Sec. 3, where we present a simple, new derivation of the time-independent CGME. We express the equation in CP (Lindblad) form, state the range of validity for different approximations made in the derivation, and thus for the equation itself. We also compare the range of applicability with other master equations, and note the spatial locality of the Lindblad generators of the CGME. At this point we are ready to address the case of time-dependent Hamiltonians. The derivation of the equation for this case (which happens to result in exactly the same form) is given in Sec. 4. In particular, this equation is well suited for the simulation of open system adiabatic quantum computation, but also for dynamical decoupling, which involves the opposite limit of very fast system dynamics. We note that while we will sometimes refer to qubits, none of the results we discuss in this work are limited to qubits, and in fact any finite multi-level system, or interacting set of such systems, is captured by the formalism. We study some applications and examples in Sec. 5, including error suppression using a dynamical decoupling protocol, and numerical results of the comparison of our master equation with the Redfield and Lindblad-Davies master equations. In the remainder of the paper we derive the range of validity of the various master equations we study. We give a detailed treatment of the Born approximation and the other approximations involved, in terms of rigorous bounds presented in Sec. 6, where we also discuss the generalization to multiple coupling terms and analyze the scaling of the error bounds for large system sizes. Conclusions are presented in Sec. 7. Various additional technical details are given in the Appendix.

2 Three master equations

This section provides a summary of our main results. We first present a brief background to define basic concepts and notation, followed by the definition and properties of the two main timescales we use later to provide bounds and ranges of applicability. We then summarize the three master equations we focus on in this work, followed by upper bounds on the distance of their solutions from the true state.

2.1 Background

Consider a system interacting with its environment as described by the total Hamiltonian

[TABLE]

Here $H$ is the time-independent system Hamiltonian (we treat the time-dependent case in Sec. 4), $H_{\text{b}}$ is the bath Hamiltonian, and $A$ and $B$ are, respectively, Hermitian system and bath operators. The more general situation, with $V=A\otimes B$ replaced by $\sum_{i}A_{i}\otimes B_{i}$, is easily treated as well (see Sec. 3.7), but we avoid it here to simplify the notation. We choose $A$ to be dimensionless with $\|A\|=1$ [the operator norm is defined in Eq. (36)] and $B$ to have energy units, but we deliberately do not set $\|B\|=1$, since we wish to account for baths (such as harmonic oscillator baths) for which $\|B\|$ can be infinite. We also set $\hbar\equiv 1$ throughout, so that energy and frequency have identical units.

Let the eigenstates of $H$ be $\{|n\rangle\}$ , and the corresponding eigenvalues $\{E_{n}\}$. Equivalently, $H=\sum_{n}E_{n}\Pi_{n}$, where $\Pi_{n}$ is a projector onto the eigenspace with eigenvalue $E_{n}$, and $\Pi_{m}\Pi_{n}=\delta_{mn}\Pi_{n}$. Note that in the interaction picture $A(t)=U^{\dagger}(t)AU(t)=\sum_{nm}A_{nm}e^{-iE_{mn}t}|n\rangle\langle m|$ where $E_{mn}=E_{m}-E_{n}$, $A_{nm}=\bra{n}A\ket{m}$, and $U$ is the solution of $\dot{U}=-iHU$. Let

[TABLE]

so that333The choice of the sign for this notation, as well as other notation choices, follow the textbook [2], p.133-134.

[TABLE]

Here $\omega$ runs over all the energy differences (Bohr frequencies) between the eigenvalues of $H$.

We assume henceforth that the bath is stationary — $[\rho_{\text{b}},H_{\text{b}}]=0$ where $\rho_{\text{b}}$ is the bath state — and introduce the bath correlation function $\mathcal{C}(t)$ and its Fourier transform $\gamma(\omega)$, known as the bath spectral density:

[TABLE]

where $B(t)=e^{iH_{\text{b}}t}Be^{-iH_{\text{b}}t}$ and

[TABLE]

Note that the real-valued functions $S(\omega)$ and $\gamma(\omega)$ are related via a Kramers-Kronig transform

[TABLE]

where $\mathcal{P}\left(\frac{1}{x}\right)[f]=\lim_{\epsilon\to 0}\int_{-\epsilon}^{\epsilon}\frac{f(x)}{x}dx$ is the Cauchy principal value.444This follows from the identity $\int_{0}^{\infty}d\tau e^{ix\tau}=\pi\delta(x)+i\mathcal{P}\left(\frac{1}{x}\right)$. When $A\otimes B$ is replaced by $\sum_{i}A_{i}\otimes B_{i}$, we note the positivity of $\gamma_{ij}(\omega)=\int dte^{i\omega t}\mathrm{Tr}[\rho_{\text{b}}B_{i}(t)B_{j}]$ as a matrix in the indices $i,j$ (of course this reduces to $\gamma(\omega)\geq 0$ in the scalar case). This is a non-trivial fact that is usually associated with Bochner’s theorem, as e.g., in the textbook [2]. In Appendix A we give two different proofs, one using Bochner’s theorem and another that is direct. If we assume not only that the bath state is stationary, but that it is also in thermal equilibrium at inverse temperature $\beta$, i.e., $\rho_{\text{b}}=e^{-\beta H_{\text{b}}}/\mathcal{Z}$, $\mathcal{Z}=\mathrm{Tr}e^{-\beta H_{\text{b}}}$, then it follows that the correlation function satisfies the time-domain Kubo-Martin-Schwinger (KMS) condition [2, 44]:

[TABLE]

If in addition the correlation function is analytic in the strip between $t=-i\beta$ and $t=0$, then it follows that the Fourier transform of the bath correlation function satisfies the frequency domain KMS condition [2, 44]:

[TABLE]

We note that the thermal equilibrium assumption is not necessary for the results we derive in this work, and we never use it in our proofs. We mention it here for completeness, and also since we use it in one of our dynamical decoupling examples later on. Finally, note that $\mathcal{C}(t)$ has units of frequency squared.

2.2 The bath correlation time $\tau_{B}$ and the fastest system decoherence timescale $\tau_{SB}$

We define the two key quantities

[TABLE]

Here $T$ is the total evolution time, used as a cutoff which can often be taken as $\infty$. We show later that in the interaction picture $\|\dot{\rho}\|_{1}\leq 4\|\rho\|_{1}/\tau_{SB}$, where $\rho$ is the system density matrix and the trace norm defined in Eq. (35), so that we can interpret $\tau_{SB}$ as the fastest system decoherence timescale, or timescale over which $\rho$ changes due to the coupling to the bath, in the interaction picture (i.e., every other system decay timescale, including the standard $T_{1}$ and $T_{2}$ relaxation and dephasing times, must be longer). The quantity $\tau_{B}$ is the characteristic timescale of the decay of $\mathcal{C}(t)$. We return in Sec. 6 to why we define these timescales in this manner, but note that Eq. (11b) becomes an identity if we choose $|\mathcal{C}(t)|\propto e^{-t/\tau_{B}}$ and take the limit $T\to\infty$.

We further note that $\tau_{SB}$ and $\tau_{B}$ are the only two parameters relevant for determining the range of applicability of the various master equations; in particular, nothing about the norm or time-dependence of $H_{\text{tot}}$ affects the accuracy of our time-dependent CGME, given below. In particular, $H_{\text{tot}}$ can be arbitrarily large or small in the operator norm (strong coupling, unbounded bath), or have an arbitrarily large derivative (non-adiabatic regime). This remark is important in light of previous approaches, so we expand on it some more.

In previous work it was common to extract a dimensionful coupling parameter $g$, i.e., to replace $B$ by $g\tilde{B}$, where $\|\tilde{B}\|=1$. One then defines $\tau_{B}=\int_{0}^{\infty}|\tilde{C}(t)|dt$ in terms of a dimensionless correlation function $\tilde{C}(t)=C(t)/g^{2}$ [31]. One problem with this approach is the arbitrariness of distributing the numerical factors between $g$ and $\tilde{B}$. Another is that it precludes unbounded baths, such as oscillator baths, for which $\|B\|=\infty$. In contrast, the formalism we present here is applicable even when $\|B\|$ diverges. Furthermore, $\|B\|$ contains an extra scale that does not carry any new information about the range of applicability of the master equations we discuss and derive. By not introducing this extra scale we highlight that there are only two free parameters in our analysis of the error: $\tau_{B}$ and $\tau_{SB}$, and no other information about the bath is needed. I.e., even though different baths will lead to different equations, our results are universal for any bath with the same values of these two parameters (in fact only their dimensionless ratio matters).

Note that we can relate $\tau_{B}$ and $\tau_{SB}$ to the spectral density. First, using $\gamma(\omega)>0$ and $C(t)=C^{*}(-t)$:

[TABLE]

The dephasing time $T_{2}$ and relaxation time $T_{1}$ for a single qubit are usually defined as $\propto 1/\gamma(0)$ and $\propto 1/\gamma(\omega)$ respectively, where $\omega$ is the qubit operating frequency (e.g., see Ref. [45]), so we see that $\tau_{SB}\lesssim T_{1},T_{2}$, which illustrates our earlier comment that $\tau_{SB}$ is the fastest system decoherence time-scale. Second, $\gamma^{\prime}(\omega)=i\int_{-\infty}^{\infty}tC(t)e^{i\omega t}dt$, so that

[TABLE]

where the limit is taken in Eq. (11b). And lastly, $\frac{d\gamma(-\omega)}{d\omega}=-\frac{d\gamma(\omega)}{d\omega}|_{-\omega}$, while differentiating the KMS condition yields $\frac{d\gamma(-\omega)}{d\omega}=-\beta e^{-\beta\omega}\gamma(\omega)+e^{-\beta\omega}\frac{d\gamma(\omega)}{d\omega}$, so that:

[TABLE]

This implies that $\gamma^{\prime}(0)>0$, so that under the KMS condition (10) we can replace Eq. (13) by

[TABLE]

The upper bounds on the approximation errors of all the master equations we present below involves the ratio $\tau_{B}/\tau_{SB}$ [see, e.g., Eqs. (26)-(28)]. The smallness of this ratio is sufficient for the CGME to be useful, while additional assumptions are required for the Davies-Lindblad equation and some versions of the Redfield equation. It follows from Eq. (15) that for our bound to be small it is necessary that the temperature $\beta^{-1}$ is not too low and/or the spectral density at $\omega=0$ (roughly the same as the dephasing rate $T_{2}^{-1}$) is not too large. The simple condition $\beta\gamma(0)\ll 1$ already rules out diverging spectral densities such as the case of $1/f$ noise, though this can be ameliorated by introducing a low-frequency cutoff.

We now present the three master equations, in the Schrödinger picture.

2.3 Redfield master equation

For the derivation see Sec. 3.2. This equation was known earlier than the Davies-Lindblad equation [19], and is often used in quantum chemistry. It is not CP and hence cannot be put in Lindblad form:

[TABLE]

where we defined the “filtered" operator:

[TABLE]

where $f(\omega)$ is defined in Eq. (7), and we used the subscript $R$ to denote that the density matrix satisfies the Redfield equation (we use a similar subscript notation below to distinguish the solutions of different master equations). There is no natural way to separate the Lamb shift. One of the benefits of the Redfield equation is that there is only one generator $A_{f}$ per interaction with the environment $A\otimes B$. We will show that as long as it preserves positivity, the Redfield equation is more accurate than the Davies-Lindblad master equation.

2.4 Davies-Lindblad master equation

For the derivation see Sec. 3.3. This is the familiar result [4]

[TABLE]

where $H_{\mathrm{LS}}$ is the Lamb shift. We note that ad hoc forms of the Lindblad equation are often written down and used without justification, e.g., with just one Lindblad term per qubit (e.g., $\sigma_{-}$ or $I-\sigma_{z}$ [46, 47]). In reality the Davies generators are derived from first principles, i.e., from the description of the total Hamiltonian of the system and bath. The number of Davies generators is unfortunately exponential in the number of qubits $n$: $\sum_{\omega}$ is over all $4^{n}$ energy differences.

2.5 Coarse-grained master equation (CGME) for time-independent or time-dependent system Hamiltonians

This is our main new set of results, generalizing Ref. [40] to the time-dependent case.

2.5.1 The time-independent case

For the derivation see Sec. 3.4. First, considering time-independent system Hamiltonians, we obtain:

[TABLE]

where the Lindblad operators are

[TABLE]

with $\operatorname{sinc}(x)\equiv\sin(x)/x$, and where $T_{a}$ is the averaging, or coarse-graining time, discussed below. The continuous family of terms labeled by $\epsilon$ is an unusual form of the Lindblad equation, but is manifestly CP. The more standard form in terms of a discrete sum over transition frequencies is equivalent to this one, and is given in Eq. (71) below. The range of applicability of the CGME is only slightly smaller than that of the Redfield ME, as we discuss in Sec. 3.5. Much like for Davies generators, there are exponentially many frequency differences, but unlike the Davies case the discretization of the continuous integral $\int d\epsilon$ makes it possible to keep the number of generators constant in the system size for each $A\otimes B$ term. These results are presented in Sec. 3.6.

The same applies for the Lamb shift, which is

[TABLE]

where $\omega_{+}=\frac{\omega+\omega^{\prime}}{2}$. Note that $F_{\omega\omega^{\prime}}$ is well defined in the limit $\omega_{+}\to 0$. The results above are generalized to the $n$-qubit setting in Sec. 3.7.

2.5.2 The time-dependent case

For the derivation see Sec. 4. In the case of a time-dependent system Hamiltonian $H(t)$ we obtain the same form, and even the same range of applicability, except that all the operators are now time-dependent:

[TABLE]

where the time-dependent Lindblad operators are

[TABLE]

where $A(t^{\prime},t)=U^{\dagger}(t^{\prime},t)AU(t^{\prime},t)$ with $U(t^{\prime},t)=\mathcal{T}\exp[-i\int_{t}^{t^{\prime}}H(s)ds]$ (the forward time-ordered exponential, from $t^{\prime}$ on the left to $t$ on the right), and the time-dependent Lamb shift is

[TABLE]

Although the coarse-graining time $T_{a}$ is a free parameter, we show that to minimize the upper bound we derive on the distance between the solutions of the Redfield and coarse-grained master equations it is nearly optimal to choose

[TABLE]

2.6 Error bounds and range of applicability

The error bound of the Redfield master equation is

[TABLE]

where $\rho_{\text{true}}(t)$ denotes the true (approximation-free) state.

The error bound of the Davies-Lindblad master equation is

[TABLE]

where $\delta E=$min${}_{i\neq j}|E_{i}-E_{j}|$ is the level spacing, with $E_{i}$ the eigenenergies of the system Hamiltonian $H$.

The error bound of both the time-independent and time-dependent coarse-grained master equation is:555Strictly, our proof is only for the time-independent case, but there do not seem to be any obstacles for its generalization to the time-dependent case.

[TABLE]

In the above expressions, $O(X)$ is understood for $X\to 0$, specifically for CGME $X=\sqrt{\tau_{B}/\tau_{SB}}e^{6t/\tau_{SB}}$. There is a subtlety in the definition of big-$O$ notation that we would like to emphasize. By definition, “$f(X)=O(X)$ at $X\to 0$” means there exist $X_{0},M$ s.t. for any $X\leq X_{0}$ we have $f(X)\leq MX$. Applying this to Eq. (28), we have the following unpacking of the above statement: there exist universal constants $x_{0},M$ such that if $\sqrt{\tau_{B}/\tau_{SB}}\leq X_{0}$ and $t\leq(\tau_{SB}/6)\text{ln}(X_{0}/\sqrt{\tau_{B}/\tau_{SB}})$, then $\|\rho_{\text{true}}(t)-\rho_{C}(t)\|_{1}\leq M\sqrt{\tau_{B}/\tau_{SB}}e^{6t/\tau_{SB}}$ holds. Note that for sufficiently small $\tau_{B}/\tau_{SB}$, a dimensionless quantity $t/\tau_{SB}$ can have an arbitrarily large constant value in this bound. The error can be made arbitrarily small at any fixed value of $t/\tau_{SB}$ just by choosing a sufficiently small $\tau_{B}/\tau_{SB}$. This is what is commonly referred to as the weak coupling limit. We emphasize that the constants $X_{0},M$ involved in this statement are universal, i.e. independent of the parameters of the equation such as the Hamiltonian and the initial state. In that sense, our bound is uniform. It does depend on time $t$ as explicitly stated, so it is not uniform in time.666Note that in both the Davies-Lindblad and CGME cases the l.h.s. is bounded above by $2$ due to the positivity of the density matrices; this is not necessarily true in the Redfield case. Thus, comparing these bounds, we observe that the Redfield bound is the tightest, which is natural given that it involves the fewest approximations. However, recall that the Redfield master equation is not CP. When comparing the bounds on the solutions of the two CP master equations, we observe that unlike the Davies-Lindblad equation the CGME does not involve the energy gap, which means that the CGME in principle has a much larger range of applicability, in particular for systems whose gap shrinks with growing system size.777By range of applicability we mean the range of parameters over which the approximation is accurate. This range can be defined formally in terms of an upper bound on the right-hand side (r.h.s.), e.g., $\sqrt{{\tau_{B}}/{\tau_{SB}}}e^{6t/\tau_{SB}}<\epsilon$ for some fixed $\epsilon<1$.

While the equations and inequalities can be derived for any bath, we made some extra assumptions to prove the error bounds as presented above. Specifically, we assumed a Gaussian bath and the convergence of our series expansion for the Born error, as defined and discussed in detail in Sec. 6.2. For a non-Gaussian bath, extra time scales relating to higher correlation functions generally appear, and the error bound can be derived by a straightforward generalization of our approach. It will contain, besides $\tau_{B}$ and $\tau_{SB}$, additional terms involving those time scales.

3 (Re-)Derivation of the coarse-grained master equation

Our goal in this section is to present a compact and simplified derivation of the CGME found in Ref. [40]. This master equation is in Lindblad form but avoids making use of the RWA, which contributes the largest approximation error. Instead, the key idea is to introduce a coarse-graining (averaging) timescale, first exploited in this context in Refs. [20, 22]. In this section we consider the case of a time-independent system Hamiltonian, and later generalize our derivation to the time-dependent case, which includes adiabatic evolution as a special case.

Along the way we also derive the Redfield and Davies-Lindblad master equations.

3.1 Setting up

We start from the total Hamiltonian given in Eq. (3), and let $V=A\otimes B$. Recall that $A$ is dimensionless and $B$ has units of energy. The first few steps in our derivation are standard [2]. In the double system-bath interaction picture $V(t)=U_{\textrm{0}}^{\dagger}(t)VU_{\textrm{0}}(t)$ [where $U_{\textrm{0}}$ is the solution of $\dot{U}_{\textrm{0}}=-i(H\otimes I_{\text{b}}+I\otimes H_{\text{b}})U_{\textrm{0}},\leavevmode\nobreak\ U_{0}(0)=1$], and the state satisfies

[TABLE]

Substituting this back into the original equation yields:

[TABLE]

The reduced system state in the interaction picture is defined as $\rho_{\textrm{true},I}(t)=\mathrm{Tr}_{\text{b}}[\rho_{\text{tot}}(t)]$, where the subscript “true" denotes that this is the correct, completely approximation-free state. The corresponding true system state in the Schrödinger picture is $\rho_{\textrm{true}}=U_{0}(t)\rho_{\textrm{true},I}U^{\dagger}_{0}(t)$. The Born approximation,

[TABLE]

together with the shift of $B$ such that $\mathrm{Tr}[\rho_{\text{b}}B]=0$ allows one to write:

[TABLE]

which can be understood as the definition of the Born approximation error $\mathcal{E}_{B}$:

[TABLE]

Henceforth we assume that the bath correlation function decays rapidly, namely that $\tau_{B}\ll\tau_{SB}$, where $\tau_{SB}$ and $\tau_{B}$ were defined in Eq. (11).

Now, let $\rho_{B,I}$ denote the solution of the master equation after the Born approximation:

[TABLE]

Later we collect other errors as letters next to $B$. The error estimate given in Eq. (33) for $\|\mathcal{E}_{B}\|_{1}$ is derived in Sec. 6.2 [Eq. (198b)].

We digress briefly to carefully explain the meaning of the norms and big-$O$ notation used here, since this is important for the remainder of this work. $\mathcal{E}_{B}$ is an operator acting on the system Hilbert space. The trace norm $\|A\|_{1}$ is:

[TABLE]

The operator norm $\|X\|$ is defined as follows:

[TABLE]

where $\lambda_{X,i}$ are the eigenvalues of $|X|\equiv\sqrt{X^{\dagger}X}$ (singular values of $X$) indexed by $i$. The big-$O$ notation $\mathcal{E}=O(x)$ is taken for times such that $t/\tau_{SB}\leq M$ and $x\to 0$, i.e., there exist an $M$-dependent number $\epsilon_{M}$ and an $M$-dependent constant $C_{M}$ such that for any $x<\epsilon_{M}$ the error $\|\mathcal{E}\|_{1}\leq C_{M}x$.

It is natural to use the trace norm for density matrices. Indeed, if we want to find the deviation in an observable $O$ relative to the difference between two states, $\delta\rho=\rho_{1}-\rho_{2}$, then

[TABLE]

The second expression gives a tighter bound, if we manage to have the same bound on $\|\delta\rho\|_{1}$ as on $\|\delta\rho\|$. Fortunately, this will turn out to be the case in this work, and was already observed for the Markov error in Ref. [31]. The key property used to prove the inequality above is submultiplicativity, valid for any unitarily invariant norm [48]:

[TABLE]

For simplicity, we assume the bath state to be stationary, such that $\mathcal{C}(t,t^{\prime})=\mathrm{Tr}[\rho_{B}B(t)B(t^{\prime})]$ is invariant with respect to shifts in time of both arguments.888The derivation can potentially be extended to the general case of a two-time correlation function, i.e., without assuming translational invariance; indeed the derivation in [40] does not. We can then replace $\mathcal{C}(t,t^{\prime})$ by $\mathcal{C}(t)=\mathrm{Tr}[\rho_{B}B(t)B]$, as in Eq. (6a). Recall that $\mathcal{C}^{*}(t)=\mathcal{C}(-t)$. After expanding the commutators, one arrives at

[TABLE]

where $K^{B,2}_{t-\tau}$ is a superoperator. Taking the trace norm, we note that:

[TABLE]

where we used our earlier choice of setting $\|A\|=1$, and the constant $c_{B}$ is chosen as an upper bound on $\|\rho_{B,I}\|_{1}$ (if this were a CP map, $c_{B}=1$ would hold). Under the condition $\forall t,\leavevmode\nobreak\ \|\rho_{\text{test}}(t)\|_{1}=1$ we have

[TABLE]

We next introduce the Markov approximation $\rho_{B,I}(\tau)\mapsto\rho_{B,I}(t)$:

[TABLE]

Equation (42b) is the Redfield equation [2], which is notoriously non-CP (though various fixes have been proposed [49, 28]). The Markov approximation error $\mathcal{E}_{M}$ is of same order as the Born approximation error given in Eq. (33) (for more detail see Sec. 6).

We note that the errors on the r.h.s. are not generally additive: if we try to write a Markovian master equation for the true (approximation-free) evolution $\rho_{\text{true},I}$, then the error on the r.h.s. will, besides $\mathcal{E}_{B}+\mathcal{E}_{M}$, contain a correction to the Markov error from the Born error:

[TABLE]

In this particular case, using methods from Section 6 the difference between bounds on $\|\mathcal{E}_{M}(\rho_{\text{true}})\|_{1}$ and $\|\mathcal{E}_{M}\|_{1}$ can be found to be subleading in $\tau_{B}/\tau_{SB}$. Since we wish to present higher orders of the error, we do not attempt to collect errors as in Eq. (43a). Instead, our derivation will present a sequence of equations, e.g., $\rho_{\text{true}},\rho_{B},\rho_{BM}$, and an error estimate for each step.

After using Eq. (5), i.e., in the frequency representation, Eq. (42b) takes the following form:

[TABLE]

We sometimes omit the explicit time dependence of $\rho$ on the r.h.s. since it is always $\rho(t)$ from hereon.

3.2 Redfield master equation

Let us digress briefly in order to establish the form of the Redfield equation presented in Sec. 2. Rotating Eq. (42b) back to the Schrödinger picture we obtain:

[TABLE]

Introducing a new integration variable $t^{\prime}=t-\tau$:

[TABLE]

We extend the upper integration limit from $t$ to $\infty$, which introduces an additional error:

[TABLE]

where for the first summand in Eq. (47b) we assumed a lower cutoff on the evolution time $t>\tau_{SB}\delta>0$ where $\delta$ is a small constant number which we absorbed in the big-$O$ notation [i.e., $1/t=O(1/\tau_{SB})$]. We analyze the transients happening for $t<\tau_{SB}\delta$ in Sec. 6.6. The constant $c_{BM}$ is chosen as an upper bound on $\|\rho_{BM,I}\|_{1}$ (again, if this were a CP map, $c_{BM}=1$ would hold). We will see in Sec. 6 that $c_{BM}=O(1)$ in terms of our big-$O$ notation, specified in Eq. (49). For the second summand, we introduced a new bath parameter

[TABLE]

For most physically motivated choices of bath correlation functions $\int_{0}^{\infty}t|\mathcal{C}(t)|dt$ converges, so we may replace the total evolution time $T\to\infty$ and then $\epsilon_{T}=0$. But, in case $\tau_{B}/\tau_{SB}$ diverges as a function of $T$ [Eq. (11b)], the above bound should be used.

The big-$O$ notation in terms of all of these parameters $\tau_{SB}\delta\leq t\leq\tau_{SB}M,\epsilon_{T},\tau_{B}/\tau_{SB},c_{BM}$ should be understood as follows: there exist constants $C(M,\delta),\epsilon(M,\delta)$ such that for $\tau_{B}/\tau_{SB}+\epsilon_{T}\leq\epsilon(M,\delta)$

[TABLE]

As seen in the derivation above, $C(\delta,\epsilon)=\max(4c_{BM}/\delta,4c_{BM})$ in fact does not depend on $\epsilon$.

Note that the change of integration limit $t\to\infty$ is optional: the Redfield equation can be used without it, but the change does simplify its form. Eq. (47a) without the error gives:

[TABLE]

Introducing the “filtered" operator $A_{f}=\int_{0}^{\infty}\mathcal{C}(-t^{\prime})A(-t^{\prime})dt^{\prime}$ immediately yields Eq. (16), so that in fact $\rho_{BMl}=\rho_{R}$. Replacing $A(-t^{\prime})$ in $A_{f}$ by its Fourier decomposition [Eq. (5)] gives Eq. (17), for which the integral can be expressed in multiple ways:

[TABLE]

3.3 Davies-Lindblad master equation

As a second digression let us establish the form of the Davies-Lindblad equation presented in Sec. 2. We start from the form given in Eq. (44), and use $\omega t+\omega^{\prime}\tau=\omega(t-\tau)+(\omega+\omega^{\prime})\tau$ and a change of variables to $t^{\prime}=\tau-t$ to write it as

[TABLE]

Next we apply the RWA, in which one considers the term $e^{-i(\omega+\omega^{\prime})(t+t^{\prime})}$ to be rapidly oscillating and hence self-cancelling, so the non-negligible contribution is only due to the terms with $\omega^{\prime}=-\omega$. This allows us to write:

[TABLE]

where we replaced $\rho_{BM,I}$ by $\rho_{D,I}$ and $\omega^{\prime}$ by $-\omega$. The error $\mathcal{E}\sim O(1/\tau_{SB})$ between the right-hand sides of Eqs. (52) and (53) is large but rapidly oscillating, so the solutions $\rho_{BM,I}(t)$ and $\rho_{D,I}(t)$ remain close (see Sec. 3.5).

We now transform back to the Schrödinger picture (recall that this requires multiplying each $A_{\omega}$ by $e^{i\omega t}$, but these cancel now due to the corresponding $A_{-\omega}$):

[TABLE]

It is convenient to separate the real and imaginary parts, using $\mathcal{C}^{*}(t)=\mathcal{C}(-t)$:

[TABLE]

Using the decomposition given by Eq. (55) in Eq. (54) directly yields the form of the Davies-Lindblad master equation presented in Eq. (18). We note that versions of the Davies-Lindblad master equation that allow a time-dependent drive have been derived before (see, e.g., Refs. [43, 31]), but such derivations must always assume that the driving is adiabatic.

3.4 Coarse-grained master equation

Let us now return to our main goal in this section, the derivation of the CGME. At this point we deviate from the standard derivations and instead follow the coarse-graining approach [20, 22, 40]. There are two main steps: (i) Time-averaging of the state, which is done in lieu of the RWA (and reduces to the RWA in the limit of large time-averaging scale; see Appendix C), (ii) removal of a small part of the integration domain, which restores complete positivity (new in this work). We describe each in turn, along with the associated error.

3.4.1 Time-averaging of the state

Before we leave the interaction picture, we coarse-grain, or time-average over $T_{a}$ such that $\tau_{B}\ll T_{a}\ll\tau_{SB}$. Specifically, we consider another equation satisfied by a new state $\rho_{BMT,I}(t)$, where the time-averaging operation $\frac{1}{T_{a}}\int_{t-T_{a}/2}^{t+T_{a}/2}dt^{\prime}$ is applied to the r.h.s. of Eq. (44) taken at time $t^{\prime}$, except that the state is still at time $t$ on the l.h.s.:

[TABLE]

Alas, the associated error $\mathcal{E}_{T}$ has $\|\mathcal{E}_{T}\|_{1}=O(1/\tau_{SB})$, which is not small. The reason is that the fast-rotating terms may have a large derivative. What can be proven is that the solutions of Eqs. (44) and (56) remain close:

[TABLE]

where $t\leq c\tau_{SB}$, and the r.h.s. $O(T_{a}/\tau_{SB})$ does not depend on $t$, only on $c$, which is considered a constant for the purposes of big-$O$ notation. The proof is given in Sec. 6.4, in Lemma 2 (this analysis first appeared in Ref. [32]). The error we track is now the error of the solution, not its derivative.

We use the following fact about time-local (Markovian) and time-nonlocal (non-Markovian) differential equations:

Lemma 1**.**

Assume that

[TABLE]

where $\mathcal{L}$ is a linear superoperator and $\Lambda$ is a positive constant such that $\sup_{\tau,x:\|x\|_{1}=1}\|\mathcal{L}_{\tau}(x)\|_{1}\leq\Lambda$. Or, assume that

[TABLE]

where $K_{t-\tau}(x)$ is a linear superoperator and $\Lambda$ is a positive constant such that $\sup_{t,x(\tau):\|x(\tau)\|_{1}=1}\int_{0}^{t}\|K_{t-\tau}(x(\tau))\|_{1}d\tau\leq\Lambda$. Then:

[TABLE]

For the proof see Sec. 6.1. We can take $\Lambda=4/\tau_{SB}$ because of Eq. (41, 42c). Using this, we rewrite the result of the Lemma as follows, for brevity: if $t\leq c/\Lambda=c\tau_{SB}/4$, then $\|x-y\|_{1}=O(\tau_{SB}\|\mathcal{E}\|_{1})$. Thus, using the previously noted results that $\|\mathcal{E}_{B}\|_{1},\|\mathcal{E}_{M}\|_{1}=O\left({\tau_{B}}/{\tau_{SB}^{2}}\right)$:

[TABLE]

3.4.2 Neglecting part of the integration domain to regain complete positivity

Upon changing variables to $s=t^{\prime}-t$, and renaming $s$ as $t^{\prime}$, Eq. (56) becomes:

[TABLE]

We rotate out of the system interaction picture into the system Schrödinger picture by defining the coarse-grained state $\rho_{BMT}=U(t)\rho_{BMT,I}U^{\dagger}(t)=e^{-iHt}\rho_{BMT,I}e^{iHt}$, leaving only the bath interaction picture, by multiplying each summand by a $e^{i(\omega+\omega^{\prime})t}$ factor:999 Note that the transformation back from the interaction picture is $U\dot{\rho}_{I}U^{\dagger}$. Thus, e.g., $A_{\omega}\rho_{I}A_{\omega^{\prime}}\mapsto UA_{\omega}\rho_{I}A_{\omega^{\prime}}U^{\dagger}=UA_{\omega}U^{\dagger}\rho UA_{\omega^{\prime}}U^{\dagger}$, where $U=\sum_{k}e^{-iE_{k}t}\Pi_{k}$. Using Eqs. (4) and (5), we have $UA_{\omega}U^{\dagger}=\sum_{kl}\sum_{mn:E_{mn}=\omega}e^{-i(E_{k}-E_{l})t}\Pi_{k}\Pi_{n}A\Pi_{m}\Pi_{l}=\sum_{mn:E_{mn}=\omega}e^{-i(E_{n}-E_{m})t}\Pi_{n}A\Pi_{m}=e^{i\omega t}A_{\omega}$, with the last equality holding due to the fact that $E_{n}-E_{m}$ is fixed at $\omega$ for all combinations of $m$ and $n$. The claim now follows.

[TABLE]

Next, we change variables to $\tau^{\prime}=\tau-t$, so that

[TABLE]

We define $\mathcal{L}^{BMT}_{t}(\rho_{BMT})$ to be the superoperator on the r.h.s. The key trick now is to replace $-t$ by $-T_{a}/2$, which lets us write (note that Lemma 1 will work for $\mathcal{E}_{p}$ either in the initial or in the resulting equation, like we do here):

[TABLE]

This replacement cuts a corner of area $\sim\tau_{B}^{2}$ out of the integration domain with non-vanishing $\mathcal{C}(t)$, itself of area $\sim T_{a}\tau_{B}$, as illustrated in Fig. 1(a) for $t>T_{a}/2$ (times $t<T_{a}/2$ introduce a transient, as discussed in more detail in Sec. 6.4). Thus, the corner being cut for $t>T_{a}/2$ is just a $\tau_{B}/T_{a}$ fraction of the whole integral in Eq. (65b) (ignoring factors of $\sqrt{2}$), which [using Eq. (11a)] is itself upper bounded by $1/\tau_{SB}$ in absolute value. Thus, the error introduced by the corner removal is

[TABLE]

[for $t<T_{a}/2$ it can be $O(1/\tau_{SB})$ for a short transient]. The constant in big-$O$ notation here is independent of $\tau_{SB},\tau_{B},T_{a}$. As we show below, $\mathcal{E}_{p}$ is the cost of recovering complete positivity, in the sense that the master equation after subtraction of $\mathcal{E}_{p}$ is in canonical Lindblad form.

Eq. (65a) is the coarse-grained master equation from Ref. [40]. To see this we need some properties of the coefficient $x_{\omega\omega^{\prime}}$. In particular, we need the properties

[TABLE]

with the first following from $\mathcal{C}^{*}(t)=\mathcal{C}(-t)$, and the second from interchanging integration limits after a reflection w.r.t. the $\tau^{\prime}=-t^{\prime}$ axis, as illustrated in Fig. 1(b). This allows us to reshuffle the h.c. part of Eq. (65a), taking the $-\omega^{\prime},-\omega$ term in the sum for the $\omega,\omega^{\prime}$ term in the original part. The resulting equation is:

[TABLE]

The $\textrm{Im}x_{\omega\omega^{\prime}}$ part cancels for $A_{\omega}\rho A_{\omega^{\prime}}$ and sends $A_{\omega^{\prime}}A_{\omega}$ into the Lamb shift:

[TABLE]

where to get the second form we used $\textrm{Im}(x_{\omega\omega^{\prime}})=\frac{1}{2i}(x_{\omega\omega^{\prime}}-x^{*}_{\omega\omega^{\prime}})$, combined Eqs. (65b), (67a) into a single integral using the $\operatorname{sgn}(t_{1}-t_{2})$ function, and used Eq. (5). The Lamb shift can be further simplified as shown in Appendix B. The simplified form is the one we presented as a part of our main result [Eq. (21)] in Sec. 2.

We define $\gamma_{\omega\omega^{\prime}}\equiv x_{\omega\omega^{\prime}}^{*}+x_{\omega\omega^{\prime}}=2\textrm{Re}x_{\omega\omega^{\prime}}$ and obtain:

[TABLE]

Therefore

[TABLE]

which is Eq. (50) from [40] up to a shift in the center of averaging and changing variables $A_{\omega^{\prime}}\to A_{\omega^{\prime}}^{\dagger},\leavevmode\nobreak\ \omega^{\prime}\to-\omega^{\prime}$.

3.4.3 Complete positivity

Complete positivity is not immediately apparent from Eq. (71). Let us demonstrate it next. Substituting $\mathcal{C}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-i\epsilon t}\gamma(\epsilon)d\epsilon$ [Eq. (6a)], where $\gamma(\epsilon)\geq 0$ (see Appendix A), into Eq. (70), we have:

[TABLE]

where

[TABLE]

is a filter function whose integrated form is given in Eq. (20). By using the above decomposition for $\gamma_{\omega\omega^{\prime}}$ in Eq. (71), we obtain the manifestly CP Lindblad form given in Eq. (19).

Note further that it was shown in Ref. [40] that the Davies-Lindblad equation [Eq. (18)] arises as the $T_{a}\to\infty$ limit of the CGME. Similar ideas have appeared in [37, 38]. We provide this proof for completeness in Appendix C. This derivation has the advantage that it arrives at the Davies-Lindblad equation without invoking the (uncontrolled) RWA, and instead shows that the Davies-Lindblad equation is the limit of infinite coarse graining time of the CGME. One of our new results is to show how the error of the Davies-Lindblad equation can be controlled for a sufficiently large (but not infinite) coarse graining time; see Sec. 6.5.

3.5 Ranges of applicability of the three master equations

We now discuss the ranges of applicability of the master equations derived above in more rigor than in Sec. 2.6, while postponing the derivations to Sec. 6. Recall that the ranges are dependent on the timescales $\tau_{SB}$ and $\tau_{B}$ defined in Eq. (11).

Errors are incurred from different approximations made along the way. The very first approximation made in the derivation of all three master equations is the Born approximation (see Sec. 3.1). Unlike the other approximations, we do not prove a rigorous error bound for the Born approximation error; we just provide a bound on the first order contribution to this error. However, under the assumption that the error converges in the first place, the Born error will turn out to be subleading. Therefore we settle for an estimate of this error for the sake of simplicity, but we note that a rigorous bound can be obtained [4].

We first collect all the errors introduced in the derivation of the CGME: the difference in solutions $\|\rho_{\text{true},I}(t)-\rho_{BMT,I}(t)\|_{1}$ [Eq. (61b), which includes the error due to time-averaging], and the error due to enforcing complete positivity, $\|\mathcal{E}_{p}\|_{1}=O(\tau_{B}/(T_{a}\tau_{SB}))$ [Eq. (66)]. The difference between the solution of the CGME $\rho_{C}(t)$ and the true time evolved state $\rho_{\text{true}}(t)$ thus satisfies (note that the norm is unitarily invariant and hence unaffected by transformation to or out of the interaction picture):

[TABLE]

Here we used Lemma 1 to multiply $\|\mathcal{E}_{p}\|_{1}$ by $\tau_{SB}$. Optimizing $T_{a}$ to minimize this error, we obtain $T_{a}=\sqrt{\tau_{SB}\tau_{B}}$, and:

[TABLE]

Note that if we define the coupling $g$ strength via $1/\tau_{SB}=g^{2}\tau_{B}$ (recall the discussion in Sec. 2.2), then the error bound of Eq. (75) is $O(g\tau_{B})$.

We can in fact improve on the result presented in Eq. (75). A careful derivation presented in Sec. 6 generalizes the above result from the $t>T_{a}/2$ case discussed here [recall Eq. (66)] to $t>0$. The big-$O$ notation in Eq. (75) means that there exist numbers $(C,\delta,c)$ such that for all $0<t<c\tau_{SB}$ and $\tau_{B}/\tau_{SB}<\delta$ the error is $\|\rho_{\text{true}}(t)-\rho_{C}(t)\|_{1}\leq C(c,\delta)\sqrt{\tau_{B}/\tau_{SB}}$. We do not know the function $C(c,\delta)$ explicitly, because we do not keep track of it in our analysis of the Born error [the details of which are given in Sec. 6.2]. The contributions from all the other errors are known explicitly. Our derivation also allows us to state a stronger result exhibiting the coefficients and the $t$-dependence of the first few leading terms of the error as a series in $\tau_{B}/\tau_{SB}$, since the unknown Born contributions start with $O(\tau_{B}^{2}/\tau_{SB}^{2})$. The improved bound is:

[TABLE]

The time $t$ can now be arbitrarily large, as long as the combination $e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}$ is small. Specifically, the big-$O$ notation in Eq. (76) means that there exist constants $C,\epsilon$ such that for all $e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\leq\epsilon$, $O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\right)\leq Ce^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}$.

The time-averaging window leading to Eq. (75) was chosen sub-optimally, as the more careful analysis yielding Eq. (76) uses $T_{a}=\sqrt{\tau_{B}\tau_{SB}/5}$. Including higher order corrections in $\tau_{B}/\tau_{SB}$ leads to a minor improvement in the subleading terms in the bound, but does not affect the leading term (see Section 6.4.4).

Let us now discuss the time-dependence of the error. From the above expression it follows that for any $t$, $\|\rho_{\text{true}}(t)-\rho_{C}(t)\|_{1}=O(\sqrt{\tau_{B}/\tau_{SB}}e^{6t/\tau_{SB}})$, which is Eq. (28). For the sake of completeness, we present the errors of the Redfield equation (for $\epsilon_{T}=0$ [recall Eq. (48)]; for the general case, see Sec. 6.6): $\|\rho_{\text{true}}(t)-\rho_{R}(t)\|_{1}=O\left(\frac{\tau_{B}}{\tau_{SB}}\text{ln}\frac{\tau_{SB}}{\tau_{B}}e^{12t/\tau_{SB}}\right)$, which is Eq. (26), and the Davies-Lindblad equation: $\|\rho_{\text{true}}(t)-\rho_{D}(t)\|_{1}=O\left(\left(\frac{\tau_{B}}{\tau_{SB}}+\sqrt{\frac{1}{\tau_{SB}\delta E}}\right)e^{12t/\tau_{SB}}\right)$, which is Eq. (27).

We note that for Eq. (46) there is no transient before the introduction of $\mathcal{E}_{l}$:

[TABLE]

thus the error grows from zero linearly at small times, and the slope of the upper bound is $O(\tau_{B}/\tau_{SB}^{2})$. In contrast, for the Redfield, CGME and Davies-Lindblad equations, there is a nonzero transient error bound at $t\to 0$. The error itself starts at zero (our bounds are not tight), but the slope is $O(1/\tau_{SB})$ instead.

It is important to emphasize that the solutions of open system master equations do not actually grow exponentially. In fact, for time-independent system Hamiltonians all of the equations presented have steady states to which they converge in the operator and trace norm at some rate. If that rate is sufficiently fast, one can prove a much stronger result [41] than our error bounds: a stability of the open system dynamics to small perturbations. The stability can be expressed in terms of a time-independent error bound on the difference in solutions.

We remark that the coefficients we presented and the time-dependence of the error are explicit, while other work hides them in big-$O$ notation. Our results allow one to put error bars on the numerical solution of the above equations.

3.6 Discretization, spatial locality, and the Lieb-Robinson bound

Different forms of the CGME are preferred depending on the application. The form given in Eq. (71) with a discrete sum $\sum_{\omega,\omega^{\prime}}$ is ready for numerical implementation but at exponential cost, since for a Hilbert space of dimension $2^{n}$ the number of terms in the sum is $2^{4n}$. Numerical solution of the equation already requires at least $2^{2n}$ numbers just to store the density matrix. It is therefore desirable to reduce the number of terms on the r.h.s. to a constant in $n$ (or polynomial in $n$ for the general case of multiple interaction terms; see Sec. 5). We first present a discretization that achieves this goal. It is also desirable to have each term locally supported, instead of acting on the whole system. This speeds up the numerics, and also allows one to meaningfully use the equation with methods that use resources that are polynomial in $n$ to store the density matrix, e.g., matrix product states (MPS) [50, 51]. For a recent implementation of these methods see [46, 47]. We present a way to write down a local equation at the end of this subsection.

In Eq. (19), the number of terms is infinite because of the integration over $\epsilon$. In Appendix D we explain how a discrete approximation can be derived. The result is:

[TABLE]

where $A_{\epsilon}$ is given by Eq. (20) and the values of $\Delta\epsilon$, $k^{*}$ and $\|\mathcal{E}_{k}\|_{1}$ are given in Appendix D, and are system-size-independent as long as $[H,A]=O(1)$. We do not present a similar discrete approximation for the Lamb shift, leaving it to future work.

We now proceed to show that storing the Lindblad operators $A_{\epsilon}$ does not require resources exponential in the system size $n$. Recall that $A_{\epsilon}=\sqrt{\frac{\gamma(\epsilon)}{2\pi T_{a}}}\int_{-T_{a}/2}^{T_{a}/2}e^{i\epsilon t}A(t)dt$. This form is interesting because $A(t)$ is evaluated over the interval $[-T_{a}/2,T_{a}/2]$. An operator $A(t)$ at those times is local (with exponentially decaying tails) with a locality radius $vT_{a}/2$, where $v$ is the Lieb-Robinson velocity ($v\sim\|H_{\textrm{local}}\|$, the coupling strength in the Hamiltonian) [52]. More rigorously, we consider a system defined on a graph. Each vertex of the graph supports a finite-dimensional Hilbert space, and the system Hilbert space is a tensor product of vertex Hilbert spaces. The edges of the graph are present if there are nontrivial interactions between the corresponding vertices in the Hamiltonian. That is, for vertices $i,j$ we can define a trace over the degrees of freedom at other vertices, denoted $\mathrm{Tr}_{\overline{i,j}}H-\mathrm{Tr}_{\overline{i}}H-\mathrm{Tr}_{\overline{j}}H$. If this operator is nonzero (for traceless $H$), then there is a nontrivial interaction and the edge $(i,j)$ is present in the graph. The Lieb-Robinson bound can be used (see, e.g., Lemma 5 in Ref. [53]) to derive the following decomposition of the operator $A_{\epsilon}$ into its local and small nonlocal part:

[TABLE]

where $c$ is a universal (system-independent and $\epsilon$-independent) constant, and

[TABLE]

is an operator strictly local within the graph distance $x$ from the original location of $A$. We define the graph distance as the length (number of edges) of the smallest connected path of edges between two vertices. Here $H^{(x)}=\mathrm{Tr}_{\overline{x}}H$ is composed of those Hamiltonian terms which involve vertices within graph distance $x$ from the original location of $A$. If we drop $\delta A$ from Eq. (79), we can control the resulting error similarly to what is described in Sec. 6 below:

[TABLE]

This master equation has local generators $A^{(x)}_{\epsilon}$ with $x=vT_{a}/2+\delta r$. We also assumed that some truncation method is applied to the Lamb shift $H_{\mathrm{LS}}\to H^{(x)}_{\mathrm{LS}}$, although we do not present it in this work. One can now use standard methods [54] to recast Eq. (81)) as a stochastic Schrödinger equation, with those operators corresponding to jumps:

[TABLE]

here $r$ represents a random number used to generate trajectories. We will not present the explicit form of $G$, but note that:

[TABLE]

where the sum is over qubits and each $G_{i}$ is a local operator of radius $x=vT_{a}/2+\delta r$ and can be computed in $O(1)$ time in the system size $n$. The time evolution is then given by a trotterization $U^{\text{Trotter}}(r)$ of

[TABLE]

$U^{\text{Trotter}}(r)$ is a geometrically local circuit of depth $O(t)$ where each gate takes $O(1)$ time to compute. Contraction of such a circuit $|\psi(t,r)\rangle=U^{\text{Trotter}}(r)|\psi(0)\rangle$ does not necessarily have a local structure, even if the initial condition $|\psi(0)\rangle$ is a product state. However, note that in 1d $|\psi(t,r)\rangle$ has the structure of a MPS [50]:

[TABLE]

Where $M_{i}^{s_{i}}$ are $\chi\times\chi$ matrices for each $i,s_{i}$, and $\chi$ is called a bond dimension. A portion of the time evolution of a MPS is illustrated in Fig. 2 for 2-local $G_{i}$, for which $U^{\text{Trotter}}=\prod_{t}(\prod_{i=2k}U_{i,i+1}\prod_{i=2k-1}U_{i,i+1})$ and $U_{i,i+1}=e^{-iG_{i}dt}$. Iterating the time step naively grows the bond dimension as $\chi=e^{O(t)}$.

The computation of any observable in $|\psi(t,r)\rangle$ requires a contraction of such a tensor network, which takes $O(n)e^{O(t)}$ computational time, and $e^{O(t)}$ memory. In practice, MPS calculations employ a truncation scheme for the bond dimension after every time step, which fixes the bond dimension to a constant and leads to $O(n,t)$ computational time, $O(n)$ memory. Such a method is generally not a controlled approximation, but it may be possible to develop a special version with an error estimate. Averaging over randomness adds an extra computational cost to either method. We leave the development of further details to a future publication.

The simultaneous locality and positivity of Eq. (81) is a new result: previously, the Davies generators were only local for commuting Hamiltonians, and the Redfield master equation leads to completely positive evolution only for flat spectral densities. Thus, this is the first method for ab initio completely positive open system evolution of MPSs.

3.7 Scaling with size

Here we discuss the system-size dependence. So far we assumed that there is only one term $V=A\otimes B$ describing the interaction with the bath. In reality, there are $n$ terms like this, one for each qubit. If qubits are coupled to the same bath $V=\sum_{i}^{n}A_{i}\otimes B$, it will lead to correlated noise (collective decoherence [55]). If the qubits are coupled to different baths $V=\sum_{i}^{n}A_{i}\otimes B_{i}$, it can lead to uncorrelated noise if $\mathrm{Tr}[\rho_{B}B_{i}(t)B_{j}]=0$. We will consider the latter case since it contains the former case if one sets $B_{i}=B$ (though distinctly different phenomena such as decoherence-free subspaces appear in this special case [56, 57]). The form of the Davies-Lindblad and Redfield equations can be found elsewhere; here we present only the form of the CGME. The sum is carried through to the frequency form of the equation — Eq. (71):

[TABLE]

The theoretically optimal $T_{a}=\sqrt{\tau_{SB}\tau_{B}/5}$ is given by the same expression in terms of time scales $\tau_{B}$ and $\tau_{SB}$, which are now defined using $\int|\mathcal{C}|\mapsto\text{max}_{ij}\int|\mathcal{C}_{ij}|$. Such a definition keeps $\tau_{B}$ and $\tau_{SB}$ consistent with the single qubit quantities. The locality of the equation is now also dependent on the decay of $\mathcal{C}_{ij}$ with the distance between $i$ and $j$. If it can be approximated by zero outside some radius (“strong locality”), the equation remains locally generated in the sense that every term on the r.h.s. contains local terms around some qubit $i$ conjugating the density matrix, and there are polynomially many of them. If we use a weaker sense of locality, when the operators on the right and left of the density matrix are allowed to be local around different $i,j$, and neglect the Lamb shift, then no assumption on $\mathcal{C}_{ij}$ is needed.

Let us present an explicitly Lindblad form. Let $\mu$ index the positive (see Appendix A) eigenvalues of $\gamma_{ij}(\omega)=\int_{-\infty}^{\infty}e^{i\omega s}\mathcal{C}_{ij}(s)ds$:

[TABLE]

This allows one to write:

[TABLE]

Next, note that if we let $\tilde{A}=\sum_{j}U_{\mu j}(\epsilon)A_{j}$, the previous derivation goes through.

[TABLE]

where the Lindblad operators are

[TABLE]

These results generalize Eqs. (19)-(20) to the $n$-qubit setting. Now the locality is more directly defined: if the operators $A_{\epsilon,\mu}$ can be approximately truncated to a ball around some point for each $\mu$, then the CGME for $n$ qubits is local in the strong sense described above; otherwise it is local in a weak sense. In any case, it has $n$ times more generator terms than the single-qubit CGME. The latter has a constant number of terms after the discretization described in Appendix D.

Let us now discuss the range of applicability of the $n$-qubit CGME. In the case of correlated noise, the norm of each $A$ effectively grows by a factor of $n$. Thus the bound on the relaxation time becomes $\Lambda\sim n^{2}/\tau_{SB}$ with $\tau_{SB}$ taken from the single qubit case. For uncorrelated noise, the bound is $\Lambda\sim n/\tau_{SB}$. In both cases, the range of applicability can be derived by directly repeating the derivation given in Sec. 6:

[TABLE]

where $\tau_{B}$ is system-independent, and thus is the same as the single qubit quantity. The apparently shrinking range with $n$ should not discourage us, since it can be mitigated under the assumption of locality. Namely, one can prove that the error in local density matrices obeys a stronger bound consistent with the isolated qubit result. Let us note that in an $n$-qubit system there are in fact nonlocal modes that relax with the rate $n/\tau_{SB}$, or $n^{2}/\tau_{SB}$ for correlated noise (e.g., superradiance), as well as modes whose relaxation is completely suppressed (e.g., subradiance, or decoherence-free subspaces, as mentioned above). These phenomena are well known [58] and experimentally documented [59, 60, 61, 62].

The benefit expressed by Eq. (91) is that the requirement on the single qubit $\tau_{B}/\tau_{SB}$ value switches from being exponential in the system size for the Davies-Lindblad case to polynomial for the CGME and the Redfield master equation. In terms of the coupling strength, we found these ranges of applicability at the relevant times of $\Lambda^{-1}$ to be between $\sqrt{n}g\tau_{B}\ll 1$ and $ng\tau_{B}\ll 1$ depending on the correlations in the noise. For local observables the earlier range $g\tau_{B}\ll 1$ for times $\sim\tau_{SB}$ is expected to hold. Here the inequality notation ($\ll$) guarantees that the error $\|\delta\rho\|_{1}$ or $\|\delta\rho_{\text{loc}}\|_{1}$ is small.

4 The CGME with a time-dependent system Hamiltonian

We start again from the total Hamiltonian (3), but this time we let the system Hamiltonian be time-dependent: $H(t)$. We denote the interaction term $A\otimes B$ by $V$. To go to the interaction picture, we now need to specify both the start time $t_{i}$ and the final time $t\leq t_{f}$. The unitary propagator $U_{tt_{i}}$ satisfies

[TABLE]

The formal solution is

[TABLE]

where $\mathcal{T}$ represents forward time-ordering. Some basic additional facts are:

[TABLE]

The relation between the Schrödinger picture and the interaction picture is:

[TABLE]

Here $B(t)$ is given by the bath Hamiltonian in the same way as in the previous derivation. We again assume the bath state to be stationary for simplicity, although the derivation will go through without it. At this point we repeat the steps involved in Eqs. (29)-(42b), with the only change being that the $A(t)$ operators now acquire a second time variable. Namely, Eq. (42b) — the Redfield equation in the interaction picture — is replaced by

[TABLE]

where we again introduced the Markov approximation $\rho_{B,I}(\tau)\mapsto\rho_{B,I}(t)$.

Let us digress briefly to derive the time-dependent Redfield equation. To do so, we transform Eq. (96) out of the interaction picture:

[TABLE]

Changing variables $t-\tau=t^{\prime}$ and switching the $dt^{\prime}$ integration limit to $t\to\infty$ we obtain:

[TABLE]

which is the time-dependent Redfield equation.

We now proceed with essentially the same series of transformations we used in the derivation for the time-independent Hamiltonian in Sec. 3.1. The main difference is that there are now additional time variables. Similarly as in the transition to Eq. (56), we time-average the r.h.s. by shifting $t\mapsto t+t_{1}$ and applying $\frac{1}{T_{a}}\int_{-T_{a}/2}^{T_{a}/2}dt_{1}$. The new equation is:

[TABLE]

and as long as $T_{a}$ is small compared to the fastest timescale over which $\rho_{I}(t)$ changes, the associated error $\|\rho_{BM}(t)-\rho_{BMT}(t)\|_{1}$ is small (see Sec. 6.4 for a detailed analysis of the time-independent case). We now rotate out of the interaction picture using Eq. (95):

[TABLE]

Next, we change variables to $t_{2}=\tau-t$:

[TABLE]

The transformations we performed after coarse graining of the Redfield equation (96) were exact. We know what to do to recover complete positivity: as in the transition to Eq. (65a) in the time-independent case, we need to neglect a part of the integral by changing the lower limit from $-t$ to $-T_{a}/2$:

[TABLE]

where we defined

[TABLE]

The same estimates on the error incurred by doing this apply, i.e., this introduces an error of order $\|\mathcal{E}_{p}(\text{const}\cdot\tau_{SB})\|=O(\tau_{B}/T_{a})$. Complete positivity yet remains to be demonstrated. To do so, we first note that we may interchange the order of integration since $\int_{-T_{a}/2}^{T_{a}/2}dt_{1}\int_{-T_{a}/2}^{t_{1}}dt_{2}=\int_{-T_{a}/2}^{T_{a}/2}dt_{2}\int_{t_{2}}^{T_{a}/2}dt_{1}$. We then swap the integration variables $t_{1}$ and $t_{2}$ and obtain:

[TABLE]

so that:

[TABLE]

and note that $X$ and $Y$ are both Hermitian. Also, $Q_{1}=\frac{1}{2}(X+iY)$ and $Q_{2}=\frac{1}{2}(X-iY)$, so that

[TABLE]

Combining Eq. (102) with the expressions in Eqs. (105)-(106), we thus find

[TABLE]

where

[TABLE]

Finally, we establish complete positivity by again introducing the Fourier transform $C(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-i\epsilon t}\gamma(\epsilon)d\epsilon$ [Eq. (6a)], where $\gamma(\epsilon)>0$ (Appendix A). In analogy with Eq. (20) for the time-independent case, let us also define

[TABLE]

We then obtain from Eq. (107) a master equation in exactly the same form as the time-independent CGME [Eq. (19)], except that all the operators now depend on time:

[TABLE]

This result, which is in a manifestly completely positive form, establishes that the CGME is valid also for arbitrary time-dependent system Hamiltonians. The time-dependent CGME, Eq. (110) [also stated earlier as Eq. (22)], is one of the main new results of this work.

5 Applications and Examples

5.1 Application: dynamical decoupling

One interesting application of our time-dependent CGME is to study methods to reduce decoherence. One such method is dynamical decoupling (DD) [63, 64] — the use of pulse sequences to cancel out the system-bath interaction (for reviews see [65, 66]). Even though the master equation has a Markovian appearance, the underlying bath can have a nonzero correlation time, and can therefore describe the effects of DD. At one level this might appear to be a counterintuitive result since it is widely accepted that in the limit $\mathcal{C}(t)=\delta(t)$ DD is useless, and hence one might naively expect that one needs to use non-Markovian methods (e.g., the Magnus expansion to second order or higher [67]) to describe the benefits of DD. However, it has recently been recognized that DD can in fact work in the Markovian setting as well, specifically the Davies-Lindblad master equation [68], and a more abstract semigroup framework [69] (see also Refs. [70, 71]). Our contribution here is to establish this in the framework of the Redfield master equation and, more importantly, the CGME. This ensures a wider range of applicability than the Davies-Lindblad master equation result. Since our result is obtained in a strictly Markovian setting, we differ with the conclusion reported in Ref. [72], that “success, however limited, of DD is a meaningful concept of non-Markovianity".

We consider two examples, both for a qubit with $H=0$ coupled to a purely dephasing environment:

[TABLE]

where $Z\equiv\sigma^{z}$ is the Pauli $z$-matrix. The noise model (111) can be decoupled with a sequence of $X$ pulses spaced by $\Delta t$. Here by a pulse we mean an instantaneous unitary being applied to the qubit, through the action of an externally controlled, time-dependent qubit Hamiltonian: instead of $H=0$ we have a time-dependent system Hamiltonian $H(t)=\frac{\pi}{2}\sum_{j}\delta(t-j\Delta t)X$, where the prefactor of $\pi/2$ is chosen so that $X$-pulses are generated.

5.1.1 Example 1: time-dependent Redfield master equation with a rectangle bath correlation function

Consider the rectangle bath correlation function

[TABLE]

where we have explicitly introduced the coupling strength $g$. To analyze this problem it will suffice to use the time-dependent version of the Redfield equation, Eq. (98). We thus need to compute

[TABLE]

and note that the time-evolution of the operator $A=Z$ under $X$ pulses results in $A(t-t^{\prime},t)=\pm Z$ with the sign alternating every $\Delta t$:

[TABLE]

Taking the integral, we find:

[TABLE]

where $c(t)$ is the remainder after the cancellation of different signs. To show why $|c(t)|\leq\Delta t$, note that there is no way to accumulate $>\Delta t$ by integrating any interval of the function $(-1)^{\lceil x/\Delta t\rceil}$

[TABLE]

To see this, we note that a translation $a,b\to a+\Delta t,b+\Delta t$ changes the sign of the integral. We then split the integration interval $[a,b]$ into three subsets: $\{[a+2n\Delta t,a+(2n+1)\Delta t]\},\leavevmode\nobreak\ \{[a+(2n+1)\Delta t,a+2(n+1)\Delta t]\},\leavevmode\nobreak\ [a+2\Delta t(n_{\text{max}}+1),b]$ where $n_{\text{max}}$ has been chosen closest to $b$: $n_{\text{max}}=\lceil(b-a)/2\Delta t\rceil-2$. The integrals over the first two subsets cancel due to the translation property, and the remaining subset $[a+2\Delta t(n_{\text{max}}+1),b]$ has length $L<2\Delta t$. Note that for $L=0$ and $L=2\Delta t$ the integral vanishes. Also note that the derivative of the integral with respect to $b$ is $(-1)^{\lceil b/\Delta t\rceil}$, and its absolute value is $1$. Any function with such a derivative, up to two inflection points and zeros at the end of the interval $[0,2L]$ stays within $[-\Delta t,+\Delta t]$, and the bound is saturated when $b=m\Delta t$.

We compare this with the case without pulses:

[TABLE]

Here $\tau_{c}$ enters the correlation function [see Eq. (112)]. Thus the decoherence rate is reduced by a factor of at least $\Delta t/\tau_{c}$. For $(g\tau_{c})^{2}\ll 1$ our analysis is accurate with the error bounds discussed above.

However, we note that the bath correlation function $\mathcal{C}(t)=g^{2}\theta(\tau_{c}-|t|)$ used here for simplicity is not physically permissible since its Fourier transform $\gamma(\omega)<0$ for some $\omega$. This leads to possibly nonpositive $c(t)$ for some $t$, and then the equation does not describe a CP map. To address this we next perform a more careful calculation with a realistic bath.

5.1.2 Example 2: time-dependent CGME with an Ohmic spectral density

Now consider a physically permissible environment with a spectral density that satisfies the KMS condition (10). We specify $\gamma(\omega)$ below. We wish to use the time-dependent CGME [Eq. (22)] to study the same dynamical decoupling protocol as in the previous example. We first need to compute $A(t+t_{1},t)=U^{\dagger}(t+t_{1},t)ZU(t+t_{1},t)$, where $U(t+t_{1},t)=\exp[-i\int_{t}^{t+t_{1}}H(s)ds]$. This unitary is a product of as many $X$-pulses as fit in the interval $[t,t+t_{1}]$, which we can express as

[TABLE]

where $m_{i}$ and $m_{f}-1$ respectively count the number of pulses in the intervals $[0,t]$ and $[0,t+t_{1}]$. Conjugating $Z$ by these pulses results in alternating signs, depending on the parity of $m_{f}-m_{i}$:

[TABLE]

This yields the time-dependent Lindblad operators $A_{\epsilon}(t)$ via Eq. (23):

[TABLE]

Note that the Lamb shift [Eq. (24)] is proportional to $Z^{2}=\mathbbm{1}$, so it drops out. The Lindblad operators without DD are simply

[TABLE]

We may thus write the time-dependent CGME as:

[TABLE]

where

[TABLE]

with $x=\text{DD}$ or $x=\text{noDD}$. The decoherence rate is thus suppressed by the factor

[TABLE]

The inner integral can be evaluated under certain simplifying assumptions. For example, let us assume that $T_{a}$ is an even integer multiple of $\Delta t$: $T_{a}=2k\Delta t$, and let us consider only times $t$ which are integer multiples of $\Delta t$: $t=\ell\Delta t$. Then $\lceil(t+\zeta T_{a})/\Delta t\rceil-\lceil t/\Delta t\rceil=\lceil 2k\zeta\rceil$, and the dependence on $t$ cancels out, i.e., the suppression factor $\xi(t)$ is time-independent. We can then rewrite the inner integral as:

[TABLE]

where in the last line we assumed for further simplicity assume that $k$ is even ($k=2k^{\prime}$), so that $\sin\left(\frac{k\pi}{2}+\frac{\epsilon T_{a}}{2}\right)=\pm\sin\left(\frac{\epsilon T_{a}}{2}\right)$. Writing $t=\ell\Delta t$ and $T_{a}={4k^{\prime}}\Delta t$, we thus have:

[TABLE]

To guarantee that DD will cause suppression [$\xi(\ell\Delta t)<1$] it is sufficient for the spectral density to have a high-frequency cutoff $\omega_{c}$ such that $|\tan\left(\omega\Delta t\right)|<1$, i.e.,

[TABLE]

in agreement with well established results [63].

An important physical example is a bath with an Ohmic spectral density

[TABLE]

where $\omega_{c}$ is a high-frequency cutoff and $\kappa$ is a positive dimensionless constant. This spectral density satisfies the KMS condition. The associated bath correlation function can be expressed in terms of the Polygamma function (see Appendix I of Ref. [31]), but unfortunately evaluating $\tau_{SB}$ and $\tau_{B}$ analytically using Eq. (11) is not possible in this case, so we simply pick a convenient value of $T_{a}$ ($=4k^{\prime}\Delta t$) rather than its optimal value $\sqrt{\tau_{B}\tau_{SB}/5}$. The suppression factor is plotted in Fig. 3 for a variety of parameter settings, including for $\Delta t>{\pi}/(4\omega_{c})$. It can be seen that the CGME correctly and consistently predicts that DD will result in the suppression of dephasing when Eq. (127) is satisfied. Thus, the CGME can be used to study DD protocols despite its Markovian appearance.

5.2 Lambless master equations

It turns out that the Lamb shift terms containing $\operatorname{sgn}(t)$, such as Eq. (18b), are more demanding to compute numerically, but at the same time the interesting dynamics is often given by the other, relaxation terms. Therefore it is convenient to have a version of the equations in the “Lambless" regime. Even though these equations do not have any correctness guarantees, one can still use them as toy models, and their numerical solution is often significantly faster than their complete counterparts. Let us list the resulting equations:

• •

Lambless Davies-Lindblad [replacing Eq. (18)]:

[TABLE]

• •

Lambless Redfield [replacing Eq. (16)]:

[TABLE]

where

[TABLE]

This simplified form of $A_{f}$ follows from Eqs. (7) and (17) by substituting $S(\omega)=0$ (vanishing Lamb shift). When the KMS condition holds we may write $A_{f}=\frac{1}{2}\sum_{\omega}e^{-\beta\omega}\gamma(\omega)A_{\omega}$.

• •

Lambless CGME [replacing Eq. (19)]:

[TABLE]

where

[TABLE]

and we choose $T_{a,\text{opt}}=\sqrt{\tau_{B}\tau_{SB}/5}$, $\Delta\epsilon$, and $k^{*}$ as prescribed in Appendix D. We study this case numerically in Sec. 5.4.

The importance of the Lamb shift in open system dynamics has been analyzed before, e.g., in Refs. [31, 73]

5.3 Toy example of a slow bath at the boundary of the range of applicability, with explicit correlation function and spectral density

Another way to speed up numerical simulations is to have a spectral density that leads to an explicit, analytically computable form of the bath correlation function $\mathcal{C}(t)$. Note that in the “Lambless" CGME case discussed in Sec. 5.2 the only direct role played by $\mathcal{C}(t)$ is in the determination of $\tau_{B}$ and $\tau_{SB}$ [via Eq. (11)], and apart from this it suffices to specify only $\gamma(\omega)$. The complete equations are more complicated: the Lamb terms contain $\mathcal{C}(t)$, or a Cauchy principal value integral of $\gamma(\omega)$ [Eq. (8)].

We present such a toy example with the property of having a bath almost at the boundary of the slowest correlation functions allowed within the range of applicability of the CGME. Baths with the same $\tau_{B}$ can have a different time decay of $|\mathcal{C}(t)|$ at large times $t\gg\tau_{B}$. A bath with $|\mathcal{C}(t)|\sim 1/t^{2}$ will have a logarithmically divergent $\tau_{B}$, which is undesirable for a toy example.101010$|\mathcal{C}(t)|\sim 1/t^{2}$ is actually the borderline case realized by an Ohmic bath, for which the divergence has a cutoff depending on the total time of the experiment $T$ in Eq. (11b).

Consider a bath with the following spectral density, which satisfies the KMS condition, such that $|\mathcal{C}(t)|\sim 1/t^{4}$—almost as slow as we are allowed:

[TABLE]

where $a>1$, $b>1/2$. Since $\frac{1}{\tau_{SB}}=\int_{0}^{\infty}|\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-i\omega t}\gamma(\omega)d\omega|dt$ [Eqs. (6b), (11a)], the normalization factor is

[TABLE]

The difference between the two exponentials in Eq. (134) is such that for small $\omega$ the term linear in $|\omega|$ cancels, and there is no inflection in $\gamma(\omega)$. Note that for a single $e^{-\beta|\omega|}$ we would have $|\mathcal{C}(t)|\sim 1/t^{2}$, but after this cancellation we get $|\mathcal{C}(t)|\sim 1/t^{4}$.

Our goal is now to use this example to compare the actual numerical errors to the theoretical bounds derived in this work. These bounds involve many relaxations and inequalities, so we do not expect them to be particularly tight. Also of interest is to compare the predicted optimal averaging time $T_{a,\text{theory opt}}$ to the numerically optimal $T_{a}$. To be specific we set $a=1.01$, $b=0.6$, $\beta=4$, and find the normalization factor to be $\mathcal{N}\approx 21.0$. This choice of parameters gives rise to a maximum of $\gamma(\omega)$ at $\omega^{*}\approx 2.0$, and the peak is wide: for $0.15\leq\omega\leq 6.08,\leavevmode\nobreak\ \gamma(\omega)\geq 0.5\gamma(\omega^{*})$. This frequency range is where we will choose system transitions to lie; see Fig. 4.

For this choice of parameters we obtain:

[TABLE]

so that indeed $\mathcal{C}(t)\sim 1/t^{4}$. Integrals of $\mathcal{C}(t)e^{i\omega t}$ over various intervals can be expressed via the Meijer G-function. Its evaluation is a part of the numerical simulation of master equations discussed in this paper, so their computational cost is sensitive to any shortcuts we can find, and the analytic form above is helpful. We will comment on the specific runtimes below.

The bath parameter $\tau_{SB}$ can be chosen freely, while $\tau_{B}\approx 0.69$ is determined by the choice of $a,b,\beta$. Note that it has dimensions of $\beta$.

We choose $\tau_{SB}=10$ so that $\tau_{B}/\tau_{SB}=0.1\ll 1$, a necessary condition for our error bounds to be small. We find $T_{a,\text{theory opt}}=\sqrt{\tau_{B}\tau_{SB}/5}=0.97$. To compare different equations, we pick an arbitrary two-qubit Hamiltonian:

[TABLE]

and choose $|11\rangle$ as the initial state. We let the coupling to the bath be given by $V=\sigma_{1}^{z}\otimes B$.

For all of our results, there will be no way to know what the true state $\rho_{\text{true}}(t)$ is. The differences found will be within the error bounds (growing at least linearly with time) of all equations. However we can simulate Eq. (46), which is presumably closest to the true solution since only the Born and Markov approximations were made:

[TABLE]

We call this the “Original Redfield" equation (ORE). It is in fact possible that the CGME returns solutions closer to the true solution than this non-CP equation, but we note that in our derivation all the equations appeared as subsequent approximations to the ORE. Therefore we may characterize how good those approximations are by studying the difference $\rho_{C}-\rho_{OR}$. We first investigate the solution $\rho_{OR}(t)$, and find that it becomes negative (and hence unphysical) at $t\approx 2.57\tau_{SB}$. Thus, we choose to compare the solutions of the CGME and ORE in the interval $0\leq t\leq 2.56\tau_{SB}$, where $\rho_{OR}\geq 0$ [see Fig. 4, and hence we can assume that $c_{BM}=1$ [recall that this is always true for a CP master equation; see the discussion following Eq. (47b)]. We note that the computational cost of the ORE is the largest (by at least on order of magnitude) among our equations.111111Precomputing the analytic form of the r.h.s. takes about $14$ minutes on a contemporary desktop computer, and the solution of the differential equation itself  — about a minute. In contrast, the CGME takes $20$ seconds total, and the Davies and Redfield equations  — even less. These numbers are for the implementation in Mathematica, where precomputing the symbolic form of the r.h.s. is possible. Faster implementations in other programming languages are certainly possible. The Mathematica code for obtaining the plots of this section are available [74].

The second observation we make is about the norm of the r.h.s. of Eq. (138). In Fig. 5 we present a probability distribution of the norms of the right hand side $\|\mathcal{L}^{BM,I}_{t}(\rho_{\text{test}})\|_{1}$, where $\rho_{\text{test}}$ is taken from the Gaussian unitary ensemble (Hermitian matrices $X$ sampled from a probability distribution $\sim e^{-2^{n-1}\mathrm{Tr}X^{2}}$, where $n=2$ is the number of qubits), and subsequently normalized such that $\|\rho_{\text{test}}\|_{1}=1$. We also combine the data for different times, as if the time was sampled uniformly in $[0,2.56\tau_{SB}]$. We find numerically that the maximum value of the norm that can be achieved is $0.44\pm 0.01$. It would suggest a higher effective $\tau_{SB}$ than our definition:

[TABLE]

Such an adjusted $\tau_{SB}^{\text{eff}}\approx 9.14=1.33\tau_{SB}$ can be used to choose an adjusted $T_{a,\text{adj}}=\sqrt{\tau_{SB}^{\text{eff}}/\tau_{SB}}T_{a,\text{theory opt}}=1.12=1.15T_{a,\text{theory opt}}$. We could non-rigorously use the typical value of the norm instead, leading to $T_{a,\text{typ.}}\approx 1.7T_{a,\text{theory opt}}$. To see if this adjustment is justified, we vary $T_{a}$ in Fig. 5 and plot the average of the trace norm of the difference in the solutions over $t\in[0,2.56\tau_{SB}]$, i.e., $\frac{1}{T}\int_{0}^{T}\|\rho_{C}(t)-\rho_{OR}(t)\|_{1}dt$. We see that the true optimum is at even higher $T_{a,\text{num. opt}}\approx 2.87\approx 3T_{a,\text{theory opt}}$, and the equation matches Davies in the limit $T_{a}\to\infty$ [Fig. 5], as formally expected (see Appendix C). Even though our adjustments of $T_{a}$ are not very large, the error is very sensitive to them and improves significantly.

We illustrate the error in more detail in Fig. 6 by plotting the differences

[TABLE]

as a function of $t$, where the CGME is taken for $T_{a,\text{adj}}$ and $T_{a,\text{num. opt}}$. Also in Fig. 6 we plot our strongest upper bound on $\|\rho_{C,\text{adj}}(t)-\rho_{OR}(t)\|_{1}$ [Eq. (240), derived in Sec. 6 below]. We use $T_{a}\to T_{a,\text{adj}}$ and $\Lambda\to 4/\tau_{SB}^{\text{eff}}$ Comparing this bound to the observed error in Fig. 6, we see that even after this adjustment, the bound is still not tight. Nonetheless, the adjustment helps to match the optimum in Fig. 5 better, and we believe the bounds can be improved by future work in this direction.

In Fig. 6 we vary $\tau_{SB}$, and evaluate

[TABLE]

for $T=2.56\tau_{SB}$. To find $T_{a,\text{num. opt}}$, we vary $T_{a}$ for each value of $\tau_{SB}$ and choose the one that gives the smallest values for the corresponding expression in Eq. (141). We again add the similar average error of Davies to the plot. Since our family of CGMEs with different $T_{a}$ includes Davies as a $T_{a}\to\infty$ limit, numerical optimization will always give an improvement over Davies. A simpler adjusted $T_{a}$ that does not require optimization still outperforms Davies at moderate $\tau_{B}/\tau_{SB}$. We also note that the $\sqrt{\tau_{B}/\tau_{SB}}$ dependence of its error at moderate $\tau_{B}/\tau_{SB}$ is consistent with the intuition from our bounds, while the linear portion suggests that our bounds are not tight in some parameter regime.

5.4 Role of the Lamb shift in the CGME

We perform another test, where we drop the Lamb shift $H_{LS}$ of the CGME and obtain a corresponding solution $\rho_{CL}(t)$, as prescribed by Eqs. (132) and (133). We then plot $\|\rho_{C,\text{adj}}(t)-\rho_{CL,\text{adj}}(t)\|_{1}$ as a function of time, as shown in Fig. 7. We again pick $T_{a,\text{adj}}=1.15T_{a,\text{opt}}$. We observe in Fig. 7 that the Lamb shift has the same role as in the Davies equation: it affects the phases of the off-diagonal matrix elements of $\rho_{C}(t)$ in the eigenbasis. Since off-diagonal elements decay to zero for the Davies case, the effect of their phase on the norm difference first grows proportional to the phase, then disappears with the magnitude of those elements.

6 Error bounds and estimates

In this section we give the details of the various error bounds and estimates mentioned without proof in our derivation of the CGME in Sec. 3. Our strategy is to derive the error terms successively, and then add them up using the triangle inequality. We work our way up through the various approximations, starting from Born in Sec. 6.2, followed by Markov in Sec. 6.3, then time-averaging (or coarse-graining) in Sec. 6.4. In the course of the latter, we bound the error due to dropping part of the integration domain to achieve complete positivity, and analyze the optimization of $T_{a}$. But we start with a proof of Lemma 1, which plays a central role in our analysis.

6.1 Bound on the difference of solutions of differential equations differing by a bounded operator

We first present a proof of Lemma 1.

Proof.

The first part of the Lemma is concerned with the pair of differential equations $\dot{x}=\mathcal{L}_{t}x+\mathcal{E}$ and $\dot{y}=\mathcal{L}_{t}y$. Their formal solution is:

[TABLE]

Let $\Lambda>0$ be a constant defined as follows: $\sup_{\tau,x}\|\mathcal{L}_{\tau}(x)\|_{1}\leq\Lambda$ for all $x$ such that $\|x\|_{1}=1$. Since $\mathcal{L}$ is linear, and writing $x=\delta/\|\delta\|_{1}$, it follows that $\forall\delta$: $\|\mathcal{L}_{\tau}(\delta)\|_{1}=\|\delta\|_{1}\|\mathcal{L}_{\tau}(x)\|_{1}\leq\|\delta\|_{1}\Lambda$. Then, assuming $x(0)=y(0)$, and taking the operator norm of the difference yields, for $\delta=x(t)-y(t)$:

[TABLE]

Many different functions $f(t)\equiv\|x(t)-y(t)\|_{1}$ satisfy this inequality for all $t$, but they are all upper-bounded by the function $b_{1}(t)=\max\|x(t)-y(t)\|_{1}$ that saturates the inequality:

[TABLE]

The solution is

[TABLE]

This proves the desired bound on $\|x(t)-y(t)\|_{1}\leq b_{1}(t)$.

The second part of the Lemma is concerned with the more general, non-Markovian pair of differential equations $\dot{x}(t)=\int_{0}^{t}K_{t-\tau}(x(\tau))d\tau+\mathcal{E}$ and $\dot{y}(t)=\int_{0}^{t}K_{t-\tau}(y(\tau))d\tau$. The first part of the Lemma in fact follows from this case as a corollary, but we presented it for clarity. We again subtract the formal solutions:

[TABLE]

Define $\Lambda_{t}=\sup_{x}\|K_{t}(x)\|_{1}$ for all $x$ such that $\|x\|_{1}=1$. By linearity of $K_{t}$ it follows as above that $\forall\delta$: $\|K_{t}(\delta)\|_{1}\leq\Lambda_{t}\|\delta\|_{1}$. Therefore:

[TABLE]

The upper bound $b_{K}(t)$ on $\|x(t)-y(t)\|_{1}$ satisfies the equation:

[TABLE]

Taking the derivative once, we obtain:

[TABLE]

Since the r.h.s. is positive, $b_{K}(t)$ is a monotonic function. So $b_{K}(\tau)\leq b_{K}(t)$, and it follows that:

[TABLE]

In other words, we can again upper-bound $b_{K}(t)\leq b_{K}^{\prime}(t)$, such that

[TABLE]

Now let $\Lambda$ be a positive constant such that $\int_{0}^{t}\Lambda_{\tau}d\tau\leq\Lambda$ for all $t$. We thus arrive at

[TABLE]

for a new upper-bound $b_{K}^{\prime\prime}(t)\geq b_{K}^{\prime}(t)$. This is exactly the same as Eq. (144), and hence the Lemma still holds for non-Markovian equations. ∎

Note that for our purposes $\Lambda=4/\tau_{SB}$ is usually taken, thanks to Eq. (41, 42c), and that the non-Markovian bound will be generalized in Sec. 6.2. The exponential on the r.h.s. of Eq. (145b) may be troublesome. Under additional assumptions we can prove a much tighter linear in $t$ bound presented in Appendix E, yet we did not find a way to recast all of our results in that tighter form.

6.2 Born approximation error

We start with an estimate of the error associated with making the Born approximation, which we presented as the estimate $\|\mathcal{E}_{B}\|_{1}=O\left(\frac{\tau_{B}}{\tau_{SB}^{2}}\right)$ in Eq. (33). I.e., we would like to bound the error $\mathcal{E}_{B}=\dot{\rho}_{\text{true},I}(t)-\int_{0}^{t}K_{t-\tau}^{2,B}(\rho_{\text{true},I}(\tau))d\tau$. It is difficult to do so directly since $\rho_{\text{true},I}(t)$ is unknown, so instead we will first use a proxy for $\rho_{\text{true},I}(t)$, which we call $\rho_{B4,I}(t)$. The latter is the solution of the master equation obtained by iterating the substitution of $\rho_{\text{tot}}=\rho_{\text{tot}}(0)-i\int_{0}^{t}[V(\tau),\rho_{\text{tot}}(\tau)]d\tau$ [Eq. (29b)] into ${\dot{\rho}}_{\text{tot}}(t)=-i[V(t),\rho_{\text{tot}}(t)]$ [Eq. (29a)] to 4th order, and only then applying the Born approximation. After developing an understanding of the 4th order we proceed to the infinite series. We make certain assumptions about its convergence, and comment below on their compatibility with the known convergence radius of the Dyson series. Let us first define the type of bath for which we can carry out this procedure in order to arrive at a bound.

We refer to a bath as exponential if its correlation function $\mathcal{C}(t)$ decays exponentially. This type of assumption is both common and convenient. For example, as we will see in more detail below, for a Gaussian bath (one which satisfies Wick’s theorem) the times $t_{1}<t_{2}<t_{3}<t_{4}$ appear in the estimate of the Born error [they enter terms such as $\mathcal{C}(t_{1}-t_{3})\mathcal{C}(t_{2}-t_{4})$]. If $\mathcal{C}(t)$ is exponentially decaying with a characteristic time $\tau_{B}$, that immediately lets us conclude that only $t_{1},t_{2},t_{3},t_{4}$ that are all within $\sim\tau_{B}$ of each other contribute to the estimate, and the bound on the Born error acquires a small factor of $\tau_{B}/\tau_{SB}$. However, an exponentially decaying $\mathcal{C}(t)$ may be an unrealistically strong assumption (e.g., it is hard to satisfy at low temperatures). Therefore it is desirable to introduce a condition weaker than exponential decay, which we refer to as a non-exponential bath. Specifically, we consider the convergence of integrals such as $\int_{0}^{\infty}|\mathcal{C}(t)|t^{n}dt$ for some integer $n\geq 0$. Such a condition was used in Ref. [31] to control the error of the subsequent Markov approximation. Here we show how these integrals for $n=0,1$ control our estimate of the Born error for a Gaussian bath. The physical meaning of these integrals was defined in Eq. (11). A more general requirement on a non-Gaussian bath can be derived in the same way. However, that requirement will involve a higher-order correlation function that does not have a straightforward and concise interpretation.

Recall the steps of the derivation where we make the Born approximation, leading to Eq. (39). For brevity, we drop the interaction picture subscript, replace the time argument by a subscript, and place the approximation level $B$ in the superscript: $\rho_{B,I}(t)=\rho_{t}^{B}$. Written out explicitly, Eq. (39) in terms of $\rho_{t}^{B}$ is:

[TABLE]

The formal solution is:

[TABLE]

We substitute this solution back into Eq. (153), a necessary step because our estimate of the Born error will come from the 4th order terms, so we need to subtract the original solution. The full expression to be subtracted is:

[TABLE]

Now, these terms constitute only a part of the expression for the derivative of the true density matrix ${\dot{\rho}}^{\text{true}}_{t}$ to the $4$th order in the system-bath interaction. The full $4$th order expression can be obtained, as mentioned above, by substituting the solution of the total evolution, $\rho_{\text{tot}}=\rho_{\text{tot}}(0)-i\int_{0}^{t}[V(\tau),\rho_{\text{tot}}(\tau)]d\tau$ into ${\dot{\rho}}_{\text{tot}}(t)=-i[V(t),\rho_{\text{tot}}(t)]$ three consecutive times (to generate a 4th order commutator), and only then applying the Born approximation. We call this $\rho_{t}^{B4}$ since it will differ from $\rho_{t}^{B}$ in the 4th order in the interaction $V(t)$. Doing so, we obtain:

[TABLE]

Defined in this way, ${\rho}^{B4}_{t}$ is close to the $4$th order of the true state ${\rho}^{\text{true}}_{t}$ up to terms of $6$th order and higher, so we use $\rho^{B4}_{t}-\rho^{B}_{t}$ to estimate the actual leading order error $\rho^{\text{true}}_{t}-\rho^{B}_{t}$. We proceed to expand all the commutators:

[TABLE]

There are 16 terms total in the last line, and the order of $B_{t},B_{\tau},B_{\theta},B_{m}$ and $\rho_{\text{b}}$ is exactly the same as the order of $A_{t},A_{\tau},A_{\theta},A_{m}$ and $\rho^{B}_{m}$.

Let us now assume that the bath is Gaussian. By definition, a Gaussian bath obeys Wick’s theorem (or Isserlis’ theorem [75]), which states that at all orders the higher order correlation functions decouple into products of two-point correlation functions. For example, for the 4-point correlation function:

[TABLE]

By definition, $\mathcal{C}(t-\tau)=\mathrm{Tr}_{\text{b}}B_{t}B_{\tau}\rho_{\text{b}}$. Thus:

[TABLE]

We see that in Eq. (155) only the first term was present. One can check that the order of $A$ is exactly the same, so upon subtracting Eq. (155) from Eq. (158) only the two terms $\mathcal{C}(t-\theta)\mathcal{C}(\tau-m)+\mathcal{C}(t-m)\mathcal{C}(\tau-\theta)$ remain. These extra terms can be interpreted in terms of the corresponding diagrammatic technique. Indeed, note that $t>\tau>\theta>m$. If in Wick’s theorem we have a pairing $t\leftrightarrow\tau,\leavevmode\nobreak\ \theta\leftrightarrow m$ and the arrows are drawn above the time axis, they do not relate to each other, and each arc can be anywhere along the time axis. But if the pairing is $t\leftrightarrow\theta,\leavevmode\nobreak\ \tau\leftrightarrow m$ or $t\leftrightarrow m,\leavevmode\nobreak\ \tau\leftrightarrow\theta$, then the two arcs are stuck together due to the condition $t>\tau>\theta>m$; see Fig. 8.

At the moment we have two equations: the Born equation in two forms Eq. (153) and Eq. (155), as well as the 4th order equation Eq. (157a):

[TABLE]

where $K(\rho)$ are corresponding linear superoperators. $K^{2,B}$ is the 2nd order (in the bath interaction) operator from Eq. (153), $K^{4,B}$ and $K^{4,B4}$ have terms up to the 4th order defined in Eq. (155) and Eq. (157a). Define

[TABLE]

for instance:

[TABLE]

Define

[TABLE]

The quantity that appeared in the derivation is $\mathcal{E}_{B}$. We will now bound $\|\mathcal{E}_{B4}(\rho^{\text{any}})\|_{1}$ and discuss its generalization $\|\mathcal{E}_{Bk}(\rho^{\text{any}})\|_{1}$. We postpone the estimate for $\|\mathcal{E}_{B}\|_{1}$ till we develop more advanced tools.

The difference $K_{t-\tau}^{4,B4}(\rho_{\tau}^{4,B4})-K_{t-\tau}^{4,B}(\rho_{\tau}^{B4})$ is given by the two second order pairings in Fig. 8, times the $16$ possible orders of $A_{t},A_{\tau},A_{\theta},A_{m}$ and $\rho_{m}$:

[TABLE]

We take the trace norm of both sides:

[TABLE]

where we used the triangle inequality $\|A+B\|_{1}\leq\|A\|_{1}+\|B\|_{1}$. We now use submultiplicativity [Eq. (38)]. Noting that the first norm in the submultiplicativity inequality is the operator norm, we obtain:

[TABLE]

The operator norm $\|A\|=1$ and we will set $c_{a}=\text{max}_{t}\|\rho^{\text{any}}_{t}\|_{1}$ such that $\|\rho^{\text{any}}\|_{1}\leq c_{a}$:

[TABLE]

In particular, $\mathcal{E}_{B4}(\rho^{B})$ is bounded with $c_{B}$ and $\mathcal{E}_{B4}(\rho^{B4})$ with $c_{B4}$. There are $16$ terms after the expansion of commutators:

[TABLE]

where we have also used that $|\mathcal{C}(t-m)|=|\mathcal{C}(m-t)^{*}|=|\mathcal{C}(m-t)|$. What is left is to bound the integrals in Eq. (164). It turns out that the following condition suffices:

[TABLE]

which is automatically satisfied for $\tau_{SB}$ and $\tau_{B}$ as defined in Eq. (11). We also change the order of integration variables using the condition $0\leq m\leq\theta\leq\tau\leq t$, e.g.:

[TABLE]

The first term in Eq. (164) is bounded as:

[TABLE]

while the second term in Eq. (164) is bounded as:

[TABLE]

Together, this yields:

[TABLE]

which is the first hint for the result given in Eq. (33). The actual error term appearing in Eq. (32) is $\mathcal{E}_{B}$ as defined in Eq. (163). We can consider repeating the construction of $\rho^{B4}$ presented here for $\rho^{B6},\rho^{B8}\dots$ as successive approximations to $\rho^{\text{true}}$. The corresponding errors can be defined:

[TABLE]

By employing the diagrammatic technique as in Fig. 8, the higher orders of perturbation theory in the system-bath coupling can be shown to give higher orders in $\tau_{B}/\tau_{SB}$. There are also extra factors of $4t/\tau_{SB}$ that appear in the higher orders. For instance, $\mathcal{E}_{B6}$ has the form:

[TABLE]

and the diagrams for the second term are shown in Fig. 9(a). We do not prove this form of $\mathcal{E}_{B6}$ here. We will bound all $t$-dependent terms of the first order in $\tau_{B}/\tau_{SB}$, by summing the diagrams in Fig. 9(b):

[TABLE]

for all $k$. Here we fix time to be a constant: $\exists\epsilon(t),C(t)$ s.t. $O\left(\frac{\tau_{B}}{\tau_{SB}}\right)\leq C(t)\frac{\tau_{B}}{\tau_{SB}}$ for all $\frac{\tau_{B}}{\tau_{SB}}\leq\epsilon(t)$. This is a weak statement. A stronger statement would include the time-dependence inside big-$O$, and have time-independent constants $\epsilon,C$:

[TABLE]

We do not know how to prove such a statement as of yet. The proof would rely on the diagrammatic calculation of subleading terms. It is straightforward, but cumbersome, and we do not attempt it here. Instead we make two assumptions:

In the $k\to\infty$ limit, there is a finite radius of convergence in $e^{4t/\tau_{SB}}\tau_{B}/\tau_{SB}$:

[TABLE] 2. 2.

The true solution is the limit of the perturbative series:

[TABLE]

To prove convergence, one needs to control the factorial number of terms appearing in Wick’s theorem. Such a proof was already presented by Davies [4]. The slight adjustments needed to include an integral condition such as Eq. (169) are a somewhat tedious problem left for a future study.

Our definition of $\mathcal{E}_{B}$ under this assumption contains $\mathcal{E}_{B\infty}$ discussed above, as well as an extra term:

[TABLE]

where $c_{a}=1$ since the true evolution is CP. The second term will turn out to be time-dependent, and quite tricky to estimate. We will return to this problem after we develop some extra machinery.

Let us first note the form of $K^{k,B}$, starting with $k=4$:

[TABLE]

The solution of

[TABLE]

is equivalent to solving a system of $4$th order equations that are our proxy for the true solution (which is $\infty$ order):

[TABLE]

We will prove that the solution of the above system remains close to $\rho^{B}$, the solution of the equation after the Born approximation:

[TABLE]

We will also present a straightforward generalization to all $k$ which will allow one to bound $\|\mathcal{E}_{B}\|_{1}$.

First write the equation in the compact form:

[TABLE]

Define a custom norm that is an infinity norm on the vector of operators, but for each operator the trace norm is taken:

[TABLE]

The triangle inequality is obeyed as required:

[TABLE]

The maximum of the custom norm of ${\bf{L}}$ over inputs of norm 1 is the same as before:

[TABLE]

so all the properties needed to prove Lemma 1 hold. We establish:

[TABLE]

in the original variables:

[TABLE]

Note that the error we are trying to bound is

[TABLE]

Our approximation to it is bounded as:

[TABLE]

where we have used the triangle inequality: $\|\rho^{(1)}_{\tau}-\rho^{B4}_{\tau}\|_{1}\leq\|\delta\rho^{(1)}_{\tau}\|_{1}+\|\delta\rho^{B4}_{\tau}\|_{1}$. Now we can repeat the above proof for arbitrary $k$:

[TABLE]

For appropriately defined ${\bf{L}}$, $\|\cdot\|_{c}$ the r.h.s. is still bounded by the same constant:

[TABLE]

which leads to

[TABLE]

and the error in the r.h.s. is bounded without change:

[TABLE]

Taking the $k\to\infty$ limit, we obtain:

[TABLE]

and, finally:

[TABLE]

where we have used that $c_{a}=1$ for $\rho^{\text{true}}$. This result is the sought-after justification of the estimate $\|\mathcal{E}_{B}\|_{1}=O\left(\frac{\tau_{B}}{\tau_{SB}^{2}}\right)$ given in Eq. (33).

An immediate consequence of Eq. (197) is that we can use it to obtain a bound on $c_{B}$, which we defined after Eq. (40) as an upper bound on $\|\rho_{B,I}\|_{1}$. Since

[TABLE]

we obtain:

[TABLE]

6.3 Markov approximation error

This subsection follows Ref. [76] initially. Before the Markov approximation but after the Born approximation, the master equation had the form given in Eq. (39). The error due the Markov approximation $\rho_{B,I}(\tau)\mapsto\rho_{B,I}(t)$, in transitioning to Eq. (42b), is:

[TABLE]

Since $\rho_{B,I}(t)-\rho_{B,I}(\tau)=\int_{\tau}^{t}{\dot{\rho}}_{B,I}(t^{\prime})dt^{\prime}$, we have

[TABLE]

where we used Eq. (40) in the second inequality. Now we can estimate the norm of the error:

[TABLE]

We found $c_{B}$ and saw that it is $O(1)$ in Eq. (200).

Let us now use the Markovian Lemma 1 with Eq. (203):

[TABLE]

A total expression for the Born and Markov error is [adding Eq. (197)]:

[TABLE]

which simplifies to:

[TABLE]

or substituting in $c_{B}$:

[TABLE]

The extra factor from $c_{B}$ can be absorbed into the big-$O$ notation:

[TABLE]

A loose bound follows from the above:

[TABLE]

Sometimes we wish to report the bound briefly as:

[TABLE]

Proof.

The original big-O notation states that for $e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\leq\epsilon$ the following inequality holds:

[TABLE]

We first note that $e^{\frac{4t}{\tau_{SB}}}\leq e^{\frac{8t}{\tau_{SB}}}$, so if $e^{\frac{8t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\leq\epsilon$, then $e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\leq\epsilon$ as well, and we are within the range of the original big-O notation. By taking $12(1+C\epsilon)$ as the new constant, we prove Eq. (210). ∎

6.4 Time averaging error and its optimization

6.4.1 General analysis of time-averaging

We consider the time-averaging operation, which leads to a more complicated bound on the difference of the solutions. We first prove a general result:

Lemma 2**.**

For any first order linear differential equation $\dot{\rho}(t)=\mathcal{L}_{t}\rho(t)$ there is a time-averaged version

[TABLE]

Note that we use a different symbol $\pi(t)$ for the solution of the time-averaged equation. Assume the same initial condition: $\rho(0)=\pi(0)$. Then the deviation between the two, $\delta\rho(t)\equiv\pi(t)-\rho(t)$, is bounded as:

[TABLE]

where $\Lambda,c_{\rho}$ are constants such that $\forall t,\rho_{0},\|\rho_{0}\|_{1}\leq 1:\leavevmode\nobreak\ \|\mathcal{L}_{t}(\rho_{0})\|_{1}\leq\Lambda$ and the solution $\rho(t)$ with the initial condition $\rho(0)=\rho_{0}$ is bounded as $\|\rho(t)\|_{1}\leq c_{\rho}$. The function $b_{2}(t)$ is:

[TABLE]

An immediate corollary is:

[TABLE]

The proof can be found in Appendix F. The $O(\Lambda T_{a})$ corollary enters the estimate for $\|\rho_{BM,I}(t)-\rho_{BMT,I}(t)\|_{1}$ we presented in Eq. (61b) in the derivation (using $\Lambda=4/\tau_{SB}$). We will also need a compact bound on:

[TABLE]

One can prove that the expression in the brackets is upper-bounded by

[TABLE]

This allows us to bound

[TABLE]

Note that the function $b_{2}(t)$ is monotonic: $b_{2}(t)\leq b_{2}(T_{a}/2)$. One can then loosen the bound in Lemma 2:

[TABLE]

6.4.2 Time averaging of the Redfield ME

We now focus on the specific equations of interest to us. We already established in Eq. (42c) that we can choose $\Lambda=4/\tau_{SB}$ (recall that $\Lambda$ needs to upper bound the r.h.s.: $\max_{t,\rho_{\text{test}},\|\rho_{\text{test}}\|_{1}=1}\|\mathcal{L}^{BM,I}_{t}(\rho_{\text{test}})\|_{1}\leq\Lambda$). The total time-averaging error (219) is now bounded as:

[TABLE]

We will use both the compact bound of Eq. (220), and the tightest bound available which is a repetition of Lemma 2:

[TABLE]

6.4.3 Bound on the error due to enforcing complete positivity

Recall that after the time-averaged equation is obtained, we need to drop an extra portion of the integral to obtain complete positivity [the transition to Eq. (65a)]. Before that, Eq. (64) was:

[TABLE]

which changes into Eq. (65a) (same terms, but the portion of the integral is subtracted):

[TABLE]

Collecting $\sum_{\omega}$ into $A(t)$ we obtain:

[TABLE]

so that:

[TABLE]

Since the resulting equation is CP, we used $\|\rho_{C}\|_{1}=1$. We would like to bound this norm by splitting the integration domain into different times. The corresponding domains are illustrated in Fig. 10. For $t<T_{a}/2$ we note that the argument of $\mathcal{C}(\tau^{\prime}-t^{\prime})$ can be both positive and negative. Thus a factor of $2$ arises when we bound $\int dt^{\prime}|\mathcal{C}(\tau^{\prime}-t^{\prime})|\leq 2/\tau_{SB}$:

[TABLE]

For $t>T_{a}/2$, we perform a change of variables: $\tau^{\prime}-t^{\prime}=-\tau,\leavevmode\nobreak\ t^{\prime}+T_{a}/2=\theta$:

[TABLE]

We now change the order of integration, noting that $0<\theta<\tau$, thus obtaining the following upper bound:

[TABLE]

Recall that the total evolution time $T$ entered the definition $\frac{\tau_{B}}{\tau_{SB}}=\int_{0}^{T}t|\mathcal{C}(t)|dt$ [Eq. (11b)]. We can usually set $T=\infty$, except for an Ohmic bath [see Eq. (128)], which is a borderline applicability case, with a logarithmic in $T$ divergence of the integral. We define an evolution time with a collar of $T_{a}/2$: $T=\max_{t}(t+T_{a}/2)$:

[TABLE]

Collecting everything, we obtain:

[TABLE]

where $\theta(x)=1,\leavevmode\nobreak\ x>0$ and zero otherwise.

Recall that $\mathcal{E}_{p}$ is defined in Eq. (65a) as the extra term appearing if the CGME is written with the same expression as Eq. (64) in the r.h.s. We can now use the same idea as in the proof of Lemma 1 (Sec. 6.1). Let us introduce $\delta\rho_{2}=\rho_{BMT}-\rho_{C}$. A function $b_{p}$ that upper-bounds $\|\delta\rho_{2}\|_{1}$ is the solution of the differential equation

[TABLE]

where we have used $\Lambda\geq\max_{t,\rho_{\text{test}},\|\rho_{\text{test}}\|_{1}=1}\|\mathcal{L}^{BM,I}_{t}(\rho_{\text{test}})\|_{1}$ from Eq. (42c) as the upper bound on the norm of the superoperator on the r.h.s. of Eq. (64). Indeed, by applying the triangle inequality to the definition of the r.h.s. of Eq. (64) we get:

[TABLE]

The solution to Eq. (232) is:

[TABLE]

Substituting Eq. (231), we obtain two cases:

[TABLE]

which we can summarize as

[TABLE]

For a concise form we change the first term using $1-e^{-x}\leq x$:

[TABLE]

Let us now focus on the difference with $\rho_{BM,I}(t)$ which is the solution of Eq. (42b). We will keep $\Lambda\geq\max_{t,\rho_{\text{test}},\|\rho_{\text{test}}\|_{1}=1}\|\mathcal{L}^{BM,I}_{t}(\rho_{\text{test}})\|_{1}$ in the bound in case it is $<4/\tau_{SB}$. Using Eqs. (220) and (238), we have:

[TABLE]

The tightest bound available that starts at zero at $t=0$ is the combination of Eqs. (221) and (237):

[TABLE]

where for $t<T_{a}/2$

[TABLE]

We may optimize $T_{a}$ to minimize the bound Eq. (239), which we proceed to do next.

6.4.4 Optimization of $T_{a}$

We have used $\Lambda\leq 4/\tau_{SB}$ to obtain the form of the bound in Eq. (239), but we kept both $\Lambda$ and $\tau_{SB}$ to demonstrate the possibility of making the bound tighter by using different values for $\Lambda$ and $4/\tau_{SB}$ evaluated independently. To that end, we denote $c_{\Lambda}=4/\Lambda\tau_{SB}\geq 1$ and treat it as a constant below. At large $t\gg\Lambda^{-1},\leavevmode\nobreak\ t>T_{a}/2$ the leading terms of Eq. (239) are:

[TABLE]

We can loosen the bound again via $e^{-\Lambda T_{a}/2}\leq 1$:

[TABLE]

Minimizing the expression in the brackets, we find that $T_{a}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/(3c_{BM}+2)}$ [$=\sqrt{\tau_{SB}\tau_{B}/5}$ if $c_{BM}=c_{\Lambda}=1$, which is Eq. (25)] for any $\tau_{B}/\tau_{SB}$ leads to an optimal bound at large enough times. The asymptotic behavior of the bound is:

[TABLE]

One may be concerned that the relaxation $e^{-\Lambda T_{a}/2}\leq 1$ reduced the quality of our bound. We checked this for $c_{BM}=c_{\Lambda}=1$, though we do not present that calculation here. The result is that after the optimization of the original Eq. (240), which is a more complicated function of $T_{a}$, the bound is unchanged to the leading order, and in fact acquires a multiplicative correction of the form $\left(1-\frac{17}{15\sqrt{5}}\frac{\sqrt{\tau_{B}}}{\sqrt{\tau_{SB}}}+O(\tau_{B}/\tau_{SB}))\right)$. We do not pursue this route since the effect is small, and the calculation is unreasonably long. Instead we use $T_{a}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/(3c_{BM}+2)}$ obtained above and find that the full expression for the error from Eq. (239) can be bounded as

[TABLE]

and it holds for all positive $\tau_{B}/\tau_{SB}$, $t/\tau_{SB}$ and $c_{BM},c_{\Lambda}\geq 1$. The derivation can be found in Appendix G. This simple form allows us to finally bound $c_{BM}$. Since for a CP coarse-grained equation $\|\rho_{C}\|_{1}=1$:

[TABLE]

This is a quadratic equation on $c_{BM}$. Its $\geq 1$ solution is:

[TABLE]

And here the big-$O$ notation means that there exist constants $x_{0},C$ such that for $0\leq x\leq x_{0},O(x)\leq Cx$. We note that the choice of $T_{a}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/(3c_{BM}+2)}$ depends on the total evolution time, if we use the above bound on $c_{BM}$. We would like to approximate it by a time-independent $T_{a}$. Note that for $\sqrt{\frac{\tau_{B}}{\tau_{SB}}}e^{\Lambda t}\ll 1$ we can approximate $c_{BM}\approx 1$ and consequently $T_{a}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/5}$. Note also that our construction was optimal for asymptotic $e^{\Lambda t}\gg 1$ (for the discussion of small times, see Appendix H). We combine the two conditions as $1\ll e^{\Lambda t}\ll\sqrt{\frac{\tau_{SB}}{\tau_{B}}}$. In units of actual time, this corresponds to

[TABLE]

This is a fairly narrow range for any realistic parameters, but it still becomes $[1,\infty]$ in the weak coupling limit. We choose to work with $T_{a}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/5}$ that is approximately optimal within that range. From a purely mathematical point of view, there was no particular reason to choose it this way. From a physical point of view, we would like to work with errors of $\leq 10\%$, which would correspond to a range like this. We have used $T_{a,\text{adj}}=\sqrt{c_{\Lambda}\tau_{SB}\tau_{B}/5}$ and $T_{a,\text{theory}}=\sqrt{\tau_{SB}\tau_{B}/5}$ for the numerical simulations in Section 5.3. For theoretical purposes of the remainder of this Section, we can always use $c_{\Lambda}=1$ which is proven. We will keep the notation $e^{\Lambda t}=e^{4t/\tau_{SB}}$ for brevity. The error bound for $T_{a}=\sqrt{\tau_{SB}\tau_{B}/5}$ is derived as Eq. (367):

[TABLE]

The constant $c_{BM}$ can be obtained numerically independently, and then one just needs to substitute it. For the theoretical use of the CGME, that route is not available, and instead one needs to substitute its upper bound Eq. (247). For the rest of this section, we will loosen the resulting bound to make the expression more compact. Since we optimized for $e^{\Lambda t}\gg 1$, dropping constants will not be significant in the parameter regime of interest. We will also update the expression for $c_{BM}$ using $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ in Eq. (247):

[TABLE]

Now we simplify Eq. (249) by relaxing the constant terms:

[TABLE]

With the relaxed form of $c_{BM}$ above this becomes:

[TABLE]

The error near $t=0$ is doubled as compared to the longer formula (249). Moreover, at each of decreasing timescales $t\sim\Lambda^{-1}$, $t\sim T_{a}$, $t=0$ the error gets progressively worse in our approximations, as compared to the tightest bound possible: the error at $t=0$ is supposed to be [math]. Nevertheless, we proceed, as this route allows us to obtain a concise form and the very small $t$ region is less important for our purposes. The bound can be further relaxed with rounding, using $6e^{2}\alpha\leq 13-2\sqrt{5}e+\frac{9e^{3}\alpha^{2}}{13-2\sqrt{5}e}$:

[TABLE]

where we plugged in $\Lambda=4/\tau_{SB}$. We will use this form while computing the total error. We still need to include the Markov and Born errors. Adding the bound in Eqs. (253) and (211), and using the triangle inequality, we finally obtain:

[TABLE]

Recall that this bound is applicable within the region $e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\leq C$ for some constant $C$ (from big-$O$), but within that region any triple $\{t,\tau_{B},\tau_{SB}\}\geq 0$ is allowed. In particular, for sufficiently small $\tau_{B}/\tau_{SB}$, the error at any value of $t/\tau_{SB}$ can be shown to be small. As a corollary, we proved the validity in the weak coupling limit, but the above inequality is more general than that. In particular, this proves Eq. (76), which uses $T_{a}=\sqrt{\tau_{B}\tau_{SB}/5}$.

We will finally show how to simplify Eq. (254) to the form $O\left(e^{\frac{6t}{\tau_{SB}}}\sqrt{\frac{\tau_{B}}{\tau_{SB}}}\right)$ presented in Eq. (28). Denoting $x=e^{\frac{6t}{\tau_{SB}}}\sqrt{\frac{\tau_{B}}{\tau_{SB}}}$, we can relax the bound as follows:

[TABLE]

We then use the property of Big O notation: for any function $f(t/\tau_{SB},\tau_{B}/\tau_{SB})$, if $f(t/\tau_{SB},\tau_{B}/\tau_{SB})=O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}\right)$, then also $f(t/\tau_{SB},\tau_{B}/\tau_{SB})=O(x^{2})$. We obtain:

[TABLE]

From this, the result of Eq. (28) immediately follows:

[TABLE]

6.5 The Davies-Lindblad error

We write the CGME with an explicit dependence on $T_{a}$:

[TABLE]

The error bound is known for arbitrary $T_{a}$: see Eq. (242). We now wish to compare the CGME with finite $T_{a}$ to the limit $T_{a}\to\infty$ of the same equation:

[TABLE]

In Appendix C, we show that this limit is just the Davies-Lindblad equation:

[TABLE]

We also show in Appendix C that the difference between the right-hand sides is bounded as:

[TABLE]

Here the system Hilbert space dimension $2^{n}$ appeared, and we note that $1/\delta E$ could contain an additional exponential dependence on $n$. The accumulated error in the solution is bounded as:

[TABLE]

In the long time limit the leading terms are:

[TABLE]

Assuming that $\frac{4^{n}}{\delta E}\gg\tau_{B}$, we find the optimal $T_{a}=2^{n}O(\sqrt{\tau_{SB}/\delta E})$. Using that together with a bound for $t\geq T_{a}/2$ we get:

[TABLE]

Assuming a constant system size $n=O(1)$, this is the form presented in Eq. (27), our second main error bound.

Then, assuming $\frac{4^{n}}{\delta E}\gg\tau_{B}$, we can neglect the second term. Recall the discussion in Sec. 2.2: $g=1/\sqrt{\tau_{B}\tau_{SB}}$, where $g$ is the coupling strength $g$. Thus, we can express $\sqrt{1/\tau_{SB}\delta E}=g\sqrt{\tau_{B}/\delta E}$. If we set $t<c\tau_{SB}$, we can write $\|\rho_{\text{true}}(t)-\rho_{\text{D}}(t)\|_{1}=O(\sqrt{1/\tau_{SB}\delta E})=O(g\sqrt{\tau_{B}/\delta E})$, as stated in Eq. (1) in the Introduction.

6.6 Redfield limits of integration error

There is one more error left to analyze: the error due to changing $\int^{t}\to\int^{\infty}$ in the Redfield equation. We already had the bound given in Eq. (47b):

[TABLE]

where $\epsilon_{T}=\tau_{SB}\int_{T}^{\infty}|\mathcal{C}(t^{\prime})|dt^{\prime}$ [Eq. (48)]. We augment Eq. (265) by a trivial bound:

[TABLE]

Since the bound on $\|\mathcal{E}_{l}\|_{1}$ is time dependent, we cannot directly use Lemma 1 for Eq. (47a). Instead, we repeat its proof and integrate the corresponding differential equations as was done in the proof of Lemma 2. We define a function $b_{R}(t)\geq\|\rho_{BM}(t)-\rho_{R}(t)\|_{1}$ (recall that $\rho_{BMl}=\rho_{R}$) that satisfies the following differential equation, analogous to Eq. (144):

[TABLE]

Integrating it, we obtain:

[TABLE]

As shown in Appendix I, this results in the error bound:

[TABLE]

In other words, the $O(\tau_{B}/\tau_{SB})$ error scaling of the master equation just after Born-Markov approximation, noted in Eq. (209) acquires a logarithmic correction.

We will now show how to simplify the total error of Redfield to obtain Eq.(26). We add the error of $\rho_{BM}$ from Eq. (211):

[TABLE]

We will set $\epsilon_{T}=0$ in this estimate and use the following bound on $c_{BM}$ from Eq. (250):

[TABLE]

where we denoted $x=e^{\frac{12t}{\tau_{SB}}}\frac{\tau_{B}}{\tau_{SB}}$. After a series of simplifications similar to the ones done after Eq. (255) we arrive at:

[TABLE]

Since big O notation is defined for $x\leq\epsilon$ for arbitrary small $\epsilon$, it follows that $\tau_{B}/\tau_{SB}\leq\epsilon$ is also small, so the logarithm is the leading term:

[TABLE]

This is Eq. (26).

7 Summary and Conclusions

Our primary goal in this work was to derive a Markovian master equation with a complete positivity guarantee that is in addition capable of incorporating arbitrary time-dependent driving, and that has a larger range of applicability than the Davies-Lindblad Markovian master equation. The master equation that achieves these desiderata is based on the coarse-graining approach introduced in Ref. [40], which only considered the time-independent (no driving) case. We rederived this result by showing that it results from the non-CP Redfield master equation after neglecting part of the time-integration domain, which restored complete positivity. We then included arbitrary time-dependent driving and obtained a new CGME that is still CP, and retains essentially the same form as its time-independent counterpart. This is our central new result, given in Eq. (22).

This time-dependent CGME depends on one free parameter, the coarse graining time $T_{a}$. We showed that it is optimal to choose it as proportional to the geometric mean of two key timescales, the bath correlation time $\tau_{B}$ and the fastest system decoherence timescale $\tau_{SB}$, given in Eq. (11). These timescales are related to the bath spectral density $\gamma(\omega)$ via the inequalities $2/\tau_{SB}\geq\gamma(\omega)$ and $2\tau_{B}/\tau_{SB}\geq|\gamma^{\prime}(\omega)|$. The $1$-parameter family of CGME equations involves a continuous integral, or equivalently, a sum over exponentially many terms in the frequency domain. Alongside this $1$-parameter family, we also derived an approximation in terms of $3$-parameter family ($T_{a}$, $\Delta\epsilon$, $k^{*}$), given in Eq. (78). We showed that if $[H,A]=O(1)$, where $H$ is the system Hamiltonian and $A\otimes B$ is the interaction with the bath [Eq. (3)], then there is no system-size dependence in $\Delta\epsilon,k^{*}$, where $k^{*}$ is the number of terms replacing the continuous integral.

Along the way we provided rigorous error bounds on three different master equations: Redfield, Davies-Lindblad, and the CGME. These bounds [Eqs. (26)-(28)] are expressed in terms of the dimensionless ratio $\tau_{B}/\tau_{SB}$, and in the Davies-Lindblad case also in terms of the quantity $1/\sqrt{\tau_{SB}\delta E}$, where $\delta E$ is the minimum level spacing. These results are valid even when $\|B\|$ diverges, so they apply in particular to, e.g., oscillator baths. When the system-bath coupling strength $g$ is finite we can write $\tau_{B}/\tau_{SB}=(g\tau_{B})^{2}$ and $1/\sqrt{\tau_{SB}\delta E}=g\sqrt{\tau_{B}/\delta E}$. These bounds provide a new perspective on the range of applicability of Markovian master equations: the Redfield and CG master equations apply in the regime where $g\tau_{B}\ll 1$, but the Davies-Lindblad master equation applies when in addition the minimum level spacing is large relative to $g^{2}\tau_{B}$. While these are all sufficient and not necessary conditions, if we informally view a large level spacing as a requirement, then it is a very onerous one indeed, since in many cases of interest the spacing shrinks rapidly (polynomially or even exponentially) with problem size, e.g., in adiabatic quantum computing [77]. Thus the guaranteed range of applicability of the Davies-Lindblad master equation is severely limited, while the CGME achieves complete positivity without imposing this cost.

To summarize the latter in the simplest possible terms, what we have shown is that the CGME and Redfield master equations are valid for fast baths, in the sense of a very simple integral condition on the bath correlation function. Ohmic and super-Ohmic baths satisfy this condition, provided the system-bath coupling strength is weak but constant in system size (for a local system). For the Davies-Lindblad master equation the coupling strength needs to be exponentially small in the system size. However, the Redfield equation and the CGME are not on equal footing: the Redfield equation provides no guarantee of complete positivity, while the CGME does.

We list a few potential applications of the CGME:

• •

The CGME opens the door to studies of open system dynamics with arbitrary time-dependent drives $H(t)$. In particular, the drive need not be adiabatic. In fact, the time-dependent CGME even correctly describes a system driven by $\delta$-function pulses, as demonstrated in Sec. 5.1 for dynamical decoupling. Moreover, the time-dependent CGME overcomes an inherent problem in applying the Davies-Lindblad ME to circuit model quantum computing, in particular in the context of fault-tolerant quantum error correction. The latter ME is incompatible with the assumption of arbitrarily fast gates needed there [43], but the time-dependent CGME can be consistently applied in this context. This opens the door to a rigorous analysis of fault-tolerance in a first-principles Markovian setting.

• •

Earlier studies of the lifetime of the toric code at finite temperature took advantage of the fact that the Davies generators are local for commuting models [78]. The CGME provides local generators of open system finite temperature dynamics for non-commuting models (see Sec. 3.6). This opens the door to studying the finite-temperature lifetime of, e.g., topological subsystem codes [79] in the thermodynamic limit.

• •

Many-body localization effects in quantum annealing — previously studied only in the closed system setting [80] — can now be extended to finite temperature open system quantum annealing.

We hope that this work, including the applications above, will stimulate further research into the uses of the CGME, the only master equation we know of that provides a complete positivity guarantee, is locally generated, and applies for arbitrarily driven open quantum systems.

Acknowledgements

The research is based upon work (partially) supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.

Appendix A Proof that $\gamma(\omega)>0$

Trivial case.

First consider $V=A\otimes B$. The function $\gamma(\omega)$ is the Fourier transform of the correlation function:

[TABLE]

Our derivations in this paper are for a stationary bath $\rho_{B}(t)=\rho_{B}$, or $[\rho_{B},H_{B}]=0$. Note that $[\sqrt{\rho_{B}},H_{B}]=0$ follows.

[TABLE]

Inserting the sum over bath eigenstates:

[TABLE]

Integration results in delta-functions:

[TABLE]

This is enough for all the results in this paper. A more complicated interaction with environment, however, requires a different proof.

General case.

We consider the case $V=\sum_{i}A_{i}\otimes B_{i}$, so that $\gamma(\omega)$ becomes a matrix. $\rho_{B}$ will still be assumed stationary: $[\rho_{B},H_{B}]=0$. Repeating the same steps, we arrive at

[TABLE]

Thus we need to prove that the matrix

[TABLE]

is positive. Let $v_{i}=\langle k|\sqrt{\rho_{B}}B_{i}|l\rangle$, then

[TABLE]

For positivity, we need to prove that for any vector $\vec{w}$ with elements $w_{i}$ it holds that $\sum_{ij}w_{i}M_{ij}w_{j}^{*}\geq 0$. Indeed:

[TABLE]

Thus

[TABLE]

as a matrix.

Bochner’s theorem

Our proof above used the eigendecomposition of the bath, and $\gamma(\omega)$ was a sum of Dirac $\delta$ functions. In practice, a smooth $\gamma(\omega)$ is used, corresponding to, e.g., a bath with a continuous spectrum. In this sense it is desirable to avoid using bath eigenstates and $\delta$-functions altogether. Here we give another proof that uses Bochner’s theorem, as suggested, e.g., in the textbook [2], and which avoids the issue mentioned above.

Lemma 3** (Bochner).**

If a function $\phi(t)$ is of positive type, i.e. $\forall n\in\mathbb{Z},\leavevmode\nobreak\ t_{i}\in\mathbb{R},\leavevmode\nobreak\ \xi_{i}\in\mathbb{C}$

[TABLE]

and also $\phi(0)$ is finite and $\phi(t)$ is continuous at $t=0$, then its Fourier transform (if it exists) is non-negative:

[TABLE]

We again consider many interaction terms $V=\sum_{i}A_{i}\otimes B_{i}$ (as in Sec. 3.7) and a stationary $\rho_{B}$ s.t. $[\rho_{B},H_{B}]=0$. Since $\gamma(\omega)$ is Hermitian we can diagonalize it using a unitary transformation:

[TABLE]

$D$ is diagonal so we only need to consider the diagonal elements (i.e., the eigenvalues of $\gamma$). Substituting $\gamma_{ij}=\int_{-\infty}^{\infty}e^{i\omega s}\mathcal{C}_{ij}(s)ds$, where $\mathcal{C}_{ij}(s)=\mathrm{Tr}[\rho_{B}B_{i}(s)B_{j}]$, gives

[TABLE]

We wish to show that $D_{\alpha}$ is non-negative for each $\alpha$. To do this we must consider the function in parenthesis. $D_{\alpha}$ is the Fourier transform of this function so if we can show that it is of positive type then $D_{\alpha}$ must be positive by Bochner’s theorem [81]. Define the following function with $\left\{t_{i}\right\}$ an arbitrary time partition:

[TABLE]

Note that

[TABLE]

Now use the property $\langle B_{\alpha}(s)B_{\beta}(0)\rangle=\langle B_{\alpha}(t)B_{\beta}(t-s)\rangle$ [Eq. (289c)] to write $f_{mn}^{\alpha}$ as

[TABLE]

We need to show that $f^{\alpha}$ is a positive matrix. For arbitrary $\ket{v}$ we have

[TABLE]

where $M_{\alpha}\equiv\sum_{i,m}v_{m}^{*}U_{\alpha i}\sqrt{\rho_{B}}B_{i}(t_{m})$.

We have established that $\bra{v}f^{\alpha}\ket{v}\geq 0$ for any time partition $\left\{t_{i}\right\}$. Therefore $D_{\alpha}$ is non-negative by Bochner’s theorem. Consequently, $\gamma\geq 0$ as a matrix since all its eigenvalues are non-negative.

Appendix B Lamb shift simplification

The Lamb shift was given in Eq. (69) as $H_{\mathrm{LS}}=\frac{i}{2T_{a}}\int_{-T_{a}/2}^{T_{a}/2}dt_{1}\int_{-T_{a}/2}^{T_{a}/2}dt_{2}\operatorname{sgn}(t_{1}-t_{2})\mathcal{C}(t_{2}-t_{1})A(t_{2})A(t_{1})$. Recall also that $A(t)=\sum_{\omega}A_{\omega}e^{-i\omega t}$ [Eq. (5)], so we can write:

[TABLE]

where

[TABLE]

where $\omega_{\pm}=\frac{1}{2}(\omega_{1}\pm\omega_{2})$ and $\theta_{\pm}=t_{1}\pm t_{2}$. This change of variables rotates the square by $\pi/4$ into a rhombus, and the integration limits correspondingly change to:

[TABLE]

The Jacobian is $1/2$. We can now perform the integral over $\theta_{+}$:

[TABLE]

where we used the identity $\mathcal{C}(\theta)=\mathcal{C}^{*}(-\theta)$ and a change of variables $\theta\mapsto-\theta$ to make all the integration domains positive. When $\omega_{+}=0$ the expression above should be understood as a lim${}_{\omega_{+}\to 0}$. Thus:

[TABLE]

which is Eq. (21b) in the main text.

Recall that $A_{\omega}^{\dagger}=A_{-\omega}$; for $H_{\mathrm{LS}}=\sum_{\omega_{1}\omega_{2}}F_{\omega_{1}\omega_{2}}A_{\omega_{2}}A_{\omega_{1}}$ to be Hermitian, it is thus sufficient if

[TABLE]

This is indeed satisfied by Eq. (296), thus ensuring the Hermiticity.

Appendix C The Davies-Lindblad equation is the infinite coarse-graining time limit of the CGME

The RWA used to derive the Davies-Lindblad equation leaves something to be desired. We simply dropped terms with different Bohr frequencies, without a rigorous mathematical justification. Here we repeat, with some extra clarifications, the derivation given in Ref. [40] which shows that the Davies-Lindblad equation is the infinite coarse-graining timescale limit of the CGME.

First we rewrite Eq. (70) in the form

[TABLE]

which we can do after changing variables $\omega^{\prime}\to-\omega^{\prime}$ as remarked below Eq. (71). Our goal is to show that $\lim_{T_{a}\to\infty}\gamma_{\omega\omega^{\prime}}(T_{a})=\gamma(\omega)\delta_{\omega\omega^{\prime}}$.

Lemma 4**.**

The following equivalent form holds for $\gamma_{\omega\omega^{\prime}}(T_{a})$:

[TABLE]

Proof.

Let us rewrite $\omega^{\prime}s-\omega s^{\prime}$ in terms of a sum and difference of Bohr frequencies:

[TABLE]

where $u=s-s^{\prime}$ and $v=s+s^{\prime}$. After this change of variables $\mathcal{C}(s-s^{\prime})=\mathcal{C}(u)$, and since $s=(v+u)/2$ and $s^{\prime}=(v-u)/2$, the Jacobian of the transformation is $1/2$. In terms of the new variables the integration domain is a rhombus, bounded between the lines $u=v$ and $u=-v$ for $v\in[0,T_{a}]$ and the lines $u=2T_{a}-v$ and $v-2T_{a}$ for $v\in[T_{a},2T_{a}]$. Thus:

[TABLE]

To get the integration limits to be the same we make a change of variables from $v$ to $2T_{a}-v$ in the second double integral:

[TABLE]

which is Eq. (298). ∎

C.1 The $\omega=\omega^{\prime}$ case

For $\omega=\omega^{\prime}$ we now have:

[TABLE]

Let $U=\int_{-v}^{v}du\ e^{i\omega u}\mathcal{C}(u)$. Then $dU=\left(e^{i\omega v}\mathcal{C}(v)+e^{-i\omega v}\mathcal{C}(-v)\right)dv$, and integrating by parts gives:

[TABLE]

Consider the second integral:

[TABLE]

where in the last step we used Eq. (11b). Since $\mathcal{C}(v)=\mathcal{C}^{*}(-v)$, the third integral in Eq. (304) satisfies the same bound and limit. We are thus left with

[TABLE]

with the last equality being due to Eq. (6b).

C.2 The $\omega\neq\omega^{\prime}$ case

For $\omega\neq\omega^{\prime}$ we also perform integration by parts of Eq. (299), but we shall see that this time the boundary terms vanish. We write $\gamma_{\omega\omega^{\prime}}(T_{a})=\frac{1}{T_{a}}e^{i\frac{\omega^{\prime}-\omega}{2}T_{a}}\int_{0}^{T_{a}}dV\ U(v)$, where now $dV=\cos\left(\frac{\omega^{\prime}-\omega}{2}(v-T_{a})\right)dv$ and $U(v)=\int_{-v}^{v}du\ e^{i\frac{\omega+\omega^{\prime}}{2}u}\mathcal{C}(u)$. Then

[TABLE]

Therefore:

[TABLE]

Changing from $v$ to $-v$ in the second term we get

[TABLE]

where we used the angle sum identity for the sine in the last equality. Thus:

[TABLE]

where $S(\omega)$ is the Lamb shift amplitude in Eq. (18b). Since nothing cancels with the overall $T_{a}^{-1}$, we find that the $\omega\neq\omega^{\prime}$ term vanishes.

A similar calculation could be done for the Lamb shift term (69), showing that $\textrm{Im}(x_{\omega\omega^{\prime}})\to\textrm{Im}(x_{\omega\omega})\delta_{\omega\omega^{\prime}}$. Therefore, the Davies-Lindblad (RWA) equation can be understood as the $T_{a}\to\infty$ limit of the CGME.

Appendix D Derivation of the discrete approximation Eq. (78)

The starting point is the time-independent CGME given in Eq. (19). We can obtain a discrete sum just by discretizing $\int d\epsilon$. Below is the estimate of the number of terms $(2k^{*}-1)$ in the sum for that case. The error is:

[TABLE]

where $A_{\epsilon}$ is given in Eq. (20a). The discretization of an integral introduces errors as follows:

[TABLE]

where $f$ is an operator or operator-valued function.

Proof.

By the triangle inequality:

[TABLE]

We apply the triangle inequality to pull the sum out:

[TABLE]

To bound this, we use $O(t)=O(0)+\int_{0}^{t}O^{\prime}(\theta)d\theta$, so that $\int_{-1/2}^{1/2}O(t)dt-O(0)=\int_{-1/2}^{1/2}\int_{0}^{t}O^{\prime}(\theta)d\theta dt$, which implies

[TABLE]

which applied to the interval $[\epsilon,\epsilon+\Delta\epsilon]$ reads:

[TABLE]

Substituting this yields:

[TABLE]

which concludes the proof. ∎

If we use the proven Eq. (315) for

[TABLE]

it will express the error $\mathcal{E}_{k}$ as the sum of:

[TABLE]

Let us start with the $\mathcal{E}_{k2}$ term corresponding to truncating the integral at $k^{*}$. As $\epsilon$ becomes far removed from all the frequencies $\omega$ in the system, this integral becomes small as $1/\epsilon$. To prove this, we integrate by parts:

[TABLE]

Taking the norm and recalling that $\|A\|=1$, we get:

[TABLE]

Note that for a local $n$ qubit Hamiltonian with $O(1)$ local terms and a local operator $A$ with $\|A\|=1$, the commutator $\|[H,A]\|=O(1)$. Combining, we find:

[TABLE]

We enforce that the error is $\sqrt{\tau_{B}/\tau_{SB}}\frac{1}{\tau_{SB}}$, so as to match the other errors in our equation:

[TABLE]

Let us now focus on the term with the derivatives $\mathcal{E}_{k1}$:

[TABLE]

The necessary norm bounds are:

[TABLE]

Using Eqs. (12) and (13) we have the bound

[TABLE]

Substituting $T_{a}=\sqrt{\tau_{B}\tau_{SB}/5}$, we get:

[TABLE]

and the error bound is then:

[TABLE]

This leads to the step-size choice:

[TABLE]

Combining this with Eq. (323), the solution is:

[TABLE]

If $[H,A]=O(1)$, there is no system-size dependence in $\Delta\epsilon,k^{*}$.

We remark that the bounds we used to derive Eq. (330) can be significantly tightened; our goal here was to show that they are $O(1)$. We do not actually recommend to use Eq. (330) for numerical calculations: one can obtain accurate results with much larger $\Delta\epsilon$ and much smaller $k^{*}$.

Appendix E An error bound that is linear in $t$

Lemma 5**.**

Let $\mathcal{L}$ be a linear superoperator and consider the pair of equations

[TABLE]

If any solution $y(t)$ possesses the property that $\forall t,y(0)$ such that $\|y(0)\|_{p}\leq 1:\|y(t)\|_{p}\leq c_{y}$, then for any finite $t$:

[TABLE]

The result holds in the same form for the operator ($p=\infty$, omitted in our notation) and trace ($p=1$) norms.

Proof.

The linear differential equation for $y$ has a unique solution in a neighborhood of any initial condition. A trivial corollary is that the evolution operator is reversible for any finite $t$:

[TABLE]

This is an equality between superoperators. Now we can express the difference as:

[TABLE]

The explicit form of the solution $x(t)$ is:

[TABLE]

which yields, using Eq. (334):

[TABLE]

We now take the norm of both sides:

[TABLE]

By linearity of the time evolution and the norm:

[TABLE]

Note that $\|y\|_{p}=1$. Since we can take $y$ as the initial condition of the equation $\dot{y}(t)=\mathcal{L}y(t)$, and since $e^{\mathcal{L}(t-\tau)}y$ will then be the solution at time $t-\tau$, we have, by assumption of the lemma:

[TABLE]

Note that we had to make these assumptions for any $y$ of bounded norm, not just density operators. It now follows from Eq. (337) that:

[TABLE]

∎

In particular, if $\mathcal{L}$ is a positive trace-preserving map then $c_{y}=1$.

The only error for which we can apply Lemma 5 in this work is the Markov error, which we present here:

[TABLE]

where $c_{y}=c_{BM}$ and we used Eq. (203), which contains $c_{B}$. We instead used a loose bound from Eq. (204) since it results in a more compact expression. Apart from this we have no occasion to use Lemma 5 in this work, since not all the equations appearing in our derivation are invertible. This is why we instead use Lemma 1 or prove the necessary results independently.

Appendix F Proof of Lemma 2

Proof.

The formal solutions allow us to write:

[TABLE]

Substituting $\pi(t)=\rho(t)+\delta\rho(t)$ gives:

[TABLE]

Making the time-averaging explicit:

[TABLE]

We would like to change the integration variables of the first term $\mathcal{L}_{\theta+\tau}\rho(\theta)$ in such a way that it becomes $\mathcal{L}_{\theta^{\prime}}\rho(\theta^{\prime}-\tau^{\prime})$. The change of variables is $\theta+\tau=\theta^{\prime},\leavevmode\nobreak\ \tau=\tau^{\prime}$. In the second term $-\mathcal{L}_{\theta}\rho(\theta)$ the change is trivial $\theta=\theta^{\prime},\leavevmode\nobreak\ \tau=\tau^{\prime}$. This introduces transient effects at the beginning and the end of the evolution, as we want to revert the integration region of the first term $\mathcal{L}_{\theta^{\prime}}\rho(\theta^{\prime}-\tau^{\prime})$ to be the same as the second term. This is illustrated in Fig. 11 and the corresponding integrals are:

[TABLE]

After this change of variables and collecting all the contributions from the $\mathcal{L}_{\theta}\rho(\theta)$ term in Eq. (344) into a single term, we get:

[TABLE]

We call the middle term transient because its contribution does not grow with $t$. One can interpret it as two instantaneous kicks to the solution at time [math] and $t$. For $t\leq T_{a}/2$ the integrals in the second term cancel over the area $(T_{a}/2-t)^{2}$, as indicated by the yellow triangles in Fig. 11. In fact by splitting the integral above differently we could get a slightly tighter (parametrically better for the transient) bound, but working with rectangular regions is easier. Now we can take the norm:

[TABLE]

where $\Lambda$ is a number such that

[TABLE]

The difference between $\rho$’s at nearby time points is bounded trivially as:

[TABLE]

which can be used to simplify:

[TABLE]

We analyzed a similar inequality in the proof of Lemma 1, and as we saw there this amounts to solving for $b_{2}(t)$ that saturates the inequality and upper bounds $\delta\rho(t)$. Note that despite the $\min$ function there is no abrupt slope change on the r.h.s. and the derivative matches, so there is no discontinuity in $\dot{b}_{2}$. The equations for $b_{2}$ such that $\|\delta\rho(t)\|_{1}\leq b_{2}(t)$ are:

[TABLE]

Solving the first equation, we find:

[TABLE]

where the last estimate is in the limit $\Lambda T_{a}\ll 1$. Now we use that as the initial condition for the second equation to find:

[TABLE]

Notice that $\dot{b}_{2}>0$ always. We can express a bound $\|\delta\rho(t)\|_{1}\leq b_{2}(t)$ for all $t$ as follows:

[TABLE]

The simplest observation is that for $t\leq c/\Lambda$:

[TABLE]

∎

Appendix G Derivation of the bound (245)

Here we discuss the derivation of the error introduced by time-averaging and restoring complete positivity, starting from Eq. (239):

[TABLE]

Similarly to Eq. (242), we loosen the bound by using $e^{-\Lambda T_{a}/2}\leq 1$. For $t>T_{a}/2$ and $c_{\Lambda}=4/\Lambda\tau_{SB}\geq 1$, we obtain:

[TABLE]

At $T_{a}=\sqrt{c_{\Lambda}\tau_{B}\tau_{SB}/(3c_{BM}+2)}$ for $t>T_{a}/2$ we have:

[TABLE]

We can also give a bound that works for all times. The simplest way to do so would be to employ the monotonicity of the error:

[TABLE]

and this holds for all $t$ and all $\tau_{B}/\tau_{SB}$. We note that the more careful optimization in App. H allows us to get a tighter bound, but the difference is most drastic for large $\tau_{B}/\tau_{SB}$ when the bound is weaker than the trivial $\|\rho_{BM,I}(t)-{\rho}_{C,I}(t)\|\leq 1+c_{BM}$.

Note that for $t\leq T_{a}/2$ the expression in Eq. (359) grows exponentially with $\sqrt{\tau_{B}/\tau_{SB}}$. But we can have a bound linear in $\sqrt{\tau_{B}/\tau_{SB}}$ if we combine Eq. (359) with the trivial bound $\|\rho_{BM,I}(t)-{\rho}_{C,I}(t)\|\leq 1+c_{BM}$:

[TABLE]

Proof.

Here we prove Eq. (361). Note that the bound is monotonic in $t$, so the lowest value of $\tau_{B}/\tau_{SB}=\epsilon^{*2}$ for which the minimum is $2$ for all $t$ is given by equating the expressions at $t=0$:

[TABLE]

Relaxing this to an inequality, we find:

[TABLE]

where we used $c_{\Lambda}\geq 1$. Thus any bound on Eq. (359) for $\tau_{B}/\tau_{SB}\leq\epsilon^{*2}$ that is monotonic in both $\tau_{B}/\tau_{SB}$ and $t$ will also be a bound for all $\tau_{B}/\tau_{SB}$. We choose to use the following in the exponent:

[TABLE]

Thus we arrive at

[TABLE]

which evaluates to the r.h.s of Eq. (361) using Eq. (363) and $c_{BM}\geq 1\Rightarrow c_{BM}/(3c_{BM}+2)\geq 1/5$. ∎

We will now redo this calculation for $T_{a}=\sqrt{\tau_{B}\tau_{SB}/5}$ and $c_{\Lambda}=1$. Using this in Eq. (357):

[TABLE]

Now we find $\epsilon^{*}$:

[TABLE]

This allows one to bound:

[TABLE]

for all $t$.

Appendix H Optimal time

We begin with Eq. (357), and use it to write down the optimum $t$-dependent coarse-graining time:

[TABLE]

In principle we could consider a family of equations with different $T_{a}(t)$ and solve each one of them just to obtain $\rho(t)$ at one point. This is not a significant overhead. We proceed with the goal of obtaining a single equation since it is conceptually simpler. To this end we now define the bound relevance time $t^{*}$, which is the largest time for which our bound is still better than the trivial bound $\|\rho_{BM,I}(t)-{\rho}_{C,I}(t)\|\leq 1+c_{BM}$:

[TABLE]

We would like to optimize $T_{a}$ at $t=t^{*}$, which would mean that $t^{*}$ is maximized. The corresponding equation for $t^{*}$ is:

[TABLE]

and solving the resulting quadratic equation yields:

[TABLE]

Thus $T_{a}$ is:

[TABLE]

And the error is:

[TABLE]

We should only use these expressions if $t^{*}>T_{a}/2$. For $\tau_{SB}\ll\tau_{B}$ (when $t^{*}$ is not likely to be $>T_{a}/2$), this changes the asymptotic form of $T_{a}$ to $T_{a}\sim\tau_{SB}$. In the regime $\tau_{SB}\gg\tau_{B}$ that we normally address with the CGME, we expect these formulas to be a minor improvement on the loose bounds presented in the main text. Their main advantage is a tighter bound on the error for $t>T_{a}/2$, $t\sim\tau_{SB}$.

Appendix I Proof of the bound (269b) for the Redfield limits of integration error

Starting from Eq. (268), we first send the integration limit to infinity:

[TABLE]

The $\min$ function evaluates to $1-\epsilon_{T}$ when $x\leq x^{*}$ and to $\frac{4\tau_{B}}{x\tau_{SB}}$ when $x>x^{*}$, where

[TABLE]

Therefore:

[TABLE]

where $E_{1}(x)$ denotes the exponential integral function: $E_{1}(x)\equiv\int_{x}^{\infty}t^{-1}e^{-t}dt$. A tight bound in terms of elementary functions is given by [82][Thm. 2]:

[TABLE]

Thus

[TABLE]

where in the second line we used that $\forall x$: $1-e^{-x}<x$.

Alternatively, and using more elementary methods, assuming $x^{*}<1$, we obtain:

[TABLE]

while when $x^{*}\geq 1$ the term with the logarithm is absent:

[TABLE]

Combining these two cases directly yields Eq. (269b), which is an upper bound on the tighter bound presented in Eq. (378b).