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:
\begin{equation}
x(t) - y(t) =\int_0^t\int_0^\tau K_{\tau-\theta}( x(\theta) - y(\theta))d\tau d\theta + \int_0^t\mathcal{E}d\tau \ .
\end{equation}
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:
\begin{equation}
\|x(t) - y(t)\|_1 \leq\int_0^t\int_0^\tau \Lambda_{\tau-\theta}\|x(\theta) - y(\theta)\|_1 d\tau d\theta + \int_0^t\|\mathcal{E}\|_1 d\tau \ .
\end{equation}
The upper bound $b_K(t)$ on $ \|x(t) - y(t)\|_1$ satisfies the equation:
\begin{equation}
b_K(t) = \int_0^t\int_0^\tau \Lambda_{\tau-\theta}b_K(\theta)d\tau d\theta + \|\mathcal{E}\|_1 t\ .
\end{equation}
Taking the derivative once, we obtain:
\begin{equation}
\dot{b}_K(t) = \int_0^t \Lambda_{t -\tau}b_K(\tau)d\tau  +  \|\mathcal{E}\|_1\ .
\end{equation}
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:
\begin{equation}
\dot{b}_K(t)\leq \int_0^t \Lambda_{t -\tau}b_K(t)d\tau  +  \|\mathcal{E}\|_1\ .
\end{equation}
In other words, we can again upper-bound $b_K(t)\leq b_K'(t)$, such that
\begin{equation}
\dot{b}_K'(t)= \int_0^t \Lambda_{t -\tau}d\tau b_K'(t)  +  \|\mathcal{E}\|_1\ .
\end{equation}
Now let $\Lambda$ be a positive constant such that $\int_0^t \Lambda_\tau d\tau \leq \Lambda$ for all $t$.
%Note that the $\Lambda$ given to us in the conditions of the Lemma satisfies:
%    \int_0^t \Lambda_t dt \leq \Lambda
We thus arrive at
\begin{equation}
\dot{b_K}''(t)= \Lambda b_K''(t)  +  \|\mathcal{E}\|_1\ ,
\end{equation}
for a new upper-bound $b_K''(t)\geq b_K'(t)$. This is exactly the same as Eq.~\eqref{eq:107}, and hence the Lemma still holds for non-Markovian equations.
\end{proof}

Note that for our purposes $\Lambda = 4/\tau_{SB}$ is usually taken, thanks to Eq.~(\ref{eq:|L|}, \ref{newLambda}), and that the non-Markovian bound will be generalized in Sec.~\ref{Born}.
The exponential on the r.h.s. of  Eq.~\eqref{eq:108} may be troublesome. Under additional assumptions we can prove a much tighter linear in $t$ bound presented in Appendix \ref{LemmaTwo}, yet we did not find a way to recast all of our results in that tighter form.

\subsection{Born approximation error}
\label{Born}

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.~\eqref{eq:Born21}. 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.~\eqref{eq:16b}]
${\dot\rho}_{\text{tot}}(t) =-i[V(t),\rho_{\text{tot}}(t)]$ [Eq.~\eqref{eq:16a}] 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.
%\DL{This can be addressed using the observations in our BC discussion}
%which can be checked in a future work. \DL{to be finalized} \EM{I've written one way to complete the sentence, but if you want to add anything about why the Dyson series is generally converging, please do}
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\ge 0$. Such a condition was used in Ref.~\cite{ABLZ:12-SI} 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.~\eqref{eq:T1tauB}. 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.~\eqref{eq:rho_Born}.
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.~\eqref{eq:rho_Born} in terms of $\rho_t^B$ is:
\begin{equation}
{\dot\rho}^B_t =\int_0^t K_{t-\tau}^{2,B}(\rho_\tau^{B}) d\tau=\int_0^t \mathcal{C}(\tau-t) (A_t \rho^B_\tau A_\tau - \rho^B_\tau  A_\tau A_t) +\mathcal{C}(t-\tau) (A_\tau \rho^B_\tau  A_t -  A_t A_\tau\rho^B_\tau )d\tau \ .
\label{eq:Born84}
\end{equation}
%where we have introduced the notation for the kernel $K_{t-\tau}^{2,B}$ of the superoperator.
The formal solution is:
\begin{equation}
{\rho}^B_\tau =\rho_0 +\int_0^{\tau}d\theta \int_0^\theta dm \mathcal{C}(m-\theta) (A_\theta \rho^B_m A_m - \rho^B_m  A_m A_\theta) +\mathcal{C}(\theta-m) (A_m \rho^B_m  A_\theta -  A_\theta A_m\rho^B_m ) .
\end{equation}
We substitute this solution back into Eq.~\eqref{eq:Born84}, 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:
\bes
\label{eq:86}
\begin{align}
{\dot\rho}^B_t &=\int_0^t K_{t-\tau}^{2,B}\left(\int_0^\tau \int_0^\theta K_{\theta-m}^{2,B}(\rho_m^{B}) dm d\theta \right) d\tau=\int_0^t K_{t-m}^{4,B}(\rho_m^{B}) d\tau \\
&=\int_0^t \mathcal{C}(\tau-t) (A_t \rho_0 A_\tau - \rho_0  A_\tau A_t) +\mathcal{C}(t-\tau) (A_\tau \rho_0  A_t -  A_t A_\tau\rho_0 )d\tau
\label{firstLine}\\
&+\int_0^td\tau\int_0^{\tau}d\theta \int_0^\theta dm  \\
[ &\mathcal{C}(\tau-t)\mathcal{C}(m-\theta) A_t (A_\theta \rho^B_m A_m - \rho^B_m  A_m A_\theta)  A_\tau + \mathcal{C}(\tau-t)\mathcal{C}(\theta-m) A_t (A_m \rho^B_m  A_\theta -  A_\theta A_m\rho^B_m )  A_\tau \\
-&\mathcal{C}(\tau-t)\mathcal{C}(m-\theta) (A_\theta \rho^B_m A_m - \rho^B_m  A_m A_\theta)   A_\tau A_t - \mathcal{C}(\tau-t)\mathcal{C}(\theta-m) (A_m \rho^B_m  A_\theta -  A_\theta A_m\rho^B_m )  A_\tau A_t \\
+&\mathcal{C}(t-\tau) \mathcal{C}(m-\theta) A_\tau (A_\theta \rho^B_m A_m - \rho^B_m  A_m A_\theta) A_t+\mathcal{C}(t-\tau)\mathcal{C}(\theta-m)  A_\tau (A_m \rho^B_m  A_\theta -  A_\theta A_m\rho^B_m )  A_t \\
-& \mathcal{C}(t-\tau)\mathcal{C}(m-\theta) A_t A_\tau(A_\theta \rho^B_m A_m - \rho^B_m  A_m A_\theta)  - \mathcal{C}(t-\tau)\mathcal{C}(\theta-m) A_t A_\tau(A_m \rho^B_m  A_\theta -  A_\theta A_m\rho^B_m )]\ . \label{lastLine}
\end{align}
\ees

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$
${\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:
\bes
\begin{align}
{\dot\rho}^{B4}_t &=\int_0^t \mathcal{C}(\tau-t) (A_t \rho_0 A_\tau - \rho_0  A_\tau A_t) +\mathcal{C}(t-\tau) (A_\tau \rho_0  A_t -  A_t A_\tau\rho_0 )d\tau \\
&+\Tr_{\text{b}}[A_t\otimes B_t,\int_0^t d\tau [A_\tau\otimes B_\tau, \int_0^\tau d\theta [A_\theta \otimes B_\theta, \int_0^\theta dm [A_m \otimes B_m, \rho^{B4}_m\otimes \rho_{\text{b}}]]]] \ .
\end{align}
\ees
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:
\bes
\begin{align}
\label{eq:B4}
&{\dot\rho}^{B4}_t =\int_0^t K_{t-\tau}^{4,B4}(\rho_\tau^{B4}) d\tau=\int_0^t K_{t-\tau}^{2,B}(\rho_0) d\tau \\
&+\int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm \Tr_{\text{b}} B_t B_\tau B_\theta B_m \rho_{\text{b}}  A_t A_\tau A_\theta A_m  \rho^{B4}_m
- \Tr_{\text{b}} B_t B_\tau B_\theta \rho_{\text{b}}  B_m A_t A_\tau A_\theta   \rho^{B4}_m A_m - \dots
\end{align}
\ees
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~\cite{Isserlis:1916}), 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:
\begin{equation}
\Tr_{\text{b}} B_t B_\tau B_\theta B_m \rho_{\text{b}}  =   \Tr_{\text{b}} B_t B_\tau  \rho_{\text{b}}    \Tr_{\text{b}} B_\theta B_m \rho_{\text{b}}   +  \Tr_{\text{b}} B_t B_\theta  \rho_{\text{b}}    \Tr_{\text{b}} B_\tau B_m \rho_{\text{b}} +\Tr_{\text{b}} B_t B_m \rho_{\text{b}}    \Tr_{\text{b}} B_\tau B_\theta \rho_{\text{b}}\ . \label{full}
\end{equation}
By definition, $\mathcal{C}(t-\tau ) =  \Tr_{\text{b}} B_t B_\tau  \rho_{\text{b}} $. Thus:
\begin{equation}
\Tr_{\text{b}} B_t B_\tau B_\theta B_m \rho_{\text{b}}  =  \mathcal{C}(t-\tau)   \mathcal{C}(\theta-m)
+ \mathcal{C}(t-\theta)   \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)
\end{equation}
We see that in Eq.~\eqref{eq:86} only the first term was present. One can check that the order of $A$ is exactly the same, so upon subtracting Eq.~\eqref{eq:86} from Eq.~\eqref{full} only the two terms $  \mathcal{C}(t-\theta)  \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)$ remain.
%That is not surprising since both expressions came from the same derivation.
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,~  \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,~  \tau\leftrightarrow m$ or $t\leftrightarrow m,~  \tau\leftrightarrow \theta$, then the two arcs are stuck together due to the condition $t>\tau>\theta>m$; see Fig.~\ref{pairIllustrated}.
\begin{figure}[h]
\centering\includegraphics[width=0.5\linewidth]{pairing.png}
\caption{Different possible pairings of $ \Tr_{\text{b}} B_t B_\tau B_\theta B_m \rho_{\text{b}}$ in Wick's theorem. Note that the first and second order are in $\tau_B/\tau_{SB}$, not in the coupling strength itself, in which both diagrams are 4th order.}
\label{pairIllustrated}
\end{figure}

At the moment we have two equations: the Born equation in two forms Eq.~\eqref{eq:Born84} and Eq.~\eqref{eq:86}, as well as the 4th order equation Eq.~\eqref{eq:B4}:
\begin{equation}
{\dot\rho}^{B}_t = \int_0^t K_{t-\tau}^{2,B}(\rho_\tau^{B}) d\tau= \int_0^t K_{t-\tau}^{4,B}(\rho_\tau^{B}) d\tau, \quad {\dot\rho}^{B4}_t = \int_0^t K_{t-\tau}^{4,B4}(\rho_\tau^{B4}) d\tau,
\label{eq:B+B4}
\end{equation}
where $K(\rho)$ are corresponding linear superoperators. $K^{2,B}$ is the 2nd order (in the bath interaction) operator from Eq.~\eqref{eq:Born84}, $K^{4,B}$ and $K^{4,B4}$ have terms up to the 4th order defined in Eq.~\eqref{eq:86} and Eq.~\eqref{eq:B4}. Define
\begin{align}
\mathcal{E}_{B4}(\rho^{\text{any}}) &= \int_0^t (K_{t-\tau}^{4,B4}(\rho_\tau^{\text{any}}) - K_{t-\tau}^{4,B}(\rho_\tau^{\text{any}})) d\tau\ ,
\end{align}
for instance:
\beq
\mathcal{E}_{B4}(\rho^{\text{B4}}) =  {\dot\rho}^{B4}_t - \int_0^t  K_{t-\tau}^{4,B}(\rho_\tau^{B4}) d\tau\ .
\eeq
\begin{align}
\mathcal{E}_{B} &= {\dot\rho}^{\text{true}}_t - \int_0^t  K_{t-\tau}^{2,B}(\rho_\tau^{\text{true}}) d\tau \ . \label{eq:EbDef}
\end{align}
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.~\ref{pairIllustrated}, times the $16$ possible orders of $A_t,A_\tau, A_\theta, A_m$ and $\rho_m$:
\begin{equation}
\mathcal{E}_{B4} =\int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm (\mathcal{C}(t-\theta)   \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)) A_t A_\tau A_\theta A_m  \rho^{\text{any}}_m
- \dots
\label{eq:92}
\end{equation}
We take the trace norm of both sides:

\begin{equation}
\|\mathcal{E}_{B4}(\rho^{\text{any}})\|_1 \leq \int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm \|(\mathcal{C}(t-\theta)   \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)) A_t A_\tau A_\theta A_m  \rho^{\text{any}}_m\|_1
+ \dots
\end{equation}
where we used the triangle inequality
$\|A+B\|_1 \leq \|A\|_1 + \|B\|_1$.
We now use submultiplicativity [Eq.~\eqref{eq:submult}].
Noting that the first norm in the submultiplicativity inequality is the operator norm, we obtain:
\begin{equation}
\|\mathcal{E}_{B4}(\rho^{\text{any}})\|_1 \leq \int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm \|(\mathcal{C}(t-\theta)   \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)) A_t A_\tau A_\theta A_m \| \|\rho^{\text{any}}_m\|_1
+ \dots
\end{equation}
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$:
\begin{equation}
\|\mathcal{E}_{B4}(\rho^{\text{any}})\|_1 \leq c_a\int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm |(\mathcal{C}(t-\theta)   \mathcal{C}(\tau-m) +\mathcal{C}(t-m)   \mathcal{C}(\tau-\theta)) |
+ \dots
\end{equation}
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:
\begin{equation}
\|\mathcal{E}_{B4}(\rho^{\text{any}})\| \leq 16c_{a}\int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm |\mathcal{C}(t-\theta)|   |\mathcal{C}(\tau-m)| +|\mathcal{C}(t-m)|   |\mathcal{C}(\tau-\theta)| \ ,\label{ineq:Eb4}
\end{equation}
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.~\eqref{eq:92}. It turns out that the following condition suffices:
\begin{equation}
\int_{0}^{\infty}|\mathcal{C}(t)| t^n dt \leq \frac{\tau_B^{n}}{\tau_{SB}} \ , \quad \text{for $n=0,1$}\ ,
\label{integrobound}
\end{equation}
which is automatically satisfied for $\tau_{SB}$ and $\tau_B$ as defined in Eq.~\eqref{eq:T1tauB}. We also change the order of integration variables using the condition $0\leq m\leq\theta\leq\tau\leq t$, e.g.:
\begin{equation}
\int_0^t d\tau \int_0^\tau d\theta = \int_0^t d\theta \int_\theta^t d\tau .
\end{equation}
The first term in Eq.~\eqref{eq:92} is bounded as:
\bes
\begin{align}
\int_0^t d\tau \int_0^\tau d\theta \left[\int_0^\theta dm |\mathcal{C}(\tau-m)|\right] |\mathcal{C}(t-\theta)|  &\leq \frac{1}{\tau_{SB}} \int_0^t d\tau \int_0^\tau d\theta|\mathcal{C}(t-\theta)|
= \frac{1}{\tau_{SB}} \int_0^t d\theta \int_\theta^t d\tau|\mathcal{C}(t-\theta)| \\
&= \frac{1}{\tau_{SB}} \int_0^t (t-\theta)d\theta |\mathcal{C}(t-\theta)| \leq \frac{\tau_B}{\tau_{SB}^2} ,
\end{align}
\ees
while the second term in Eq.~\eqref{eq:92} is bounded as:
\bes
\begin{align}
\int_0^t d\tau \int_0^\tau d\theta \int_0^\theta dm |\mathcal{C}(t-m)|   |\mathcal{C}(\tau-\theta)| &= \int_0^t dm \int_m^t d\theta \left[\int_\theta^t d\tau |\mathcal{C}(\tau-\theta)|\right] |\mathcal{C}(t-m)|  \\
&\leq  \frac{1}{\tau_{SB}} \int_0^t dm \int_m^t d\theta|\mathcal{C}(t-m)|  = \frac{1}{\tau_{SB}} \int_0^t (t-m)dm |\mathcal{C}(t-m)| \\
&\leq \frac{\tau_B}{\tau_{SB}^2} \ .
\end{align}
\ees
Together, this yields:
\begin{equation}
\|\mathcal{E}_{B4}(\rho^{\text{any}})\|_1 \leq 32c_a\frac{\tau_B}{\tau_{SB}^2}\ ,
\end{equation}
which is the first hint for the result given in Eq.~\eqref{eq:Born21}.
%%DL: Commenting this out since I don't understand its relevance:
%We note that in Eq.~\eqref{eq:|L|} we had $\|{\dot\rho}^B_t\| \leq 4/\tau_{SB}$, so the proxy $\tau_{SB}\|{\dot\rho}_{4th}(t) - {\dot\rho}_{2nd}(t)\|$ for the relative error is $O(\tau_B/\tau_{SB})$.
The actual error term appearing in Eq.~\eqref{eq:Born20} is $ \mathcal{E}_B$ as defined in Eq.~\eqref{eq:EbDef}. 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:
\begin{equation}
\mathcal{E}_{Bk}(\rho^{\text{any}}) = \int_0^t (K_{t-\tau}^{k,Bk}(\rho_\tau^{\text{any}}) - K_{t-\tau}^{k,B}(\rho_\tau^{\text{any}})) d\tau\ ,
\end{equation}
By employing the diagrammatic technique as in Fig.~\ref{pairIllustrated}, 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:
\begin{equation}
\|\mathcal{E}_{B6}(\rho^{\text{any}})\|_1 \leq 32c_a\frac{\tau_B}{\tau_{SB}^2}\left(1+ \frac{8t}{\tau_{SB}}+ O\left(\frac{\tau_B}{\tau_{SB}}\right)\right) ,
\end{equation}
and the diagrams for the second term are shown in Fig.~\ref{EBkfig}(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.~\ref{EBkfig}(b):
\begin{figure}[t]
\includegraphics[width=0.8\linewidth]{MoreDiagrams.png}
\caption{ Diagrams contributing to the second order of (a) $\mathcal{E}_{B6}$, (b) $\mathcal{E}_{Bk}$ for arbitrary $k$.}
\label{EBkfig}
\end{figure}
\bes
\begin{align}
\|\mathcal{E}_{Bk}(\rho^{\text{any}})\|_1 &\leq 32c_a\frac{\tau_B}{\tau_{SB}^2}\left(\sum_{m=1}^{k/2-1} \frac{m}{(m-1)!}\left(\frac{4t}{\tau_{SB}}\right)^{m-1}+ O\left(\frac{\tau_B}{\tau_{SB}}\right)\right) \\
&\leq 32c_a\frac{\tau_B}{\tau_{SB}^2}\left((1+4 t/\tau_{SB})e^{4 t/\tau_{SB}}+ O\left(\frac{\tau_B}{\tau_{SB}}\right)\right)\ ,
\end{align}
\ees
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$:
\begin{equation}
\|\mathcal{E}_{Bk}(\rho^{\text{any}})\|_1 \leq 32c_a e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}^2}  (1+4 t/\tau_{SB}+O(e^{4 t/\tau_{SB}}\tau_B/\tau_{SB}))
\end{equation}
We do not know how to prove such a statement as of yet.
%\EM{I actually did prove the above with diverging in $k$ coefficient, but it requires even more infinite summations and combinatorics.}
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 \emph{assumptions}:
\begin{enumerate}
\item In the $k\to \infty$ limit, there is a finite radius of convergence in
%\EM{I did not see the convergence in my back of the envelope estimate with loose combinatoric bounds}
\begin{equation}
\|\mathcal{E}_{B\infty}(\rho^{\text{any}})\|_1 \leq 32c_a e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}^2}  (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB}))\ .
\label{eq:128}
\end{equation}
\item The true solution is the limit of the perturbative series:
\begin{equation}
{\dot\rho}^{\text{true}}_t = \int_0^t  K_{t-\tau}^{\infty,B\infty}(\rho_\tau^{\text{true}}) d\tau\ .
\end{equation}
\end{enumerate}

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~\cite{Davies:74}. The slight adjustments needed to include an integral condition such as Eq.~\eqref{integrobound} 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:
\bes
\begin{align}
\mathcal{E}_{B} &= {\dot\rho}^{\text{true}}_t - \int_0^t  K_{t-\tau}^{2,B}(\rho_\tau^{\text{true}}) d\tau  = \int_0^t  (K_{t-\tau}^{\infty,B\infty}(\rho_\tau^{\text{true}}) - K_{t-\tau}^{2,B}(\rho_\tau^{\text{true}}))  d\tau \\
\| \mathcal{E}_{B}\|_1 &\leq\left\|\int_0^t  (K_{t-\tau}^{\infty,B\infty}(\rho_\tau^{\text{true}}) - K_{t-\tau}^{\infty,B}(\rho_\tau^{\text{true}}))  d\tau \right\|_1 + \left\|\int_0^t  (K_{t-\tau}^{\infty,B}(\rho_\tau^{\text{true}}) - K_{t-\tau}^{2,B}(\rho_\tau^{\text{true}}))  d\tau \right\|_1 \\
32\frac{\tau_B}{\tau_{SB}^2}  (1+O(\tau_B/\tau_{SB})) + \left\|\int_0^t  (K_{t-\tau}^{\infty,B}(\rho_\tau^{\text{true}}) - K_{t-\tau}^{2,B}(\rho_\tau^{\text{true}}))  d\tau \right\|_1\ ,
\end{align}
\ees
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.

%Now we generalize the methods developed above to finally bound the Born approximation error, picking up the thread from where we left off at the end of Sec.~\ref{Born}.

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

\begin{equation}
\int_0^t K^{4,B}_{t-\tau}(\rho^{\text{any}}_\tau) d\tau =\int_0^t K^{2,B}_{t-\tau}\left(\rho_0 + \int_0^\tau\int_0^m K^{2,B}_{m-\theta}(\rho^{\text{any}}_m)dm d\theta\right)d\tau\ .
\end{equation}
The solution of
\begin{equation}
\dot{\rho}^{B4} =\int_0^t K^{4,B4}_{t-\tau}(\rho_\tau^{B4}) d\tau= \int_0^t K^{4,B}_{t-\tau}(\rho^{B4}_\tau) d\tau + \mathcal{E}_{B4}
\end{equation}
is equivalent to solving a system of $4$th order equations that are our proxy for the true solution (which is $\infty$ order):
\begin{equation}
\left\{ \begin{array}{ll}
\dot{\rho}^{B4} = \int_0^t K^{2,B}_{t-\tau}(\rho^{(1)}_\tau) d\tau  + \mathcal{E}_{B4},& \rho^{B4}(0) = \rho_0  \\
\dot{\rho}^{(1)} = \int_0^t K^{2,B}_{t-\tau}(\rho^{B4}_\tau) d\tau, & \rho^{(1)}(0) = \rho_0 \ .\end{array} \right.
\end{equation}
We will prove that the solution of the above system remains close to $\rho^{B}$, the solution of the equation after the Born approximation:
\begin{equation}
\dot{\rho}^{B} = \int_0^t K^{2,B}_{t-\tau}(\rho^{B}_\tau) d\tau\ .
\end{equation}
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:
\begin{equation}
\dot{\bf{p}} = \bf{Lp} + \bf{E}, \quad \bf{p} = \left(~^{\rho^{B4}}_{\rho^{(1)}}\right), \quad \bf{L} =  \left(~^{0~ \quad\quad ~ \int K_{t-\tau}d\tau}_{\int K_{t-\tau}d\tau ~ \quad\quad ~0}\right), \quad \bf{E} =  \left(~^{ \mathcal{E}_{B4}(\rho^{B4})}_{0 }\right)\ .
\end{equation}
Define a custom norm that is an infinity norm on the vector of operators, but for each operator the trace norm is taken:
\begin{equation}
\|{\bf{v}}\|_c =\left\| ~^{v^1}_{v^2}\right\|_c = \text{max}(\|v^1\|_1, \|v^2\|_1)
\end{equation}
The triangle inequality is obeyed as required:
\bes
\begin{align}
\|{\bf{v} + \bf{w}}\|_c  &= \text{max}(\|v^1 + w^1\|_1, \|v^2 +w^2\|_1) \leq \text{max}(\|v^1\|_1 + \|w^1\|_1, \|v^2\|_1 +\|w^2\|_1)  \\
&\leq \text{max}(\|v^1\|_1, \|v^2\|_1) +\text{max}( \|w^1\|_1, \|w^2\|_1) = \|{\bf{v}}\|_c + \|{\bf{w}}\|_c\ .
\end{align}
\ees
The maximum of the custom norm of ${\bf{L}}$ over inputs of norm 1 is the same as before:
\begin{equation}
\text{max}_{{\bf{p}}: \|{\bf{p}}\|_c =1}\|{\bf{Lp}}\|_c = \Lambda = \frac{4}{\tau_{SB}} \ ,
\end{equation}
so all the properties needed to prove Lemma~\ref{factDif} hold. We establish:
\begin{equation}
\|\delta{\bf{p}}_t\|_c \leq \frac{e^{\Lambda t}-1}{\Lambda}\|{\bf{E}}\|_c
\end{equation}
in the original variables:
\begin{equation}
\text{max}(\|\delta\rho^{(1)}_t\|_1,  \|\delta\rho^{B4}_t\|_1) \leq \frac{e^{\Lambda t}-1}{\Lambda} \|\mathcal{E}_{B4}  \|_1\ .
\end{equation}
Note that the error we are trying to bound is
\begin{equation}
\| \mathcal{E}_B\|_1 =\left\|\dot{\rho}^{\text{true}} - \int_0^t K^{2,B}_{t-\tau}(\rho^{\text{true}}_\tau) d\tau \right\|_1 = \left\| \int_0^t \left(K^{\infty,B\infty}_{t-\tau}(\rho^{\text{true}}_\tau)-K^{2,B}_{t-\tau}(\rho^{\text{true}}_\tau\right) d\tau \right\|_1\ .
\end{equation}
Our approximation to it is bounded as:
\begin{equation}
\left\|\dot{\rho}^{B4} - \int_0^t K^{2,B}_{t-\tau}(\rho^{B4}_\tau) d\tau \right\|_1 \leq  \left\| \int_0^t K^{2,B}_{t-\tau}(\rho^{(1)}_\tau -\rho^{B4}_\tau  ) d\tau \right\|_1 + \|\mathcal{E}_{B4}\|_1 \leq (2e^{\Lambda t}-1) \|\mathcal{E}_{B4}  \|_1 \ ,
\end{equation}
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$:
\begin{equation}
\left\{ \begin{array}{ll}
\dot{\rho}^{Bk} = \int_0^t K^{2,B}_{t-\tau}(\rho^{(1)}_\tau) d\tau  + \mathcal{E}_{Bk}(\rho^{Bk})& \rho^{Bk}(0) = \rho_0 \\
\dot{\rho}^{(1)} = \int_0^t K^{2,B}_{t-\tau}(\rho^{(2)}_\tau) d\tau & \rho^{(1)}(0) = \rho_0\\ \dots & \dots\\
\dot{\rho}^{(k/2-1)} = \int_0^t K^{2,B}_{t-\tau}(\rho^{B4}_\tau) d\tau & \rho^{(k/2-1)}(0) = \rho_0\end{array} \right.
\end{equation}
For appropriately defined ${\bf{L}}$, $\| \cdot \|_c$ the r.h.s. is still bounded by the same constant:
\begin{equation}
\text{max}_{{\bf{p}}: \|{\bf{p}}\|_c =1}\|{\bf{Lp}}\|_c = \Lambda = \frac{4}{\tau_{SB}}\ ,
\end{equation}
which leads to
\begin{equation}
\text{max}(\|\delta\rho^{(1)}_t\|_1, \dots, \|\delta\rho^{Bk}_t\|_1) \leq \frac{e^{\Lambda t}-1}{\Lambda} \|\mathcal{E}_{Bk}  \|_1\ ,
\end{equation}
and the error in the r.h.s. is bounded without change:
\begin{equation}
\left\|\dot{\rho}^{Bk} - \int_0^t K^{2,B}_{t-\tau}(\rho^{Bk}_\tau) d\tau \right\|_1 \leq  \left\| \int_0^t K^{2,B}_{t-\tau}(\rho^{(1)}_\tau -\rho^{Bk}_\tau  ) d\tau \right\|_1 + \|\mathcal{E}_{Bk}\|_1 \leq (2e^{\Lambda t}-1) \|\mathcal{E}_{Bk}  \|_1 \ .
\end{equation}
Taking the $k\to \infty$ limit, we obtain:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{B}_t\|_1 \leq \frac{e^{\Lambda t}-1}{\Lambda} \|\mathcal{E}_{B\infty} (\rho^{\text{true}}) \|_1 \leq  8(e^{4t/\tau_{SB}}-1) e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}}  (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB})) \ ,
\label{BorErrFirst}
\end{equation}
and, finally:
\bes
\begin{align}
\|\mathcal{E}_B\|_1 &\leq (2e^{\Lambda t}-1) \|\mathcal{E}_{B\infty} (\rho^{\text{true}})\|_1  \\
&\leq 32(2e^{4t/\tau_{SB}}-1) e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}^2}  (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB}))\ ,
\label{eq:193b}
\end{align}
\ees
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.~\eqref{eq:Born21}.

An immediate consequence of Eq.~\eqref{BorErrFirst} is that we can use it to obtain a bound on $c_B$, which we defined after Eq.~\eqref{eq:|L|B} as an upper bound on $\|\rho_{B,I}\|_1$. Since
\begin{equation}
\|\rho^{B}_t\|_1
%= \|\rho^{B}_t-\rho^{\text{true}}_t+\rho^{\text{true}}_t\|_1
\leq \|\rho^{\text{true}}_t\|_1 + \|\rho^{B}_t-\rho^{\text{true}}_t\|_1  = 1 +\|\rho^{\text{true}}_t-\rho^{B}_t\|_1 \ ,
\end{equation}
we obtain:
\begin{equation}
c_B \leq 1 + 8(e^{4t/\tau_{SB}}-1) e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}}  (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB})) \ .
\label{cBbound}
\end{equation}

\subsection{Markov approximation error}
\label{sec:Markov-error}

This subsection follows Ref.~\cite{ALMZ:12} initially. Before the Markov approximation but after the Born approximation, the master equation had the form given in Eq.~\eqref{eq:rho_Born}.
%    {\dot\rho}_{B,I}(t) =\mathcal{L}\rho_{B,I}  = \int_0^t \mathcal{C}(\tau-t) [A(t) \rho_I(\tau) A(\tau) - \rho_I(\tau) A(\tau) A(t)] d\tau +\textrm{h.c.}
%where $\|\mathcal{L}\| \leq 4/\tau_{SB} $. Then $\|{\dot\rho}_{B,I}\| \leq \frac{4}{\tau_{SB}}$.
The error due the Markov approximation $\rho_{B,I}(\tau) \mapsto \rho_{B,I}(t)$, in transitioning to Eq.~\eqref{RedStart}, is:
\begin{equation}
\mathcal{E}_M =   \int_0^t \mathcal{C}(\tau-t) [A(t) (\rho_{B,I}(\tau) - \rho_{B,I}(t)) A(\tau) - (\rho_{B,I}(\tau) - \rho_{B,I}(t)) A(\tau) A(t)] d\tau +\textrm{h.c.}
\end{equation}
Since $\rho_{B,I}(t) - \rho_{B,I}(\tau) = \int_\tau^t {\dot\rho}_{B,I}(t') dt'$, we have
\beq
\|\rho_{B,I}(t) - \rho_{B,I}(\tau)\|_1 \leq \|{\dot\rho}_{B,I}\|_1(t-\tau) \leq\frac{4(t-\tau)c_B}{\tau_{SB}}\  ,
\end{equation}
where we used Eq.~\eqref{eq:|L|B} in the second inequality.
Now we can estimate the norm of the error:
\begin{equation}
\| \mathcal{E}_M\|_1 \le 4\int_0^t |\mathcal{C}(\tau-t)| \frac{4(t-\tau)}{\tau_{SB}}   d\tau  \leq 16 \frac{\tau_B c_B}{\tau_{SB}^2}\ .
\label{eq:E_M}
\end{equation}
We found $c_B$ and saw that it is $O(1)$ in Eq.~\eqref{cBbound}.

Let us now use the Markovian Lemma  \ref{factDif} with Eq.~\eqref{eq:E_M}:
\begin{equation}
\|\rho^{BM}_t-\rho^{B}_t\|_1  \leq 4 \frac{\tau_B  c_{B} }{\tau_{SB}} \left(e^{4t/\tau_{SB}}-1\right)\ . \label{cBMarkov}
\end{equation}
A total expression for the Born and Markov error is [adding Eq.~\eqref{BorErrFirst}]:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq 4 \frac{\tau_B  c_{B} }{\tau_{SB}} \left(e^{4t/\tau_{SB}}-1\right)
+ 8(e^{4t/\tau_{SB}}-1) e^{4t/\tau_{SB}}\frac{\tau_B}{\tau_{SB}}  (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB}))\ ,
\end{equation}
which simplifies to:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq 4(e^{4t/\tau_{SB}}-1)\frac{\tau_B}{\tau_{SB}} \left( c_{B}+ 2 e^{4t/\tau_{SB}} (1+4t/\tau_{SB}+O(e^{4t/\tau_{SB}}\tau_B/\tau_{SB}))\right)\ , \label{BMerrorFirst}
\end{equation}
or substituting in $c_B$:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq \frac{\left(e^{\frac{4t}{\tau_{SB}}} -1\right)\tau_B \left(4 + 8\left(1+\frac{4\left(e^{\frac{4t}{\tau_{SB}}} -1\right)\tau_B}{\tau_{SB}}\right) e^{\frac{4t}{\tau_{SB}}} \left(1+\frac{4t}{\tau_{SB}}+O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right)\right)}{\tau_{SB}}  \ .
\end{equation}
The extra factor from $c_{B}$ can be absorbed into the big-$O$ notation:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq \frac{\left(e^{\frac{4t}{\tau_{SB}}} -1\right)\tau_B \left(4 + 8 e^{\frac{4t}{\tau_{SB}}} \left(1+\frac{4t}{\tau_{SB}}\right)\left(1+O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right)\right)}{\tau_{SB}}  \ .
\end{equation}
A loose bound follows from the above:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq \frac{12\left(e^{\frac{4t}{\tau_{SB}}} -1\right)e^{\frac{8t}{\tau_{SB}}}\tau_B \left(1+O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right)}{\tau_{SB}} \ . \label{BMloose}
\end{equation}
Sometimes we wish to report the bound briefly as:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  = \left(e^{\frac{4t}{\tau_{SB}}} -1\right)O\left(e^{\frac{8t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right) \label{bigOexplain}
\end{equation}
\begin{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:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{BM}_t\|_1  \leq \frac{12\left(e^{\frac{4t}{\tau_{SB}}} -1\right)e^{\frac{8t}{\tau_{SB}}}\tau_B \left(1+Ce^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)}{\tau_{SB}} \leq  \frac{12\left(e^{\frac{4t}{\tau_{SB}}} -1\right)e^{\frac{8t}{\tau_{SB}}}\tau_B \left(1+C\epsilon\right)}{\tau_{SB}}\ . \label{BornObtained}
\end{equation}
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.~\eqref{bigOexplain}.
\end{proof}

\subsection{Time averaging error and its optimization}
\label{sec:t-ave-err2}

\subsubsection{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:
\begin{mylemma}
\label{LemTave}
For any first order linear differential equation $\dot{\rho}(t) = \mathcal{L}_t\rho(t)$ there is a time-averaged version
\begin{equation}
\dot{\pi}(t) = \overline{\mathcal{L}_t}\pi(t)  ,\quad \overline{\mathcal{L}_t} = \frac{1}{T_a}\int_{-T_a/2}^{T_a/2}\mathcal{L}_{t+\tau}d\tau\ .
\end{equation}
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:
\begin{align}
\|\delta\rho(t)\|_1 \leq     [b_2(\min(t,T_a/2)) + c_\rho\Lambda  T_a/4]\max(e^{\Lambda (t-T_a/2)},1) - c_\rho\Lambda  T_a/4\ ,
\label{eq:125}
\end{align}
where $\Lambda,c_\rho$ are constants such that $\forall t,\rho_0, \|\rho_0\|_1\leq1: ~ \|\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:
\beq
b_2(t) \equiv 4c_\rho\frac{t}{T_a}+
c_\rho\frac{4-2\Lambda T_a -(\Lambda T_a/2)^2}{\Lambda T_a} (1 -e^{\Lambda t})   \ .
\eeq
An immediate corollary is:
\beq
\|\delta \rho (t)\|_1 = O(\Lambda T_a)  \text{    for    } t = O(\Lambda^{-1})\ .
\eeq
\end{mylemma}
The proof can be found in Appendix~\ref{LemProof}. The $O(\Lambda T_a)$ corollary enters the estimate for $\|\rho_{BM,I}(t) -\rho_{BMT,I}(t)\|_1$ we presented in Eq.~\eqref{solErr} in the derivation (using $\Lambda = 4/\tau_{SB}$).
We will also need a compact bound on:
\begin{equation}
\frac{b_2(T_a/2)}{c_\rho} =   e^{\Lambda T_a/2}\frac{(\Lambda T_a/2)^2 +4(e^{-\Lambda T_a/2}   -1 +\Lambda T_a/2)}{\Lambda T_a} -\frac{ \Lambda T_a}{4} \ .
\end{equation}
One can prove that the expression in the brackets is upper-bounded by
\begin{equation}
e^{-\Lambda T_a/2}   -1 +\Lambda T_a/2 \leq (\Lambda T_a)^2/8\ .
\end{equation}
This allows us to bound
\begin{equation}
b_2(T_a/2) \leq c_\rho\left(  e^{\Lambda T_a/2} \frac{3\Lambda T_a}{4} -\frac{ \Lambda T_a}{4}\right)\ .
\label{eq:143}
\end{equation}
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 \ref{LemTave}:
\begin{equation}
\|\delta\rho(t)\|_1 \leq   \left\{ \begin{array}{ll} ~ c_\rho\left( e^{\Lambda T_a/2} \frac{3\Lambda T_a}{4} -\frac{ \Lambda T_a}{4} \right), &t\leq T_a/2
\\ ~ c_\rho\left( e^{\Lambda t} \frac{3\Lambda T_a}{4} -\frac{ \Lambda T_a}{4} \right), &t> T_a/2
\end{array} \right. = c_\rho\left(  e^{\Lambda \max (t,T_a/2)} \frac{3\Lambda T_a}{4} -\frac{ \Lambda T_a}{4}\right) \ .
\label{Lambded}
\end{equation}

\subsubsection{Time averaging of the Redfield ME}

We now focus on the specific equations of interest to us.
We already established in Eq.~\eqref{newLambda} 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~\eqref{Lambded} is now bounded as:
\begin{align}
\|\rho_{BM}-\rho_{BMT}\|_1 =\|\delta\rho(t)\|_1\leq c_{BM}\left(  e^{\Lambda \max (t,T_a/2)} \frac{3 T_a}{\tau_{SB}} -\frac{  T_a}{\tau_{SB}}\right) \ ,
\label{TAforAdj}
\end{align}
We will use both the compact bound of Eq.~(\ref{TAforAdj}), and the tightest bound available which is a repetition of Lemma \ref{LemTave}:
\begin{align}
\|\delta\rho(t)\|_1 \leq b_2(t)=     [b_2(\min(t,T_a/2)) + c_{BM}\Lambda  T_a/4]\max(e^{\Lambda (t-T_a/2)},1) - c_{BM}\Lambda  T_a/4\ , \label{repeatTAer}\\  b_2(t) = 4c_{BM}\frac{t}{T_a} +
c_{BM}\frac{4-2\Lambda T_a -(\Lambda T_a/2)^2}{\Lambda T_a} (1 -e^{\Lambda t}), \quad t< T_a/2 \ .
\end{align}

\subsubsection{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.~\eqref{eq:26}]. Before that, Eq.~\eqref{eq:11} was:
\begin{equation}
{\dot\rho}_{BMT}(t) =-i[H,\rho_{BMT}] + \sum_{\omega, \omega'}\frac{1}{T_a}\int_{-T_a/2}^{T_a/2} dt' \int_{-t}^{t'} d\tau'\mathcal{C} (\tau' -t')e^{-i(\omega t'+\omega' \tau')}   (A_{\omega} \rho_{BMT} A_{\omega'} - \rho_{BMT} A_{\omega'} A_{\omega}) +\textrm{h.c.}\ ,
\end{equation}
which changes into Eq.~\eqref{eq:26} (same terms, but the portion of the integral is subtracted):
\bes
\begin{align}
{\dot\rho}_C(t) &=-i[H,\rho_C] + \sum_{\omega, \omega'}
\frac{1}{T_a}\int_{-T_a/2}^{T_a/2} dt' \int_{-t}^{t'} d\tau'\ \mathcal{C} (\tau'-t')e^{-i(\omega t'+\omega' \tau')}(A_{\omega} \rho_C A_{\omega'} - \rho_C A_{\omega'} A_{\omega}) +\textrm{h.c.} - \mathcal{E}_p \\
\mathcal{E}_p &= \sum_{\omega, \omega'}
\frac{1}{T_a}\int_{-T_a/2}^{T_a/2} dt' \int_{-t}^{-T_a/2} d\tau'\ \mathcal{C} (\tau'-t')e^{-i(\omega t'+\omega' \tau')}(A_{\omega} \rho_C A_{\omega'} - \rho_C A_{\omega'} A_{\omega}) +\textrm{h.c.}
\end{align}
\ees
Collecting $\sum_{\omega}$ into $A(t)$ we obtain:
\begin{equation}
\mathcal{E}_p =
\frac{1}{T_a}\int_{-T_a/2}^{T_a/2} dt' \int_{-t}^{-T_a/2} d\tau'\ \mathcal{C} (\tau'-t')(A(t') \rho_C A(\tau') - \rho_C A(\tau') A(t')) +\textrm{h.c.}\ ,
\end{equation}
so that:
\begin{equation}
\|\mathcal{E}_p\|_1 \leq\frac{4}{T_a}\int_{-T_a/2}^{T_a/2} dt' \left|\int_{-t}^{-T_a/2} d\tau'|\mathcal{C}(\tau' -t')|\right| \ .
\end{equation}
\begin{figure}[t]
\includegraphics[width=0.8\linewidth]{regionTa.png}
\caption{Times contributing to $\mathcal{E}_p$ for (a) $t>T_a/2$ (b) $t<T_a/2$}
\label{regionTa}
\end{figure}
\noindent 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.~\ref{regionTa}.
For $t<T_a/2$ we note that the argument of $\mathcal{C}(\tau' - t')$ can be both positive and negative. Thus a factor of $2$ arises when we bound $\int dt'|\mathcal{C}(\tau' - t')| \leq 2/\tau_{SB}$:
\begin{equation}
\|\mathcal{E}_p\|_1 \leq\frac{8}{\tau_{SB}} \frac{T_a/2 -t}{T_a} \leq \frac{4}{\tau_{SB}}\ .
\end{equation}

For $t>T_a/2$, we perform a change of variables: $\tau' -t' = -\tau, ~ t' + T_a/2 = \theta$:
\begin{equation}
\|\mathcal{E}_p\|_1 \leq\frac{4}{T_a}\int_{-T_a/2}^{-T_a/2} dt' \int_{-t}^{-T_a/2} d\tau'|\mathcal{C}(\tau' -t')| = \frac{4}{T_a}\int_{0}^{T_a} d\theta \int_{\theta}^{t + \theta -T_a/2 } d\tau|\mathcal{C}(\tau)|\ .
\end{equation}
We now change the order of integration, noting that $0<\theta<\tau$, thus obtaining the following upper bound:
\begin{equation}
\|\mathcal{E}_p\|_1 \leq \frac{4}{T_a}\int_{0}^{T_a} d\theta \int_{\theta}^{t + \theta -T_a/2 } d\tau|\mathcal{C}(\tau)| \leq\frac{4}{T_a} \int_{0}^{t  +T_a/2 } \tau|\mathcal{C}(\tau)|d\tau\ .
\end{equation}
Recall that the total evolution time  $T$ entered the definition $\frac{\tau_B}{\tau_{SB}} = \int_0^T t |\mathcal{C}(t)| dt$ [Eq.~\eqref{eq:T1tauB-b}]. We can usually set $T=\infty$, except for an Ohmic bath [see Eq.~\eqref{Ohmm}], 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)$:
\begin{equation}
\|\mathcal{E}_p\|_1  \leq\frac{4}{T_a} \int_{0}^{T} \tau|\mathcal{C}(\tau)|d\tau = \frac{4\tau_B}{T_a\tau_{SB}}\ .
\end{equation}
Collecting everything, we obtain:
\begin{equation}
\|\mathcal{E}_p\|_1(t) \leq \frac{4}{\tau_{SB}}\left(\theta(t-T_a/2)\frac{\tau_B}{T_a} +\theta(T_a/2 -t)\right) \ ,
\label{eq:155}
\end{equation}
where $\theta(x) = 1, ~ x>0$ and zero otherwise.

Recall that $\mathcal{E}_p$ is defined in Eq.~\eqref{eq:26} as the extra term appearing if the CGME is written with the same expression as Eq.~\eqref{eq:11} in the r.h.s.
We can now use the same idea as in the proof of Lemma~\ref{factDif} (Sec.~\ref{sec:t-ave-err}). 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
\beq
\dot{b}_p = \Lambda b_p +\|\mathcal{E}_p\|_1(t), \quad b_p(0) =0 \ ,
\label{eq:156}
\eeq
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.~(\ref{newLambda}) as the upper bound on the norm of the superoperator on the r.h.s. of Eq.~\eqref{eq:11}. Indeed, by applying the triangle inequality to the definition of the r.h.s. of Eq.~\eqref{eq:11} we get:
\begin{equation}
\mathcal{L}^{BMT,I}_t(\rho_{\text{test}}) = \frac{1}{T_a}\int_{-T_a/2}^{T_a/2} \mathcal{L}^{BM,I}_{t+\tau}(\rho_{\text{test}})d\tau, \quad \|\mathcal{L}^{BMT,I}_t(\rho_{\text{test}})\|_1 \leq \|\mathcal{L}^{BM,I}_t(\rho_{\text{test}})\|_1 \leq \Lambda \ .
\end{equation}
The solution to Eq.~(\ref{eq:156}) is:
\beq
b_p(t) = e^{\Lambda t} \int_0^t e^{-\Lambda \tau}\|\mathcal{E}_p\|_1(\tau) d\tau \ .
\eeq
Substituting Eq.~\eqref{eq:155}, we obtain two cases:
\begin{align}
b_p(t) &= \frac{4}{\Lambda\tau_{SB}}e^{\Lambda t} \left(1-e^{-\Lambda t } \right),\quad  t< T_a/2\\
b_p(t) &= \frac{4}{\Lambda\tau_{SB}}e^{\Lambda t} \left((1-e^{-\Lambda T_a/2 } )  +\frac{\tau_B}{T_a}(e^{-\Lambda T_a/2 }-e^{-\Lambda t} ) \right),\quad  t> T_a/2
\ ,
\end{align}
which we can summarize as
\beq
b_p(t) \leq \frac{4}{\Lambda\tau_{SB}}[(e^{\Lambda t}  - \max(1,e^{\Lambda (t-T_a/2)})) + \frac{\tau_B}{T_a}\max(e^{\Lambda (t-T_a/2 )}-1,0 ) ] \ . \label{zeroStart}
\eeq
For a concise form we change the first term using  $1-e^{-x}\leq x$:
\beq
b_p(t) \leq e^{\Lambda t}\frac{2T_a}{\tau_{SB}} + \frac{4\tau_B}{\Lambda T_a\tau_{SB}}\max(e^{\Lambda (t-T_a/2 )}-1,0 )  \ . \label{LambdaStay}
\eeq

Let us now focus on the difference with $\rho_{BM,I}(t)$ which is the solution of Eq.~\eqref{RedStart}.
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.~\eqref{TAforAdj} and~\eqref{LambdaStay}, we have:
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq c_{BM}\left(  e^{\Lambda \max (t,T_a/2)} \frac{3 T_a}{\tau_{SB}} -\frac{  T_a}{\tau_{SB}}\right) +  e^{\Lambda t}\frac{2T_a}{\tau_{SB}} + \frac{4\tau_B}{\Lambda T_a\tau_{SB}}\max(e^{\Lambda (t-T_a/2 )}-1,0 )\ . \label{fullWLam}
\end{align}
The tightest bound available that starts at zero at $t=0$ is the combination of Eqs.~\eqref{repeatTAer} and~\eqref{zeroStart}:
\bes
\label{strongestBound}
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq [b_2(\min(t,T_a/2)) + c_{BM}\Lambda  T_a/4]\max(e^{\Lambda (t-T_a/2)},1) - c_{BM}\Lambda  T_a/4+\\  +\frac{4}{\Lambda\tau_{SB}}[(e^{\Lambda t}  - \max(1,e^{\Lambda (t-T_a/2)})) + \frac{\tau_B}{T_a}\max(e^{\Lambda (t-T_a/2 )}-1,0 ) ]\ . \end{align}
\ees
where for $  t< T_a/2$
\begin{equation}
b_2(t) = 4c_{BM}\frac{t}{T_a} +
c_{BM}\frac{4-2\Lambda T_a -(\Lambda T_a/2)^2}{\Lambda T_a} (1 -e^{\Lambda t})
\end{equation}
We may optimize $T_a$ to minimize the bound Eq.~(\ref{fullWLam}), which we proceed to do next.

\subsubsection{Optimization of $T_a$}
\label{meatOptim}
We have used $\Lambda \leq 4/\tau_{SB}$ to obtain the form of the bound in Eq.~(\ref{fullWLam}), 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}, ~ t>T_a/2$ the leading terms of Eq.~(\ref{fullWLam})  are:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq c_{BM}e^{\Lambda t}\frac{3T_a}{\tau_{SB}} + e^{\Lambda t}  \frac{2T_a}{\tau_{SB}}   +c_{\Lambda}e^{\Lambda(t-T_a/2 )}\frac{\tau_B}{T_a} +o(e^{\Lambda t}) \ ,
\label{weakBound}
\end{equation}
We can loosen the bound again via $e^{-\Lambda T_a/2} \leq 1$:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq e^{\Lambda t} \left(\frac{(3c_{BM}+2)T_a}{\tau_{SB}}    +c_{\Lambda}\frac{\tau_B}{T_a}\right) +o(e^{\Lambda t}) \ .
\end{equation}
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.~\eqref{eq:T_a-opt}] for any $\tau_B/\tau_{SB}$ leads to an optimal bound at large enough times. The asymptotic behavior of the bound is:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq 2e^{\Lambda t}\sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}  + o(e^{\Lambda t}) \ . \label{startApp0}
\end{equation}
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.~(\ref{strongestBound}), 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.~(\ref{fullWLam}) can be bounded as
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq 2\sqrt{\frac{c_{\Lambda}(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(e^{\Lambda t+1}-3/5\right)\ ,  \label{looseBound}
\end{equation}
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~\ref{boundingDetails}. This simple form allows us to finally bound $c_{BM}$. Since for a CP coarse-grained equation $\|\rho_C\|_1 =1$:
\begin{equation}
\|\rho_{BM,I}(t)\|_1 \leq 1 +  2\sqrt{\frac{c_{\Lambda}(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(e^{\Lambda t+1}-3/5\right), \quad c_{BM}= 1 + 2\sqrt{\frac{c_{\Lambda}(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(e^{\Lambda t+1}-3/5\right)
\end{equation}
This is a quadratic equation on $c_{BM}$. Its $\geq 1$ solution is:
\begin{equation}
c_{BM} = 1 + \frac12 \left( \sqrt{20 X + 9 X^2} +3 X \right) = 1 + O\left(\sqrt{\frac{\tau_B}{\tau_{SB}}}e^{\Lambda t}\right), \quad X  = 4\frac{c_{\Lambda}\tau_B}{\tau_{SB}}\left(e^{\Lambda t+1}-3/5\right)^2 \label{cBMexact}
\end{equation}
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 C x$. 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 \ref{optAside}). 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
\begin{equation}
1\ll \Lambda t \ll \text{ln}\sqrt{\frac{\tau_{SB}}{\tau_B}}\ .
\end{equation}
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 \ref{numeric}. 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.~\eqref{sqrt5allT}:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq \sqrt{\frac{\tau_B}{5\tau_{SB}}} \left( \left(3e^{\Lambda t+1} -1\right)c_{BM} + \left(7e^{\Lambda t+1} -5\right)\right)\ . \label{simplifySqrt5}
\end{equation}
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.~\eqref{cBMexact}. 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.~\eqref{cBMexact}:
\begin{equation}
c_{BM} \leq 1 + \sqrt{5X} +3X, \quad \sqrt{X} \leq  2\sqrt{\frac{\tau_B}{\tau_{SB}}}e^{\Lambda t+1}\ . \label{cBMX}
\end{equation}
Now we simplify Eq.~\eqref{simplifySqrt5} by relaxing the constant terms:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq \sqrt{\frac{\tau_B}{5\tau_{SB}}}e^{\Lambda t
+1} \left( 3 c_{BM} + 7 \right)  \ .
\end{equation}
With the relaxed form of $c_{BM}$ above this becomes:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq \alpha(t) e (2\sqrt{5} + 6e \alpha(t) +\frac{36e^2}{\sqrt{5}}\alpha^2(t)) ,\quad \alpha(t) =\sqrt{\frac{\tau_B}{\tau_{SB}}}e^{\Lambda t}\ .
\end{equation}
The error near $t=0$ is doubled as compared to the longer formula~\eqref{simplifySqrt5}. 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 $0$. 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}$:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq 13\sqrt{\frac{\tau_B}{\tau_{SB}}}e^{4t/\tau_{SB}}\left(1 + 29 \frac{\tau_B}{\tau_{SB}}e^{8t/\tau_{SB}}\right)\ .
\label{compactErr}
\end{equation}
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.~\eqref{compactErr} and~\eqref{BornObtained}, and using the triangle inequality, we finally obtain:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{C\sqrt{5}}_t\|_1  \leq \frac{\left(e^{\frac{4t}{\tau_{SB}}} -1\right)e^{\frac{8t}{\tau_{SB}}}\tau_B \left(12+  O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right)}{\tau_{SB}}  +\frac{13\sqrt{\tau_B}e^{\frac{4t}{\tau_{SB}}}\left(1 +  \frac{29\tau_B e^{\frac{8t}{\tau_{SB}}}}{\tau_{SB}}\right)}{\sqrt{\tau_{SB}}}\ . \label{BoundObtained}
\end{equation}
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.~\eqref{eq:58}, which uses $T_a= \sqrt{\tau_B \tau_{SB}/5}$.

We will finally show how to simplify Eq.~(\ref{BoundObtained}) to the form $O\left(e^{\frac{6t}{\tau_{SB}}}\sqrt{\frac{\tau_B}{\tau_{SB}}}\right)$ presented in Eq.~\eqref{eq:59b}. Denoting $x =e^{\frac{6t}{\tau_{SB}}}\sqrt{\frac{\tau_B}{\tau_{SB}}} $, we can relax the bound as follows:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{C\sqrt{5}}_t\|_1  \leq x^2 \left(12+  O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right) +13x\left(1 +  29x^2\right)\ . \label{bigOwork}
\end{equation}
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:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{C\sqrt{5}}_t\|_1  \leq x^2 \left(12+  O\left(x^2\right)\right) +13x\left(1 +  29x^2\right)\ .
\end{equation}
From this, the result of Eq.~\eqref{eq:59b} immediately follows:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{C\sqrt{5}}_t\|_1  \leq O(x) = O\left(e^{\frac{6t}{\tau_{SB}}}\sqrt{\frac{\tau_B}{\tau_{SB}}}\right)\ .
\end{equation}

\subsection{The Davies-Lindblad error}
\label{DavSec}

We write the CGME with an explicit dependence on $T_a$:
\begin{equation}
{\dot\rho_C}(t) =-i[H +H_{\mathrm{LS}},\rho_C] + \sum_{\omega, \omega'}\gamma_{\omega\omega'}^{T_a} (A_{\omega} \rho_C A_{\omega'} - \frac{1}{2}\left\{ \rho_C , A_{\omega'} A_{\omega}\right\} )\ .
\end{equation}
The error bound is known for arbitrary $T_a$: see Eq.~\eqref{weakBound}. We now wish to compare the CGME with finite $T_a$ to the limit $T_a \to \infty$ of the same equation:
\begin{equation}
{\dot\rho_C}(t) =-i[H +H_{\mathrm{LS}},\rho_C] + \sum_{\omega, \omega'}\gamma_{\omega\omega'}^{T_a =\infty} (A_{\omega} \rho_C A_{\omega'} - \frac{1}{2}\left\{ \rho_C , A_{\omega'} A_{\omega}\right\} )  + \mathcal{E}_R\ .
\end{equation}
In Appendix \ref{app:D}, we show that this limit is just the Davies-Lindblad equation:
\begin{equation}
{\dot\rho_D}(t) =-i[H +H_{\mathrm{LS}},\rho_D] + \sum_{\omega, \omega'}\gamma_{\omega\omega'}^{T_a =\infty} (A_{\omega} \rho_D A_{\omega'} - \frac{1}{2}\left\{ \rho_D , A_{\omega'} A_{\omega}\right\} )\ .
\end{equation}
We also show in Appendix \ref{app:D} that the difference between the right-hand sides is bounded as:
\begin{equation}
\|\mathcal{E}_R\|_1 \leq 4\sum_{\omega\omega'} |\gamma_{\omega \omega'}^{T_a} - \gamma_{\omega \omega'}^{T_a=\infty}|\leq \frac{4^n}{\tau_{SB}} O\left(\frac{1}{T_a\delta E}\right)\ .
\end{equation}
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:
\begin{equation}
\|\rho_{\text{true}}(t) -\rho_{\text{D}}(t)\|_1 \leq (e^{\Lambda t} -1) \frac{\|\mathcal{E}_R\|_1}{\Lambda} + \|\rho_{\text{true}} -\rho_{\text{C}}\|_1\ .
\end{equation}
In the long time limit the leading terms are:
\begin{equation}
\|\rho_{\text{true}}(t) -\rho_{\text{D}}(t)\|_1 \leq e^{4t/\tau_{SB}} 4^n O\left(\frac{1}{T_a\delta E}\right) + \frac{5T_a}{\tau_{SB}}e^{4t/\tau_{SB}} +e^{4t/\tau_{SB}}\frac{\tau_B}{T_a} +\|\rho_{\text{true}} -\rho_{\text{BM}}\|_1+o(e^{4t/\tau_{SB}}) \ .
\end{equation}
Assuming that $\frac{4^n}{\delta E} \gg \tau_B $, we find the optimal $T_a =2^nO(\sqrt{\tau_{SB}/\delta E})$. Using that together with a bound for $t\geq T_a/2$ we get:
\begin{equation}
\|\rho_{\text{true}}(t) -\rho_{\text{D}}(t)\|_1\leq  2^n O\left( \frac{e^{4t/\tau_{SB}}}{\sqrt{\tau_{SB}\delta E}} \left(1 +\frac{\tau_B}{\tau_{SB}} \right) \right) +\|\rho_{\text{true}} -\rho_{\text{BM}}\|_1 \leq  2^n O\left( e^{12t/\tau_{SB}} \left(\frac{1}{\sqrt{\tau_{SB}\delta E}} +\frac{\tau_B}{\tau_{SB}} \right) \right) \ .
\label{LblBound}
\end{equation}
%but $\sim \tau_B$ error is still there as well.
Assuming a constant system size $n=O(1)$, this is the form presented in Eq.~\eqref{eq:72}, 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.~\ref{sec:times}: $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.~\eqref{eq:1} in the Introduction.

\subsection{Redfield limits of integration error}
\label{errLimits}

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.~\eqref{eq:31c}:
\begin{equation}
\|\mathcal{E}_l\|_1 \leq 4c_{BM}\int_t^{\infty} |\mathcal{C} (t')|dt' \leq \frac{4c_{BM}}{t}\int_t^{T}t' |\mathcal{C} (t')|dt' + 4c_{BM}\int_T^{\infty} |\mathcal{C} (t')|dt'  \leq \frac{4c_{BM}\tau_B}{t\tau_{SB}} + \frac{4c_{BM}\epsilon_T}{\tau_{SB}}\ ,
\label{eq:184}
\end{equation}
where $\epsilon_T= \tau_{SB}\int_T^{\infty} |\mathcal{C} (t')|dt'$ [Eq.~\eqref{eq:eps_T}]. We augment Eq.~\eqref{eq:184} by a trivial bound:
\begin{equation}
\|\mathcal{E}_l\|_1 \leq 4c_{BM}\int_t^{\infty} |\mathcal{C} (t')|dt' \leq \frac{4c_{BM}}{\tau_{SB}} \min (1, \tau_B/t + \epsilon_T)\ .
\end{equation}
Since the bound on $\|\mathcal{E}_l\|_1$ is time dependent, we cannot directly use Lemma~\ref{factDif} for Eq.~\eqref{eq:31a}. Instead, we repeat its proof and integrate the corresponding differential equations as was done in the proof of Lemma~\ref{LemTave}. 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.~\eqref{eq:107}:
\begin{equation}
\dot{b}_R = \frac{4b_R}{\tau_{SB}} + \frac{4c_{BM}}{\tau_{SB}} \min (1, \tau_B/t + \epsilon_T)
\end{equation}
Integrating it, we obtain:
\begin{equation}
\|\rho_{BM}(t) - \rho_R(t)\|_1 { \leq \delta(t)} = c_{BM} e^{4t/\tau_{SB}} \int_0^{4t/\tau_{SB}}\min (1, \frac{4\tau_B}{x\tau_{SB}} + \epsilon_T)e^{-x}dx\ .
\label{eq:182}
\end{equation}
As shown in Appendix~\ref{app:Red-details}, this results in the error bound:
\bes
\begin{align}
\|\rho_{BM}(t) - \rho_R(t)\|_1  &\leq c_{BM}e^{4t/\tau_{SB}}  \left[  \frac{4\tau_B}{\tau_{SB}}\left(1 - \ln\left(1-e^{-\frac{4\tau_B}{\tau_{SB} (1-\epsilon_T)}}\right)\right)+\epsilon_T \right]\\
&\leq c_{BM}e^{4t/\tau_{SB}} \left[\frac{4\tau_B}{{\tau_{SB}}} \left({1 + \frac{1}{e}} +\max\left(\text{ln}\frac{\tau_{SB}(1-\epsilon_T)}{4\tau_B},0\right) \right) + \epsilon_T\right]\ .
\label{eq:183}
\end{align}
\ees
In other words, the $O(\tau_B/\tau_{SB})$ error scaling of the master equation just after Born-Markov approximation, noted in Eq.~\eqref{BMloose} acquires a logarithmic correction.

We will now show how to simplify the total error of Redfield to obtain Eq.\eqref{err:Redfield}. We add the error of $\rho_{BM}$ from Eq.~(\ref{BornObtained}):
\begin{align}
\|\rho^{\text{true}}_t-\rho^{R}_t\|_1  \leq &\frac{\left(e^{\frac{4t}{\tau_{SB}}} -1\right)e^{\frac{8t}{\tau_{SB}}}\tau_B \left(12+  O\left(e^{\frac{4t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\right)}{\tau_{SB}}  +\\&+c_{BM}e^{4t/\tau_{SB}} \left[\frac{4\tau_B}{{\tau_{SB}}} \left({1 + \frac{1}{e}} +\max\left(\text{ln}\frac{\tau_{SB}(1-\epsilon_T)}{4\tau_B},0\right) \right) + \epsilon_T\right]\ .
\end{align}
We will set $\epsilon_T =0$ in this estimate and use the following bound on $c_{BM}$ from Eq.~(\ref{cBMX}):
\begin{equation}
c_{BM} \leq 1 + \sqrt{5X} +3X, \quad \sqrt{X} \leq  2\sqrt{\frac{\tau_B}{\tau_{SB}}}e^{\frac{4t}{\tau_{SB}}+1}\ .
\end{equation}
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.~(\ref{bigOwork}) we arrive at:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{R}_t\|_1  \leq O(x)(1 +\max\left(\text{ln}\frac{\tau_{SB}}{4\tau_B},0\right))\ .
\end{equation}
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:
\begin{equation}
\|\rho^{\text{true}}_t-\rho^{R}_t\|_1  \leq O(x)\text{ln}\frac{\tau_{SB}}{\tau_B} = O\left(e^{\frac{12t}{\tau_{SB}}}\frac{\tau_B}{\tau_{SB}}\right)\text{ln}\frac{\tau_{SB}}{\tau_B}\ .
\end{equation}
This is Eq.~\eqref{err:Redfield}.

\section{Summary and Conclusions}
\label{sec:conc}

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.~\cite{Majenz:2013qw}, 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.~\eqref{tdCoarse}.

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.~\eqref{eq:T1tauB}. These timescales are related to the bath spectral density $\g(\o)$ via the inequalities $2/\tau_{SB} \geq \g(\o)$ and $2\tau_B/\tau_{SB} \geq |\g'(\o)|$. 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.~\eqref{eq:discrete-approx}. 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.~\eqref{eq:Htot}], 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.~\eqref{err:Redfield}-\eqref{eq:59b}] 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~\cite{Albash-Lidar:RMP}. 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:
\begin{itemize}
\item 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.~\ref{ddSec} 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~\cite{PhysRevA.73.052311}, 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.}
\item 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~\cite{Alicki:07a}. The CGME  provides local generators of open system finite temperature dynamics for non-commuting models (see Sec.~\ref{posLoc}). This opens the door to studying the finite-temperature lifetime of, e.g., topological subsystem codes~\cite{Bombin:2010aa} in the thermodynamic limit.
\item Many-body localization effects in quantum annealing --- previously studied only in the closed system setting~\cite{Altshuler2010} --- can now be extended to finite temperature open system quantum annealing.
\end{itemize}

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.

\newpage

\appendix

\section{Proof that $\g(\o) > 0$}
\label{app:A}

\paragraph{Trivial case.} First consider $V = A \otimes B$. The function $\gamma(\omega)$ is the Fourier transform of the correlation function:
\begin{equation}
\gamma(\omega) = \int_{-\infty}^\infty e^{i\omega t}\mathcal{C}(t) dt = \int_{-\infty}^\infty e^{i\omega t}\Tr[\rho_B B(t) B] dt
\end{equation}
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.
\begin{equation}
\gamma(\omega)= \sum_{kl} \int_{-\infty}^\infty e^{i\omega t}\Tr[\sqrt{\rho_B} e^{iH_Bt}\sqrt{\rho_B} B e^{-iH_Bt} B] dt
\end{equation}
Inserting the sum over bath eigenstates:
\begin{equation}
\gamma(\omega)= \sum_{kl} \int_{-\infty}^\infty e^{i(\omega+ E_{kl}) t }\langle k|\sqrt{\rho_B} B |l\rangle \langle l| B\sqrt{\rho_B}|k\rangle dt
\end{equation}
Integration results in delta-functions:
\begin{equation}
\gamma(\omega)= 2\pi\sum_{kl} \delta(\omega -E_{lk})|\langle l| B\sqrt{\rho_B}|k\rangle|^2 \geq 0
\end{equation}
This is enough for all the results in this paper. A more complicated interaction with environment, however, requires a different proof.
\paragraph{General case.}  We consider the case $V = \sum_i A_i\otimes B_i$, so that $\g(\o)$ becomes a matrix. $\rho_B$ will still be assumed stationary: $[\rho_B,H_B] =0$. Repeating the same steps, we arrive at
\begin{equation}
\gamma_{ij}(\omega)= 2\pi\sum_{kl} \delta(\omega -E_{lk})\langle k|\sqrt{\rho_B} B_i |l\rangle \langle l| B_j\sqrt{\rho_B}|k\rangle
\end{equation}
Thus we need to prove that the matrix
\begin{equation}
M_{ij} = \langle k|\sqrt{\rho_B} B_i |l\rangle \langle l| B_j\sqrt{\rho_B}|k\rangle
\end{equation}
is positive. Let $v_i = \langle k|\sqrt{\rho_B} B_i |l\rangle$, then
\begin{equation}
M_{ij} = v_i v_j^*\ .
\end{equation}
For positivity, we need to prove that for any vector  $\vec{w}$ with elements $w_i$ it holds that $\sum_{ij}w_iM_{ij}w_j^* \geq 0$. Indeed:
\begin{equation}
\forall \vec{w}, \quad \sum_{ij}w_iM_{ij}w_j^* = \sum_{ij} w_iv_i v_j^* w_j = \left|\sum_i v_i w_i\right|^2\geq 0\ .
\end{equation}
\begin{equation}
\gamma(\omega) \geq 0
\end{equation}
as a matrix.
\paragraph{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.
%The summation over bath states then has to be replaced by an integration. Placing that integral next to a $\delta$-function should remove it, but how to treat that operation mathematically is not perfectly clear. So our previous proof contains objects that are not perfectly well defined for the continuous spectrum baths.
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~\cite{Breuer:book}, and which avoids the issue mentioned above.

\begin{mylemma}[Bochner]
If a function $\phi(t)$ is of positive type, i.e. $\forall n\in\mathbb{Z}, ~ t_i\in \mathbb{R}, ~\xi_i\in\mathbb{C}$
\begin{equation}
\sum_{i,j=1}^n \phi(t_i-t_j)\xi_i\xi_j^* \geq 0,
\end{equation}
and also $\phi(0)$ is finite and $\phi(t)$ is continuous at $t=0$, then its Fourier transform (if it exists) is non-negative:
\begin{equation}
\int_{-\infty}^\infty\phi(t) e^{i\omega t} dt \geq 0
\end{equation}
\end{mylemma}

We again consider many interaction terms $V = \sum_i A_i\otimes B_i$ (as in Sec.~\ref{nTerms}) and a stationary $\rho_B$ s.t. $[\rho_B,H_B]=0$.
Since $\g(\o)$ is Hermitian we can diagonalize it using a unitary transformation:
\beq
D \equiv U\gamma U^\dagger \Rightarrow D_{\alpha\beta} = \sum_{i,j}U_{\alpha i}\gamma_{ij}U_{\beta j}^* \ .
\eeq
$ 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) = \Tr[\rho_B B_i(s) B_j]$, gives
\beq
D_\alpha = \int_{-\infty}^\infty e^{i\omega s}\pen{\sum_{i,j}U_{\alpha i}\mathcal{C}_{ij}(s)U_{\alpha j}^*}ds\ .
\eeq
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 \emph{positive type} then $ D_\alpha $ must be positive by Bochner's theorem \cite{Reed:1975aa}. Define the following function with $ \cen{t_i} $ an arbitrary time partition:
\beq
f_{mn}^\alpha \equiv \sum_{i,j}U_{\alpha i}\mathcal{C}_{ij}(t_m-t_n)U_{\alpha j}^* \ .
\eeq
Note that \bes
\begin{align}
\mc{B}_{\a\b}(t,t-\tau) &\equiv
%\langle B^\dagger(s) B(s') \rangle = \Tr[\r_B(0) B^\dagger(s) B(s')] = (\Tr[B(s')^\dagger B(s) \r_B(0) ]^\dagger)^* = (\Tr[\r_B(0)B(s')^\dagger B(s)]^\dagger)^* = \langle B^\dagger(s') B(s) \rangle^*\ .
\langle B_{\alpha}(t)B_{\beta}(t-\tau)\rangle_{B}=\Tr\bigl(e^{iH_{B}t}B_{\alpha}e^{-iH_{B}t}e^{iH_{B}(t-\tau)}B_{\beta}e^{-iH_{B}(t-\tau)}\rho_{B}\bigr)\\
\Tr\bigl(e^{iH_{B}\tau}B_{\alpha}e^{-iH_{B}\tau}B_{\beta}\rho_{B}\bigr)=\langle B_{\alpha}(\tau) B_{\beta}\,\rangle_{B} \\
& = \mc{B}_{\a\b}(\tau,0) \equiv \mc{B}_{\a\b}(\tau)\ ,
%= \mc{B}_{\b\a}^*(0,\tau) = \mc{B}_{\b\a}^*(-\tau,0)\ ,
\label{eq:477c}
\end{align}
\ees
Now use the property $ \expv{B_\alpha(s)B_\beta(0)} = \expv{B_\alpha(t)B_\beta(t-s)} $ [Eq.~\eqref{eq:477c}] to write $ f_{mn}^\alpha $ as
\beq
f_{mn}^\alpha = \sum_{i,j}U_{\alpha i}\Tr\ben{\rho_BB_i(t_m)B_j(t_n)}U_{\alpha j}^* = \Tr\pen{\rho_B\sum_iU_{\alpha i}B_i(t_m)\sum_jU_{\alpha j}^*B_j(t_n)} \ .
\eeq
We need to show that $ f^\alpha $ is a positive matrix. For arbitrary $ \ket{v} $ we have
\bes
\begin{align}
\bra{v}f^\alpha\ket{v} &= \sum_{m,n}v_m^*v_nf_{mn}^\alpha = \Tr\ben{\pen{\sum_{i,m}v_m^*U_{\alpha i}\sqrt{\rho_B}B_i(t_m)}\pen{\sum_{j,n}v_nU_{\alpha j}^*B_j(t_n)\sqrt{\rho_B}}} \\
& =
%\Tr\ben{\pen{\sum_{i,m}v_m^*U_{\alpha i}\sqrt{\rho_B}B_i(t_m)}\pen{\sum_{i,m}v_mU_{\alpha i}^*B_i(t_m)\sqrt{\rho_B}}} \\
\Tr\pen{M_\alpha^\dagger M_\alpha} \geq 0 \ ,
\end{align}
\ees
where $M_\alpha \equiv \sum_{i,m}v_m^*U_{\alpha i}\sqrt{\rho_B}B_i(t_m)$.
%	the final inequality follows from the fact that $ M_\alpha^\dagger M_\alpha $ is non-negative which follows immediately from right polar decomposing $ M_\alpha $ (then $ M_\alpha^\dagger M_\alpha = RU^\dagger UR = R^2 \geq 0 $).

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

\section{Lamb shift simplification}
\label{app:LS-simp}
The Lamb shift was given in Eq.~\eqref{eq:HLS1} 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 \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.~\eqref{eq:A(t)}],
so we can write:
\begin{equation}
H_{\mathrm{LS}}=\sum_{\omega_1\omega_2}F_{\omega_1\omega_2}A_{\omega_2}A_{\omega_1}\ ,
\label{eq:HLS2}
\end{equation}
\begin{subequations}
\begin{align}
F_{\omega_1\omega_2}&= \frac{i}{2T_a}\int_{-T_a/2}^{T_a/2}dt_1\int_{-T_a/2}^{T_a/2}  dt_2 \sgn(t_1-t_2)\mathcal{C}(t_2-t_1)  e^{-i(\omega_2t_2 + \omega_1 t_1)}\\
&=\frac{i}{2T_a}\int_{-T_a/2}^{T_a/2}dt_1\int_{-T_a/2}^{T_a/2}  dt_2 \sgn(\theta_-)\mathcal{C}(-\theta_-)  e^{-i(\omega_-\theta_- + \omega_+ \theta_+)}\ ,
\end{align}
\end{subequations}
%Note that many of the terms in the sum in Eq.~\eqref{eq:HLS2} vanish because of the product of two hop operators $A_{\omega_1}A_{\omega_2}$ that force the state in between to be the same. One can use that to speed up the calculation of the sum, but one needs to be wary of accidental degeneracies.
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:
\bes
\begin{align}
&\theta_->0: ~ -T_a+\theta_-\leq   \theta_+\leq T_a -\theta_-\\
&\theta_-<0: ~ -T_a-\theta_-\leq   \theta_+\leq T_a +\theta_- \ .
\end{align}
\ees
The Jacobian is $1/2$.
We can now perform the integral over $\theta_+$:
\bes
\begin{align}
F_{\omega_1\omega_2}  &= \frac{i}{4T_a}\left( -\int_{-T_a}^0 d\theta_- \int_{ -T_a-\theta_-}^{T_a +\theta_- }d\theta_+  +  \int_{0}^{T_a} d\theta_- \int_{-T_a+\theta_-}^{T_a -\theta_-}d\theta_+\right) \mathcal{C}(-\theta_-)  e^{-i(\omega_-\theta_- + \omega_+ \theta_+)}\\
&=\frac{1}{4T_a\o_+}\left( \int_{-T_a}^0 d\theta_- ( e^{-i\omega_+(T_a+\theta_-)} -e^{i\omega_+(T_a+\theta_-)} )  -  \int_{0}^{T_a} d\theta_-( e^{-i\omega_+(T_a-\theta_-)} -e^{i\omega_+(T_a-\theta_-)}) \right) \mathcal{C}(-\theta_-)  e^{-i\omega_-\theta_- }\\
%So it is the same as the Redfield integral, but with smaller bounds.
%Let us rewrite it again:
%&= \frac{1}{4T_a\o_+}\left( \int_{-T_a}^0 d\theta (e^{-iT_a\o_+-i\omega_1\theta } -e^{iT_a\o_++i\omega_2\theta} )  -  \int_{0}^{T_a} d\theta( e^{-iT_a\o_++i\omega_2\theta} -e^{iT_a\o_+-i\omega_1\theta}) \right) \mathcal{C}(-\theta)  \\
&= \frac{1}{4T_a\o_+}\left( \int_{-T_a}^0 d\theta (e^{-iT_a\o_+-i\omega_1\theta } -e^{iT_a\o_++i\omega_2\theta} )\mathcal{C}(-\theta)  -  \int_{0}^{T_a} d\theta( e^{-iT_a\o_++i\omega_2\theta} -e^{iT_a\o_+-i\omega_1\theta}) \mathcal{C}(-\theta)\right)   \\
&= \frac{1}{4T_a\o_+}\left( \int_0^{T_a} d\theta (e^{-iT_a\o_++i\omega_1\theta } -e^{iT_a\o_+-i\omega_2\theta} )\mathcal{C}(\theta)  -  \int_{0}^{T_a} d\theta( e^{-iT_a\o_++i\omega_2\theta} -e^{iT_a\o_+-i\omega_1\theta}) \mathcal{C}^*(\theta)\right)   \\
&= \frac{1}{4T_a\o_+} \int_0^{T_a} d\theta (e^{-iT_a\o_++i\omega_1\theta } -e^{iT_a\o_+-i\omega_2\theta} )\mathcal{C}(\theta)  +  (  e^{-iT_a\o_++i\omega_1\theta}-e^{iT_a\o_+-i\omega_2\theta})^* \mathcal{C}^*(\theta)
%&= \frac{1}{4T_a\o_+}\left(e^{i T_a\o_+}  \int_{-T_a}^0 d\theta e^{i\o_1\theta}  + e^{-i T_a\o_+}  \int_{0}^{T_a} d\theta e^{i\o_1\theta}  - e^{-i T_a\o_+}  \int_{-T_a}^0 d\theta e^{-i\o_2\theta}  - e^{i T_a\o_+}  \int_{0}^{T_a} d\theta e^{-i\o_2\theta}  \right)\mathcal{C}(\theta)\\
%&=  \frac{1}{2T_a\o_+}\left( \Re \int_{0}^{T_a} d\theta e^{i(\o_1\theta-T_a\o_+)}  - \Re \int_{0}^{T_a} d\theta e^{-i(\o_2\theta-T_a\o_+)} \right)\mathcal{C}(\theta)\ ,
\end{align}
\ees
where %in the last equality
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}$.
%This can be further separated into a real and imaginary part:
%     F_{\omega_1\omega_2}  &=\frac{\cos T_a\o_+}{4T_a\o_+}\left( \int_{-T_a}^0 d\theta_- (  e^{-i\omega_1\theta_- } -e^{i\omega_2\theta_-} )  -  \int_{0}^{T_a} d\theta_-( e^{i\omega_2\theta_-} -e^{-i\omega_1\theta_-}) \right) \mathcal{C}(-\theta_-) + \\
%     &\quad -\frac{i\sin T_a\o_+}{4T_a\o_+}\left( \int_{-T_a}^0 d\theta_- (e^{-i\omega_1\theta_- } +e^{i\omega_2\theta_-} )  -  \int_{0}^{T_a} d\theta_-( e^{i\omega_2\theta_-} +e^{-i\omega_1\theta_-}) \right) \mathcal{C}(-\theta_-)  \\
%       &=\left(\frac{\cos T_a\o_+}{4T_a\o_+} \int_{-T_a}^{T_a} d\theta_- (e^{-i\omega_1\theta_- } -e^{i\omega_2\theta_-} )  +\frac{\sin T_a\o_+}{4iT_a\o_+} \int_{-T_a}^{T_a} d\theta_- \sgn(\theta_-) (  e^{-i\omega_1\theta_- } +e^{i\omega_2\theta_-} ) \right) \mathcal{C}(-\theta_-)\ .
%In other words, the result depends on the integrals:
%    R(\omega) \equiv \int_{-T_a}^{T_a} d\theta   e^{-i\omega\theta }\mathcal{C}(\theta), \quad  I(\omega) \equiv -\int_{-T_a}^{T_a} d\theta  e^{-i\omega\theta }\sgn(\theta)\mathcal{C}(\theta)\ ,
%and then the coefficient is:
%    F_{\omega_1\omega_2}  =\frac{1}{4T_a\o_+} \left(\cos T_a\o_+ \left[R(-\omega_1)-R(\omega_2)\right] -i\sin T_a\o_+ \left[I(-\omega_1) +  I(\omega_2)\right]\right)\ .
%    \label{eq:C9}
%Note that $\o_+=0$ when $\o_1=-\o_2$, so the divergence in the first term of Eq.~\eqref{eq:C9} is then cancelled by $R(-\omega_1)-R(\omega_2)=0$. Also note that, using $\mathcal{C}^*(\theta) = \mathcal{C}(-\theta)$, $R$ purely real and $I$ is purely imaginary.
%Eq.~\eqref{eq:C9} can be further simplified by splitting the integrals into their negative and positive integration domains:
%\cos T_a\o_+ R(-\omega_1) -i\sin T_a\o_+ I(-\omega_1) &= e^{-i T_a\o_+}  \int_{-T_a}^0 d\theta e^{i\o_1\theta} \mathcal{C}(\theta) + e^{i T_a\o_+}  \int_{0}^{T_a} d\theta e^{i\o_1\theta} \mathcal{C}(\theta) \\
%& = 2 \Re \int_{0}^{T_a} d\theta e^{i(\o_1\theta+T_a\o_+)} \mathcal{C}(\theta)  \\
%\cos T_a\o_+ R(\omega_2) +i\sin T_a\o_+ I(\omega_2) &= e^{i T_a\o_+}  \int_{-T_a}^0 d\theta e^{-i\o_2\theta} \mathcal{C}(\theta) + e^{-i T_a\o_+}  \int_{0}^{T_a} d\theta e^{-i\o_2\theta} \mathcal{C}(\theta) \\
%& = 2 \Re \int_{0}^{T_a} d\theta e^{-i(\o_2\theta+T_a\o_+)} \mathcal{C}(\theta)\ .
\begin{equation}
F_{\omega_1\omega_2}  =\frac{1}{2T_a\o_+} \Re  \int_{0}^{T_a} d\theta \left( e^{i(\o_1\theta-T_a\o_+)} - e^{-i(\o_2\theta-T_a\o_+)} \right) \mathcal{C}(\theta)\ ,
\label{eq:Fw1w2-copy}
\end{equation}
which is Eq.~\eqref{eq:Fw1w2} in the main text.

Recall that $A_\omega^\dag = 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
\begin{equation}
F_{\omega_1\omega_2} = F_{\omega_1\omega_2}^* = F_{-\omega_2,-\omega_1}\ .
\end{equation}
This is indeed satisfied by Eq.~\eqref{eq:Fw1w2-copy}, thus ensuring the Hermiticity.

\section{The Davies-Lindblad equation is the infinite coarse-graining time limit of the CGME}
\label{app:D}

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.~\cite{Majenz:2013qw} which shows that the Davies-Lindblad equation is the infinite coarse-graining timescale limit of the CGME.

%\subsection{Quick summary}
%For convenience, let us collect the main results of each of the two approaches. The Davies-Lindblad equation is:
%\dot{\tilde{\rho}}(t) = -i \left[ H_{\mathrm{LS}}, \tilde{\rho}(t) \right] + \sum_{\omega} \gamma({\omega}) \left( A_{\omega} \tilde{\rho}(t) A_{\omega}^\dagger - \frac{1}{2} \left\{A_{\omega}^{\dagger} A_{\omega}, \tilde{\rho}(t) \right\} \right)
%\gamma(\omega) = \int_{-\infty}^{\infty} d s e^{i \omega s} C(s,0)\ .
%%H_{\mathrm{LS}} &= \sum_{\omega} Q(\omega) A_{\omega}^{\dagger} A_{\omega} \ , \quad  Q(\omega) = \int_{-\infty}^{\infty} d \omega' \gamma(\omega') \mathcal{P} \left(\frac{1}{\omega - \omega'} \right)\ .
%The CGME is
%\beq \label{eqt:SME_Average}
%\dot{\tilde{\rho}}(t) = -i \left[ H'_{\mathrm{LS}}, \tilde{\rho}(t) \right] + \sum_{\omega, \omega'} \gamma_{\omega \omega'}(T_a) \left( A_{\omega} \tilde{\rho}(t) A_{\omega'}^\dagger - \frac{1}{2} \left\{A_{\omega'}^{\dagger} A_{\omega}, \tilde{\rho}(t) \right\} \right) \ ,
%where the rates $\gamma$ keep a dependence on two different Bohr frequencies $\omega$ and $\o'$:
First we rewrite Eq.~\eqref{eq:18} in the form
\begin{align}
\label{eq:gammaww'}
\gamma_{\omega \omega'} (T_a)  =  \frac{1}{T_a}  \int_0^{T_a} d s \int_0^{T_a} ds'  e^{i ( \omega' s - \omega s' )}  \mathcal{C}(s-s')  \ ,
\end{align}
which we can do after changing variables $\omega' \to -\omega'$ as remarked below Eq.~\eqref{eq:CGME}.
Our goal is to show that $\lim_{T_a\to\infty} \gamma_{\omega \omega'} (T_a) = \gamma(\omega)\delta_{\o\o'}$.

\begin{mylemma}
The following equivalent form holds for $\g_{\omega \omega'} (T_a)$:
\beq
\g_{\omega \omega'} (T_a) = \frac{1}{T_a}e^{i\frac{\omega'-\omega}{2}T_a}\int_0^{T_a}dv\ \cos\left(\frac{\omega'-\omega}{2}(v-T_a)\right)\int_{-v}^{v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u) \ .
\label{eq:g-simplified}
\eeq
\end{mylemma}

\begin{proof}
Let us rewrite $\omega' s - \omega s'$ in terms of a sum and difference of Bohr frequencies:
\beq
\omega' s - \omega s' = \frac{1}{2}(\o'-\o)v+\frac{1}{2}(\o'+\o)u\ ,
\eeq
where $u=s-s'$ and $v=s+s'$. After this change of variables $\mathcal{C}(s-s') = \mathcal{C}(u)$, and since $s=(v+u)/2$ and $s'=(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]$.
\begin{align}
\g_{\omega \omega'}(T_a)= \frac{1}{2T_a}\int_0^{T_a}dv\ e^{i\frac{\omega'-\omega}{2} v}\int_{-v}^{v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u)+\frac{1}{2} \int_{T_a}^{2 T_a}dv\ e^{i\frac{\omega'-\omega}{2} v}\int_{-(2T_a-v)}^{2 T_a-v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u) \ .
\end{align}
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:
\bes
\begin{align}
\g_{\omega \omega'}(T_a)&=\frac{1}{2T_a}\int_0^{T_a}dv\ e^{i\frac{\omega'-\omega}{2} [(v-T_a)+T_a]}\int_{-v}^{v}du\  e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u)+ \frac{1}{2} \int_{0}^{T_a}dv\ e^{-i\frac{\omega'-\omega}{2} [(v-T_a)-T_a]}\int_{-v}^{v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u)\\
&=\frac{1}{T_a}e^{i\frac{\omega'-\omega}{2}T_a}\int_0^{T_a}dv\ \cos\left(\frac{\omega'-\omega}{2}(v-T_a)\right)\int_{-v}^{v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u) \ ,
\end{align}
\ees
which is Eq.~\eqref{eq:gammaww'}.
\end{proof}

\subsection{The $\o=\o'$ case}

For $\omega=\omega'$ we now have:

\begin{equation}
\g_{\omega \omega}(t)=\frac{1}{T_a}\int_0^{T_a} dv\ \int_{-v}^{v}du\ e^{i\omega u}\mathcal{C}(u) \ .
\end{equation}
Let $U = \int_{-v}^{v}du\ e^{i\omega u}\mathcal{C}(u)$.
%Recall the Leibnitz rule for differentiating a definite integral:
%\partial_z \int_{a(z)}^{b(z)} f(x,z)dx = \int_{a(z)}^{b(z)}\partial_z f(x,z)dx + f(b(z),z)b' - f(a(z),z)a'\ .
Then $dU = \left(e^{i\o v} \mathcal{C}(v)+e^{-i\o v} \mathcal{C}(-v)\right)dv$, and integrating by parts
%($\int_0^{T_a} Udv = \left[Uv\right]_0^{T_a} - \int_0^{T_a} vdU$)
\begin{align}
\g_{\omega\omega}(T_a)&=\frac{1}{T_a}\left[v \int_{-v}^{v}du\ e^{i\omega u}\mathcal{C}(u)\right]_0^{T_a}-\frac{1}{T_a}\int_0^{T_a}dv\ v\left( e^{i\omega v}\mathcal{C}(v)+e^{-i\omega v}\mathcal{C}(-v)\right)\ .
% &=\  T_a\int_{-T_a}^{T_a}du\left(1-\frac{u}{2 T_a}\right)e^{i\omega u}C(u) \ .
\label{eq:539}
\end{align}
Consider the second integral:
\bes
\begin{align}
\left| \frac{1}{T_a} \int_0^{T_a}dv\ v e^{i\omega v}\mathcal{C}(v) \right| &\leq \frac{1}{T_a} \int_0^{T_a}dv\ v \left| \mathcal{C}(v) \right| \leq \frac{1}{T_a}\int_0^{\infty}dv\ v \left| \mathcal{C}(v) \right| \\
&= \frac{\tau_B}{\tau_{SB} T_a}  \stackrel{T_a\to\infty}{\longrightarrow} 0\ ,
\end{align}
\ees
where in the last step we used Eq.~\eqref{eq:T1tauB-b}. Since $\mathcal{C}(v) =\mathcal{C}^*(-v)$, the third integral in Eq.~\eqref{eq:539} satisfies the same bound and limit. We are thus left with
\beq
\lim_{T_a \to \infty} \g_{\o\o}(T_a) =\int_{-\infty}^{\infty}du\ e^{i\omega u}\mathcal{C}(u)=\gamma(\omega) \  ,
\eeq
with the last equality being due to Eq.~\eqref{eq:gamma31}.

\subsection{The $\o\neq\o'$ case}
For $\omega\neq\omega'$ we also perform integration by parts of Eq.~\eqref{eq:g-simplified}, but we shall see that this time the boundary terms vanish.
We write
$\g_{\omega \omega'} (T_a) = \frac{1}{T_a}e^{i\frac{\omega'-\omega}{2}T_a} \int_0^{T_a}dV\  U(v)$,
where now $dV = \cos\left(\frac{\omega'-\omega}{2}(v-T_a)\right)dv$ and $U(v) = \int_{-v}^{v}du\ e^{i\frac{\omega+\omega'}{2}u}\mathcal{C}(u)$. Then
\bes
\begin{align}
V(v) &= \frac{2}{\omega'-\omega}\sin\left(\frac{\omega'-\omega}{2}(v-T_a)\right) \\
dU/dv &= e^{i\frac{\omega+\omega'}{2}v} \mathcal{C}(v)+e^{-i\frac{\omega+\omega'}{2}v} \mathcal{C}(-v)\\
\left[U(v) V(v) \right]_0^{T_a} &= U(T_a) V(T_a) - U(0)V(0) = 0 \ .
\label{eq:695c}
\end{align}
\ees
\begin{equation}
\g_{\omega \omega'}(T_a)=-\int_{0}^{T_a}V dU = -\frac{2 e^{i\frac{\omega'-\omega}{2}T_a}}{(\omega'-\omega)T_a}\int_0^{T_a}dv \sin\left( \frac{(\omega'-\omega)}{2}(v-T_a)\right)\left[e^{i\frac{\omega+\omega'}{2}v}\mathcal{C}(v)+e^{-i\frac{\omega+\omega'}{2}v}\mathcal{C}(-v)\right] \ .
\end{equation}
Changing from $v$ to $-v$ in the second term we get
\bes
\begin{align}
\g_{\omega \omega'}(T_a)=&-\frac{2 e^{i\frac{\omega'-\omega}{2}T_a}}{(\omega'-\omega)T_a}\Bigg[\int_0^{T_a}dv\sin\left(\frac{(\omega'-\omega)}{2}(v-T_a)\right)e^{i\frac{\omega+\omega'}{2}v}\mathcal{C}(v)+\int_{-T_a}^{0}dv\sin\left(\frac{(\omega'-\omega)}{2} (-v-T_a)\right)e^{i\frac{\omega+\omega'}{2}v}\mathcal{C}(v)\Bigg]\\
=&\frac{e^{i\frac{\omega'-\omega}{2}T_a}}{(\omega'-\omega)T_a}\int_{-T_a}^{T_a}dv\left[\sin\left(\frac{\omega'-\omega}{2}T_a\right)\left(e^{i \omega v}+e^{i \omega' v}\right)+\frac{\mathrm{sgn(v)}}{i}\cos\left(\frac{\omega'-\omega}{2}T_a\right)\left(e^{i \omega v}-e^{i \omega' v}\right)\right]\mathcal{C}(v) \ ,
\end{align}
\ees
where we used the angle sum identity for the sine in the last equality. Thus:
\beq
\lim_{T_a \to \infty} \g_{\omega \omega'}(T_a) = \lim_{T_a \to \infty} \frac{e^{i\frac{\omega'-\omega}{2}T_a}}{(\omega'-\omega)T_a}\left[\sin\left(\frac{\omega'-\omega}{2}T_a\right)(\gamma(\omega)+\gamma(\omega'))+2\cos\left(\frac{\omega'-\omega}{2}T_a\right)(S(\omega)-S(\omega'))\right] \ .
\eeq
where $S(\omega)$ is the Lamb shift amplitude in Eq.~\eqref{eq:LS-final}.  Since nothing cancels with the overall $T_a^{-1}$, we find that the $\omega \neq \omega'$ term vanishes.

A similar calculation could be done for the Lamb shift term~\eqref{eq:HLS1}, showing that $\Im(x_{\o\o'}) \to \Im(x_{\o\o})\delta_{\o\o'}$.  Therefore, the Davies-Lindblad (RWA) equation can be understood as the $T_a\to\infty$ limit of the CGME.

\section{Derivation of the discrete approximation Eq.~\eqref{eq:discrete-approx}}
\label{app:C}

The starting point is the time-independent CGME given in Eq.~\eqref{eq:Lind}.
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:
\begin{equation}
\mathcal{E}_k = \int_{-\infty}^\infty d\epsilon (A_{\epsilon} \rho_C A_{\epsilon}^\dag - \frac{1}{2}\left\{ \rho_C , A_{\epsilon}^\dag A_{\epsilon}\right\} )-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*} (A_{\epsilon} \rho_C A_{\epsilon}^\dag - \frac{1}{2}\left\{ \rho_C , A_{\epsilon}^\dag A_{\epsilon}\right\} )\ ,
\end{equation}
where $A_\epsilon$ is given in Eq.~\eqref{Atime}.
%    A_\epsilon =  \sqrt{\frac{\gamma(\epsilon)}{2\pi T_a}}\int_{-T_a/2}^{T_a/2}e^{i\epsilon t}A(t)dt
The discretization of an integral introduces errors as follows:
%    \left|\int d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ k<k^*}f(\epsilon) \right| \leq k^* \max|f'|\Delta \epsilon^2 +\left(\int_{-\infty}^{-\Delta \epsilon k^* } + \int_{\Delta \epsilon k^*}^\infty\right) d\epsilon |f(\epsilon)|
%Enhancing this to operators, we get:
\begin{equation}
\left\|\int d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*}f(\epsilon) \right\|_1 \leq \frac{k^*-0.5}{2} \max\|f'\|_1\Delta \epsilon^2 +\left(\int_{-\infty}^{-\Delta \epsilon (k^*-0.5) } + \int_{\Delta \epsilon (k^*-0.5)}^\infty\right) d\epsilon \|f(\epsilon)\|_1 \ ,
\label{eq:C2}
\end{equation}
where $f$ is an operator or operator-valued function.
\begin{proof}
By the triangle inequality:
\bes
\begin{align}
\left\|\int d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*}f(\epsilon) \right\|_1 \leq \left\|\int_{-\Delta \epsilon (k^*-0.5) }^{\Delta \epsilon (k^*-0.5) } d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*}f(\epsilon) \right\|_1
\\+\left(\int_{-\infty}^{-\Delta \epsilon (k^*-0.5) } + \int_{\Delta \epsilon (k^*-0.5)}^\infty\right) d\epsilon \|f(\epsilon)\|_1 \ .
\end{align}
\ees
We apply the triangle inequality to pull the sum out:
\begin{equation}
\left\|\int_{-\Delta \epsilon (k^*-0.5) }^{\Delta \epsilon(k^*-0.5) } d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*}f(\epsilon) \right\|_1 \leq \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*} \left\|\int_{\epsilon-\Delta \epsilon/2 }^{\epsilon +\Delta \epsilon/2 } d\mu f(\mu)-\Delta \epsilon f(\epsilon) \right\|_1\ .
\end{equation}
To bound this, we use $O(t) = O(0) + \int_0^t O'(\theta)d\theta$, so that $\int_{-1/2}^{1/2} O(t)dt - O(0)  = \int_{-1/2}^{1/2} \int_0^t O'(\theta)d\theta dt$, which implies
\begin{equation}
\left\|\int_{-1/2}^{1/2} O(t)dt - O(0)\right\|_1\leq \frac{1}{4}\max_{\theta\in [0,1]}\|O'\|_1\ ,
\label{eq:C3}
\end{equation}
which applied to the interval $[\epsilon, \epsilon +\Delta \epsilon]$ reads:
\begin{equation}
\left\|\int_{\epsilon-\Delta \epsilon/2 }^{\epsilon +\Delta \epsilon/2 } d\mu f(\mu) - f(\epsilon)\right\|_1\leq \frac{\Delta \epsilon^2}{4}\max_{\mu\in [\epsilon, \epsilon +\Delta \epsilon]}\|f'\|_1\ .
\end{equation}
Substituting this yields:
\begin{equation}
\left\|\int_{-\Delta \epsilon (k^*-0.5) }^{\Delta \epsilon (k^*-0.5) } d\epsilon f(\epsilon)-\Delta \epsilon \sum_{k, \epsilon = \Delta \epsilon k, ~ |k|<k^*}f(\epsilon) \right\|_1 \leq (2k^* -1) \frac{\Delta \epsilon^2}{4}\max_{\mu\in [-\Delta \epsilon k^*, \Delta \epsilon k^*]}\|f'\|_1\ ,
\end{equation}
which concludes the proof.
\end{proof}
If we use the proven Eq.~\eqref{eq:C3} for
\begin{equation}
f_\epsilon = (A_{\epsilon} \rho_C A_{\epsilon}^\dag - \frac{1}{2}\left\{ \rho_C , A_{\epsilon}^\dag A_{\epsilon}\right\} )\ , \label{fDef}
\end{equation}
it will express the error $\mathcal{E}_k$ as the sum of:
\begin{equation}
\mathcal{E}_k =\mathcal{E}_{k1} +\mathcal{E}_{k2} ,\quad \mathcal{E}_{k1} = \frac{k^*-0.5}{2} \max\|f'\|_1\Delta \epsilon^2, \quad  \mathcal{E}_{k1} =\left(\int_{-\infty}^{-\Delta \epsilon (k^*-0.5) } + \int_{\Delta \epsilon (k^*-0.5)}^\infty\right) d\epsilon \|f(\epsilon)\|_1\ .
\end{equation}
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:
\begin{equation}
A_\epsilon =  \sqrt{\frac{\gamma(\epsilon)}{2\pi T_a}}\int_{-T_a/2}^{T_a/2}A(t) e^{i\epsilon t} dt = \sqrt{\frac{\gamma(\epsilon)}{2\pi T_a}}\frac{1}{i\epsilon} \left(A(T_a/2) e^{i\epsilon T_a/2} - A(-T_a/2) e^{-i\epsilon T_a/2} +i\int_{-T_a/2}^{T_a/2}[H,A(t)] e^{i\epsilon t}dt  \right)\ .
\end{equation}
Taking the norm and recalling that $\| A\|=1$, we get:
\begin{equation}
\|A_\epsilon\| \leq \sqrt{\frac{\gamma(\epsilon)}{2\pi T_a}} \frac{2 + \|[H,A]\| T_a}{|\epsilon|} \leq \sqrt{\frac{1}{\pi \tau_{SB} T_a}}\frac{2 + \|[H,A]\|T_a}{|\epsilon|}\ .
\end{equation}
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:
\bes
\begin{align}
\|\mathcal{E}_{k2}\|_1 &\leq \left(\int^{-\Delta\epsilon (k^*-0.5)}_{-\infty}+\int_{\Delta\epsilon (k^*-0.5)}^{\infty}\right) \|A_\epsilon\|^2 d\epsilon \leq 2 \int_{\Delta\epsilon (k^*-0.5)}^{\infty}\frac{1}{\pi T_a\tau_{SB}}\frac{(2 + \|[H,A]\| T_a)^2}{|\epsilon|^2}d\epsilon \\
&= 2 \frac{1}{\pi T_a\tau_{SB}}\frac{(2 + \|[H,A]\| T_a)^2}{\Delta\epsilon (k^*-0.5)}\ .
\end{align}
\ees
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:
\begin{equation}
\sqrt{\frac{\tau_B}{\tau_{SB}}} \frac{1}{\tau_{SB}} =   2 \frac{1}{\pi T_a\tau_{SB}}\frac{(2 + \|[H,A]\| T_a)^2}{\Delta\epsilon (k^*-0.5)} \quad \Rightarrow\quad  \Delta\epsilon (k^*-0.5) = 2 \frac{1}{\pi T_a}\sqrt{\frac{\tau_{SB}}{\tau_B}}(2 + \|[H,A]\| T_a)^2\ .
\label{eq:C8}
\end{equation}

Let us now focus on the term with the derivatives $\mathcal{E}_{k1}$:
\begin{equation}
\|\mathcal{E}_{k1}\|_1 \leq 2(k^*-0.5) \max\|A_\epsilon'\|\|A_\epsilon\| \Delta \epsilon^2\ .
\end{equation}
The necessary norm bounds are:
\begin{equation}
\|A_\epsilon\| \leq \sqrt{\frac{\gamma(\epsilon) T_a}{2\pi}} ,\quad \|A_\epsilon' \| \leq\frac{1}{2}\sqrt{\frac{|\gamma'(\epsilon)| T_a}{2\pi\gamma(\epsilon)}} + \frac{T_a}{4}\sqrt{\frac{\gamma(\epsilon) T_a}{2\pi}}\ .
\end{equation}
Using Eqs.~\eqref{eq:22new} and~\eqref{eq:23} we have the bound
\begin{equation}
\|A_\epsilon\| \|A_\epsilon' \| \leq\frac{1}{2}\frac{\sqrt{|\gamma'(\epsilon)|} T_a}{2\pi} + \frac{T_a}{4}\frac{\gamma(\epsilon) T_a}{2\pi}  \leq \frac{1}{2}\frac{\sqrt{2\tau_B/\tau_{SB}} T_a}{2\pi} + \frac{1}{4}\frac{ T_a^2}{2\pi \tau_{SB}}\ .
\end{equation}
Substituting $T_a = \sqrt{\tau_B \tau_{SB}/5}$, we get:
\begin{equation}
\|A_\epsilon\| \|A_\epsilon' \| \leq\frac{\sqrt{2/5} \tau_B}{4\pi} + \frac{ \tau_B}{40\pi} = \frac{(10\sqrt{2/5}+1) \tau_B}{40\pi }\ ,
\end{equation}
and the error bound is then:
\begin{equation}
\|\mathcal{E}_{k1}\|_1 \leq 2(k^*-0.5) \frac{(10\sqrt{2/5}+1) \tau_B}{10\pi }\Delta \epsilon^2\ .
\end{equation}
This leads to the step-size choice:
\begin{equation}
2(k^*-0.5) \frac{(10\sqrt{2/5}+1) \tau_B}{10\pi }\Delta \epsilon^2 = \frac{1}{\tau_{SB}}\sqrt{\frac{\tau_B}{\tau_{SB}}}, \quad \Delta\epsilon = \frac{1}{\tau_{SB}\sqrt{\tau_B \tau_{SB}} (k^*-0.5)\Delta \epsilon}\frac{5\pi }{(10\sqrt{2/5}+1) }\ .
\end{equation}
Combining this with Eq.~\eqref{eq:C8}, the solution is:
\begin{equation}
\Delta\epsilon =\frac{1}{\tau_{SB}}\sqrt{\frac{\tau_B}{\tau_{SB}}}\frac{1}{(2 + \|[H,A]\| T_a)^2}\frac{\sqrt{5}\pi^2 }{2(10\sqrt{2/5}+1) }, \quad k^* = \text{ceil}\left( \sqrt{\frac{\tau_{SB}}{\tau_B}}^3(2 + \|[H,A]\| T_a)^4\frac{4(10\sqrt{2/5}+1) }{\pi^3 } +0.5\right) \ .
\label{eq:C15}
\end{equation}
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.~\eqref{eq:C15} can be significantly tightened; our goal here was to show that they are $O(1)$. We do not actually recommend to use Eq.~\eqref{eq:C15}  for numerical calculations: one can obtain accurate results with much larger $\Delta \epsilon$ and much smaller $k^*$.

%\section{The induced superoperator norm}
%Consider a superoperator $X$ mapping between matrices $A\in \mathcal{M}_n$, the space of $n\times n$ matrices. We define the induced superoperator norm~\cite{Perez-Garcia:2006aa}:
%\|X\|_{p,q} = \sup_{A\in \mathcal{M}_n}\frac{\|X(A)\|_p}{\|A\|_q} ,
%where the $p$-norm of $A$ is $\|A\|_p = (\Tr|A|^p)^{1/p} = (\sum_i \lambda_i^p)^{1/p}$, and $\lambda_i$ are the singular values of $A$ (eigenvalues of $|A|= \sqrt{A^\dagger A}$). The case $p=1$ is the trace-norm and $p=\infty$ is the operator norm. Clearly,
%\| X(A)\|_p  \leq  \| X \|_{q,p} \|A\|_q \ .
%For our purposes $A$ is usually a density matrix or a difference of density matrices, $n$ is the Hilbert space dimension, and we are interested in the case $p=q=1$, where we bound the trace-norm $\|\cdot\|_1$, or $p=q=\infty$, where we bound the operator-norm $\|\cdot\|$, or $p=\infty$ and $q=1$, when $A$ is a density matrix and hence $\| A\|_1=1$.
%In them main text we dropped the $p,q$ subscripts, but they are always implied by context when we take the norm of a superoperator.

\section{An error bound that is linear in $t$}
\label{LemmaTwo}
\begin{mylemma}
\label{LemLinear}
%Consider the positive trace-preserving (PTP) maps $\Lambda_i(t) = e^{\mathcal{L}_i t}$, generated by the time-independent superoperators $\mathcal{L}_1$ and $\mathcal{L}_2 = \mathcal{L}_1 + \Delta\mathcal{L}$ (we do not need complete positivity). The associated states are $\rho_i(t) = \Lambda_i(t) \rho(0)$.
%Let $n$ be the Hilbert space dimension.
%If $\Lambda_1(t)$ is invertible then
%\|\rho_2(t) - \rho_1(t)\|_1 \leq t \|\Delta\mathcal{L} \| \ .
Let $\mathcal{L}$ be a linear superoperator and consider the pair of equations
\begin{equation}
\dot{x}(t) = \mathcal{L}x(t) +\mathcal{E}, \quad  \dot{y}(t) = \mathcal{L}y(t)\ .
\end{equation}
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$:
\begin{equation}
\|x(t)-y(t)\|_p \leq c_y t\|\mathcal{E}\|_p\ .
\end{equation}
The result holds in the same form for the operator ($p=\infty$, omitted in our notation) and trace ($p=1$) norms.
\end{mylemma}

\begin{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$:
\begin{equation}
e^{\mathcal{L}t}e^{-\mathcal{L}t} =\openone \ .
\end{equation}
This is an equality between superoperators. Now we can express the difference as:
\begin{equation}
x(t) -y(t) =e^{\mathcal{L}t}(e^{-\mathcal{L}t}x(t) - y(0)) \ .
\label{eq:119}
\end{equation}
The explicit form of the solution $x(t)$ is:
\begin{equation}
x(t) = e^{\mathcal{L}t}x(0) + e^{\mathcal{L}t}\int_0^t e^{-\mathcal{L}\tau}\mathcal{E} d\tau\ ,
\end{equation}
which yields, using Eq.~\eqref{eq:119}:
\begin{equation}
x(t) - y(t) =e^{\mathcal{L}t} \int_0^t e^{-\mathcal{L}\tau}\mathcal{E} d\tau = \int_0^t e^{\mathcal{L}(t- \tau)}\mathcal{E} d\tau  \ .
\end{equation}
We now take the norm of both sides:
\begin{equation}
\|x(t) - y(t)\|_p \leq  t \|e^{\mathcal{L}(t- \tau)}\mathcal{E} \|_p \ .
\label{eq:122}
\end{equation}
By linearity of the time evolution and the norm:
\begin{equation}
\|e^{\mathcal{L}(t- \tau)}\mathcal{E} \|_p  = \|\mathcal{E}\|_p \|e^{\mathcal{L}(t- \tau)}y \|_p, \quad y = \frac{\mathcal{E}}{\|\mathcal{E}\|_p}
\end{equation}
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:
\begin{equation}
\|e^{\mathcal{L}(t- \tau)}y \|_p \leq c_y\ .
\end{equation}
Note that we had to make these assumptions for any $y$ of bounded norm, not just density operators. It now follows from Eq.~\eqref{eq:122} that:
\begin{equation}
\|x(t) - y(t)\|_p \leq c_y t \|\mathcal{E}\|_p\ .
\end{equation}
\end{proof}
In particular, if $\mc{L}$ is a positive trace-preserving map then $c_y=1$.

The only error for which we can apply Lemma~\ref{LemLinear} in this work is the Markov error, which we present here:
\begin{equation}
\|\rho^{BM}_t-\rho^{B}_t\|_1  \leq 16 \frac{\tau_B c_B c_{BM} t}{\tau_{SB}^2} \ ,
\end{equation}
where $c_y =c_{BM}$ and we used Eq.~\eqref{eq:E_M}, which contains $c_{B}$. We instead used a loose bound from Eq.~\eqref{cBMarkov} since it results in a more compact expression.
Apart from this we have no occasion to use Lemma~\ref{LemLinear} in this work, since not all the equations appearing in our derivation are invertible. This is why we instead use Lemma~\ref{factDif} or prove the necessary results independently.

\section{Proof of Lemma \ref{LemTave}}
\label{LemProof}
\begin{proof}
The formal solutions allow us to write:
\begin{equation}
\delta\rho(t)=   \int_0^t \overline{\mathcal{L}_\theta}\pi(\theta)  d\theta - \int_0^t  \mathcal{L}_{\theta}\rho(\theta)d\theta \ .
\end{equation}
Substituting $\pi(t) = \rho(t) + \delta\rho(t)$ gives:
\begin{equation}
\delta\rho(t)= \int_0^t (\overline{\mathcal{L}_\theta}\rho(\theta) -\mathcal{L}_{\theta}\rho(\theta))  d\theta + \int_0^t \overline{\mathcal{L}_\theta}\delta\rho(\theta)  d\theta\ .
\end{equation}
Making the time-averaging explicit:
\begin{equation}
\delta\rho(t)= \frac{1}{T_a}\int_0^t d\theta \int_{-T_a/2}^{T_a/2} d\tau [\mathcal{L}_{\theta+\tau}\rho(\theta) -\mathcal{L}_{\theta}\rho(\theta)]   + \int_0^t \overline{\mathcal{L}_\theta}\delta\rho(\theta)  d\theta\ .
\label{eq:113}
\end{equation}
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'}\rho(\theta' -\tau')$. The change of variables is $\theta + \tau = \theta', ~ \tau = \tau'$. In the second term $-\mathcal{L}_{\theta}\rho(\theta)$ the change is trivial $\theta =\theta', ~ \tau = \tau'$. 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'}\rho(\theta' -\tau')$ to be the same as the second term. This is illustrated in Fig.~\ref{timeAvtrick} and the corresponding integrals are:
\begin{figure}[t]
\includegraphics[width=0.8\linewidth]{timeAvtrick.png}
\caption{Times contributing to $\delta \rho(t)$ for (a) $t>T_a/2$ (b) $t<T_a/2$.  }
\label{timeAvtrick}
\end{figure}
\begin{equation}
\int_0^t d\theta \int_{-T_a/2}^{T_a/2}d\tau =\int_{-T_a/2}^{T_a/2}d\tau'\int_{\tau'}^{t+\tau'}d\theta'  = \int_{-T_a/2}^{T_a/2}d\tau' \left(\int_0^{t}d\theta' + \int_t^{t+\tau'}d\theta' -\int_0^{\tau'}d\theta' \right)
\end{equation}
After this change of variables and collecting all the contributions from the $\mathcal{L}_{\theta}\rho(\theta)$ term in Eq.~\eqref{eq:113} into a single term, we get:
\begin{equation}
\delta\rho(t)= \frac{1}{T_a}\int_{-T_a/2}^{T_a/2}d\tau' \int_0^{t}d\theta'\mathcal{L}_{\theta'}(\rho(\theta' -\tau')-\rho(\theta')) +{\frac{1}{T_a}\int_{-T_a/2}^{T_a/2}d\tau' \left(\int_t^{t+\tau'}d\theta' -\int_0^{\tau'}d\theta' \right)\mathcal{L}_{\theta'}\rho(\theta' -\tau')} + \int_0^t \overline{\mathcal{L}_\theta}\delta\rho(\theta)  d\theta .
\end{equation}
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 $0$ 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.~\ref{timeAvtrick}. 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: \begin{equation}
\|\delta\rho(t)\|_1\leq \frac{\Lambda }{T_a}\int_0^t d\theta' \int_{-T_a/2}^{T_a/2} d\tau' \|\rho(\theta'-\tau') -\rho(\theta')\|_1   + {\frac{c_\rho \Lambda }{2T_a}\left(T_a^2 -\left(\max\left(T_a-2t,0\right)\right)^2\right)} + \Lambda \int_{0}^t\|\delta\rho(\theta)\|_1d\theta \ ,
\end{equation}
where $\Lambda$ is a number such that
\beq
\Lambda  \geq \sup_{\theta \in [0,t],\|\rho\|_1= 1} \|\mathcal{L}_{\theta}(\rho)\|_1\ .
\eeq
The difference between $\rho$'s at nearby time points is bounded trivially as:
\begin{equation}
\dot\rho = \mathcal{L}\rho\quad \Rightarrow \quad \|\rho(\theta' -\tau') - \rho(\theta')\|_1 \leq c_\rho\Lambda |\tau'| \ ,
\end{equation}
which can be used to simplify:
\bes
\begin{align}
\|\delta\rho(t)\|_1 &\leq \frac{c_\rho\Lambda }{T_a}\int_0^t d\theta' \int_{-T_a/2}^{T_a/2} d\tau' \Lambda |\tau'|   + \frac{c_\rho \Lambda \min(T_a,4t - 4t^2/T_a)}{2} + \Lambda \int_{0}^t\|\delta\rho(\theta)\|_1d\theta   \\
&\leq \frac{c_\rho\Lambda ^2 t T_a}{4}   + \frac{c_\rho\Lambda \min(T_a,4t - 4t^2/T_a)}{2} + \Lambda \int_{0}^t\|\delta\rho(\theta)\|_1d\theta\ .
\end{align}
\ees
We analyzed a similar inequality in the proof of Lemma~\ref{factDif}, 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:
\bes
\begin{align}
\dot{b}_2(t) &= \frac{c_\rho\Lambda ^2  T_a}{4}   + c_\rho\Lambda (2 - 4t/T_a) + \Lambda b_2(t), \quad b_2(0)=0,\quad t\leq T_a/2 \\
\dot{b}_2(t) &= \frac{c_\rho\Lambda ^2  T_a}{4}   + \Lambda b_2(t), \quad t\geq T_a/2
\end{align}
\ees
Solving the first equation, we find:
\bes
\begin{align}
b_2(t) &= 4c_\rho t/T_a + c_\rho\frac{4-2\Lambda T_a -(\Lambda T_a/2)^2}{\Lambda T_a} (1 -e^{\Lambda t}) \label{dFull}\\
b_2(T_a/2) &= c_\rho\frac{4- (\Lambda T_a/2)^2}{\Lambda T_a} + c_\rho e^{\Lambda T_a/2}\frac{(\Lambda T_a/2)^2 + 2\Lambda T_a -4}{\Lambda T_a}  = \frac{c_\rho\Lambda T_a }{2}+ O(\Lambda ^2T_a^2)
\end{align}
\ees
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:
\begin{equation}
b_2(t) = [b_2(T_a/2) + c_\rho\Lambda  T_a/4]e^{\Lambda (t-T_a/2)} - c_\rho\Lambda  T_a/4,   \quad t>T_a/2 \ .
\end{equation}
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:
\begin{equation}
\|\delta\rho(t)\|_1\leq (b_2(\min(t,T_a/2)) + c_\rho\Lambda  T_a/4)\max(e^{\Lambda (t-T_a/2)},1) - c_\rho\Lambda  T_a/4\ .
\end{equation}
The simplest observation is that for $t\leq c/\Lambda $:
\begin{equation}
\|\delta\rho(t)\|_1 \leq O(\Lambda T_a)\ .
\end{equation}
\end{proof}

\section{Derivation of the bound~\eqref{looseBound}}
\label{boundingDetails}

Here we discuss the derivation of the error introduced by time-averaging and restoring complete positivity, starting from Eq.~(\ref{fullWLam}):
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq c_{BM}\left(  e^{\Lambda \max (t,T_a/2)} \frac{3 T_a}{\tau_{SB}} -\frac{  T_a}{\tau_{SB}}\right) +  e^{\Lambda t}\frac{2T_a}{\tau_{SB}} + \frac{4\tau_B}{\Lambda T_a\tau_{SB}}\max(e^{\Lambda (t-T_a/2 )}-1,0 )\ .
\end{align}
Similarly to Eq.~\eqref{weakBound}, 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:
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq ((3c_{BM}+2)e^{\Lambda t}-c_{BM}) \frac{ T_a}{\tau_{SB}}   +(e^{\Lambda t}-1 )\frac{c_\Lambda\tau_B}{ T_a}\ .
\label{startApp1}
\end{align}
At $T_a =\sqrt{c_\Lambda\tau_B\tau_{SB}/(3c_{BM}+2)}$ for $t>T_a/2$ we have:
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq \sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{\Lambda t}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right)\ .
\end{align}
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:
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq\sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{\Lambda \max(t,T_a/2)}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right), \quad T_a = \sqrt{\frac{c_\Lambda\tau_B \tau_{SB}}{3c_{BM}+2}} \ ,
\label{allT}
\end{align}
and this holds for all $t$ and all $\tau_B/\tau_{SB}$. We note that the more careful optimization in App.~\ref{optAside} 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}$.
% it does not really matter since the error is so big so pretty much any bound that is concise enough will do our result justice.

Note that for $t\leq T_a/2$ the expression in Eq.~\eqref{allT} 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.~\eqref{allT} with the trivial bound $ \|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|\leq 1 +c_{BM}$:
\begin{align}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|&\leq \text{min}\left(\sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{\Lambda \max(t,T_a/2)}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right), 1 +c_{BM}\right) \\
&\leq  \sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{\Lambda t+1}-6/5\right)\ .
\label{eq:F4}
\end{align}
\begin{proof} Here we prove Eq.~\eqref{eq:F4}. 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$:
\bes
\begin{align}
\sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{2T_a/c_\Lambda\tau_{SB}}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right)= 1 +c_{BM},  \\
\epsilon^*\sqrt{c_\Lambda(3c_{BM}+2)}\left(2e^{2\epsilon^*/\sqrt{c_\Lambda(3c_{BM}+2)}}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right)= 1 +c_{BM}
\label{maxFrac}
\end{align}
\ees
Relaxing this to an inequality, we find:
\begin{equation}
\epsilon^* \sqrt{c_\Lambda(3c_{BM}+2)}\left(2-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right)\leq 1 +c_{BM}, \quad  \frac{2\epsilon^* \sqrt{c_\Lambda}}{\sqrt{3c_{BM}+2}}\leq 1, \quad e^{2\epsilon^*/\sqrt{c_\Lambda(c_{BM}+2})}\leq e^{1/c_\Lambda}\leq e \ .\label{eIneq}
\end{equation}
where we used $c_\Lambda\geq 1$. Thus any bound on Eq.~\eqref{allT} 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:
\begin{equation}
\frac{\max(t,T_a/2)}{\tau_{SB}} \leq \frac{t}{\tau_{SB}} + \frac{1}{2}\sqrt{\frac{c_\Lambda\tau_B}{(3c_{BM}+2)\tau_{SB}}} \leq \frac{t}{\tau_{SB}} + \frac{\epsilon^* \sqrt{c_\Lambda}}{2\sqrt{(3c_{BM}+2)}}, \quad \text{for} \quad \frac{\tau_B }{\tau_{SB}}\leq \epsilon^{*2}\ .
\end{equation}
Thus we arrive at
\begin{align*}
\|\rho_{BM,I}(t) - {\rho}_{C,I}(t)\|_1\leq\sqrt{\frac{c_\Lambda(3c_{BM}+2)\tau_B}{\tau_{SB}}}\left(2e^{2\epsilon^*/\sqrt{c_\Lambda(3c_{BM}+2)}}e^{\Lambda t}-\left(1 +\frac{c_{BM}}{(3c_{BM}+2)}\right)\right)\ ,
\end{align*}
which evaluates to the r.h.s of Eq.~\eqref{eq:F4} using Eq.~\eqref{eIneq} and $c_{BM}\geq 1  \Rightarrow c_{BM}/(3c_{BM}+2) \geq 1/5$.
\end{proof}
We will now redo this calculation for $T_a = \sqrt{\tau_B\tau_{SB}/5}$ and $c_\Lambda =1$. Using this in Eq.~\eqref{startApp1}:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq \sqrt{\frac{\tau_B}{5\tau_{SB}}} \left( \left(3e^{\Lambda t} -1\right)c_{BM} + \left(7e^{\Lambda t} -5\right)\right)\ .
\end{equation}
Now we find $\epsilon^*$:
\begin{equation}
1+c_{BM}= \frac{\epsilon^*}{\sqrt{5}} \left( \left(3e^{2\epsilon^*/\sqrt{5}} -1\right)c_{BM} + \left(7e^{2\epsilon^*/\sqrt{5}} -5\right)\right) \geq \frac{\epsilon^*}{\sqrt{5}} (2c_{BM}+2), \quad \Rightarrow \quad \frac{2\epsilon^*}{\sqrt{5}} \leq 1\ .
\end{equation}
This allows one to bound:
\begin{equation}
\|\rho_{BM,I}(t) - {\rho}_{C\sqrt{5},I}(t)\|_1\leq \sqrt{\frac{\tau_B}{5\tau_{SB}}} \left( \left(3e^{\Lambda t+1} -1\right)c_{BM} + \left(7e^{\Lambda t+1} -5\right)\right) \label{sqrt5allT}
\end{equation}
for all $t$.

\section{Optimal time}
\label{optAside}

We begin with Eq.~\eqref{startApp1}, and use it to write down the optimum $t$-dependent coarse-graining time:
\begin{equation}
T_a(t) = \sqrt{c_\Lambda\tau_B \tau_{SB}\frac{(e^{\Lambda t}-1 )}{(3c_{BM}+2)e^{\Lambda t}-c_{BM}}}\ .
\end{equation}
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}$:
\begin{equation}
1+c_{BM}(t^*)=((3c_{BM}(t^*)+2)e^{\Lambda t^*}-c_{BM}(t^*)) \frac{ T_a}{\tau_{SB}}   +(e^{\Lambda t^*}-1 )\frac{c_\Lambda\tau_B}{T_a}\ .
\end{equation}
We would like to optimize $T_a$ at $t=t^*$, which would mean that $t^*$ is maximized. The corresponding equation for $t^*$ is:
\begin{equation}
1+c_{BM}(t^*) =2 \sqrt{c_\Lambda\frac{\tau_B}{ \tau_{SB}}(e^{\Lambda t^*}-1 )((3c_{BM}(t^*)+2)e^{\Lambda t^*}-c_{BM}(t^*))}\ ,
\end{equation}
and solving the resulting quadratic equation yields:
\begin{equation}
e^{\Lambda t^*} = \frac{4c_{BM}(t^*)+2 + Y}{2(3c_{BM}(t^*)+2)}, \quad Y =\sqrt{ 4( c_{BM}(t^*) +1)^2 +2(1+c_{BM}(t^*))^2 (3c_{BM}(t^*)+2)(\tau_{SB}/\tau_B c_\Lambda)}\ .
\end{equation}
Thus $T_a$ is:
\begin{equation}
T_a = \sqrt{c_\Lambda\tau_B \tau_{SB}\frac{-2c_{BM}(t^*)-2 + Y}{(3c_{BM}(t^*)+2)(2c_{BM}(t^*) +2 + Y)}}\ .
\end{equation}
And the error is:
\bes
\begin{align}
\|\rho_{BM}(t) - {\rho}_{C*}(t)\|_1&\leq  ((3c_{BM}(t^*)+2)e^{\Lambda\text{max}(t,T_a/2)}-c_{BM}(t^*)) \sqrt{\frac{c_\Lambda\tau_B }{\tau_{SB}} \frac{-2c_{BM}(t^*)-2 + Y}{(3c_{BM}(t^*)+2)(2c_{BM}(t^*) +2 + Y)}}   \\
&+(e^{\Lambda\text{max}(t,T_a/2)}-1 )\sqrt{\frac{c_\Lambda\tau_B}{\tau_{SB}} \frac{(3c_{BM}(t^*)+2)(2c_{BM}(t^*) +2 + Y)}{-2c_{BM}(t^*)-2 + Y}}
\end{align}
\ees
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}$.

\section{Proof of the bound~\eqref{eq:183} for the Redfield limits of integration error}
\label{app:Red-details}

Starting from Eq.~\eqref{eq:182}, we first send the integration limit to infinity:
\beq
\int_0^{4t/\tau_{SB}}\min (1, \frac{4\tau_B}{x\tau_{SB}} + \epsilon_T)e^{-x}dx \leq\int_0^{\infty}\min (1, \frac{4\tau_B}{x\tau_{SB}} + \epsilon_T)e^{-x}dx = \epsilon_T + \int_0^{\infty} \min (1-\epsilon_T, \frac{4\tau_B}{x\tau_{SB}}  )e^{-x}dx  \ .
\eeq
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 \beq
x^* = \frac{4\tau_B}{\tau_{SB} (1-\epsilon_T)}\ .
\eeq
\bes
\begin{align}
\int_0^{\infty} \min (1-\epsilon_T, \frac{4\tau_B}{x\tau_{SB}}  )e^{-x}dx
&=  \int_0^{x^*}(1-\epsilon_T) e^{-x}dx + \int_{x^*}^\infty \frac{4\tau_B}{x\tau_{SB}} e^{-x}dx \\
&= (1-\epsilon_T)(1-e^{-x^*}) +  \frac{4\tau_B}{\tau_{SB}}E_1(x^*) \ ,
\end{align}
\ees
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~\cite{Alzer:1997}[Thm.~2]:
\beq
E_1(x) \leq -\ln(1-e^{-x}) \quad (0<x<\infty) \ .
\eeq
\bes
\begin{align}
\|\rho_{BM}(t) - \rho_R(t)\|_1  &\leq  c_{BM}e^{4t/\tau_{SB}} \left[ \epsilon_T+(1-\epsilon_T)(1-e^{-x^*}) - \frac{4\tau_B}{\tau_{SB}}\ln(1-e^{-x^*}) \right] \\
&\leq c_{BM}e^{4t/\tau_{SB}}  \left[ \epsilon_T + \frac{4\tau_B}{\tau_{SB}}\left(1 - \ln\left(1-e^{-\frac{4\tau_B}{\tau_{SB} (1-\epsilon_T)}}\right)\right)\right] \ ,
\label{eq:183new}
\end{align}
\ees
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:
\bes
\begin{align}
\int_0^{\infty}\min (1, \frac{4\tau_B}{x\tau_{SB}} + \epsilon_T)e^{-x}dx &=\epsilon_T +(1-\epsilon_T) (1-e^{-x^*}) + \int_{x^*}^1 \frac{4\tau_B}{x\tau_{SB}} e^{-x}dx  +\int_{1}^\infty \frac{4\tau_B}{x\tau_{SB}} e^{-x}dx \\
&\leq \epsilon_T +(1-\epsilon_T) x^* + \int_{x^*}^1 \frac{4\tau_B}{x\tau_{SB}} dx  + \int_1^\infty\frac{4\tau_B}{\tau_{SB}}e^{-x}dx  \\
&=  \epsilon_T +(1-\epsilon_T) x^* -  \frac{4\tau_B}{\tau_{SB}} \text{ln}x^* + \frac{4\tau_B}{e\tau_{SB}}    \\
&= \epsilon_T + \frac{4\tau_B}{\tau_{SB} } +  \frac{4\tau_B}{\tau_{SB}} \text{ln}\frac{\tau_{SB} (1-\epsilon_T)}{4\tau_B} + \frac{4\tau_B}{e\tau_{SB}} \quad (x^*< 1 ) \ ,
\end{align}
\ees
while when $x^*\geq 1$ the term with the logarithm is absent:
\begin{align}
{\int_0^{\infty}\min (1, \frac{4\tau_B}{x\tau_{SB}} + \epsilon_T)e^{-x}dx \leq\epsilon_T + \frac{4\tau_B}{\tau_{SB} } + \frac{4\tau_B}{e\tau_{SB}} \quad (x^*\geq 1 ) } \ .
\end{align}
Combining these two cases directly yields Eq.~\eqref{eq:183}, which is an upper bound on the tighter bound presented in Eq.~\eqref{eq:183new}.

%%merlin.mbs apsrev4-1.bst 2010-07-25 4.21a (PWD, AO, DPC) hacked
%Control: key (0)
%Control: author (72) initials jnrlst
%Control: editor formatted (1) identically to author
%Control: production of article title (0) allowed
%Control: page (1) range
%Control: year (0) verbatim
%Control: production of eprint (0) enabled
\begin{thebibliography}{82}%
\makeatletter
\providecommand \@ifxundefined [1]{%
\@ifx{#1\undefined}
\providecommand \@ifnum [1]{%
\ifnum #1\expandafter \@firstoftwo
\else \expandafter \@secondoftwo
\fi
\providecommand \@ifx [1]{%
\ifx #1\expandafter \@firstoftwo
\else \expandafter \@secondoftwo
\fi
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
`\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\providecommand \doibase [0]{http://dx.doi.org/}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
\let\auto@bib@innerbib\@empty
\bibitem [{\citenamefont {Alicki}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
{Alicki}}\ and\ \bibinfo {author} {\bibfnamefont {K.}~\bibnamefont {Lendi}},\
}\href {https://doi.org/10.1007/3-540-18276-4} {\emph {\bibinfo {title}
{Quantum {Dynamical} {Semigroups} and {Applications}}}}\ (\bibinfo
{publisher} {Springer Science \& Business Media},\ \bibinfo {year}
{2007})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Breuer}\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.-P.}\ \bibnamefont
{Breuer}}\ and\ \bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
{Petruccione}},\ }\href
{https://doi.org/10.1093/acprof:oso/9780199213900.001.0001} {\emph {\bibinfo
{title} {The Theory of Open Quantum Systems}}}\ (\bibinfo  {publisher}
{Oxford University Press},\ \bibinfo {address} {Oxford},\ \bibinfo {year}
{2002})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {{C.W. Gardiner and P.
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibnamefont {{C.W. Gardiner and P.
Zoller}}},\ }\href {https://www.springer.com/gp/book/9783540223016} {\emph
{\bibinfo {title} {{Quantum Noise}}}},\ \bibinfo {series} {Springer Series in
Synergetics}, Vol.~\bibinfo {volume} {56}\ (\bibinfo  {publisher}
{Springer},\ \bibinfo {address} {Berlin},\ \bibinfo {year}
{2000})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Davies}(1974)}]{Davies:74}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.~B.}\ \bibnamefont
{Davies}},\ }\bibfield  {title} {\enquote {\bibinfo {title} {Markovian master
equations},}\ }\href {http://dx.doi.org/10.1007/BF01608389} {\bibfield
{journal} {\bibinfo  {journal} {Communications in Mathematical Physics}\
}\textbf {\bibinfo {volume} {39}},\ \bibinfo {pages} {91--110} (\bibinfo
{year} {1974})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lindblad}(1976)}]{Lindblad:76}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
{Lindblad}},\ }\bibfield  {title} {\enquote {\bibinfo {title} {On the
generators of quantum dynamical semigroups},}\ }\href
{http://dx.doi.org/10.1007/BF01608499} {\bibfield  {journal} {\bibinfo
{journal} {Comm. Math. Phys.}\ }\textbf {\bibinfo {volume} {48}},\ \bibinfo
{pages} {119--130} (\bibinfo {year} {1976})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Albash}\ \emph
{et~al.}(2015{\natexlab{a}})\citenamefont {Albash}, \citenamefont {Vinci},
\citenamefont {Mishra}, \citenamefont {Warburton},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
{Albash}}, \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Vinci}},
\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Mishra}}, \bibinfo
{author} {\bibfnamefont {P.~A.}\ \bibnamefont {Warburton}}, \ and\ \bibinfo
{author} {\bibfnamefont {D.~A.}\ \bibnamefont {Lidar}},\ }\bibfield  {title}
{\enquote {\bibinfo {title} {Consistency tests of classical and quantum
models for a quantum annealer},}\ }\href
{https://doi.org/10.1103/PhysRevA.91.042314} {\bibfield  {journal} {\bibinfo
{journal} {Phys. Rev. A}\ }\textbf {\bibinfo {volume} {91}},\ \bibinfo
{pages} {042314--} (\bibinfo {year} {2015}{\natexlab{a}})}\BibitemShut
\bibitem [{\citenamefont {Boixo}\ \emph {et~al.}(2016)\citenamefont {Boixo},
\citenamefont {Smelyanskiy}, \citenamefont {Shabani}, \citenamefont {Isakov},
\citenamefont {Dykman}, \citenamefont {Denchev}, \citenamefont {Amin},
\citenamefont {Smirnov}, \citenamefont {Mohseni},\ and\ \citenamefont
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
{Boixo}}, \bibinfo {author} {\bibfnamefont {V.~N.}\ \bibnamefont
{Smelyanskiy}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
{Shabani}}, \bibinfo {author} {\bibfnamefont {S.~V.}\ \bibnamefont {Isakov}},
\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Dykman}}, \bibinfo
{author} {\bibfnamefont {V.~S.}\ \bibnamefont {Denchev}}, \bibinfo {author}
{\bibfnamefont {M.~H.}\ \bibnamefont {Amin}}, \bibinfo {author}
{\bibfnamefont {A.~Y.}\ \bibnamefont {Smirnov}}, \bibinfo {author}
{\bibfnamefont {M.}~\bibnamefont {Mohseni}}, \ and\ \bibinfo {author}
{\bibfnamefont {H.}~\bibnamefont {Neven}},\ }\bibfield  {title} {\enquote
{\bibinfo {title} {Computational multiqubit tunnelling in programmable
quantum annealers},}\ }\href {http://dx.doi.org/10.1038/ncomms10327}
{\bibfield  {journal} {\bibinfo  {journal} {Nat Commun}\ }\textbf {\bibinfo
{volume} {7}} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Albash}\ \emph
{et~al.}(2015{\natexlab{b}})\citenamefont {Albash}, \citenamefont {Hen},
\citenamefont {Spedalieri},\ and\ \citenamefont {Lidar}}]{Albash:2015pd}%
\BibitemOpen
\bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
{Albash}}, \bibinfo {author} {\bibfnamefont {I.}~\bibnamefont {Hen}},
\bibinfo {author} {\bibfnamefont {F.~M.}\ \bibnamefont {Spedalieri}}, \ and\
\bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont {Lidar}},\ }\bibfield
{title} {\enquote {\bibinfo {title} {Reexamination of the evidence for
entanglement in a quantum annealer},}\ }\href
{https://doi.org/10.1103/PhysRevA.92.062328} {\bibfield  {journal} {\bibinfo
{journal} {Physical Review A}\ }\textbf {\bibinfo {volume} {92}},\ \bibinfo
{pages} {062328--} (\bibinfo {year} {2015}{\natexlab{b}})}\BibitemShut
\bibitem [{\citenamefont {Mishra}\ \emph {et~al.}(2018)\citenamefont {Mishra},
\citenamefont {Albash},\ and\ \citenamefont {Lidar}}]{Mishra:2018}%
\BibitemOpen