Completely positive master equation for arbitrary driving and
small level spacing
Evgeny Mozgunov
1
and Daniel Lidar
1,2,3,4
1
Center for Quantum Information Science & Technology, University of Southern California, Los Angeles, California 90089, USA
2
Department of Electrical and Computer Engineering, University of Southern California, Los Angeles, California 90089, USA
3
Department of Chemistry, University of Southern California, Los Angeles, California 90089, USA
4
Department of Physics and Astronomy, University of Southern California, Los Angeles, California 90089, USA
Markovian master equations are a ubiquitous tool in the study of open quantum systems,
but deriving them from first principles involves a series of compromises. On the one hand, the
Redfield equation is valid for fast environments (whose correlation function decays much faster
than the system relaxation time) regardless of the relative strength of the coupling to the system
Hamiltonian, but is notoriously non-completely-positive. On the other hand, the Davies equation
preserves complete positivity but is valid only in the ultra-weak coupling limit and for systems with
a finite level spacing, which makes it incompatible with arbitrarily fast time-dependent driving.
Here we show that a recently derived Markovian coarse-grained master equation (CGME), already
known to be completely positive, has a much expanded range of applicability compared to the
Davies equation, and moreover, is locally generated and can be generalized to accommodate
arbitrarily fast driving. This generalization, which we refer to as the time-dependent CGME,
is thus suitable for the analysis of fast operations in gate-model quantum computing, such as
quantum error correction and dynamical decoupling. Our derivation proceeds directly from the
Redfield equation and allows us to place rigorous error bounds on all three equations: Redfield,
Davies, and coarse-grained. Our main result is thus a completely positive Markovian master
equation that is a controlled approximation to the true evolution for any time-dependence of the
system Hamiltonian, and works for systems with arbitrarily small level spacing. We illustrate
this with an analysis showing that dynamical decoupling can extend coherence times even in a
strictly Markovian setting.
1 Introduction
Modeling experiments requires taking into account that physical systems are open, i.e., not ideally isolated
from their environments. Usually an environment contains many more degrees of freedom than the system
itself, but is not in any interesting phase of matter where the detailed modeling of each of those degrees of
freedom is required. This makes the problem of modeling open systems tractable. More precisely, the problem
is tractable under a list of assumptions on the parameters of the bath (environment) and its interaction with
the system [1–3].
One natural approach is to consider the case when this interaction is weak compared to all other energy
scales in the problem. In this limit, Davies showed [4] that the dynamics of the system are given by a
Markovian master equation with what is now called Davies generators. Rigorous bounds on the error between
the solution of this equation and the true evolution can be obtained [presented here in Eq. (264)]. The largest
contribution to the error comes from making the rotating wave approximation (RWA).
1
After the system
evolves for a time set by a relevant relaxation timescale, the error in its density matrix is given by:
system-bath coupling strength ×
q
bath correlation time/system level-spacing . (1)
1
This approximation is stated explicitly in Ref. [3], page 87, around Eq. (3.6.67).
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 1
arXiv:1908.01095v2 [quant-ph] 30 Jan 2020
This quantity, in general, becomes exponentially large in the system size due its denominator, so that the
coupling strength or the bath correlation time need to be sufficiently small for the use of Davies generators to
be justified, up to some upper system size limit. In practice, for convenience and also because it may apply
in a range that is larger than the rigorous but conservative bounds imply, the Lindblad master equation [5]
with Davies generators is often used outside its formal range of applicability in modeling of experiments
(e.g., Ref. [6–9]). In such cases it should be understood as a semi-empirical model of open system dynamics,
that possesses all the physical properties of the dynamics without reproducing them exactly. Unfortunately,
since we naturally wish to apply modeling tools to problems for which we do not know a priori what the
correct answer is, this approach cannot provide correctness guarantees. However, other methods are even
more susceptible to this problem, e.g., a path-integral (non-master equation) approach based on integrating
out the bath [10–12]: this involves integration over a very high-dimensional space, so methods [13] that bring
it back into the realm of numerical feasibility usually involve uncontrollable approximations.
One of the important requirements of open system dynamics is complete positivity [14], under which
density matrices whose eigenvalues are by definition non-negative are mapped to other such matrices, even
when applied to a subsystem of a larger system.
2
This property is equivalent to the Markovian master
equation being in Lindblad canonical form [5, 16]. There are numerous other master equations (e.g., Refs. [17–
39]), some in Lindblad form and some not, whose domains of validity and ranges of applicability were not
always studied in detail. One promising candidate is a coarse-grained master equation (CGME) [40], which
is in Lindblad form and thus automatically completely positive (CP). In this work we show that a suitably
generalized variant of this equation, that includes an arbitrarily time-dependent system Hamiltonian (hence
called the “time-dependent CGME”), has a much milder error than (1):
system-bath coupling strength ×bath correlation time , (2)
which does not diverge with the system size n for local observables, and has only a polynomial [O(n
α
) with
1/2 ≤ α ≤ 1] prefactor for nonlocal observables. We note that the full form of the error estimate (2) is not
uniform over the evolution time, unlike some of the error bounds developed elsewhere [39, 41]. We also show
that this new master equation is an improvement over the Lindblad equation with Davies generators in that
it has a local structure, compatible with approximate simulation methods with asymptotically smaller costs
[O(n) vs. O(exp(n))], and that the error estimate (2) is valid for an arbitrary time-dependence of the system
Hamiltonian. In particular, the time-dependent CGME is compatible with the assumption of arbitrarily fast
gates, often made in the circuit model of quantum computing, e.g., in the analysis of fault-tolerant quantum
error correction [42]. This assumption was shown to be incompatible with the derivation of the Davies master
equation [43].
This paper is structured as follows. We summarize our main results in Sec. 2, where we present the main
master equations involved in this work: Redfield and Davies-Lindblad (well known), and coarse-grained, both
time-independent and time-dependent (new). We also present bounds that provide the range of applicability
of each type of master equation (new). We show that these bounds can be expressed entirely in terms of
only two timescales, namely the fastest system decoherence timescale, and the characteristic timescale of the
decay of the bath correlation function. The reader interested only in the results and not in the details of
derivations can safely read only this section and then skip to the conclusions in Sec. 7. Derivations begin in
Sec. 3, where we present a simple, new derivation of the time-independent CGME. We express the equation
in CP (Lindblad) form, state the range of validity for different approximations made in the derivation, and
thus for the equation itself. We also compare the range of applicability with other master equations, and note
the spatial locality of the Lindblad generators of the CGME. At this point we are ready to address the case
of time-dependent Hamiltonians. The derivation of the equation for this case (which happens to result in
exactly the same form) is given in Sec. 4. In particular, this equation is well suited for the simulation of open
system adiabatic quantum computation, but also for dynamical decoupling, which involves the opposite limit
of very fast system dynamics. We note that while we will sometimes refer to qubits, none of the results we
discuss in this work are limited to qubits, and in fact any finite multi-level system, or interacting set of such
systems, is captured by the formalism. We study some applications and examples in Sec. 5, including error
suppression using a dynamical decoupling protocol, and numerical results of the comparison of our master
2
It should be noted that complete positivity can be relaxed without losing physicality; see, e.g., Ref. [15].
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 2
equation with the Redfield and Lindblad-Davies master equations. In the remainder of the paper we derive
the range of validity of the various master equations we study. We give a detailed treatment of the Born
approximation and the other approximations involved, in terms of rigorous bounds presented in Sec. 6, where
we also discuss the generalization to multiple coupling terms and analyze the scaling of the error bounds for
large system sizes. Conclusions are presented in Sec. 7. Various additional technical details are given in the
Appendix.
2 Three master equations
This section provides a summary of our main results. We first present a brief background to define basic
concepts and notation, followed by the definition and properties of the two main timescales we use later to
provide bounds and ranges of applicability. We then summarize the three master equations we focus on in
this work, followed by upper bounds on the distance of their solutions from the true state.
2.1 Background
Consider a system interacting with its environment as described by the total Hamiltonian
H
tot
= H ⊗ I
b
+ A ⊗ B + I ⊗ H
b
. (3)
Here H is the time-independent system Hamiltonian (we treat the time-dependent case in Sec. 4), H
b
is the
bath Hamiltonian, and A and B are, respectively, Hermitian system and bath operators. The more general
situation, with V = A ⊗ B replaced by
P
i
A
i
⊗ B
i
, is easily treated as well (see Sec. 3.7), but we avoid it
here to simplify the notation. We choose A to be dimensionless with kAk = 1 [the operator norm is defined
in Eq. (36)] and B to have energy units, but we deliberately do not set kBk = 1, since we wish to account
for baths (such as harmonic oscillator baths) for which kBk can be infinite. We also set ~ ≡ 1 throughout,
so that energy and frequency have identical units.
Let the eigenstates of H be {|ni} , and the corresponding eigenvalues {E
n
}. Equivalently, H =
P
n
E
n
Π
n
,
where Π
n
is a projector onto the eigenspace with eigenvalue E
n
, and Π
m
Π
n
= δ
mn
Π
n
. Note that in the
interaction picture A(t) = U
†
(t)AU(t) =
P
nm
A
nm
e
−iE
mn
t
|nihm| where E
mn
= E
m
− E
n
, A
nm
= hn|A |mi,
and U is the solution of
˙
U = −iHU. Let
A
ω
=
X
mn:E
mn
=ω
A
nm
|nihm| =
X
mn:E
mn
=ω
Π
n
AΠ
m
= A
†
−ω
, (4)
so that
3
A(t) =
X
ω
A
ω
e
−iωt
. (5)
Here ω runs over all the energy differences (Bohr frequencies) between the eigenvalues of H.
We assume henceforth that the bath is stationary — [ρ
b
, H
b
] = 0 where ρ
b
is the bath state — and
introduce the bath correlation function C(t) and its Fourier transform γ(ω), known as the bath spectral
density:
C(t) = Tr[ρ
b
B(t)B] = C
∗
(−t) =
1
2π
Z
∞
−∞
e
−iωt
γ(ω)dω (6a)
γ(ω) =
Z
∞
−∞
C(t)e
iωt
dt = f(ω) + f
∗
(ω) , (6b)
where B(t) = e
iH
b
t
Be
−iH
b
t
and
f(ω) =
Z
∞
0
C(t)e
iωt
dt =
1
2
γ(ω) + iS(ω) , S(ω) =
1
2i
[f(ω) −f
∗
(ω)] = S
∗
(ω) . (7)
3
The choice of the sign for this notation, as well as other notation choices, follow the textbook [2], p.133-134.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 3
Note that the real-valued functions S(ω) and γ(ω) are related via a Kramers-Kronig transform
S(ω) =
1
2π
Z
∞
−∞
P
1
ω − ω
0
[γ(ω
0
)]dω
0
, (8)
where P
1
x
[f] = lim
→0
R
−
f(x)
x
dx is the Cauchy principal value.
4
When A ⊗B is replaced by
P
i
A
i
⊗B
i
,
we note the positivity of γ
ij
(ω) =
R
dte
iωt
Tr[ρ
b
B
i
(t)B
j
] as a matrix in the indices i, j (of course this reduces
to γ(ω) ≥ 0 in the scalar case). This is a non-trivial fact that is usually associated with Bochner’s theorem,
as e.g., in the textbook [2]. In Appendix A we give two different proofs, one using Bochner’s theorem and
another that is direct. If we assume not only that the bath state is stationary, but that it is also in thermal
equilibrium at inverse temperature β, i.e., ρ
b
= e
−βH
b
/Z, Z = Tre
−βH
b
, then it follows that the correlation
function satisfies the time-domain Kubo-Martin-Schwinger (KMS) condition [2, 44]:
C(t) = Tr[ρ
b
B(0)B(t + iβ)] . (9)
If in addition the correlation function is analytic in the strip between t = −iβ and t = 0, then it follows that
the Fourier transform of the bath correlation function satisfies the frequency domain KMS condition [2, 44]:
γ(−ω) = e
−βω
γ(ω) ∀ω ≥ 0 . (10)
We note that the thermal equilibrium assumption is not necessary for the results we derive in this work, and
we never use it in our proofs. We mention it here for completeness, and also since we use it in one of our
dynamical decoupling examples later on. Finally, note that C(t) has units of frequency squared.
2.2 The bath correlation time τ
B
and the fastest system decoherence timescale τ
SB
We define the two key quantities
1
τ
SB
=
Z
∞
0
|C(t)|dt (11a)
τ
B
=
R
T
0
t|C(t)|dt
R
∞
0
|C(t)|dt
(11b)
Here T is the total evolution time, used as a cutoff which can often be taken as ∞. We show later that in the
interaction picture k ˙ρk
1
≤ 4kρk
1
/τ
SB
, where ρ is the system density matrix and the trace norm defined in
Eq. (35), so that we can interpret τ
SB
as the fastest system decoherence timescale, or timescale over which ρ
changes due to the coupling to the bath, in the interaction picture (i.e., every other system decay timescale,
including the standard T
1
and T
2
relaxation and dephasing times, must be longer). The quantity τ
B
is the
characteristic timescale of the decay of C(t). We return in Sec. 6 to why we define these timescales in this
manner, but note that Eq. (11b) becomes an identity if we choose |C(t)| ∝ e
−t/τ
B
and take the limit T → ∞.
We further note that τ
SB
and τ
B
are the only two parameters relevant for determining the range of
applicability of the various master equations; in particular, nothing about the norm or time-dependence of
H
tot
affects the accuracy of our time-dependent CGME, given below. In particular, H
tot
can be arbitrarily
large or small in the operator norm (strong coupling, unbounded bath), or have an arbitrarily large derivative
(non-adiabatic regime). This remark is important in light of previous approaches, so we expand on it some
more.
In previous work it was common to extract a dimensionful coupling parameter g, i.e., to replace B by
g
˜
B, where k
˜
Bk = 1. One then defines τ
B
=
R
∞
0
|
˜
C(t)|dt in terms of a dimensionless correlation function
˜
C(t) = C(t)/g
2
[31]. One problem with this approach is the arbitrariness of distributing the numerical
factors between g and
˜
B. Another is that it precludes unbounded baths, such as oscillator baths, for which
kBk = ∞. In contrast, the formalism we present here is applicable even when kBk diverges. Furthermore,
kBk contains an extra scale that does not carry any new information about the range of applicability of the
master equations we discuss and derive. By not introducing this extra scale we highlight that there are only
two free parameters in our analysis of the error: τ
B
and τ
SB
, and no other information about the bath is
4
This follows from the identity
R
∞
0
dτe
ixτ
= πδ(x) + iP
1
x
.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 4
needed. I.e., even though different baths will lead to different equations, our results are universal for any
bath with the same values of these two parameters (in fact only their dimensionless ratio matters).
Note that we can relate τ
B
and τ
SB
to the spectral density. First, using γ(ω) > 0 and C(t) = C
∗
(−t):
γ(ω) = |γ(ω)| ≤
Z
∞
−∞
|C(t)|dt = 2
Z
∞
0
|C(t)|dt =
2
τ
SB
. (12)
The dephasing time T
2
and relaxation time T
1
for a single qubit are usually defined as ∝ 1/γ(0) and ∝ 1/γ(ω)
respectively, where ω is the qubit operating frequency (e.g., see Ref. [45]), so we see that τ
SB
. T
1
, T
2
, which
illustrates our earlier comment that τ
SB
is the fastest system decoherence time-scale. Second, γ
0
(ω) =
i
R
∞
−∞
tC(t)e
iωt
dt, so that
|γ
0
(ω)| ≤
Z
∞
−∞
|t||C(t)|dt = 2
Z
∞
0
t|C(t)|dt
T →∞
= 2
τ
B
τ
SB
, (13)
where the limit is taken in Eq. (11b). And lastly,
dγ(−ω)
dω
= −
dγ(ω)
dω
|
−ω
, while differentiating the KMS condition
yields
dγ(−ω)
dω
= −βe
−βω
γ(ω) + e
−βω
dγ(ω)
dω
, so that:
γ
0
(0) =
1
2
βγ(0) . (14)
This implies that γ
0
(0) > 0, so that under the KMS condition (10) we can replace Eq. (13) by
βγ(0) ≤ 4
τ
B
τ
SB
. (15)
The upper bounds on the approximation errors of all the master equations we present below involves the
ratio τ
B
/τ
SB
[see, e.g., Eqs. (26)-(28)]. The smallness of this ratio is sufficient for the CGME to be useful,
while additional assumptions are required for the Davies-Lindblad equation and some versions of the Redfield
equation. It follows from Eq. (15) that for our bound to be small it is necessary that the temperature β
−1
is not too low and/or the spectral density at ω = 0 (roughly the same as the dephasing rate T
−1
2
) is not too
large. The simple condition βγ(0) 1 already rules out diverging spectral densities such as the case of 1/f
noise, though this can be ameliorated by introducing a low-frequency cutoff.
We now present the three master equations, in the Schrödinger picture.
2.3 Redfield master equation
For the derivation see Sec. 3.2. This equation was known earlier than the Davies-Lindblad equation [19], and
is often used in quantum chemistry. It is not CP and hence cannot be put in Lindblad form:
˙ρ
R
(t) = −i[H, ρ
R
(t)] + (Aρ
R
A
f
− ρ
R
A
f
A + h.c.) , (16)
where we defined the “filtered" operator:
A
f
=
Z
∞
0
C(−t)A(−t)dt =
X
ω
A
ω
Z
∞
0
C(−t)e
iωt
dt =
X
ω
f
∗
(−ω)A
ω
, (17)
where f(ω) is defined in Eq. (7), and we used the subscript R to denote that the density matrix satisfies the
Redfield equation (we use a similar subscript notation below to distinguish the solutions of different master
equations). There is no natural way to separate the Lamb shift. One of the benefits of the Redfield equation
is that there is only one generator A
f
per interaction with the environment A ⊗B. We will show that as long
as it preserves positivity, the Redfield equation is more accurate than the Davies-Lindblad master equation.
2.4 Davies-Lindblad master equation
For the derivation see Sec. 3.3. This is the familiar result [4]
˙ρ
D
(t) = −i[H + H
LS
, ρ
D
] +
X
ω
γ(ω)(A
ω
ρ
D
A
†
ω
−
1
2
{ρ
D
, A
†
ω
A
ω
}) , (18a)
H
LS
= −
X
ω
S(ω)A
†
ω
A
ω
, S(ω) =
1
2i
Z
∞
−∞
sgn(t)C(t)e
iωt
dt , (18b)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 5
where H
LS
is the Lamb shift. We note that ad hoc forms of the Lindblad equation are often written down and
used without justification, e.g., with just one Lindblad term per qubit (e.g., σ
−
or I −σ
z
[46, 47]). In reality
the Davies generators are derived from first principles, i.e., from the description of the total Hamiltonian of
the system and bath. The number of Davies generators is unfortunately exponential in the number of qubits
n:
P
ω
is over all 4
n
energy differences.
2.5 Coarse-grained master equation (CGME) for time-independent or time-dependent system Hamilto-
nians
This is our main new set of results, generalizing Ref. [40] to the time-dependent case.
2.5.1 The time-independent case
For the derivation see Sec. 3.4. First, considering time-independent system Hamiltonians, we obtain:
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
(t)] +
Z
∞
−∞
d
A
ρ
C
(t)A
†
−
1
2
n
ρ
C
(t), A
†
A
o
, (19)
where the Lindblad operators are
A
=
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
e
it
A(t)dt (20a)
=
X
ω
f(, ω)A
ω
, f(, ω) =
s
γ()T
a
2π
sinc[T
a
( − ω)/2] , (20b)
with sinc(x) ≡ sin(x)/x, and where T
a
is the averaging, or coarse-graining time, discussed below. The
continuous family of terms labeled by is an unusual form of the Lindblad equation, but is manifestly CP.
The more standard form in terms of a discrete sum over transition frequencies is equivalent to this one, and
is given in Eq. (71) below. The range of applicability of the CGME is only slightly smaller than that of
the Redfield ME, as we discuss in Sec. 3.5. Much like for Davies generators, there are exponentially many
frequency differences, but unlike the Davies case the discretization of the continuous integral
R
d makes it
possible to keep the number of generators constant in the system size for each A ⊗B term. These results are
presented in Sec. 3.6.
The same applies for the Lamb shift, which is
H
LS
=
X
ωω
0
F
ωω
0
A
ω
0
A
ω
(21a)
F
ωω
0
=
1
2T
a
ω
+
Re
Z
T
a
0
dθ
e
i(ωθ−T
a
ω
+
)
− e
−i(ω
0
θ−T
a
ω
+
)
C(θ) = F
∗
ωω
0
= F
−ω
0
,−ω
, (21b)
where ω
+
=
ω+ω
0
2
. Note that F
ωω
0
is well defined in the limit ω
+
→ 0. The results above are generalized to
the n-qubit setting in Sec. 3.7.
2.5.2 The time-dependent case
For the derivation see Sec. 4. In the case of a time-dependent system Hamiltonian H(t) we obtain the same
form, and even the same range of applicability, except that all the operators are now time-dependent:
˙ρ
C
(t) = −i[H(t) + H
LS
(t), ρ
C
] +
Z
∞
−∞
d
A
(t)ρ
C
A
†
(t) −
1
2
n
ρ
C
, A
†
(t)A
(t)
o
, (22)
where the time-dependent Lindblad operators are
A
(t) =
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
e
it
1
A(t + t
1
, t)dt
1
, (23)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 6
where A(t
0
, t) = U
†
(t
0
, t)AU(t
0
, t) with U(t
0
, t) = T exp[−i
R
t
0
t
H(s)ds] (the forward time-ordered exponential,
from t
0
on the left to t on the right), and the time-dependent Lamb shift is
H
LS
(t) =
i
2T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)A(t + t
2
, t)A(t + t
1
, t) . (24)
Although the coarse-graining time T
a
is a free parameter, we show that to minimize the upper bound we
derive on the distance between the solutions of the Redfield and coarse-grained master equations it is nearly
optimal to choose
T
a
=
q
τ
SB
τ
B
/5 . (25)
2.6 Error bounds and range of applicability
The error bound of the Redfield master equation is
kρ
true
(t) − ρ
R
(t)k
1
≤ O
τ
B
τ
SB
e
12t/τ
SB
ln
τ
SB
τ
B
, (26)
where ρ
true
(t) denotes the true (approximation-free) state.
The error bound of the Davies-Lindblad master equation is
kρ
true
(t) − ρ
D
(t)k
1
≤ O
τ
B
τ
SB
+
1
√
τ
SB
δE
e
12t/τ
SB
, (27)
where δE =min
i6=j
|E
i
− E
j
| is the level spacing, with E
i
the eigenenergies of the system Hamiltonian H.
The error bound of both the time-independent and time-dependent coarse-grained master equation is:
5
kρ
true
(t) − ρ
C
(t)k
1
≤ O
r
τ
B
τ
SB
e
6t/τ
SB
. (28)
In the above expressions, O(X) is understood for X → 0, specifically for CGME X =
p
τ
B
/τ
SB
e
6t/τ
SB
. There
is a subtlety in the definition of big-O notation that we would like to emphasize. By definition, “f (X) = O(X)
at X → 0” means there exist X
0
, M s.t. for any X ≤ X
0
we have f(X) ≤ MX. Applying this to Eq. (28),
we have the following unpacking of the above statement: there exist universal constants x
0
, M such that
if
p
τ
B
/τ
SB
≤ X
0
and t ≤ (τ
SB
/6)ln(X
0
/
p
τ
B
/τ
SB
), then kρ
true
(t) − ρ
C
(t)k
1
≤ M
p
τ
B
/τ
SB
e
6t/τ
SB
holds.
Note that for sufficiently small τ
B
/τ
SB
, a dimensionless quantity t/τ
SB
can have an arbitrarily large constant
value in this bound. The error can be made arbitrarily small at any fixed value of t/τ
SB
just by choosing a
sufficiently small τ
B
/τ
SB
. This is what is commonly referred to as the weak coupling limit. We emphasize
that the constants X
0
, M involved in this statement are universal, i.e. independent of the parameters of the
equation such as the Hamiltonian and the initial state. In that sense, our bound is uniform. It does depend
on time t as explicitly stated, so it is not uniform in time.
6
Thus, comparing these bounds, we observe
that the Redfield bound is the tightest, which is natural given that it involves the fewest approximations.
However, recall that the Redfield master equation is not CP. When comparing the bounds on the solutions
of the two CP master equations, we observe that unlike the Davies-Lindblad equation the CGME does not
involve the energy gap, which means that the CGME in principle has a much larger range of applicability,
in particular for systems whose gap shrinks with growing system size.
7
While the equations and inequalities can be derived for any bath, we made some extra assumptions to
prove the error bounds as presented above. Specifically, we assumed a Gaussian bath and the convergence
of our series expansion for the Born error, as defined and discussed in detail in Sec. 6.2. For a non-Gaussian
bath, extra time scales relating to higher correlation functions generally appear, and the error bound can be
derived by a straightforward generalization of our approach. It will contain, besides τ
B
and τ
SB
, additional
terms involving those time scales.
5
Strictly, our proof is only for the time-independent case, but there do not seem to be any obstacles for its generalization to
the time-dependent case.
6
Note that in both the Davies-Lindblad and CGME cases the l.h.s. is bounded above by 2 due to the positivity of the density
matrices; this is not necessarily true in the Redfield case.
7
By range of applicability we mean the range of parameters over which the approximation is accurate. This range can be
defined formally in terms of an upper bound on the right-hand side (r.h.s.), e.g.,
p
τ
B
/τ
SB
e
6t/τ
S B
< for some fixed < 1.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 7
3 (Re-)Derivation of the coarse-grained master equation
Our goal in this section is to present a compact and simplified derivation of the CGME found in Ref. [40].
This master equation is in Lindblad form but avoids making use of the RWA, which contributes the largest
approximation error. Instead, the key idea is to introduce a coarse-graining (averaging) timescale, first
exploited in this context in Refs. [20, 22]. In this section we consider the case of a time-independent sys-
tem Hamiltonian, and later generalize our derivation to the time-dependent case, which includes adiabatic
evolution as a special case.
Along the way we also derive the Redfield and Davies-Lindblad master equations.
3.1 Setting up
We start from the total Hamiltonian given in Eq. (3), and let V = A ⊗ B. Recall that A is dimensionless
and B has units of energy. The first few steps in our derivation are standard [2]. In the double system-bath
interaction picture V (t) = U
†
0
(t)V U
0
(t) [where U
0
is the solution of
˙
U
0
= −i(H ⊗I
b
+ I ⊗H
b
)U
0
, U
0
(0) = 1],
and the state satisfies
˙ρ
tot
(t) = −i[V (t), ρ
tot
(t)] (29a)
ρ
tot
(t) = ρ
tot
(0) − i
Z
t
0
[V (τ), ρ
tot
(τ)]dτ. (29b)
Substituting this back into the original equation yields:
˙ρ
tot
(t) = −i[V (t), ρ
tot
(0)] − [V (t),
Z
t
0
[V (τ), ρ
tot
(τ)]]dτ. (30)
The reduced system state in the interaction picture is defined as ρ
true,I
(t) = Tr
b
[ρ
tot
(t)], where the subscript
“true" denotes that this is the correct, completely approximation-free state. The corresponding true system
state in the Schrödinger picture is ρ
true
= U
0
(t)ρ
true,I
U
†
0
(t). The Born approximation,
ρ
tot
(t) = ρ
true,I
(t) ⊗ ρ
b
+ δρ
tot
, (31)
together with the shift of B such that Tr[ρ
b
B] = 0 allows one to write:
˙ρ
true,I
(t) = −Tr
b
[A(t) ⊗B(t),
Z
t
0
[A(τ) ⊗B(τ), ρ
true,I
(τ) ⊗ρ
b
]]dτ + E
B
, (32)
which can be understood as the definition of the Born approximation error E
B
:
E
B
≡ ˙ρ
true,I
(t) + Tr
b
[A(t) ⊗B(t),
Z
t
0
[A(τ) ⊗B(τ), ρ
true,I
(τ) ⊗ρ
b
]]dτ, kE
B
k
1
= O
τ
B
τ
2
SB
!
. (33)
Henceforth we assume that the bath correlation function decays rapidly, namely that τ
B
τ
SB
, where τ
SB
and τ
B
were defined in Eq. (11).
Now, let ρ
B,I
denote the solution of the master equation after the Born approximation:
˙ρ
B,I
(t) = −Tr
b
[A(t) ⊗B(t),
Z
t
0
[A(τ) ⊗B(τ), ρ
B
(τ) ⊗ρ
b
]]dτ . (34)
Later we collect other errors as letters next to B. The error estimate given in Eq. (33) for kE
B
k
1
is derived
in Sec. 6.2 [Eq. (198b)].
We digress briefly to carefully explain the meaning of the norms and big-O notation used here, since this
is important for the remainder of this work. E
B
is an operator acting on the system Hilbert space. The trace
norm kAk
1
is:
kAk
1
≡ Tr
√
A
†
A (35)
The operator norm kXk is defined as follows:
kXk = max
i
λ
X,i
(36)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 8
where λ
X,i
are the eigenvalues of |X| ≡
√
X
†
X (singular values of X) indexed by i. The big-O notation
E = O(x) is taken for times such that t/τ
SB
≤ M and x → 0, i.e., there exist an M-dependent number
M
and an M-dependent constant C
M
such that for any x <
M
the error kEk
1
≤ C
M
x.
It is natural to use the trace norm for density matrices. Indeed, if we want to find the deviation in an
observable O relative to the difference between two states, δρ = ρ
1
− ρ
2
, then
|TrOδρ| ≤ 2
n
kOkkδρk, |TrOδρ| ≤ kOkkδρk
1
. (37)
The second expression gives a tighter bound, if we manage to have the same bound on kδρk
1
as on kδρk.
Fortunately, this will turn out to be the case in this work, and was already observed for the Markov error in
Ref. [31]. The key property used to prove the inequality above is submultiplicativity, valid for any unitarily
invariant norm [48]:
kABk
1
≤ kAkkBk
1
. (38)
For simplicity, we assume the bath state to be stationary, such that C(t, t
0
) = Tr[ρ
B
B(t)B(t
0
)] is invariant
with respect to shifts in time of both arguments.
8
We can then replace C(t, t
0
) by C(t) = Tr[ρ
B
B(t)B], as in
Eq. (6a). Recall that C
∗
(t) = C(−t). After expanding the commutators, one arrives at
˙ρ
B,I
(t) =
Z
t
0
C(τ − t)[A(t)ρ
B,I
(τ)A(τ) −ρ
B,I
(τ)A(τ)A(t)]dτ + h.c. ≡
Z
t
0
K
B,2
t−τ
(ρ
B,I
(τ))dτ , (39)
where K
B,2
t−τ
is a superoperator. Taking the trace norm, we note that:
k˙ρ
B,I
k
1
≤ 4c
B
Z
∞
0
|C(t)|dt = 4c
B
/τ
SB
, (40)
where we used our earlier choice of setting kAk = 1, and the constant c
B
is chosen as an upper bound on
kρ
B,I
k
1
(if this were a CP map, c
B
= 1 would hold). Under the condition ∀t, kρ
test
(t)k
1
= 1 we have
Z
t
0
K
B,2
t−τ
(ρ
test
(t))dτ
1
≤ 4/τ
SB
. (41)
We next introduce the Markov approximation ρ
B,I
(τ) 7→ ρ
B,I
(t):
˙ρ
B,I
(t) =
Z
t
0
C(τ − t)[A(t)ρ
B,I
(t)A(τ) −ρ
B,I
(t)A(τ)A(t)]dτ + h.c. + E
M
, kE
M
k
1
= O
τ
B
τ
2
SB
!
(42a)
˙ρ
BM,I
(t) =
Z
t
0
C(τ − t)[A(t)ρ
BM,I
(t)A(τ) −ρ
BM,I
(t)A(τ)A(t)]dτ + h.c. = L
BM,I
t
(ρ
BM,I
(t)) (42b)
∀t, kρ
test
k
1
= 1 : kL
BM,I
t
(ρ
test
)k
1
≤ 4/τ
SB
(42c)
Equation (42b) is the Redfield equation [2], which is notoriously non-CP (though various fixes have been
proposed [28, 49]). The Markov approximation error E
M
is of same order as the Born approximation error
given in Eq. (33) (for more detail see Sec. 6).
We note that the errors on the r.h.s. are not generally additive: if we try to write a Markovian master
equation for the true (approximation-free) evolution ρ
true,I
, then the error on the r.h.s. will, besides E
B
+E
M
,
contain a correction to the Markov error from the Born error:
˙ρ
true,I
(t) =
Z
t
0
C(τ − t)[A(t)ρ
true,I
(t)A(τ) −ρ
true,I
(t)A(τ)A(t)]dτ + h.c. + E
B
+ E
M
(ρ
true
), (43a)
E
M
(ρ
true
) 6= E
M
(43b)
In this particular case, using methods from Section 6 the difference between bounds on kE
M
(ρ
true
)k
1
and
kE
M
k
1
can be found to be subleading in τ
B
/τ
SB
. Since we wish to present higher orders of the error, we do
not attempt to collect errors as in Eq. (43a). Instead, our derivation will present a sequence of equations,
e.g., ρ
true
, ρ
B
, ρ
BM
, and an error estimate for each step.
8
The derivation can potentially be extended to the general case of a two-time correlation function, i.e., without assuming
translational invariance; indeed the derivation in [40] does not.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 9
After using Eq. (5), i.e., in the frequency representation, Eq. (42b) takes the following form:
˙ρ
BM,I
(t) =
X
ω,ω
0
Z
t
0
dτC(τ − t)e
−i(ωt+ω
0
τ)
(A
ω
ρ
BM,I
A
ω
0
− ρ
BM,I
A
ω
0
A
ω
) + h.c. (44)
We sometimes omit the explicit time dependence of ρ on the r.h.s. since it is always ρ(t) from hereon.
3.2 Redfield master equation
Let us digress briefly in order to establish the form of the Redfield equation presented in Sec. 2. Rotating
Eq. (42b) back to the Schrödinger picture we obtain:
˙ρ
BM
(t) = −i[H, ρ
BM
(t)] +
Z
t
0
C(τ − t)[Aρ
BM
(t)A(τ − t) − ρ
BM
(t)A(τ − t)A]dτ + h.c. (45)
Introducing a new integration variable t
0
= t − τ:
˙ρ
BM
(t) = −i[H, ρ
BM
(t)] +
Z
t
0
C(−t
0
)[Aρ
BM
(t)A(−t
0
) − ρ
BM
(t)A(−t
0
)A]dt
0
+ h.c. (46)
We extend the upper integration limit from t to ∞, which introduces an additional error:
˙ρ
BM
(t) = −i[H, ρ
BM
(t)] +
Z
∞
0
C(−t
0
)[Aρ
BM
(t)A(−t
0
) − ρ
BM
(t)A(−t
0
)A]dt
0
+ h.c. + E
l
(47a)
kE
l
k
1
≤ 4c
BM
Z
∞
t
|C(t
0
)|dt
0
≤
4c
BM
t
Z
T
t
t
0
|C(t
0
)|dt
0
+ 4c
BM
Z
∞
T
|C(t
0
)|dt
0
= O
τ
B
τ
2
SB
+
T
τ
SB
!
, (47b)
where for the first summand in Eq. (47b) we assumed a lower cutoff on the evolution time t > τ
SB
δ > 0
where δ is a small constant number which we absorbed in the big-O notation [i.e., 1/t = O(1/τ
SB
)]. We
analyze the transients happening for t < τ
SB
δ in Sec. 6.6. The constant c
BM
is chosen as an upper bound
on kρ
BM,I
k
1
(again, if this were a CP map, c
BM
= 1 would hold). We will see in Sec. 6 that c
BM
= O(1)
in terms of our big-O notation, specified in Eq. (49). For the second summand, we introduced a new bath
parameter
T
≡ τ
SB
Z
∞
T
|C(t)|dt . (48)
For most physically motivated choices of bath correlation functions
R
∞
0
t|C(t)|dt converges, so we may replace
the total evolution time T → ∞ and then
T
= 0. But, in case τ
B
/τ
SB
diverges as a function of T [Eq. (11b)],
the above bound should be used.
The big-O notation in terms of all of these parameters τ
SB
δ ≤ t ≤ τ
SB
M,
T
, τ
B
/τ
SB
, c
BM
should be
understood as follows: there exist constants C(M, δ), (M, δ) such that for τ
B
/τ
SB
+
T
≤ (M, δ)
kE
l
k
1
≤ C(M, δ, )
τ
B
τ
2
SB
+
T
τ
SB
!
. (49)
As seen in the derivation above, C(δ, ) = max(4c
BM
/δ, 4c
BM
) in fact does not depend on .
Note that the change of integration limit t → ∞ is optional: the Redfield equation can be used without
it, but the change does simplify its form. Eq. (47a) without the error gives:
˙ρ
BMl
(t) = −i[H, ρ
BMl
(t)] +
Z
∞
0
C(−t
0
)[Aρ
BMl
(t)A(−t
0
) − ρ
BMl
(t)A(−t
0
)A]dt
0
+ h.c. (50)
Introducing the “filtered" operator A
f
=
R
∞
0
C(−t
0
)A(−t
0
)dt
0
immediately yields Eq. (16), so that in fact
ρ
BMl
= ρ
R
. Replacing A(−t
0
) in A
f
by its Fourier decomposition [Eq. (5)] gives Eq. (17), for which the
integral can be expressed in multiple ways:
Z
∞
0
C(−t
0
)e
iωt
0
dt
0
θ=-t
0
=
1
2
Z
∞
−∞
(C(θ) −C(θ) sgn(θ))e
−iωθ
dθ =
1
2
γ(−ω) −iS(−ω) = f
∗
(−ω) . (51)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 10
3.3 Davies-Lindblad master equation
As a second digression let us establish the form of the Davies-Lindblad equation presented in Sec. 2. We
start from the form given in Eq. (44), and use ωt + ω
0
τ = ω(t − τ ) + (ω + ω
0
)τ and a change of variables to
t
0
= τ − t to write it as
˙ρ
BM,I
(t) =
X
ω,ω
0
Z
0
−t
dt
0
C(t
0
)e
i[ωt
0
−(ω+ω
0
)(t+t
0
)]
(A
ω
ρ
BM,I
A
ω
0
− ρ
BM,I
A
ω
0
A
ω
) + h.c. (52)
Next we apply the RWA, in which one considers the term e
−i(ω+ω
0
)(t+t
0
)
to be rapidly oscillating and hence
self-cancelling, so the non-negligible contribution is only due to the terms with ω
0
= −ω. This allows us to
write:
˙ρ
D,I
(t) =
X
ω
Z
0
−∞
dt
0
C(t
0
)e
iωt
0
(A
ω
ρ
D,I
A
−ω
− ρ
D,I
A
−ω
A
ω
) + h.c., (53)
where we replaced ρ
BM,I
by ρ
D,I
and ω
0
by −ω. The error E ∼ O(1/τ
SB
) between the right-hand sides of
Eqs. (52) and (53) is large but rapidly oscillating, so the solutions ρ
BM,I
(t) and ρ
D,I
(t) remain close (see
Sec. 3.5).
We now transform back to the Schrödinger picture (recall that this requires multiplying each A
ω
by e
iωt
,
but these cancel now due to the corresponding A
−ω
):
˙ρ
D
(t) = −i[H, ρ
D
] +
X
ω
Z
0
−∞
C(t
0
)e
iωt
0
(A
ω
ρ
D
A
†
ω
− ρ
D
A
†
ω
A
ω
) + C
∗
(t
0
)e
−iωt
0
(A
ω
ρ
D
A
†
ω
− A
†
ω
A
ω
ρ
D
)
dt
0
.
(54)
It is convenient to separate the real and imaginary parts, using C
∗
(t) = C(−t):
Z
0
−∞
dt
0
C(t
0
)e
iωt
0
=
1
2
γ(ω) −
1
2
Z
∞
−∞
dt
0
sgn(t
0
)C(t
0
)e
iωt
0
(55a)
Z
0
−∞
dt
0
C
∗
(t
0
)e
−iωt
0
=
1
2
γ(ω) +
1
2
Z
∞
−∞
dt
0
sgn(t
0
)C(t
0
)e
iωt
0
. (55b)
Using the decomposition given by Eq. (55) in Eq. (54) directly yields the form of the Davies-Lindblad master
equation presented in Eq. (18). We note that versions of the Davies-Lindblad master equation that allow a
time-dependent drive have been derived before (see, e.g., Refs. [31, 43]), but such derivations must always
assume that the driving is adiabatic.
3.4 Coarse-grained master equation
Let us now return to our main goal in this section, the derivation of the CGME. At this point we deviate
from the standard derivations and instead follow the coarse-graining approach [20, 22, 40]. There are two
main steps: (i) Time-averaging of the state, which is done in lieu of the RWA (and reduces to the RWA in the
limit of large time-averaging scale; see Appendix C), (ii) removal of a small part of the integration domain,
which restores complete positivity (new in this work). We describe each in turn, along with the associated
error.
3.4.1 Time-averaging of the state
Before we leave the interaction picture, we coarse-grain, or time-average over T
a
such that τ
B
T
a
τ
SB
.
Specifically, we consider another equation satisfied by a new state ρ
BMT,I
(t), where the time-averaging
operation
1
T
a
R
t+T
a
/2
t−T
a
/2
dt
0
is applied to the r.h.s. of Eq. (44) taken at time t
0
, except that the state is still at
time t on the l.h.s.:
˙ρ
BMT,I
(t) =
X
ω,ω
0
1
T
a
Z
t+T
a
/2
t−T
a
/2
dt
0
Z
t
0
0
dτC(τ − t
0
)e
−i(ωt
0
+ω
0
τ)
(A
ω
ρ
BMT,I
(t)A
ω
0
− ρ
BMT,I
(t)A
ω
0
A
ω
) + h.c. (56)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 11
Alas, the associated error E
T
has kE
T
k
1
= O(1/τ
SB
), which is not small. The reason is that the fast-rotating
terms may have a large derivative. What can be proven is that the solutions of Eqs. (44) and (56) remain
close:
kρ
BM,I
(t) − ρ
BMT,I
(t)k
1
= O(T
a
/τ
SB
) , (57)
where t ≤ cτ
SB
, and the r.h.s. O(T
a
/τ
SB
) does not depend on t, only on c, which is considered a constant
for the purposes of big-O notation. The proof is given in Sec. 6.4, in Lemma 2 (this analysis first appeared
in Ref. [32]). The error we track is now the error of the solution, not its derivative.
We use the following fact about time-local (Markovian) and time-nonlocal (non-Markovian) differential
equations:
Lemma 1. Assume that
˙x(t) = L(x(t)) + E, ˙y(t) = L(y(t)), x(0) = y(0) , (58)
where L is a linear superoperator and Λ is a positive constant such that sup
τ,x:kxk
1
=1
kL
τ
(x)k
1
≤ Λ. Or,
assume that
˙x(t) =
Z
t
0
K
t−τ
(x(τ))dτ + E, ˙y(t) =
Z
t
0
K
t−τ
(y(τ))dτ, x(0) = y(0) , (59)
where K
t−τ
(x) is a linear superoperator and Λ is a positive constant such that sup
t,x(τ):kx(τ)k
1
=1
R
t
0
kK
t−τ
(x(τ))k
1
dτ ≤
Λ. Then:
∀t ≤
c
Λ
: kx(t) − y(t)k
1
≤ (e
c
− 1)
kEk
1
Λ
. (60)
For the proof see Sec. 6.1. We can take Λ = 4/τ
SB
because of Eq. (41, 42c). Using this, we rewrite the
result of the Lemma as follows, for brevity: if t ≤ c/Λ = cτ
SB
/4, then kx −yk
1
= O(τ
SB
kEk
1
). Thus, using
the previously noted results that kE
B
k
1
, kE
M
k
1
= O
τ
B
/τ
2
SB
:
kρ
true,I
(t) − ρ
BMT,I
(t)k
1
≤ kρ
true,I
(t) − ρ
B,I
(t)k
1
+ kρ
B,I
(t) − ρ
BM,I
(t)k
1
+ kρ
BM,I
(t) − ρ
BMT,I
(t)k
1
(61a)
= O(τ
SB
kE
B
k
1
) + O(τ
SB
kE
M
k
1
) + O(T
a
/τ
SB
) = O
τ
B
+ T
a
τ
SB
. (61b)
3.4.2 Neglecting part of the integration domain to regain complete positivity
Upon changing variables to s = t
0
− t, and renaming s as t
0
, Eq. (56) becomes:
˙ρ
BMT,I
(t) =
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t+t
0
0
dτC(τ −t −t
0
)e
−iω(t+t
0
)−iω
0
τ
(A
ω
ρ
BMT,I
A
ω
0
−ρ
BMT,I
A
ω
0
A
ω
) + h.c. (62)
We rotate out of the system interaction picture into the system Schrödinger picture by defining the coarse-
grained state ρ
BMT
= U (t)ρ
BMT,I
U
†
(t) = e
−iHt
ρ
BMT,I
e
iHt
, leaving only the bath interaction picture, by
multiplying each summand by a e
i(ω+ω
0
)t
factor:
9
˙ρ
BMT
(t) = −i[H, ρ
BMT
]+
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t+t
0
0
dτC(τ −t−t
0
)e
−iωt
0
−iω
0
(τ−t)
(A
ω
ρ
BMT
A
ω
0
−ρ
BMT
A
ω
0
A
ω
)+h.c.
(63)
Next, we change variables to τ
0
= τ − t, so that
˙ρ
BMT
(t) = −i[H, ρ
BMT
] +
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t
0
−t
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
(A
ω
ρ
BMT
A
ω
0
− ρ
BMT
A
ω
0
A
ω
) + h.c.
(64)
We define L
BMT
t
(ρ
BMT
) to be the superoperator on the r.h.s. The key trick now is to replace −t by −T
a
/2,
9
Note that the transformation back from the interaction picture is U ˙ρ
I
U
†
. Thus, e.g., A
ω
ρ
I
A
ω
0
7→
UA
ω
ρ
I
A
ω
0
U
†
= UA
ω
U
†
ρUA
ω
0
U
†
, where U =
P
k
e
−iE
k
t
Π
k
. Using Eqs. (4) and (5), we have UA
ω
U
†
=
P
kl
P
mn:E
mn
=ω
e
−i(E
k
−E
l
)t
Π
k
Π
n
AΠ
m
Π
l
=
P
mn:E
mn
=ω
e
−i(E
n
−E
m
)t
Π
n
AΠ
m
= e
iωt
A
ω
, with the last equality holding due
to the fact that E
n
− E
m
is fixed at ω for all combinations of m and n. The claim now follows.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 12
Figure 1: (a) The box outside the green line is neglected in the integration, which restores complete positivity. The area of
the non-negligible part (yellow) is ∼ τ
B
T
a
while the area we nevertheless neglect is ∼ τ
2
B
. (b) Illustration of the interchange
of the order of integration limits used to prove Eq. (67b).
which lets us write (note that Lemma 1 will work for E
p
either in the initial or in the resulting equation, like
we do here):
˙ρ
C
(t) = L
BMT
t
(ρ
C
) − E
p
= −i[H, ρ
C
] +
X
ω,ω
0
x
ωω
0
(A
ω
ρ
C
A
ω
0
− ρ
C
A
ω
0
A
ω
) + h.c. (65a)
x
ωω
0
≡
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t
0
−T
a
/2
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
, (65b)
This replacement cuts a corner of area ∼ τ
2
B
out of the integration domain with non-vanishing C(t), itself of
area ∼ T
a
τ
B
, as illustrated in Fig. 1(a) for t > T
a
/2 (times t < T
a
/2 introduce a transient, as discussed in
more detail in Sec. 6.4). Thus, the corner being cut for t > T
a
/2 is just a τ
B
/T
a
fraction of the whole integral
in Eq. (65b) (ignoring factors of
√
2), which [using Eq. (11a)] is itself upper bounded by 1/τ
SB
in absolute
value. Thus, the error introduced by the corner removal is
kE
p
k
1
= O(τ
B
/(T
a
τ
SB
)) for t > T
a
/2 (66)
[for t < T
a
/2 it can be O(1/τ
SB
) for a short transient]. The constant in big-O notation here is independent of
τ
SB
, τ
B
, T
a
. As we show below, E
p
is the cost of recovering complete positivity, in the sense that the master
equation after subtraction of E
p
is in canonical Lindblad form.
Eq. (65a) is the coarse-grained master equation from Ref. [40]. To see this we need some properties of
the coefficient x
ωω
0
. In particular, we need the properties
x
∗
ωω
0
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
T
a
/2
t
0
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
(67a)
x
ωω
0
= x
−ω
0
−ω
, (67b)
with the first following from C
∗
(t) = C(−t), and the second from interchanging integration limits after a
reflection w.r.t. the τ
0
= −t
0
axis, as illustrated in Fig. 1(b). This allows us to reshuffle the h.c. part of
Eq. (65a), taking the −ω
0
, −ω term in the sum for the ω, ω
0
term in the original part. The resulting equation
is:
˙ρ
C
(t) = −i[H, ρ
C
] +
X
ω,ω
0
x
ωω
0
(A
ω
ρ
C
A
ω
0
− ρ
C
A
ω
0
A
ω
) + x
∗
ωω
0
(A
ω
ρ
C
A
ω
0
− A
ω
0
A
ω
ρ
C
) . (68)
The Imx
ωω
0
part cancels for A
ω
ρA
ω
0
and sends A
ω
0
A
ω
into the Lamb shift:
H
LS
≡ −
X
ωω
0
Im(x
ωω
0
)A
ω
0
A
ω
=
i
2T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)A(t
2
)A(t
1
) , (69)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 13
where to get the second form we used Im(x
ωω
0
) =
1
2i
(x
ωω
0
− x
∗
ωω
0
), combined Eqs. (65b), (67a) into a single
integral using the sgn(t
1
−t
2
) function, and used Eq. (5). The Lamb shift can be further simplified as shown
in Appendix B. The simplified form is the one we presented as a part of our main result [Eq. (21)] in Sec. 2.
We define γ
ωω
0
≡ x
∗
ωω
0
+ x
ωω
0
= 2Rex
ωω
0
and obtain:
γ
ωω
0
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
T
a
/2
−T
a
/2
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
. (70)
Therefore
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
] +
X
ω,ω
0
γ
ωω
0
(A
ω
ρ
C
A
ω
0
−
1
2
{ρ
C
, A
ω
0
A
ω
}) , (71)
which is Eq. (50) from [40] up to a shift in the center of averaging and changing variables A
ω
0
→ A
†
ω
0
, ω
0
→
−ω
0
.
3.4.3 Complete positivity
Complete positivity is not immediately apparent from Eq. (71). Let us demonstrate it next. Substituting
C(t) =
1
2π
R
∞
−∞
e
−it
γ()d [Eq. (6a)], where γ() ≥ 0 (see Appendix A), into Eq. (70), we have:
γ
ωω
0
=
Z
∞
−∞
d
γ()
2πT
a
Z
T
a
/2
−T
a
/2
dt
0
Z
T
a
/2
−T
a
/2
dτ
0
e
−i(τ
0
−t
0
)
e
−iωt
0
−iω
0
τ
0
=
Z
∞
−∞
df(, ω)f
∗
(, −ω
0
) , (72)
where
f(, ω) ≡
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
dte
i(−ω)t
(73)
is a filter function whose integrated form is given in Eq. (20). By using the above decomposition for γ
ωω
0
in
Eq. (71), we obtain the manifestly CP Lindblad form given in Eq. (19).
Note further that it was shown in Ref. [40] that the Davies-Lindblad equation [Eq. (18)] arises as the
T
a
→ ∞ limit of the CGME. Similar ideas have appeared in [37, 38]. We provide this proof for completeness
in Appendix C. This derivation has the advantage that it arrives at the Davies-Lindblad equation without
invoking the (uncontrolled) RWA, and instead shows that the Davies-Lindblad equation is the limit of infinite
coarse graining time of the CGME. One of our new results is to show how the error of the Davies-Lindblad
equation can be controlled for a sufficiently large (but not infinite) coarse graining time; see Sec. 6.5.
3.5 Ranges of applicability of the three master equations
We now discuss the ranges of applicability of the master equations derived above in more rigor than in
Sec. 2.6, while postponing the derivations to Sec. 6. Recall that the ranges are dependent on the timescales
τ
SB
and τ
B
defined in Eq. (11).
Errors are incurred from different approximations made along the way. The very first approximation
made in the derivation of all three master equations is the Born approximation (see Sec. 3.1). Unlike the
other approximations, we do not prove a rigorous error bound for the Born approximation error; we just
provide a bound on the first order contribution to this error. However, under the assumption that the error
converges in the first place, the Born error will turn out to be subleading. Therefore we settle for an estimate
of this error for the sake of simplicity, but we note that a rigorous bound can be obtained [4].
We first collect all the errors introduced in the derivation of the CGME: the difference in solutions
kρ
true,I
(t) − ρ
BMT,I
(t)k
1
[Eq. (61b), which includes the error due to time-averaging], and the error due to
enforcing complete positivity, kE
p
k
1
= O(τ
B
/(T
a
τ
SB
)) [Eq. (66)]. The difference between the solution of the
CGME ρ
C
(t) and the true time evolved state ρ
true
(t) thus satisfies (note that the norm is unitarily invariant
and hence unaffected by transformation to or out of the interaction picture):
kρ
true
(t) − ρ
C
(t)k
1
≤ kρ
true,I
(t) − ρ
BMT,I
(t)k
1
+ kρ
BMT
(t) − ρ
C
(t)k
1
= O
T
a
+ τ
B
τ
SB
+ O(τ
SB
kE
p
k
1
)
(74a)
= O
τ
B
T
a
+
T
a
+ τ
B
τ
SB
(74b)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 14
Here we used Lemma 1 to multiply kE
p
k
1
by τ
SB
. Optimizing T
a
to minimize this error, we obtain T
a
=
√
τ
SB
τ
B
, and:
kρ
true
(t) − ρ
C
(t)k
1
= O
2
r
τ
B
τ
SB
+
τ
B
τ
SB
= O
r
τ
B
τ
SB
for t > T
a
/2 . (75)
Note that if we define the coupling g strength via 1/τ
SB
= g
2
τ
B
(recall the discussion in Sec. 2.2), then the
error bound of Eq. (75) is O(gτ
B
).
We can in fact improve on the result presented in Eq. (75). A careful derivation presented in Sec. 6
generalizes the above result from the t > T
a
/2 case discussed here [recall Eq. (66)] to t > 0. The big-O
notation in Eq. (75) means that there exist numbers (C, δ, c) such that for all 0 < t < cτ
SB
and τ
B
/τ
SB
< δ
the error is kρ
true
(t) − ρ
C
(t)k
1
≤ C(c, δ)
p
τ
B
/τ
SB
. We do not know the function C(c, δ) explicitly, because
we do not keep track of it in our analysis of the Born error [the details of which are given in Sec. 6.2]. The
contributions from all the other errors are known explicitly. Our derivation also allows us to state a stronger
result exhibiting the coefficients and the t-dependence of the first few leading terms of the error as a series
in τ
B
/τ
SB
, since the unknown Born contributions start with O(τ
2
B
/τ
2
SB
). The improved bound is:
kρ
true
(t) − ρ
C
(t)k
1
≤
13e
4t
τ
SB
√
τ
B
1 +
29τ
B
e
8t
τ
SB
τ
SB
!
√
τ
SB
+
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
12 + O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
. (76)
The time t can now be arbitrarily large, as long as the combination e
4t
τ
SB
τ
B
τ
SB
is small. Specifically, the big-O
notation in Eq. (76) means that there exist constants C, such that for all e
4t
τ
SB
τ
B
τ
SB
≤ , O
e
4t
τ
SB
τ
B
τ
SB
≤
Ce
4t
τ
SB
τ
B
τ
SB
.
The time-averaging window leading to Eq. (75) was chosen sub-optimally, as the more careful analysis
yielding Eq. (76) uses T
a
=
p
τ
B
τ
SB
/5. Including higher order corrections in τ
B
/τ
SB
leads to a minor
improvement in the subleading terms in the bound, but does not affect the leading term (see Section 6.4.4).
Let us now discuss the time-dependence of the error. From the above expression it follows that for
any t, kρ
true
(t) − ρ
C
(t)k
1
= O(
p
τ
B
/τ
SB
e
6t/τ
SB
), which is Eq. (28). For the sake of completeness, we
present the errors of the Redfield equation (for
T
= 0 [recall Eq. (48)]; for the general case, see Sec. 6.6):
kρ
true
(t) −ρ
R
(t)k
1
= O
τ
B
τ
SB
ln
τ
SB
τ
B
e
12t/τ
SB
, which is Eq. (26), and the Davies-Lindblad equation: kρ
true
(t) −
ρ
D
(t)k
1
= O
τ
B
τ
SB
+
q
1
τ
SB
δE
e
12t/τ
SB
, which is Eq. (27).
We note that for Eq. (46) there is no transient before the introduction of E
l
:
kρ
true
(t) − ρ
BM
(t)k
1
= (e
4t/τ
SB
− 1)O
τ
B
τ
SB
e
8t/τ
SB
, (77)
thus the error grows from zero linearly at small times, and the slope of the upper bound is O(τ
B
/τ
2
SB
). In
contrast, for the Redfield, CGME and Davies-Lindblad equations, there is a nonzero transient error bound
at t → 0. The error itself starts at zero (our bounds are not tight), but the slope is O(1/τ
SB
) instead.
It is important to emphasize that the solutions of open system master equations do not actually grow
exponentially. In fact, for time-independent system Hamiltonians all of the equations presented have steady
states to which they converge in the operator and trace norm at some rate. If that rate is sufficiently fast,
one can prove a much stronger result [41] than our error bounds: a stability of the open system dynamics
to small perturbations. The stability can be expressed in terms of a time-independent error bound on the
difference in solutions.
We remark that the coefficients we presented and the time-dependence of the error are explicit, while
other work hides them in big-O notation. Our results allow one to put error bars on the numerical solution
of the above equations.
3.6 Discretization, spatial locality, and the Lieb-Robinson bound
Different forms of the CGME are preferred depending on the application. The form given in Eq. (71) with a
discrete sum
P
ω,ω
0
is ready for numerical implementation but at exponential cost, since for a Hilbert space
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 15
of dimension 2
n
the number of terms in the sum is 2
4n
. Numerical solution of the equation already requires
at least 2
2n
numbers just to store the density matrix. It is therefore desirable to reduce the number of terms
on the r.h.s. to a constant in n (or polynomial in n for the general case of multiple interaction terms; see
Sec. 5). We first present a discretization that achieves this goal. It is also desirable to have each term locally
supported, instead of acting on the whole system. This speeds up the numerics, and also allows one to
meaningfully use the equation with methods that use resources that are polynomial in n to store the density
matrix, e.g., matrix product states (MPS) [50, 51]. For a recent implementation of these methods see [46, 47].
We present a way to write down a local equation at the end of this subsection.
In Eq. (19), the number of terms is infinite because of the integration over . In Appendix D we explain
how a discrete approximation can be derived. The result is:
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
] + ∆
X
−k
∗
<k<k
∗
X
=k∆
(A
ρ
C
A
†
−
1
2
n
ρ
C
, A
†
A
o
) + E
k
(78)
where A
is given by Eq. (20) and the values of ∆, k
∗
and kE
k
k
1
are given in Appendix D, and are system-
size-independent as long as [H, A] = O(1). We do not present a similar discrete approximation for the Lamb
shift, leaving it to future work.
We now proceed to show that storing the Lindblad operators A
does not require resources exponential in
the system size n. Recall that A
=
q
γ()
2πT
a
R
T
a
/2
−T
a
/2
e
it
A(t)dt. This form is interesting because A(t) is evaluated
over the interval [−T
a
/2, T
a
/2]. An operator A(t) at those times is local (with exponentially decaying tails)
with a locality radius vT
a
/2, where v is the Lieb-Robinson velocity (v ∼ kH
local
k, the coupling strength in
the Hamiltonian) [52]. More rigorously, we consider a system defined on a graph. Each vertex of the graph
supports a finite-dimensional Hilbert space, and the system Hilbert space is a tensor product of vertex Hilbert
spaces. The edges of the graph are present if there are nontrivial interactions between the corresponding
vertices in the Hamiltonian. That is, for vertices i, j we can define a trace over the degrees of freedom at
other vertices, denoted Tr
i,j
H − Tr
i
H − Tr
j
H. If this operator is nonzero (for traceless H), then there is a
nontrivial interaction and the edge (i, j) is present in the graph. The Lieb-Robinson bound can be used (see,
e.g., Lemma 5 in Ref. [53]) to derive the following decomposition of the operator A
into its local and small
nonlocal part:
A
= A
(vT
a
/2+δr)
+ δA, kδAk ≤
s
γ()T
a
2π
kAke
−cδr
, (79)
where c is a universal (system-independent and -independent) constant, and
A
(x)
=
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
e
it
e
iH
(x)
t
Ae
−iH
(x)
t
dt (80)
is an operator strictly local within the graph distance x from the original location of A. We define the graph
distance as the length (number of edges) of the smallest connected path of edges between two vertices. Here
H
(x)
= Tr
x
H is composed of those Hamiltonian terms which involve vertices within graph distance x from
the original location of A. If we drop δA from Eq. (79), we can control the resulting error similarly to what
is described in Sec. 6 below:
˙ρ
(x)
C
(t) = −i[H + H
(x)
LS
, ρ
(x)
C
] + ∆
X
−k
∗
<k<k
∗
X
=k∆
A
(x)
ρ
(x)
C
A
(x)†
−
1
2
n
ρ
(x)
C
, A
(x)†
A
(x)
o
. (81)
This master equation has local generators A
(x)
with x = vT
a
/2 + δr. We also assumed that some truncation
method is applied to the Lamb shift H
LS
→ H
(x)
LS
, although we do not present it in this work. One can now
use standard methods [54] to recast Eq. (81)) as a stochastic Schrödinger equation, with those operators
corresponding to jumps:
d
dt
|ψ(t, r)i = G(t, r)|ψ(t)i, ρ
(x)
C
(t) = Av
r
|ψ(t, r)ihψ(t, r)| (82)
here r represents a random number used to generate trajectories. We will not present the explicit form of G,
but note that:
G =
n
X
i=1
G
i
, (83)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 16
Figure 2: Definition of a MPS |ψ(t, r)i = U
Trotter
(r)|ψ(0)i.
where the sum is over qubits and each G
i
is a local operator of radius x = vT
a
/2 + δr and can be computed
in O(1) time in the system size n. The time evolution is then given by a trotterization U
Trotter
(r) of
U(r) = T e
−i
R
t
0
G(τ,r)dτ
. (84)
U
Trotter
(r) is a geometrically local circuit of depth O(t) where each gate takes O(1) time to compute. Con-
traction of such a circuit |ψ(t, r)i = U
Trotter
(r)|ψ(0)i does not necessarily have a local structure, even if the
initial condition |ψ(0)i is a product state. However, note that in 1d |ψ(t, r)i has the structure of a MPS [50]:
|ψ(t, r)i =
X
s
1
...s
n
Y
i
M
s
i
i
|s
1
. . . s
n
i (85)
Where M
s
i
i
are χ × χ matrices for each i, s
i
, and χ is called a bond dimension. A portion of the time
evolution of a MPS is illustrated in Fig. 2 for 2-local G
i
, for which U
Trotter
=
Q
t
(
Q
i=2k
U
i,i+1
Q
i=2k−1
U
i,i+1
)
and U
i,i+1
= e
−iG
i
dt
. Iterating the time step naively grows the bond dimension as χ = e
O(t)
. The compu-
tation of any observable in |ψ(t, r)i requires a contraction of such a tensor network, which takes O(n)e
O(t)
computational time, and e
O(t)
memory. In practice, MPS calculations employ a truncation scheme for the
bond dimension after every time step, which fixes the bond dimension to a constant and leads to O(n, t)
computational time, O(n) memory. Such a method is generally not a controlled approximation, but it may
be possible to develop a special version with an error estimate. Averaging over randomness adds an extra
computational cost to either method. We leave the development of further details to a future publication.
The simultaneous locality and positivity of Eq. (81) is a new result: previously, the Davies generators
were only local for commuting Hamiltonians, and the Redfield master equation leads to completely positive
evolution only for flat spectral densities. Thus, this is the first method for ab initio completely positive open
system evolution of MPSs.
3.7 Scaling with size
Here we discuss the system-size dependence. So far we assumed that there is only one term V = A ⊗ B
describing the interaction with the bath. In reality, there are n terms like this, one for each qubit. If qubits are
coupled to the same bath V =
P
n
i
A
i
⊗B, it will lead to correlated noise (collective decoherence [55]). If the
qubits are coupled to different baths V =
P
n
i
A
i
⊗B
i
, it can lead to uncorrelated noise if Tr[ρ
B
B
i
(t)B
j
] = 0.
We will consider the latter case since it contains the former case if one sets B
i
= B (though distinctly
different phenomena such as decoherence-free subspaces appear in this special case [56, 57]). The form of
the Davies-Lindblad and Redfield equations can be found elsewhere; here we present only the form of the
CGME. The sum is carried through to the frequency form of the equation — Eq. (71):
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
] +
X
ω,ω
0
,i,j
γ
ij
ωω
0
(A
j
ω
ρ
C
A
i
ω
0
−
1
2
{ρ
C
, A
i
ω
0
A
j
ω
}) , (86a)
γ
ij
ωω
0
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
T
a
/2
−T
a
/2
dτ
0
C
ij
(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
, (86b)
H
LS
=
i
2T
a
X
ij
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C
ij
(t
2
− t
1
)A
i
(t
2
)A
j
(t
1
) , (86c)
C
ij
(t) = Tr[ρ
B
B
i
(t)B
j
] . (86d)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 17
The theoretically optimal T
a
=
p
τ
SB
τ
B
/5 is given by the same expression in terms of time scales τ
B
and
τ
SB
, which are now defined using
R
|C| 7→ max
ij
R
|C
ij
|. Such a definition keeps τ
B
and τ
SB
consistent with
the single qubit quantities. The locality of the equation is now also dependent on the decay of C
ij
with the
distance between i and j. If it can be approximated by zero outside some radius (“strong locality”), the
equation remains locally generated in the sense that every term on the r.h.s. contains local terms around
some qubit i conjugating the density matrix, and there are polynomially many of them. If we use a weaker
sense of locality, when the operators on the right and left of the density matrix are allowed to be local around
different i, j, and neglect the Lamb shift, then no assumption on C
ij
is needed.
Let us present an explicitly Lindblad form. Let µ index the positive (see Appendix A) eigenvalues of
γ
ij
(ω) =
R
∞
−∞
e
iωs
C
ij
(s)ds:
γ
ij
(ω) =
X
µ
U
∗
µi
(ω)D
µ
(ω)U
µj
(ω) . (87)
This allows one to write:
γ
ij
ωω
0
=
Z
d
X
µ
U
∗
µi
()˜γ
ωω
0
(D
µ
(), )U
µj
() (88a)
˜γ
ωω
0
(D
µ
(), ) =
1
2πT
a
Z
T
a
/2
−T
a
/2
dt
0
Z
T
a
/2
−T
a
/2
dτ
0
D
µ
()e
−i((τ
0
−t
0
)+ωt
0
+ω
0
τ
0
)
. (88b)
Next, note that if we let
˜
A =
P
j
U
µj
()A
j
, the previous derivation goes through.
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
(t)] +
Z
∞
−∞
d
n
X
µ=1
A
,µ
ρ
C
(t)A
†
,µ
−
1
2
n
ρ
C
(t), A
†
,µ
A
,µ
o
, (89)
where the Lindblad operators are
A
,µ
=
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
X
j
e
it
U
µj
()A
j
(t)dt (90a)
=
X
ω,j
f(, µ, j, ω)A
j
ω
, f(, µ, j, ω) =
s
D
µ
()T
a
2π
sinc[T
a
( − ω)/2]U
µj
() . (90b)
These results generalize Eqs. (19)-(20) to the n-qubit setting. Now the locality is more directly defined: if
the operators A
,µ
can be approximately truncated to a ball around some point for each µ, then the CGME
for n qubits is local in the strong sense described above; otherwise it is local in a weak sense. In any case, it
has n times more generator terms than the single-qubit CGME. The latter has a constant number of terms
after the discretization described in Appendix D.
Let us now discuss the range of applicability of the n-qubit CGME. In the case of correlated noise, the
norm of each A effectively grows by a factor of n. Thus the bound on the relaxation time becomes Λ ∼ n
2
/τ
SB
with τ
SB
taken from the single qubit case. For uncorrelated noise, the bound is Λ ∼ n/τ
SB
. In both cases,
the range of applicability can be derived by directly repeating the derivation given in Sec. 6:
kρ
C
(t) − ρ
true
(t)k
1
≤ O(
p
Λτ
B
e
1.5Λt
) , (91)
where τ
B
is system-independent, and thus is the same as the single qubit quantity. The apparently shrinking
range with n should not discourage us, since it can be mitigated under the assumption of locality. Namely,
one can prove that the error in local density matrices obeys a stronger bound consistent with the isolated
qubit result. Let us note that in an n-qubit system there are in fact nonlocal modes that relax with the rate
n/τ
SB
, or n
2
/τ
SB
for correlated noise (e.g., superradiance), as well as modes whose relaxation is completely
suppressed (e.g., subradiance, or decoherence-free subspaces, as mentioned above). These phenomena are
well known [58] and experimentally documented [59–62].
The benefit expressed by Eq. (91) is that the requirement on the single qubit τ
B
/τ
SB
value switches
from being exponential in the system size for the Davies-Lindblad case to polynomial for the CGME and the
Redfield master equation. In terms of the coupling strength, we found these ranges of applicability at the
relevant times of Λ
−1
to be between
√
ngτ
B
1 and ngτ
B
1 depending on the correlations in the noise.
For local observables the earlier range gτ
B
1 for times ∼ τ
SB
is expected to hold. Here the inequality
notation () guarantees that the error kδρk
1
or kδρ
loc
k
1
is small.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 18
4 The CGME with a time-dependent system Hamiltonian
We start again from the total Hamiltonian (3), but this time we let the system Hamiltonian be time-
dependent: H(t). We denote the interaction term A ⊗ B by V . To go to the interaction picture, we now
need to specify both the start time t
i
and the final time t ≤ t
f
. The unitary propagator U
tt
i
satisfies
dU
tt
i
dt
= −iH(t)U
tt
i
, U
t
i
t
i
= I . (92)
The formal solution is
U
t
f
t
i
= T exp
−i
Z
t
f
t
i
H(t)dt
, (93)
where T represents forward time-ordering. Some basic additional facts are:
U
t
3
t
1
= U
t
3
t
2
U
t
2
t
1
,
dU
tt
i
dt
i
= iU
tt
i
H(t) , U
†
t
1
t
2
= U
t
2
t
1
. (94)
The relation between the Schrödinger picture and the interaction picture is:
ρ(t) = U
t0
ρ
I
(t)U
0t
, V
I
= V (t) = A(t, 0) ⊗B(t) = U
0t
AU
t0
⊗ B(t) . (95)
Here B(t) is given by the bath Hamiltonian in the same way as in the previous derivation. We again assume
the bath state to be stationary for simplicity, although the derivation will go through without it. At this
point we repeat the steps involved in Eqs. (29)-(42b), with the only change being that the A(t) operators
now acquire a second time variable. Namely, Eq. (42b) — the Redfield equation in the interaction picture
— is replaced by
˙ρ
B,I
(t) =
Z
t
0
C(τ − t)(A(t, 0)ρ
B,I
(t)A(τ, 0) − ρ
B,I
(t)A(τ, 0)A(t, 0))dτ + h.c. + E
M
(96a)
=
Z
t
0
dτ [C(τ − t)(A(t, 0)ρ
B,I
(t)A(τ, 0) − ρ
B,I
(t)A(τ, 0)A(t, 0))
+C(t −τ)(A(τ, 0)ρ
B,I
(t)A(t, 0) − A(t, 0)A(τ, 0)ρ
B,I
(t))] + E
M
, (96b)
where we again introduced the Markov approximation ρ
B,I
(τ) 7→ ρ
B,I
(t).
Let us digress briefly to derive the time-dependent Redfield equation. To do so, we transform Eq. (96)
out of the interaction picture:
˙ρ
BM
(t) = −i[H(t), ρ
BM
] +
Z
t
0
C(τ − t)(Aρ
B
(t)A(τ, t) − ρ
B
(t)A(τ, t)A + h.c.)dτ . (97)
Changing variables t − τ = t
0
and switching the dt
0
integration limit to t → ∞ we obtain:
˙ρ
R
(t) = −i[H(t), ρ
R
] + (Aρ
B
(t)A
f
(t) − ρ
B
(t)A
f
(t)Adτ + h.c.), A
f
(t) =
Z
∞
0
C(−t
0
)A(t − t
0
, t)dt
0
, (98)
which is the time-dependent Redfield equation.
We now proceed with essentially the same series of transformations we used in the derivation for the time-
independent Hamiltonian in Sec. 3.1. The main difference is that there are now additional time variables.
Similarly as in the transition to Eq. (56), we time-average the r.h.s. by shifting t 7→ t + t
1
and applying
1
T
a
R
T
a
/2
−T
a
/2
dt
1
. The new equation is:
˙ρ
BMT,I
(t) =
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t+t
1
0
dτ [C(τ − t − t
1
)(A(t + t
1
, 0)ρ
BMT,I
(t)A(τ, 0) − ρ
BMT,I
(t)A(τ, 0)A(t + t
1
, 0))
+C(t + t
1
− τ )(A(τ, 0)ρ
BMT,I
(t)A(t + t
1
, 0) − A(t + t
1
, 0)A(τ, 0)ρ
BMT,I
(t))] ,
(99)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 19
and as long as T
a
is small compared to the fastest timescale over which ρ
I
(t) changes, the associated error
kρ
BM
(t) − ρ
BMT
(t)k
1
is small (see Sec. 6.4 for a detailed analysis of the time-independent case). We now
rotate out of the interaction picture using Eq. (95):
˙ρ
BMT
(t) = −i[H(t), ρ
BMT
(t)]
+
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t+t
1
0
dτ [C(τ − t − t
1
)(A(t + t
1
, t)ρ
BMT
(t)A(τ, t) − ρ
BMT
(t)A(τ, t)A(t + t
1
, t))
+C(t + t
1
− τ )(A(τ, t)ρ
BMT
(t)A(t + t
1
, t) − A(t + t
1
, t)A(τ, t)ρ
BMT
(t))] . (100)
Next, we change variables to t
2
= τ − t:
˙ρ
BMT
(t) = −i[H(t), ρ(t)]
+
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t
1
−t
dt
2
[C(t
2
− t
1
)(A(t + t
1
, t)ρ
BMT
(t)A(t + t
2
, t) − ρ
BMT
(t)A(t + t
2
, t)A(t + t
1
, t))
+C(t
1
− t
2
)(A(t + t
2
, t)ρ
BMT
(t)A(t + t
1
, t) − A(t + t
1
, t)A(t + t
2
, t)ρ
BMT
(t))] . (101)
The transformations we performed after coarse graining of the Redfield equation (96) were exact. We know
what to do to recover complete positivity: as in the transition to Eq. (65a) in the time-independent case, we
need to neglect a part of the integral by changing the lower limit from −t to −T
a
/2:
˙ρ
BMT
(t) = −i[H(t), ρ
BMT
(t)] + P
1
+ P
2
− ρ
BMT
(t)Q
1
− Q
2
ρ
BMT
(t) + E
p
, (102)
where we defined
P
1
≡
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t
1
−T
a
/2
dt
2
C(t
2
− t
1
)A(t + t
1
, t)ρ
BMT
(t)A(t + t
2
, t) (103a)
P
2
≡
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t
1
−T
a
/2
dt
2
C(t
1
− t
2
)A(t + t
2
, t)ρ
BMT
(t)A(t + t
1
, t) (103b)
Q
1
≡
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t
1
−T
a
/2
dt
2
C(t
2
− t
1
)A(t + t
2
, t)A(t + t
1
, t) (103c)
Q
2
≡
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
t
1
−T
a
/2
dt
2
C(t
1
− t
2
)A(t + t
1
, t)A(t + t
2
, t) . (103d)
The same estimates on the error incurred by doing this apply, i.e., this introduces an error of order kE
p
(const·
τ
SB
)k = O(τ
B
/T
a
). Complete positivity yet remains to be demonstrated. To do so, we first note that we
may interchange the order of integration since
R
T
a
/2
−T
a
/2
dt
1
R
t
1
−T
a
/2
dt
2
=
R
T
a
/2
−T
a
/2
dt
2
R
T
a
/2
t
2
dt
1
. We then swap the
integration variables t
1
and t
2
and obtain:
P
2
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
t
1
dt
2
C(t
2
− t
1
)A(t + t
1
, t)ρ
BMT
(t)A(t + t
2
, t) (104a)
Q
2
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
t
1
dt
2
C(t
2
− t
1
)A(t + t
2
, t)ρ
BMT
(t)A(t + t
1
, t) , (104b)
so that:
P
1
+ P
2
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
C(t
2
− t
1
)A(t + t
1
, t)ρ
BMT
(t)A(t + t
2
, t) (105a)
X ≡ Q
1
+ Q
2
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
C(t
2
− t
1
)A(t + t
2
, t)A(t + t
1
, t) (105b)
Y ≡ −i(Q
1
− Q
2
) =
1
iT
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)A(t + t
2
, t)A(t + t
1
, t) , (105c)
and note that X and Y are both Hermitian. Also, Q
1
=
1
2
(X + iY ) and Q
2
=
1
2
(X − iY ), so that
ρ
BMT
Q
1
+ Q
2
ρ
BMT
= −
i
2
[Y, ρ
BMT
] +
1
2
{X, ρ
BMT
} . (106)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 20
Combining Eq. (102) with the expressions in Eqs. (105)-(106), we thus find
˙ρ
C
(t) = −i[H(t) + H
LS
(t), ρ
C
(t)] +
1
T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
C(t
2
− t
1
)[A(t + t
1
, t)ρ
C
(t)A(t + t
2
, t)
−
1
2
{A(t + t
2
, t)A(t + t
1
, t), ρ
C
(t)}] , (107)
where
H
LS
(t) ≡ −
1
2
Y =
i
2T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)A(t + t
2
, t)A(t + t
1
, t) . (108)
Finally, we establish complete positivity by again introducing the Fourier transform C(t) =
1
2π
R
∞
−∞
e
−it
γ()d
[Eq. (6a)], where γ() > 0 (Appendix A). In analogy with Eq. (20) for the time-independent case, let us also
define
A
(t) ≡
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
e
it
1
A(t + t
1
, t)dt
1
. (109)
We then obtain from Eq. (107) a master equation in exactly the same form as the time-independent CGME
[Eq. (19)], except that all the operators now depend on time:
˙ρ
C
(t) = −i[H(t) + H
LS
(t), ρ
C
] +
Z
∞
−∞
d
A
(t)ρ
C
A
†
(t) −
1
2
n
ρ
C
, A
†
(t)A
(t)
o
. (110)
This result, which is in a manifestly completely positive form, establishes that the CGME is valid also for
arbitrary time-dependent system Hamiltonians. The time-dependent CGME, Eq. (110) [also stated earlier
as Eq. (22)], is one of the main new results of this work.
5 Applications and Examples
5.1 Application: dynamical decoupling
One interesting application of our time-dependent CGME is to study methods to reduce decoherence. One
such method is dynamical decoupling (DD) [63, 64] — the use of pulse sequences to cancel out the system-
bath interaction (for reviews see [65, 66]). Even though the master equation has a Markovian appearance,
the underlying bath can have a nonzero correlation time, and can therefore describe the effects of DD. At one
level this might appear to be a counterintuitive result since it is widely accepted that in the limit C(t) = δ(t)
DD is useless, and hence one might naively expect that one needs to use non-Markovian methods (e.g., the
Magnus expansion to second order or higher [67]) to describe the benefits of DD. However, it has recently
been recognized that DD can in fact work in the Markovian setting as well, specifically the Davies-Lindblad
master equation [68], and a more abstract semigroup framework [69] (see also Refs. [70, 71]). Our contribution
here is to establish this in the framework of the Redfield master equation and, more importantly, the CGME.
This ensures a wider range of applicability than the Davies-Lindblad master equation result. Since our result
is obtained in a strictly Markovian setting, we differ with the conclusion reported in Ref. [72], that “success,
however limited, of DD is a meaningful concept of non-Markovianity".
We consider two examples, both for a qubit with H = 0 coupled to a purely dephasing environment:
H
tot
= Z ⊗ B + I ⊗ H
b
, (111)
where Z ≡ σ
z
is the Pauli z-matrix. The noise model (111) can be decoupled with a sequence of X pulses
spaced by ∆t. Here by a pulse we mean an instantaneous unitary being applied to the qubit, through the
action of an externally controlled, time-dependent qubit Hamiltonian: instead of H = 0 we have a time-
dependent system Hamiltonian H(t) =
π
2
P
j
δ(t − j∆t)X, where the prefactor of π/2 is chosen so that
X-pulses are generated.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 21
5.1.1 Example 1: time-dependent Redfield master equation with a rectangle bath correlation function
Consider the rectangle bath correlation function
C(t) = g
2
θ(τ
c
− |t|) , (112)
where we have explicitly introduced the coupling strength g. To analyze this problem it will suffice to use
the time-dependent version of the Redfield equation, Eq. (98). We thus need to compute
A
f
(t) =
Z
∞
0
C(t
0
)A(t − t
0
, t)dt
0
, (113)
and note that the time-evolution of the operator A = Z under X pulses results in A(t −t
0
, t) = ±Z with the
sign alternating every ∆t:
A(t − t
0
, t) = (−1)
dt/∆te−d(t−t
0
)/∆te
Z , (114)
Taking the integral, we find:
A
f
(t) = g
2
Z
Z
τ
c
0
(−1)
dt/∆te−d(t−t
0
)/∆te
dt
0
= g
2
c(t)Z, |c(t)| ≤ ∆t , (115)
where c(t) is the remainder after the cancellation of different signs. To show why |c(t)| ≤ ∆t, note that there
is no way to accumulate > ∆t by integrating any interval of the function (−1)
dx/∆te
∀ a, b :
Z
b
a
(−1)
dx/∆te
dx
≤ ∆t (116)
To see this, we note that a translation a, b → a + ∆t, b + ∆t changes the sign of the integral. We then
split the integration interval [a, b] into three subsets: {[a + 2n∆t, a + (2n + 1)∆t]}, {[a + (2n + 1)∆t, a +
2(n + 1)∆t]}, [a + 2∆t(n
max
+ 1), b] where n
max
has been chosen closest to b: n
max
= d(b − a)/2∆te − 2.
The integrals over the first two subsets cancel due to the translation property, and the remaining subset
[a + 2∆t(n
max
+ 1), b] has length L < 2∆t. Note that for L = 0 and L = 2∆t the integral vanishes. Also note
that the derivative of the integral with respect to b is (−1)
db/∆te
, and its absolute value is 1. Any function
with such a derivative, up to two inflection points and zeros at the end of the interval [0, 2L] stays within
[−∆t, +∆t], and the bound is saturated when b = m∆t.
We compare this with the case without pulses:
A
f
(t) = g
2
τ
c
Z , (117)
Here τ
c
enters the correlation function [see Eq. (112)]. Thus the decoherence rate is reduced by a factor of
at least ∆t/τ
c
. For (gτ
c
)
2
1 our analysis is accurate with the error bounds discussed above.
However, we note that the bath correlation function C(t) = g
2
θ(τ
c
− |t|) used here for simplicity is not
physically permissible since its Fourier transform γ(ω) < 0 for some ω. This leads to possibly nonpositive
c(t) for some t, and then the equation does not describe a CP map. To address this we next perform a more
careful calculation with a realistic bath.
5.1.2 Example 2: time-dependent CGME with an Ohmic spectral density
Now consider a physically permissible environment with a spectral density that satisfies the KMS condi-
tion (10). We specify γ(ω) below. We wish to use the time-dependent CGME [Eq. (22)] to study the
same dynamical decoupling protocol as in the previous example. We first need to compute A(t + t
1
, t) =
U
†
(t + t
1
, t)ZU(t + t
1
, t), where U(t + t
1
, t) = exp[−i
R
t+t
1
t
H(s)ds]. This unitary is a product of as many
X-pulses as fit in the interval [t, t + t
1
], which we can express as
U(t + t
1
, t) = X
m
f
−m
i
, m
i
= dt/∆te , m
f
= d(t + t
1
)/∆te , (118)
where m
i
and m
f
−1 respectively count the number of pulses in the intervals [0, t] and [0, t+t
1
]. Conjugating
Z by these pulses results in alternating signs, depending on the parity of m
f
− m
i
:
A(t + t
1
, t) = (−1)
∆m
Z , ∆m = m
f
− m
i
. (119)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 22
This yields the time-dependent Lindblad operators A
(t) via Eq. (23):
A
DD
(t) = f
DD
(t)Z , f
DD
=
s
γ()T
a
2π
Z
1/2
−1/2
e
iζT
a
(−1)
d(t+ζT
a
)/∆te−dt/∆te
dζ . (120)
Note that the Lamb shift [Eq. (24)] is proportional to Z
2
= , so it drops out. The Lindblad operators
without DD are simply
A
noDD
(t) = f
noDD
(t)Z , f
noDD
(t) =
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
e
it
1
dt
1
=
s
γ()T
a
2π
sinc(T
a
/2) . (121)
We may thus write the time-dependent CGME as:
˙ρ
C
(t) = −i[H(t), ρ
C
] + r
x
(t) (Zρ
C
Z − ρ
C
) , (122)
where
r
x
(t) =
Z
∞
−∞
d |f
x
|
2
, (123)
with x = DD or x = noDD. The decoherence rate is thus suppressed by the factor
ξ(t) ≡
r
DD
(t)
r
noDD
(t)
=
R
∞
−∞
dγ()
R
1/2
−1/2
e
iζT
a
(−1)
d(t+ζT
a
)/∆te−dt/∆te
dζ
2
R
∞
−∞
dγ()
sinc
T
a
2
2
. (124)
The inner integral can be evaluated under certain simplifying assumptions. For example, let us assume that
T
a
is an even integer multiple of ∆t: T
a
= 2k∆t, and let us consider only times t which are integer multiples
of ∆t: t = `∆t. Then d(t + ζT
a
)/∆te − dt/∆te = d2kζe, and the dependence on t cancels out, i.e., the
suppression factor ξ(t) is time-independent. We can then rewrite the inner integral as:
Z
1/2
−1/2
e
iζT
a
(−1)
d(t+ζT
a
)/∆te−dt/∆te
dζ
2
=
k−1
X
j=0
(−1)
j
Z
−
1
2
+
j+1
k
−
1
2
+
j
k
e
iζT
a
dζ
2
(T
a
= 2k∆t, t = `∆t) (125a)
=
k−1
X
j=0
(−1)
j
e
iT
a
[2(j−k)+1]/(2k)
k
sinc
T
a
2k
2
(125b)
=
2
T
a
sin
kπ
2
+
T
a
2
tan
T
a
2k
2
(125c)
=
sinc
T
a
2
tan
T
a
2k
2
(k even) , (125d)
where in the last line we assumed for further simplicity assume that k is even (k = 2k
0
), so that sin
kπ
2
+
T
a
2
=
±sin
T
a
2
. Writing t = `∆t and T
a
= 4k
0
∆t, we thus have:
ξ(`∆t) =
R
∞
−∞
dωγ(ω) |sinc (2k
0
ω∆t) tan (ω∆t)|
2
R
∞
−∞
dωγ(ω) |sinc (2k
0
ω∆t)|
2
(k
0
∈ N) . (126)
To guarantee that DD will cause suppression [ξ(`∆t) < 1] it is sufficient for the spectral density to have a
high-frequency cutoff ω
c
such that |tan (ω∆t) | < 1, i.e.,
ω
c
∆t <
π
4
, (127)
in agreement with well established results [63].
An important physical example is a bath with an Ohmic spectral density
γ(ω) = 2πκ
ωe
−|ω|/ω
c
1 − e
−βω
, (128)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 23
-6 -5 -4 -3 -2 -1
log
10
Δt
-14
-12
-10
-8
-6
-4
-2
log
10
ξ
(a) β = 0.2
-6 -5 -4 -3 -2 -1
log
10
Δt
-14
-12
-10
-8
-6
-4
-2
log
10
ξ
(b) β = 5
ω
c
=0.125×
π
4
ω
c
=0.25×
π
4
ω
c
=0.5×
π
4
ω
c
=1×
π
4
ω
c
=2×
π
4
Figure 3: The DD suppression factor ξ [Eq. (126)] as a function of the pulse interval ∆t for different high-frequency cutoffs
ω
c
, at two different inverse temperatures (a) β = 0.2 and (b) β = 5. We have set T
a
= 4∆t and κ = 1. The sufficient
condition (127) is satisfied everywhere, except for the case ω
c
= 2 × π/4 (purple line) for ∆t . 1, where indeed for the
β = 5 case the suppression factor is slightly > 1.
where ω
c
is a high-frequency cutoff and κ is a positive dimensionless constant. This spectral density satisfies
the KMS condition. The associated bath correlation function can be expressed in terms of the Polygamma
function (see Appendix I of Ref. [31]), but unfortunately evaluating τ
SB
and τ
B
analytically using Eq. (11)
is not possible in this case, so we simply pick a convenient value of T
a
(= 4k
0
∆t) rather than its optimal
value
p
τ
B
τ
SB
/5. The suppression factor is plotted in Fig. 3 for a variety of parameter settings, including for
∆t > π/(4ω
c
). It can be seen that the CGME correctly and consistently predicts that DD will result in the
suppression of dephasing when Eq. (127) is satisfied. Thus, the CGME can be used to study DD protocols
despite its Markovian appearance.
5.2 Lambless master equations
It turns out that the Lamb shift terms containing sgn(t), such as Eq. (18b), are more demanding to compute
numerically, but at the same time the interesting dynamics is often given by the other, relaxation terms.
Therefore it is convenient to have a version of the equations in the “Lambless" regime. Even though these
equations do not have any correctness guarantees, one can still use them as toy models, and their numerical
solution is often significantly faster than their complete counterparts. Let us list the resulting equations:
• Lambless Davies-Lindblad [replacing Eq. (18)]:
˙ρ
LD
(t) = −i[H, ρ] +
1
2
X
ω
γ(ω)(A
ω
ρ
LD
A
−ω
− ρ
LD
A
−ω
A
ω
) + h.c. (129)
• Lambless Redfield [replacing Eq. (16)]:
˙ρ
LR
(t) = −i[H, ρ
LR
(t)] + (Aρ
LR
A
f
− ρ
LR
A
f
A + h.c.) , (130)
where
A
f
=
1
2
X
ω
γ(−ω)A
ω
. (131)
This simplified form of A
f
follows from Eqs. (7) and (17) by substituting S(ω) = 0 (vanishing Lamb
shift). When the KMS condition holds we may write A
f
=
1
2
P
ω
e
−βω
γ(ω)A
ω
.
• Lambless CGME [replacing Eq. (19)]:
˙ρ
LC
(t) = −i[H, ρ
LC
] + ∆
X
k,=∆k,|k|<k
∗
(A
ρ
LC
A
†
−
1
2
n
ρ
LC
, A
†
A
o
) (132)
where
A
=
X
ω
A
ω
s
γ()T
a
2π
sinc [T
a
( − ω)/2] , (133)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 24
and we choose T
a,opt
=
p
τ
B
τ
SB
/5, ∆, and k
∗
as prescribed in Appendix D. We study this case
numerically in Sec. 5.4.
The importance of the Lamb shift in open system dynamics has been analyzed before, e.g., in Refs. [31, 73]
5.3 Toy example of a slow bath at the boundary of the range of applicability, with explicit correlation
function and spectral density
Another way to speed up numerical simulations is to have a spectral density that leads to an explicit,
analytically computable form of the bath correlation function C(t). Note that in the “Lambless" CGME case
discussed in Sec. 5.2 the only direct role played by C(t) is in the determination of τ
B
and τ
SB
[via Eq. (11)],
and apart from this it suffices to specify only γ(ω). The complete equations are more complicated: the Lamb
terms contain C(t), or a Cauchy principal value integral of γ(ω) [Eq. (8)].
We present such a toy example with the property of having a bath almost at the boundary of the slowest
correlation functions allowed within the range of applicability of the CGME. Baths with the same τ
B
can have
a different time decay of |C(t)| at large times t τ
B
. A bath with |C(t)| ∼ 1/t
2
will have a logarithmically
divergent τ
B
, which is undesirable for a toy example.
10
Consider a bath with the following spectral density,
which satisfies the KMS condition, such that |C(t)| ∼ 1/t
4
—almost as slow as we are allowed:
γ(ω) = N
e
βω/2
τ
SB
(e
−bβ|ω|
− a
−1
e
−abβ|ω|
) , (134)
where a > 1, b > 1/2. Since
1
τ
SB
=
R
∞
0
|
1
2π
R
∞
−∞
e
−iωt
γ(ω)dω|dt [Eqs. (6b), (11a)], the normalization factor is
N
−1
=
1
2π
Z
∞
0
Z
∞
−∞
e
−iωt+βω/2
(e
−b|βω|
− a
−1
e
−ab|βω|
)dω
dt . (135)
The difference between the two exponentials in Eq. (134) is such that for small ω the term linear in |ω|
cancels, and there is no inflection in γ(ω). Note that for a single e
−β|ω|
we would have |C(t)| ∼ 1/t
2
, but
after this cancellation we get |C(t)| ∼ 1/t
4
.
Our goal is now to use this example to compare the actual numerical errors to the theoretical bounds
derived in this work. These bounds involve many relaxations and inequalities, so we do not expect them to
be particularly tight. Also of interest is to compare the predicted optimal averaging time T
a,theory opt
to the
numerically optimal T
a
. To be specific we set a = 1.01, b = 0.6, β = 4, and find the normalization factor to
be N ≈ 21.0. This choice of parameters gives rise to a maximum of γ(ω) at ω
∗
≈ 2.0, and the peak is wide:
for 0.15 ≤ ω ≤ 6.08, γ(ω) ≥ 0.5γ(ω
∗
). This frequency range is where we will choose system transitions to
lie; see Fig. 4(a).
For this choice of parameters we obtain:
C(t) =
325620N
τ
SB
π(22i + 5t)(553i + 125t)(−106 −515it + 625t
2
)
, (136)
so that indeed C(t) ∼ 1/t
4
. Integrals of C(t)e
iωt
over various intervals can be expressed via the Meijer G-
function. Its evaluation is a part of the numerical simulation of master equations discussed in this paper, so
their computational cost is sensitive to any shortcuts we can find, and the analytic form above is helpful. We
will comment on the specific runtimes below.
The bath parameter τ
SB
can be chosen freely, while τ
B
≈ 0.69 is determined by the choice of a, b, β. Note
that it has dimensions of β.
We choose τ
SB
= 10 so that τ
B
/τ
SB
= 0.1 1, a necessary condition for our error bounds to be small.
We find T
a,theory opt
=
p
τ
B
τ
SB
/5 = 0.97. To compare different equations, we pick an arbitrary two-qubit
Hamiltonian:
H = 0.5σ
z
1
− 0.7σ
z
2
+ 0.3σ
z
1
σ
z
2
+ σ
x
1
+ σ
x
2
, (137)
and choose |11i as the initial state. We let the coupling to the bath be given by V = σ
z
1
⊗ B.
10
|C(t)| ∼ 1/t
2
is actually the borderline case realized by an Ohmic bath, for which the divergence has a cutoff depending on
the total time of the experiment T in Eq. (11b).
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 25
(a) (b)
Figure 4: (a) The spectral density γ(ω) for a = 1.01, b = 0.6, β = 4 and τ
SB
= 10.0. The position of the transition
frequencies of the Hamiltonian defined by Eq. (137) is shown by the red tickmarks. (b) Populations, i.e., the probabilities of
eigenstates (labeled by energy in the figure) p
n
(t) for the solution ρ
OR
(t) in the region of its positivity: 0 ≤ t ≤ 2.56τ
SB
.
For all of our results, there will be no way to know what the true state ρ
true
(t) is. The differences found
will be within the error bounds (growing at least linearly with time) of all equations. However we can simulate
Eq. (46), which is presumably closest to the true solution since only the Born and Markov approximations
were made:
˙ρ
OR
(t) = −i[H, ρ
OR
(t)] +
Z
t
0
C(τ − t)[Aρ
OR
(t)A(τ − t) − ρ
OR
(t)A(τ − t)A]dτ + h.c. (138)
We call this the “Original Redfield" equation (ORE). It is in fact possible that the CGME returns solutions
(a) (b)
Figure 5: (a) Histogram of the norm of the right-hand side kL
BM,I
t
(ρ
test
)k
1
sampled uniformly over time and over normalized
matrices from the Gaussian unitary ensemble. The upper bound 4/τ
SB
is shown as a dashed green line. Dotted: the
numerically computed maximum of the norm. Solid vertical line: typical value of the norm. Dash-dotted: a value of the
norm that would have explained the optimum, for comparison. (b) Thick: the average error
1
T
R
T
0
kρ
C
(t) − ρ
OR
(t)k
1
dt.
Thick dotted: same for Davies
1
T
R
T
0
kρ
D
(t) − ρ
OR
(t)k
1
dt. Dashed: position of T
a,theory opt
, which does not lead to an
advantage over Davies, in contrast to the following two. Dotted: T
a,adj
. Dash-Dotted: T
a,num. opt
. Solid vertical line: T
a
adjusted by the typical value of the norm.
closer to the true solution than this non-CP equation, but we note that in our derivation all the equations
appeared as subsequent approximations to the ORE. Therefore we may characterize how good those approx-
imations are by studying the difference ρ
C
− ρ
OR
. We first investigate the solution ρ
OR
(t), and find that it
becomes negative (and hence unphysical) at t ≈ 2.57τ
SB
. Thus, we choose to compare the solutions of the
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 26
CGME and ORE in the interval 0 ≤ t ≤ 2.56τ
SB
, where ρ
OR
≥ 0 [see Fig. 4(b), and hence we can assume
that c
BM
= 1 [recall that this is always true for a CP master equation; see the discussion following Eq. (47b)].
We note that the computational cost of the ORE is the largest (by at least on order of magnitude) among
our equations.
11
The second observation we make is about the norm of the r.h.s. of Eq. (138). In Fig. 5(a) we present
a probability distribution of the norms of the right hand side kL
BM,I
t
(ρ
test
)k
1
, where ρ
test
is taken from the
Gaussian unitary ensemble (Hermitian matrices X sampled from a probability distribution ∼ e
−2
n−1
TrX
2
,
where n = 2 is the number of qubits), and subsequently normalized such that kρ
test
k
1
= 1. We also combine
the data for different times, as if the time was sampled uniformly in [0, 2.56τ
SB
]. We find numerically that
the maximum value of the norm that can be achieved is 0.44 ± 0.01. It would suggest a higher effective τ
SB
than our definition:
max
t,ρ
test
,kρ
test
k
1
=1
kL
BM,I
t
(ρ
test
)k
1
=
4
τ
eff
SB
≈ 0.44 < 0.58 ≈
4
τ
SB
. (139)
Such an adjusted τ
eff
SB
≈ 9.14 = 1.33τ
SB
can be used to choose an adjusted T
a,adj
=
q
τ
eff
SB
/τ
SB
T
a,theory opt
=
1.12 = 1.15T
a,theory opt
. We could non-rigorously use the typical value of the norm instead, leading to
T
a,typ.
≈ 1.7T
a,theory opt
. To see if this adjustment is justified, we vary T
a
in Fig. 5(b) and plot the average of
the trace norm of the difference in the solutions over t ∈ [0, 2.56τ
SB
], i.e.,
1
T
R
T
0
kρ
C
(t) − ρ
OR
(t)k
1
dt. We see
that the true optimum is at even higher T
a,num. opt
≈ 2.87 ≈ 3T
a,theory opt
, and the equation matches Davies
in the limit T
a
→ ∞ [Fig. 5(b)], as formally expected (see Appendix C). Even though our adjustments of T
a
are not very large, the error is very sensitive to them and improves significantly.
(a) (b)
Figure 6: (a) Thick: the error kρ
C
(t) − ρ
OR
(t)k
1
for T
a,adj
. Thin: same for Davies kρ
D
(t) − ρ
OR
(t)k
1
. Dashed: The
change in the solution kρ
OR,I
(t) − ρ
OR,I
(0)k
1
induced by relaxation (our error should be small w.r.t. it). Dotted: the
tightest upper bound derived in this paper [Eq. (240)], bounding the thick line. Dash-Dotted: the error kρ
C
(t) − ρ
OR
(t)k
1
for T
a,num. opt
. (b) The average error
1
T
R
T
0
kρ
C
(t) − ρ
OR
(t)k
1
dt as a function of varying τ
SB
, holding τ
B
fixed, for T
a,adj
(thick) and T
a,num. opt
(dash-dotted). The thin line shows the average error
1
T
R
T
0
kρ
D
(t) − ρ
OR
(t)k
1
dt for Davies.
11
Precomputing the analytic form of the r.h.s. takes about 14 minutes on a contemporary desktop computer, and the solution
of the differential equation itself — about a minute. In contrast, the CGME takes 20 seconds total, and the Davies and Redfield
equations — even less. These numbers are for the implementation in Mathematica, where precomputing the symbolic form of
the r.h.s. is possible. Faster implementations in other programming languages are certainly possible. The Mathematica code
for obtaining the plots of this section are available [74].
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 27
We illustrate the error in more detail in Fig. 6(a) by plotting the differences
kρ
C,adj
(t) − ρ
OR
(t)k
1
, kρ
C,num. opt
(t) − ρ
OR
(t)k
1
, kρ
D
(t) − ρ
OR
(t)k
1
, kρ
OR,I
(t) − ρ
OR,I
(0)k
1
, (140)
as a function of t, where the CGME is taken for T
a,adj
and T
a,num. opt
. Also in Fig. 6(a) we plot our strongest
upper bound on kρ
C,adj
(t)−ρ
OR
(t)k
1
[Eq. (240), derived in Sec. 6 below]. We use T
a
→ T
a,adj
and Λ → 4/τ
eff
SB
Comparing this bound to the observed error in Fig. 6(a), we see that even after this adjustment, the bound is
still not tight. Nonetheless, the adjustment helps to match the optimum in Fig. 5(b) better, and we believe
the bounds can be improved by future work in this direction.
In Fig. 6(b) we vary τ
SB
, and evaluate
1
T
Z
T
0
kρ
C,adj
(t) − ρ
OR
(t)k
1
dt,
1
T
Z
T
0
kρ
C,num. opt
(t) − ρ
OR
(t)k
1
dt,
1
T
Z
T
0
kρ
D
(t) − ρ
OR
(t)k
1
dt , (141)
for T = 2.56τ
SB
. To find T
a,num. opt
, we vary T
a
for each value of τ
SB
and choose the one that gives the
smallest values for the corresponding expression in Eq. (141). We again add the similar average error of Davies
to the plot. Since our family of CGMEs with different T
a
includes Davies as a T
a
→ ∞ limit, numerical
optimization will always give an improvement over Davies. A simpler adjusted T
a
that does not require
optimization still outperforms Davies at moderate τ
B
/τ
SB
. We also note that the
p
τ
B
/τ
SB
dependence
of its error at moderate τ
B
/τ
SB
is consistent with the intuition from our bounds, while the linear portion
suggests that our bounds are not tight in some parameter regime.
5.4 Role of the Lamb shift in the CGME
We perform another test, where we drop the Lamb shift H
LS
of the CGME and obtain a corresponding
solution ρ
CL
(t), as prescribed by Eqs. (132) and (133). We then plot kρ
C,adj
(t) − ρ
CL,adj
(t)k
1
as a function
of time, as shown in Fig. 7. We again pick T
a,adj
= 1.15T
a,opt
. We observe in Fig. 7 that the Lamb shift has
the same role as in the Davies equation: it affects the phases of the off-diagonal matrix elements of ρ
C
(t) in
the eigenbasis. Since off-diagonal elements decay to zero for the Davies case, the effect of their phase on the
norm difference first grows proportional to the phase, then disappears with the magnitude of those elements.
Figure 7: Thin: the error kρ
C
(t) − ρ
CL
(t)k
1
of dropping the Lamb shift H
LS
for T
a,adj
. Thick: the maximum difference in
the magnitude of the matrix elements max
ij
δ|ρ
C,ij
|. Dashed: the maximum difference in the phase of the matrix elements.
Dotted: the decay of the off-diagonal matrix elements.
6 Error bounds and estimates
In this section we give the details of the various error bounds and estimates mentioned without proof in our
derivation of the CGME in Sec. 3. Our strategy is to derive the error terms successively, and then add them
up using the triangle inequality. We work our way up through the various approximations, starting from
Born in Sec. 6.2, followed by Markov in Sec. 6.3, then time-averaging (or coarse-graining) in Sec. 6.4. In the
course of the latter, we bound the error due to dropping part of the integration domain to achieve complete
positivity, and analyze the optimization of T
a
. But we start with a proof of Lemma 1, which plays a central
role in our analysis.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 28
6.1 Bound on the difference of solutions of differential equations differing by a bounded operator
We first present a proof of Lemma 1.
Proof. The first part of the Lemma is concerned with the pair of differential equations ˙x = L
t
x + E and
˙y = L
t
y. Their formal solution is:
x(t) = x(0) +
Z
t
0
L
τ
(x(τ))dτ +
Z
t
0
Edτ, y(t) = y(0) +
Z
t
0
L
τ
(y(τ))dτ . (142)
Let Λ > 0 be a constant defined as follows: sup
τ,x
kL
τ
(x)k
1
≤ Λ for all x such that kxk
1
= 1. Since L
is linear, and writing x = δ/kδk
1
, it follows that ∀δ: kL
τ
(δ)k
1
= kδk
1
kL
τ
(x)k
1
≤ kδk
1
Λ. Then, assuming
x(0) = y(0), and taking the operator norm of the difference yields, for δ = x(t) −y(t):
kx(t) − y(t)k
1
≤ Λ
Z
t
0
kx(τ) −y(τ)k
1
dτ + tkEk
1
, (143)
Many different functions f (t) ≡ kx(t) −y(t)k
1
satisfy this inequality for all t, but they are all upper-bounded
by the function b
1
(t) = max kx(t) −y(t)k
1
that saturates the inequality:
˙
b
1
= Λb
1
+ kEk
1
, b
1
(0) = 0 . (144)
The solution is
b
1
(t) = (e
Λt
− 1)
kEk
1
Λ
(145a)
≤ (e
c
− 1)
kEk
1
Λ
∀t ≤
c
Λ
. (145b)
This proves the desired bound on kx(t) −y(t)k
1
≤ b
1
(t).
The second part of the Lemma is concerned with the more general, non-Markovian pair of differential
equations ˙x(t) =
R
t
0
K
t−τ
(x(τ))dτ + E and ˙y(t) =
R
t
0
K
t−τ
(y(τ))dτ . 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:
x(t) − y(t) =
Z
t
0
Z
τ
0
K
τ−θ
(x(θ) − y(θ))dτdθ +
Z
t
0
Edτ . (146)
Define Λ
t
= sup
x
kK
t
(x)k
1
for all x such that kxk
1
= 1. By linearity of K
t
it follows as above that ∀δ:
kK
t
(δ)k
1
≤ Λ
t
kδk
1
. Therefore:
kx(t) − y(t)k
1
≤
Z
t
0
Z
τ
0
Λ
τ−θ
kx(θ) − y(θ)k
1
dτdθ +
Z
t
0
kEk
1
dτ . (147)
The upper bound b
K
(t) on kx(t) −y(t)k
1
satisfies the equation:
b
K
(t) =
Z
t
0
Z
τ
0
Λ
τ−θ
b
K
(θ)dτdθ + kEk
1
t . (148)
Taking the derivative once, we obtain:
˙
b
K
(t) =
Z
t
0
Λ
t−τ
b
K
(τ)dτ + kEk
1
. (149)
Since the r.h.s. is positive, b
K
(t) is a monotonic function. So b
K
(τ) ≤ b
K
(t), and it follows that:
˙
b
K
(t) ≤
Z
t
0
Λ
t−τ
b
K
(t)dτ + kEk
1
. (150)
In other words, we can again upper-bound b
K
(t) ≤ b
0
K
(t), such that
˙
b
0
K
(t) =
Z
t
0
Λ
t−τ
dτb
0
K
(t) + kEk
1
. (151)
Now let Λ be a positive constant such that
R
t
0
Λ
τ
dτ ≤ Λ for all t. We thus arrive at
˙
b
K
00
(t) = Λb
00
K
(t) + kEk
1
, (152)
for a new upper-bound b
00
K
(t) ≥ b
0
K
(t). This is exactly the same as Eq. (144), and hence the Lemma still
holds for non-Markovian equations.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 29
Note that for our purposes Λ = 4/τ
SB
is usually taken, thanks to Eq. (41, 42c), and that the non-
Markovian bound will be generalized in Sec. 6.2. The exponential on the r.h.s. of Eq. (145b) may be
troublesome. Under additional assumptions we can prove a much tighter linear in t bound presented in
Appendix E, yet we did not find a way to recast all of our results in that tighter form.
6.2 Born approximation error
We start with an estimate of the error associated with making the Born approximation, which we presented
as the estimate kE
B
k
1
= O
τ
B
τ
2
SB
in Eq. (33). I.e., we would like to bound the error E
B
= ˙ρ
true,I
(t) −
R
t
0
K
2,B
t−τ
(ρ
true,I
(τ))dτ. It is difficult to do so directly since ρ
true,I
(t) is unknown, so instead we will first use
a proxy for ρ
true,I
(t), which we call ρ
B4,I
(t). The latter is the solution of the master equation obtained by
iterating the substitution of ρ
tot
= ρ
tot
(0) − i
R
t
0
[V (τ), ρ
tot
(τ)]dτ [Eq. (29b)] into ˙ρ
tot
(t) = −i[V (t), ρ
tot
(t)]
[Eq. (29a)] to 4th order, and only then applying the Born approximation. After developing an understanding
of the 4th order we proceed to the infinite series. We make certain assumptions about its convergence, and
comment below on their compatibility with the known convergence radius of the Dyson series. Let us first
define the type of bath for which we can carry out this procedure in order to arrive at a bound.
We refer to a bath as exponential if its correlation function 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 C(t
1
−t
3
)C(t
2
−t
4
)]. If C(t) is exponentially decaying with a characteristic time τ
B
,
that immediately lets us conclude that only t
1
, t
2
, t
3
, t
4
that are all within ∼ τ
B
of each other contribute to
the estimate, and the bound on the Born error acquires a small factor of τ
B
/τ
SB
. However, an exponentially
decaying 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
R
∞
0
|C(t)|t
n
dt for some integer
n ≥ 0. Such a condition was used in Ref. [31] to control the error of the subsequent Markov approximation.
Here we show how these integrals for n = 0, 1 control our estimate of the Born error for a Gaussian bath. The
physical meaning of these integrals was defined in Eq. (11). A more general requirement on a non-Gaussian
bath can be derived in the same way. However, that requirement will involve a higher-order correlation
function that does not have a straightforward and concise interpretation.
Recall the steps of the derivation where we make the Born approximation, leading to Eq. (39). For
brevity, we drop the interaction picture subscript, replace the time argument by a subscript, and place the
approximation level B in the superscript: ρ
B,I
(t) = ρ
B
t
. Written out explicitly, Eq. (39) in terms of ρ
B
t
is:
˙ρ
B
t
=
Z
t
0
K
2,B
t−τ
(ρ
B
τ
)dτ =
Z
t
0
C(τ − t)(A
t
ρ
B
τ
A
τ
− ρ
B
τ
A
τ
A
t
) + C(t −τ)(A
τ
ρ
B
τ
A
t
− A
t
A
τ
ρ
B
τ
)dτ . (153)
The formal solution is:
ρ
B
τ
= ρ
0
+
Z
τ
0
dθ
Z
θ
0
dmC(m −θ)(A
θ
ρ
B
m
A
m
− ρ
B
m
A
m
A
θ
) + C(θ − m)(A
m
ρ
B
m
A
θ
− A
θ
A
m
ρ
B
m
). (154)
We substitute this solution back into Eq. (153), a necessary step because our estimate of the Born error
will come from the 4th order terms, so we need to subtract the original solution. The full expression to be
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 30
subtracted is:
˙ρ
B
t
=
Z
t
0
K
2,B
t−τ
Z
τ
0
Z
θ
0
K
2,B
θ−m
(ρ
B
m
)dmdθ
!
dτ =
Z
t
0
K
4,B
t−m
(ρ
B
m
)dτ (155a)
=
Z
t
0
C(τ − t)(A
t
ρ
0
A
τ
− ρ
0
A
τ
A
t
) + C(t −τ)(A
τ
ρ
0
A
t
− A
t
A
τ
ρ
0
)dτ (155b)
+
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dm (155c)
[C(τ − t)C(m −θ)A
t
(A
θ
ρ
B
m
A
m
− ρ
B
m
A
m
A
θ
)A
τ
+ C(τ − t)C(θ − m)A
t
(A
m
ρ
B
m
A
θ
− A
θ
A
m
ρ
B
m
)A
τ
(155d)
−C(τ − t)C(m − θ)(A
θ
ρ
B
m
A
m
− ρ
B
m
A
m
A
θ
)A
τ
A
t
− C(τ − t)C(θ − m)(A
m
ρ
B
m
A
θ
− A
θ
A
m
ρ
B
m
)A
τ
A
t
(155e)
+C(t −τ)C(m − θ)A
τ
(A
θ
ρ
B
m
A
m
− ρ
B
m
A
m
A
θ
)A
t
+ C(t − τ)C(θ −m)A
τ
(A
m
ρ
B
m
A
θ
− A
θ
A
m
ρ
B
m
)A
t
(155f)
−C(t −τ)C(m − θ)A
t
A
τ
(A
θ
ρ
B
m
A
m
− ρ
B
m
A
m
A
θ
) − C(t −τ)C(θ − m)A
t
A
τ
(A
m
ρ
B
m
A
θ
− A
θ
A
m
ρ
B
m
)] . (155g)
Now, these terms constitute only a part of the expression for the derivative of the true density matrix
˙ρ
true
t
to the 4th order in the system-bath interaction. The full 4th order expression can be obtained, as
mentioned above, by substituting the solution of the total evolution, ρ
tot
= ρ
tot
(0) − i
R
t
0
[V (τ), ρ
tot
(τ)]dτ
into ˙ρ
tot
(t) = −i[V (t), ρ
tot
(t)] three consecutive times (to generate a 4th order commutator), and only then
applying the Born approximation. We call this ρ
B4
t
since it will differ from ρ
B
t
in the 4th order in the
interaction V (t). Doing so, we obtain:
˙ρ
B4
t
=
Z
t
0
C(τ − t)(A
t
ρ
0
A
τ
− ρ
0
A
τ
A
t
) + C(t −τ)(A
τ
ρ
0
A
t
− A
t
A
τ
ρ
0
)dτ (156a)
+ Tr
b
[A
t
⊗ B
t
,
Z
t
0
dτ[A
τ
⊗ B
τ
,
Z
τ
0
dθ[A
θ
⊗ B
θ
,
Z
θ
0
dm[A
m
⊗ B
m
, ρ
B4
m
⊗ ρ
b
]]]] . (156b)
Defined in this way, ρ
B4
t
is close to the 4th order of the true state ρ
true
t
up to terms of 6th order and higher,
so we use ρ
B4
t
− ρ
B
t
to estimate the actual leading order error ρ
true
t
− ρ
B
t
. We proceed to expand all the
commutators:
˙ρ
B4
t
=
Z
t
0
K
4,B4
t−τ
(ρ
B4
τ
)dτ =
Z
t
0
K
2,B
t−τ
(ρ
0
)dτ (157a)
+
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dmTr
b
B
t
B
τ
B
θ
B
m
ρ
b
A
t
A
τ
A
θ
A
m
ρ
B4
m
− Tr
b
B
t
B
τ
B
θ
ρ
b
B
m
A
t
A
τ
A
θ
ρ
B4
m
A
m
− . . . (157b)
There are 16 terms total in the last line, and the order of B
t
, B
τ
, B
θ
, B
m
and ρ
b
is exactly the same as the
order of A
t
, A
τ
, A
θ
, A
m
and ρ
B
m
.
Let us now assume that the bath is Gaussian. By definition, a Gaussian bath obeys Wick’s theorem (or
Isserlis’ theorem [75]), which states that at all orders the higher order correlation functions decouple into
products of two-point correlation functions. For example, for the 4-point correlation function:
Tr
b
B
t
B
τ
B
θ
B
m
ρ
b
= Tr
b
B
t
B
τ
ρ
b
Tr
b
B
θ
B
m
ρ
b
+ Tr
b
B
t
B
θ
ρ
b
Tr
b
B
τ
B
m
ρ
b
+ Tr
b
B
t
B
m
ρ
b
Tr
b
B
τ
B
θ
ρ
b
. (158)
By definition, C(t −τ) = Tr
b
B
t
B
τ
ρ
b
. Thus:
Tr
b
B
t
B
τ
B
θ
B
m
ρ
b
= C(t −τ)C(θ −m) + C(t −θ)C(τ − m) + C(t − m)C(τ − θ) (159)
We see that in Eq. (155) only the first term was present. One can check that the order of A is exactly the
same, so upon subtracting Eq. (155) from Eq. (158) only the two terms C(t −θ)C(τ −m) + C(t −m)C(τ −θ)
remain. These extra terms can be interpreted in terms of the corresponding diagrammatic technique. Indeed,
note that t > τ > θ > m. If in Wick’s theorem we have a pairing t ↔ τ, θ ↔ 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 ↔ θ, τ ↔ m or t ↔ m, τ ↔ θ, then the two arcs are stuck together due to the condition
t > τ > θ > m; see Fig. 8.
At the moment we have two equations: the Born equation in two forms Eq. (153) and Eq. (155), as well
as the 4th order equation Eq. (157a):
˙ρ
B
t
=
Z
t
0
K
2,B
t−τ
(ρ
B
τ
)dτ =
Z
t
0
K
4,B
t−τ
(ρ
B
τ
)dτ, ˙ρ
B4
t
=
Z
t
0
K
4,B4
t−τ
(ρ
B4
τ
)dτ, (160)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 31
Figure 8: Different possible pairings of Tr
b
B
t
B
τ
B
θ
B
m
ρ
b
in Wick’s theorem. Note that the first and second order are in
τ
B
/τ
SB
, not in the coupling strength itself, in which both diagrams are 4th order.
where K(ρ) are corresponding linear superoperators. K
2,B
is the 2nd order (in the bath interaction) operator
from Eq. (153), K
4,B
and K
4,B4
have terms up to the 4th order defined in Eq. (155) and Eq. (157a). Define
E
B4
(ρ
any
) =
Z
t
0
(K
4,B4
t−τ
(ρ
any
τ
) − K
4,B
t−τ
(ρ
any
τ
))dτ , (161)
for instance:
E
B4
(ρ
B4
) = ˙ρ
B4
t
−
Z
t
0
K
4,B
t−τ
(ρ
B4
τ
)dτ . (162)
Define
E
B
= ˙ρ
true
t
−
Z
t
0
K
2,B
t−τ
(ρ
true
τ
)dτ . (163)
The quantity that appeared in the derivation is E
B
. We will now bound kE
B4
(ρ
any
)k
1
and discuss its gener-
alization kE
Bk
(ρ
any
)k
1
. We postpone the estimate for kE
B
k
1
till we develop more advanced tools.
The difference K
4,B4
t−τ
(ρ
4,B4
τ
) −K
4,B
t−τ
(ρ
B4
τ
) is given by the two second order pairings in Fig. 8, times the 16
possible orders of A
t
, A
τ
, A
θ
, A
m
and ρ
m
:
E
B4
=
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dm(C(t −θ)C(τ − m) + C(t −m)C(τ − θ))A
t
A
τ
A
θ
A
m
ρ
any
m
− . . . (164)
We take the trace norm of both sides:
kE
B4
(ρ
any
)k
1
≤
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dmk(C(t −θ)C(τ − m) + C(t −m)C(τ − θ))A
t
A
τ
A
θ
A
m
ρ
any
m
k
1
+ . . . (165)
where we used the triangle inequality kA + Bk
1
≤ kAk
1
+ kBk
1
. We now use submultiplicativity [Eq. (38)].
Noting that the first norm in the submultiplicativity inequality is the operator norm, we obtain:
kE
B4
(ρ
any
)k
1
≤
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dmk(C(t −θ)C(τ − m) + C(t − m)C(τ − θ))A
t
A
τ
A
θ
A
m
kkρ
any
m
k
1
+ . . . (166)
The operator norm kAk = 1 and we will set c
a
= max
t
kρ
any
t
k
1
such that kρ
any
k
1
≤ c
a
:
kE
B4
(ρ
any
)k
1
≤ c
a
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dm|(C(t −θ)C(τ − m) + C(t −m)C(τ − θ))| + . . . (167)
In particular, E
B4
(ρ
B
) is bounded with c
B
and E
B4
(ρ
B4
) with c
B4
. There are 16 terms after the expansion
of commutators:
kE
B4
(ρ
any
)k ≤ 16c
a
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dm|C(t −θ)||C(τ − m)| + |C(t −m)||C(τ − θ)| , (168)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 32
where we have also used that |C(t −m)| = |C(m −t)
∗
| = |C(m −t)|. What is left is to bound the integrals in
Eq. (164). It turns out that the following condition suffices:
Z
∞
0
|C(t)|t
n
dt ≤
τ
n
B
τ
SB
, for n = 0, 1 , (169)
which is automatically satisfied for τ
SB
and τ
B
as defined in Eq. (11). We also change the order of integration
variables using the condition 0 ≤ m ≤ θ ≤ τ ≤ t, e.g.:
Z
t
0
dτ
Z
τ
0
dθ =
Z
t
0
dθ
Z
t
θ
dτ. (170)
The first term in Eq. (164) is bounded as:
Z
t
0
dτ
Z
τ
0
dθ
"
Z
θ
0
dm|C(τ − m)|
#
|C(t −θ)| ≤
1
τ
SB
Z
t
0
dτ
Z
τ
0
dθ|C(t −θ)| =
1
τ
SB
Z
t
0
dθ
Z
t
θ
dτ|C(t −θ)| (171a)
=
1
τ
SB
Z
t
0
(t − θ)dθ|C(t −θ)| ≤
τ
B
τ
2
SB
, (171b)
while the second term in Eq. (164) is bounded as:
Z
t
0
dτ
Z
τ
0
dθ
Z
θ
0
dm|C(t −m)||C(τ − θ)| =
Z
t
0
dm
Z
t
m
dθ
Z
t
θ
dτ|C(τ − θ)|
|C(t −m)| (172a)
≤
1
τ
SB
Z
t
0
dm
Z
t
m
dθ|C(t −m)| =
1
τ
SB
Z
t
0
(t − m)dm|C(t −m)|
(172b)
≤
τ
B
τ
2
SB
. (172c)
Together, this yields:
kE
B4
(ρ
any
)k
1
≤ 32c
a
τ
B
τ
2
SB
, (173)
which is the first hint for the result given in Eq. (33). The actual error term appearing in Eq. (32) is E
B
as
defined in Eq. (163). We can consider repeating the construction of ρ
B4
presented here for ρ
B6
, ρ
B8
. . . as
successive approximations to ρ
true
. The corresponding errors can be defined:
E
Bk
(ρ
any
) =
Z
t
0
(K
k,Bk
t−τ
(ρ
any
τ
) − K
k,B
t−τ
(ρ
any
τ
))dτ , (174)
By employing the diagrammatic technique as in Fig. 8, the higher orders of perturbation theory in the
system-bath coupling can be shown to give higher orders in τ
B
/τ
SB
. There are also extra factors of 4t/τ
SB
that appear in the higher orders. For instance, E
B6
has the form:
kE
B6
(ρ
any
)k
1
≤ 32c
a
τ
B
τ
2
SB
1 +
8t
τ
SB
+ O
τ
B
τ
SB
, (175)
and the diagrams for the second term are shown in Fig. 9(a). We do not prove this form of E
B6
here. We
will bound all t-dependent terms of the first order in τ
B
/τ
SB
, by summing the diagrams in Fig. 9(b):
kE
Bk
(ρ
any
)k
1
≤ 32c
a
τ
B
τ
2
SB
k/2−1
X
m=1
m
(m − 1)!
4t
τ
SB
m−1
+ O
τ
B
τ
SB
(176a)
≤ 32c
a
τ
B
τ
2
SB
(1 + 4t/τ
SB
)e
4t/τ
SB
+ O
τ
B
τ
SB
, (176b)
for all k. Here we fix time to be a constant: ∃(t), C(t) s.t. O
τ
B
τ
SB
≤ C(t)
τ
B
τ
SB
for all
τ
B
τ
SB
≤ (t). This is
a weak statement. A stronger statement would include the time-dependence inside big-O, and have time-
independent constants , C:
kE
Bk
(ρ
any
)k
1
≤ 32c
a
e
4t/τ
SB
τ
B
τ
2
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) (177)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 33
Figure 9: Diagrams contributing to the second order of (a) E
B6
, (b) E
Bk
for arbitrary k.
We do not know how to prove such a statement as of yet. The proof would rely on the diagrammatic
calculation of subleading terms. It is straightforward, but cumbersome, and we do not attempt it here.
Instead we make two assumptions:
1. In the k → ∞ limit, there is a finite radius of convergence in e
4t/τ
SB
τ
B
/τ
SB
:
kE
B∞
(ρ
any
)k
1
≤ 32c
a
e
4t/τ
SB
τ
B
τ
2
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) . (178)
2. The true solution is the limit of the perturbative series:
˙ρ
true
t
=
Z
t
0
K
∞,B∞
t−τ
(ρ
true
τ
)dτ . (179)
To prove convergence, one needs to control the factorial number of terms appearing in Wick’s theorem.
Such a proof was already presented by Davies [4]. The slight adjustments needed to include an integral
condition such as Eq. (169) are a somewhat tedious problem left for a future study.
Our definition of E
B
under this assumption contains E
B∞
discussed above, as well as an extra term:
E
B
= ˙ρ
true
t
−
Z
t
0
K
2,B
t−τ
(ρ
true
τ
)dτ =
Z
t
0
(K
∞,B∞
t−τ
(ρ
true
τ
) − K
2,B
t−τ
(ρ
true
τ
))dτ (180a)
kE
B
k
1
≤
Z
t
0
(K
∞,B∞
t−τ
(ρ
true
τ
) − K
∞,B
t−τ
(ρ
true
τ
))dτ
1
+
Z
t
0
(K
∞,B
t−τ
(ρ
true
τ
) − K
2,B
t−τ
(ρ
true
τ
))dτ
1
(180b)
≤ 32
τ
B
τ
2
SB
(1 + O(τ
B
/τ
SB
)) +
Z
t
0
(K
∞,B
t−τ
(ρ
true
τ
) − K
2,B
t−τ
(ρ
true
τ
))dτ
1
, (180c)
where c
a
= 1 since the true evolution is CP. The second term will turn out to be time-dependent, and quite
tricky to estimate. We will return to this problem after we develop some extra machinery.
Let us first note the form of K
k,B
, starting with k = 4:
Z
t
0
K
4,B
t−τ
(ρ
any
τ
)dτ =
Z
t
0
K
2,B
t−τ
ρ
0
+
Z
τ
0
Z
m
0
K
2,B
m−θ
(ρ
any
m
)dmdθ
dτ . (181)
The solution of
˙ρ
B4
=
Z
t
0
K
4,B4
t−τ
(ρ
B4
τ
)dτ =
Z
t
0
K
4,B
t−τ
(ρ
B4
τ
)dτ + E
B4
(182)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 34
is equivalent to solving a system of 4th order equations that are our proxy for the true solution (which is ∞
order):
(
˙ρ
B4
=
R
t
0
K
2,B
t−τ
(ρ
(1)
τ
)dτ + E
B4
, ρ
B4
(0) = ρ
0
˙ρ
(1)
=
R
t
0
K
2,B
t−τ
(ρ
B4
τ
)dτ, ρ
(1)
(0) = ρ
0
.
(183)
We will prove that the solution of the above system remains close to ρ
B
, the solution of the equation after
the Born approximation:
˙ρ
B
=
Z
t
0
K
2,B
t−τ
(ρ
B
τ
)dτ . (184)
We will also present a straightforward generalization to all k which will allow one to bound kE
B
k
1
.
First write the equation in the compact form:
˙
p = Lp + E, p =
ρ
B4
ρ
(1)
, L =
0
R
K
t−τ
dτ
R
K
t−τ
dτ 0
, E =
E
B4
(ρ
B4
)
0
. (185)
Define a custom norm that is an infinity norm on the vector of operators, but for each operator the trace
norm is taken:
kvk
c
=
v
1
v
2
c
= max(kv
1
k
1
, kv
2
k
1
) (186)
The triangle inequality is obeyed as required:
kv + wk
c
= max(kv
1
+ w
1
k
1
, kv
2
+ w
2
k
1
) ≤ max(kv
1
k
1
+ kw
1
k
1
, kv
2
k
1
+ kw
2
k
1
) (187a)
≤ max(kv
1
k
1
, kv
2
k
1
) + max(kw
1
k
1
, kw
2
k
1
) = kvk
c
+ kwk
c
. (187b)
The maximum of the custom norm of L over inputs of norm 1 is the same as before:
max
p:kpk
c
=1
kLpk
c
= Λ =
4
τ
SB
, (188)
so all the properties needed to prove Lemma 1 hold. We establish:
kδp
t
k
c
≤
e
Λt
− 1
Λ
kEk
c
(189)
in the original variables:
max(kδρ
(1)
t
k
1
, kδρ
B4
t
k
1
) ≤
e
Λt
− 1
Λ
kE
B4
k
1
. (190)
Note that the error we are trying to bound is
kE
B
k
1
=
˙ρ
true
−
Z
t
0
K
2,B
t−τ
(ρ
true
τ
)dτ
1
=
Z
t
0
K
∞,B∞
t−τ
(ρ
true
τ
) − K
2,B
t−τ
(ρ
true
τ
dτ
1
. (191)
Our approximation to it is bounded as:
˙ρ
B4
−
Z
t
0
K
2,B
t−τ
(ρ
B4
τ
)dτ
1
≤
Z
t
0
K
2,B
t−τ
(ρ
(1)
τ
− ρ
B4
τ
)dτ
1
+ kE
B4
k
1
≤ (2e
Λt
− 1)kE
B4
k
1
, (192)
where we have used the triangle inequality: kρ
(1)
τ
− ρ
B4
τ
k
1
≤ kδρ
(1)
τ
k
1
+ kδρ
B4
τ
k
1
. Now we can repeat the
above proof for arbitrary k:
˙ρ
Bk
=
R
t
0
K
2,B
t−τ
(ρ
(1)
τ
)dτ + E
Bk
(ρ
Bk
) ρ
Bk
(0) = ρ
0
˙ρ
(1)
=
R
t
0
K
2,B
t−τ
(ρ
(2)
τ
)dτ ρ
(1)
(0) = ρ
0
. . . . . .
˙ρ
(k/2−1)
=
R
t
0
K
2,B
t−τ
(ρ
B4
τ
)dτ ρ
(k/2−1)
(0) = ρ
0
(193)
For appropriately defined L, k · k
c
the r.h.s. is still bounded by the same constant:
max
p:kpk
c
=1
kLpk
c
= Λ =
4
τ
SB
, (194)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 35
which leads to
max(kδρ
(1)
t
k
1
, . . . , kδρ
Bk
t
k
1
) ≤
e
Λt
− 1
Λ
kE
Bk
k
1
, (195)
and the error in the r.h.s. is bounded without change:
˙ρ
Bk
−
Z
t
0
K
2,B
t−τ
(ρ
Bk
τ
)dτ
1
≤
Z
t
0
K
2,B
t−τ
(ρ
(1)
τ
− ρ
Bk
τ
)dτ
1
+ kE
Bk
k
1
≤ (2e
Λt
− 1)kE
Bk
k
1
. (196)
Taking the k → ∞ limit, we obtain:
kρ
true
t
− ρ
B
t
k
1
≤
e
Λt
− 1
Λ
kE
B∞
(ρ
true
)k
1
≤ 8(e
4t/τ
SB
− 1)e
4t/τ
SB
τ
B
τ
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) , (197)
and, finally:
kE
B
k
1
≤ (2e
Λt
− 1)kE
B∞
(ρ
true
)k
1
(198a)
≤ 32(2e
4t/τ
SB
− 1)e
4t/τ
SB
τ
B
τ
2
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) , (198b)
where we have used that c
a
= 1 for ρ
true
. This result is the sought-after justification of the estimate
kE
B
k
1
= O
τ
B
τ
2
SB
given in Eq. (33).
An immediate consequence of Eq. (197) is that we can use it to obtain a bound on c
B
, which we defined
after Eq. (40) as an upper bound on kρ
B,I
k
1
. Since
kρ
B
t
k
1
≤ kρ
true
t
k
1
+ kρ
B
t
− ρ
true
t
k
1
= 1 + kρ
true
t
− ρ
B
t
k
1
, (199)
we obtain:
c
B
≤ 1 + 8(e
4t/τ
SB
− 1)e
4t/τ
SB
τ
B
τ
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) . (200)
6.3 Markov approximation error
This subsection follows Ref. [76] initially. Before the Markov approximation but after the Born approximation,
the master equation had the form given in Eq. (39). The error due the Markov approximation ρ
B,I
(τ) 7→
ρ
B,I
(t), in transitioning to Eq. (42b), is:
E
M
=
Z
t
0
C(τ − t)[A(t)(ρ
B,I
(τ) −ρ
B,I
(t))A(τ) −(ρ
B,I
(τ) −ρ
B,I
(t))A(τ)A(t)]dτ + h.c. (201)
Since ρ
B,I
(t) − ρ
B,I
(τ) =
R
t
τ
˙ρ
B,I
(t
0
)dt
0
, we have
kρ
B,I
(t) − ρ
B,I
(τ)k
1
≤ k˙ρ
B,I
k
1
(t − τ) ≤
4(t − τ)c
B
τ
SB
, (202)
where we used Eq. (40) in the second inequality. Now we can estimate the norm of the error:
kE
M
k
1
≤ 4
Z
t
0
|C(τ − t)|
4(t − τ)
τ
SB
dτ ≤ 16
τ
B
c
B
τ
2
SB
. (203)
We found c
B
and saw that it is O(1) in Eq. (200).
Let us now use the Markovian Lemma 1 with Eq. (203):
kρ
BM
t
− ρ
B
t
k
1
≤ 4
τ
B
c
B
τ
SB
e
4t/τ
SB
− 1
. (204)
A total expression for the Born and Markov error is [adding Eq. (197)]:
kρ
true
t
−ρ
BM
t
k
1
≤ 4
τ
B
c
B
τ
SB
e
4t/τ
SB
− 1
+ 8(e
4t/τ
SB
−1)e
4t/τ
SB
τ
B
τ
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
)) , (205)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 36
which simplifies to:
kρ
true
t
− ρ
BM
t
k
1
≤ 4(e
4t/τ
SB
− 1)
τ
B
τ
SB
c
B
+ 2e
4t/τ
SB
(1 + 4t/τ
SB
+ O(e
4t/τ
SB
τ
B
/τ
SB
))
, (206)
or substituting in c
B
:
kρ
true
t
− ρ
BM
t
k
1
≤
e
4t
τ
SB
− 1
τ
B
4 + 8
1 +
4
e
4t
τ
SB
−1
τ
B
τ
SB
e
4t
τ
SB
1 +
4t
τ
SB
+ O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
. (207)
The extra factor from c
B
can be absorbed into the big-O notation:
kρ
true
t
− ρ
BM
t
k
1
≤
e
4t
τ
SB
− 1
τ
B
4 + 8e
4t
τ
SB
1 +
4t
τ
SB
1 + O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
. (208)
A loose bound follows from the above:
kρ
true
t
− ρ
BM
t
k
1
≤
12
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
1 + O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
. (209)
Sometimes we wish to report the bound briefly as:
kρ
true
t
− ρ
BM
t
k
1
=
e
4t
τ
SB
− 1
O
e
8t
τ
SB
τ
B
τ
SB
(210)
Proof. The original big-O notation states that for e
4t
τ
SB
τ
B
τ
SB
≤ the following inequality holds:
kρ
true
t
− ρ
BM
t
k
1
≤
12
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
1 + Ce
4t
τ
SB
τ
B
τ
SB
τ
SB
≤
12
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
(1 + C)
τ
SB
. (211)
We first note that e
4t
τ
SB
≤ e
8t
τ
SB
, so if e
8t
τ
SB
τ
B
τ
SB
≤ , then e
4t
τ
SB
τ
B
τ
SB
≤ as well, and we are within the range
of the original big-O notation. By taking 12(1 + C) as the new constant, we prove Eq. (210).
6.4 Time averaging error and its optimization
6.4.1 General analysis of time-averaging
We consider the time-averaging operation, which leads to a more complicated bound on the difference of the
solutions. We first prove a general result:
Lemma 2. For any first order linear differential equation ˙ρ(t) = L
t
ρ(t) there is a time-averaged version
˙π(t) = L
t
π(t), L
t
=
1
T
a
Z
T
a
/2
−T
a
/2
L
t+τ
dτ . (212)
Note that we use a different symbol π(t) for the solution of the time-averaged equation. Assume the same
initial condition: ρ(0) = π(0). Then the deviation between the two, δρ(t) ≡ π(t) −ρ(t), is bounded as:
kδρ(t)k
1
≤ [b
2
(min(t, T
a
/2)) + c
ρ
ΛT
a
/4] max(e
Λ(t−T
a
/2)
, 1) − c
ρ
ΛT
a
/4 , (213)
where Λ, c
ρ
are constants such that ∀t, ρ
0
, kρ
0
k
1
≤ 1 : kL
t
(ρ
0
)k
1
≤ Λ and the solution ρ(t) with the initial
condition ρ(0) = ρ
0
is bounded as kρ(t)k
1
≤ c
ρ
. The function b
2
(t) is:
b
2
(t) ≡ 4c
ρ
t
T
a
+ c
ρ
4 − 2ΛT
a
− (ΛT
a
/2)
2
ΛT
a
(1 − e
Λt
) . (214)
An immediate corollary is:
kδρ(t)k
1
= O(ΛT
a
) for t = O(Λ
−1
) . (215)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 37
The proof can be found in Appendix F. The O(ΛT
a
) corollary enters the estimate for kρ
BM,I
(t) −
ρ
BMT,I
(t)k
1
we presented in Eq. (61b) in the derivation (using Λ = 4/τ
SB
). We will also need a compact
bound on:
b
2
(T
a
/2)
c
ρ
= e
ΛT
a
/2
(ΛT
a
/2)
2
+ 4(e
−ΛT
a
/2
− 1 + ΛT
a
/2)
ΛT
a
−
ΛT
a
4
. (216)
One can prove that the expression in the brackets is upper-bounded by
e
−ΛT
a
/2
− 1 + ΛT
a
/2 ≤ (ΛT
a
)
2
/8 . (217)
This allows us to bound
b
2
(T
a
/2) ≤ c
ρ
e
ΛT
a
/2
3ΛT
a
4
−
ΛT
a
4
. (218)
Note that the function b
2
(t) is monotonic: b
2
(t) ≤ b
2
(T
a
/2). One can then loosen the bound in Lemma 2:
kδρ(t)k
1
≤
c
ρ
e
ΛT
a
/2
3ΛT
a
4
−
ΛT
a
4
, t ≤ T
a
/2
c
ρ
e
Λt
3ΛT
a
4
−
ΛT
a
4
, t > T
a
/2
= c
ρ
e
Λ max(t,T
a
/2)
3ΛT
a
4
−
ΛT
a
4
. (219)
6.4.2 Time averaging of the Redfield ME
We now focus on the specific equations of interest to us. We already established in Eq. (42c) that we can
choose Λ = 4/τ
SB
(recall that Λ needs to upper bound the r.h.s.: max
t,ρ
test
,kρ
test
k
1
=1
kL
BM,I
t
(ρ
test
)k
1
≤ Λ).
The total time-averaging error (219) is now bounded as:
kρ
BM
− ρ
BMT
k
1
= kδρ(t)k
1
≤ c
BM
e
Λ max(t,T
a
/2)
3T
a
τ
SB
−
T
a
τ
SB
, (220)
We will use both the compact bound of Eq. (220), and the tightest bound available which is a repetition of
Lemma 2:
kδρ(t)k
1
≤ b
2
(t) = [b
2
(min(t, T
a
/2)) + c
BM
ΛT
a
/4] max(e
Λ(t−T
a
/2)
, 1) − c
BM
ΛT
a
/4 , (221)
b
2
(t) = 4c
BM
t
T
a
+ c
BM
4 − 2ΛT
a
− (ΛT
a
/2)
2
ΛT
a
(1 − e
Λt
), t < T
a
/2 . (222)
6.4.3 Bound on the error due to enforcing complete positivity
Recall that after the time-averaged equation is obtained, we need to drop an extra portion of the integral to
obtain complete positivity [the transition to Eq. (65a)]. Before that, Eq. (64) was:
˙ρ
BMT
(t) = −i[H, ρ
BMT
] +
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t
0
−t
dτ
0
C(τ
0
−t
0
)e
−i(ωt
0
+ω
0
τ
0
)
(A
ω
ρ
BMT
A
ω
0
−ρ
BMT
A
ω
0
A
ω
) + h.c. ,
(223)
which changes into Eq. (65a) (same terms, but the portion of the integral is subtracted):
˙ρ
C
(t) = −i[H, ρ
C
] +
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
t
0
−t
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
(A
ω
ρ
C
A
ω
0
− ρ
C
A
ω
0
A
ω
) + h.c. − E
p
(224a)
E
p
=
X
ω,ω
0
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
−T
a
/2
−t
dτ
0
C(τ
0
− t
0
)e
−i(ωt
0
+ω
0
τ
0
)
(A
ω
ρ
C
A
ω
0
− ρ
C
A
ω
0
A
ω
) + h.c. (224b)
Collecting
P
ω
into A(t) we obtain:
E
p
=
1
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
−T
a
/2
−t
dτ
0
C(τ
0
− t
0
)(A(t
0
)ρ
C
A(τ
0
) − ρ
C
A(τ
0
)A(t
0
)) + h.c. , (225)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 38
Figure 10: Times contributing to E
p
for (a) t > T
a
/2 (b) t < T
a
/2
so that:
kE
p
k
1
≤
4
T
a
Z
T
a
/2
−T
a
/2
dt
0
Z
−T
a
/2
−t
dτ
0
|C(τ
0
− t
0
)|
. (226)
Since the resulting equation is CP, we used kρ
C
k
1
= 1. We would like to bound this norm by splitting the
integration domain into different times. The corresponding domains are illustrated in Fig. 10. For t < T
a
/2
we note that the argument of C(τ
0
−t
0
) can be both positive and negative. Thus a factor of 2 arises when we
bound
R
dt
0
|C(τ
0
− t
0
)| ≤ 2/τ
SB
:
kE
p
k
1
≤
8
τ
SB
T
a
/2 − t
T
a
≤
4
τ
SB
. (227)
For t > T
a
/2, we perform a change of variables: τ
0
− t
0
= −τ, t
0
+ T
a
/2 = θ:
kE
p
k
1
≤
4
T
a
Z
−T
a
/2
−T
a
/2
dt
0
Z
−T
a
/2
−t
dτ
0
|C(τ
0
− t
0
)| =
4
T
a
Z
T
a
0
dθ
Z
t+θ−T
a
/2
θ
dτ|C(τ )| . (228)
We now change the order of integration, noting that 0 < θ < τ, thus obtaining the following upper bound:
kE
p
k
1
≤
4
T
a
Z
T
a
0
dθ
Z
t+θ−T
a
/2
θ
dτ|C(τ )| ≤
4
T
a
Z
t+T
a
/2
0
τ|C(τ )|dτ . (229)
Recall that the total evolution time T entered the definition
τ
B
τ
SB
=
R
T
0
t|C(t)|dt [Eq. (11b)]. We can usually set
T = ∞, except for an Ohmic bath [see Eq. (128)], which is a borderline applicability case, with a logarithmic
in T divergence of the integral. We define an evolution time with a collar of T
a
/2: T = max
t
(t + T
a
/2):
kE
p
k
1
≤
4
T
a
Z
T
0
τ|C(τ )|dτ =
4τ
B
T
a
τ
SB
. (230)
Collecting everything, we obtain:
kE
p
k
1
(t) ≤
4
τ
SB
θ(t − T
a
/2)
τ
B
T
a
+ θ(T
a
/2 − t)
, (231)
where θ(x) = 1, x > 0 and zero otherwise.
Recall that E
p
is defined in Eq. (65a) as the extra term appearing if the CGME is written with the same
expression as Eq. (64) in the r.h.s. We can now use the same idea as in the proof of Lemma 1 (Sec. 6.1). Let
us introduce δρ
2
= ρ
BMT
− ρ
C
. A function b
p
that upper-bounds kδρ
2
k
1
is the solution of the differential
equation
˙
b
p
= Λb
p
+ kE
p
k
1
(t), b
p
(0) = 0 , (232)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 39
where we have used Λ ≥ max
t,ρ
test
,kρ
test
k
1
=1
kL
BM,I
t
(ρ
test
)k
1
from Eq. (42c) as the upper bound on the norm
of the superoperator on the r.h.s. of Eq. (64). Indeed, by applying the triangle inequality to the definition
of the r.h.s. of Eq. (64) we get:
L
BMT,I
t
(ρ
test
) =
1
T
a
Z
T
a
/2
−T
a
/2
L
BM,I
t+τ
(ρ
test
)dτ, kL
BMT,I
t
(ρ
test
)k
1
≤ kL
BM,I
t
(ρ
test
)k
1
≤ Λ . (233)
The solution to Eq. (232) is:
b
p
(t) = e
Λt
Z
t
0
e
−Λτ
kE
p
k
1
(τ)dτ . (234)
Substituting Eq. (231), we obtain two cases:
b
p
(t) =
4
Λτ
SB
e
Λt
1 − e
−Λt
, t < T
a
/2 (235)
b
p
(t) =
4
Λτ
SB
e
Λt
(1 − e
−ΛT
a
/2
) +
τ
B
T
a
(e
−ΛT
a
/2
− e
−Λt
)
, t > T
a
/2 , (236)
which we can summarize as
b
p
(t) ≤
4
Λτ
SB
[(e
Λt
− max(1, e
Λ(t−T
a
/2)
)) +
τ
B
T
a
max(e
Λ(t−T
a
/2)
− 1, 0)] . (237)
For a concise form we change the first term using 1 − e
−x
≤ x:
b
p
(t) ≤ e
Λt
2T
a
τ
SB
+
4τ
B
ΛT
a
τ
SB
max(e
Λ(t−T
a
/2)
− 1, 0) . (238)
Let us now focus on the difference with ρ
BM,I
(t) which is the solution of Eq. (42b). We will keep
Λ ≥ max
t,ρ
test
,kρ
test
k
1
=1
kL
BM,I
t
(ρ
test
)k
1
in the bound in case it is < 4/τ
SB
. Using Eqs. (220) and (238), we
have:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ c
BM
e
Λ max(t,T
a
/2)
3T
a
τ
SB
−
T
a
τ
SB
+ e
Λt
2T
a
τ
SB
+
4τ
B
ΛT
a
τ
SB
max(e
Λ(t−T
a
/2)
− 1, 0) .
(239)
The tightest bound available that starts at zero at t = 0 is the combination of Eqs. (221) and (237):
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ [b
2
(min(t, T
a
/2)) + c
BM
ΛT
a
/4] max(e
Λ(t−T
a
/2)
, 1) − c
BM
ΛT
a
/4+ (240a)
+
4
Λτ
SB
[(e
Λt
− max(1, e
Λ(t−T
a
/2)
)) +
τ
B
T
a
max(e
Λ(t−T
a
/2)
− 1, 0)] . (240b)
where for t < T
a
/2
b
2
(t) = 4c
BM
t
T
a
+ c
BM
4 − 2ΛT
a
− (ΛT
a
/2)
2
ΛT
a
(1 − e
Λt
) (241)
We may optimize T
a
to minimize the bound Eq. (239), which we proceed to do next.
6.4.4 Optimization of T
a
We have used Λ ≤ 4/τ
SB
to obtain the form of the bound in Eq. (239), but we kept both Λ and τ
SB
to
demonstrate the possibility of making the bound tighter by using different values for Λ and 4/τ
SB
evaluated
independently. To that end, we denote c
Λ
= 4/Λτ
SB
≥ 1 and treat it as a constant below. At large
t Λ
−1
, t > T
a
/2 the leading terms of Eq. (239) are:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ c
BM
e
Λt
3T
a
τ
SB
+ e
Λt
2T
a
τ
SB
+ c
Λ
e
Λ(t−T
a
/2)
τ
B
T
a
+ o(e
Λt
) , (242)
We can loosen the bound again via e
−ΛT
a
/2
≤ 1:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ e
Λt
(3c
BM
+ 2)T
a
τ
SB
+ c
Λ
τ
B
T
a
+ o(e
Λt
) . (243)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 40
Minimizing the expression in the brackets, we find that T
a
=
p
c
Λ
τ
SB
τ
B
/(3c
BM
+ 2) [=
p
τ
SB
τ
B
/5 if c
BM
=
c
Λ
= 1, which is Eq. (25)] for any τ
B
/τ
SB
leads to an optimal bound at large enough times. The asymptotic
behavior of the bound is:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ 2e
Λt
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
+ o(e
Λt
) . (244)
One may be concerned that the relaxation e
−ΛT
a
/2
≤ 1 reduced the quality of our bound. We checked this
for c
BM
= c
Λ
= 1, though we do not present that calculation here. The result is that after the optimization
of the original Eq. (240), which is a more complicated function of T
a
, the bound is unchanged to the leading
order, and in fact acquires a multiplicative correction of the form
1 −
17
15
√
5
√
τ
B
√
τ
SB
+ O(τ
B
/τ
SB
))
. We do
not pursue this route since the effect is small, and the calculation is unreasonably long. Instead we use
T
a
=
p
c
Λ
τ
SB
τ
B
/(3c
BM
+ 2) obtained above and find that the full expression for the error from Eq. (239)
can be bounded as
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ 2
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
e
Λt+1
− 3/5
, (245)
and it holds for all positive τ
B
/τ
SB
, t/τ
SB
and c
BM
, c
Λ
≥ 1. The derivation can be found in Appendix G.
This simple form allows us to finally bound c
BM
. Since for a CP coarse-grained equation kρ
C
k
1
= 1:
kρ
BM,I
(t)k
1
≤ 1 + 2
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
e
Λt+1
− 3/5
, c
BM
= 1 + 2
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
e
Λt+1
− 3/5
(246)
This is a quadratic equation on c
BM
. Its ≥ 1 solution is:
c
BM
= 1 +
1
2
p
20X + 9X
2
+ 3X
= 1 + O
r
τ
B
τ
SB
e
Λt
, X = 4
c
Λ
τ
B
τ
SB
e
Λt+1
− 3/5
2
(247)
And here the big-O notation means that there exist constants x
0
, C such that for 0 ≤ x ≤ x
0
, O(x) ≤ Cx. We
note that the choice of T
a
=
p
c
Λ
τ
SB
τ
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
q
τ
B
τ
SB
e
Λt
1 we
can approximate c
BM
≈ 1 and consequently T
a
=
p
c
Λ
τ
SB
τ
B
/5. Note also that our construction was optimal
for asymptotic e
Λt
1 (for the discussion of small times, see Appendix H). We combine the two conditions
as 1 e
Λt
q
τ
SB
τ
B
. In units of actual time, this corresponds to
1 Λt ln
r
τ
SB
τ
B
. (248)
This is a fairly narrow range for any realistic parameters, but it still becomes [1, ∞] in the weak coupling
limit. We choose to work with T
a
=
p
c
Λ
τ
SB
τ
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 ≤ 10%, which would correspond to a range like this. We
have used T
a,adj
=
p
c
Λ
τ
SB
τ
B
/5 and T
a,theory
=
p
τ
SB
τ
B
/5 for the numerical simulations in Section 5.3. For
theoretical purposes of the remainder of this Section, we can always use c
Λ
= 1 which is proven. We will
keep the notation e
Λt
= e
4t/τ
SB
for brevity. The error bound for T
a
=
p
τ
SB
τ
B
/5 is derived as Eq. (367):
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤
r
τ
B
5τ
SB
3e
Λt+1
− 1
c
BM
+
7e
Λt+1
− 5
. (249)
The constant c
BM
can be obtained numerically independently, and then one just needs to substitute it. For
the theoretical use of the CGME, that route is not available, and instead one needs to substitute its upper
bound Eq. (247). For the rest of this section, we will loosen the resulting bound to make the expression more
compact. Since we optimized for e
Λt
1, dropping constants will not be significant in the parameter regime
of interest. We will also update the expression for c
BM
using
√
a + b ≤
√
a +
√
b in Eq. (247):
c
BM
≤ 1 +
√
5X + 3X,
√
X ≤ 2
r
τ
B
τ
SB
e
Λt+1
. (250)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 41
Now we simplify Eq. (249) by relaxing the constant terms:
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤
r
τ
B
5τ
SB
e
Λt+1
(3c
BM
+ 7) . (251)
With the relaxed form of c
BM
above this becomes:
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤ α(t)e(2
√
5 + 6eα(t) +
36e
2
√
5
α
2
(t)), α(t) =
r
τ
B
τ
SB
e
Λt
. (252)
The error near t = 0 is doubled as compared to the longer formula (249). Moreover, at each of decreasing
timescales t ∼ Λ
−1
, t ∼ 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
α ≤ 13 − 2
√
5e +
9e
3
α
2
13−2
√
5e
:
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤ 13
r
τ
B
τ
SB
e
4t/τ
SB
1 + 29
τ
B
τ
SB
e
8t/τ
SB
. (253)
where we plugged in Λ = 4/τ
SB
. We will use this form while computing the total error. We still need to
include the Markov and Born errors. Adding the bound in Eqs. (253) and (211), and using the triangle
inequality, we finally obtain:
kρ
true
t
− ρ
C
√
5
t
k
1
≤
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
12 + O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
+
13
√
τ
B
e
4t
τ
SB
1 +
29τ
B
e
8t
τ
SB
τ
SB
!
√
τ
SB
. (254)
Recall that this bound is applicable within the region e
4t
τ
SB
τ
B
τ
SB
≤ C for some constant C (from big-O), but
within that region any triple {t, τ
B
, τ
SB
} ≥ 0 is allowed. In particular, for sufficiently small τ
B
/τ
SB
, the
error at any value of t/τ
SB
can be shown to be small. As a corollary, we proved the validity in the weak
coupling limit, but the above inequality is more general than that. In particular, this proves Eq. (76), which
uses T
a
=
p
τ
B
τ
SB
/5.
We will finally show how to simplify Eq. (254) to the form O
e
6t
τ
SB
q
τ
B
τ
SB
presented in Eq. (28). Denoting
x = e
6t
τ
SB
q
τ
B
τ
SB
, we can relax the bound as follows:
kρ
true
t
− ρ
C
√
5
t
k
1
≤ x
2
12 + O
e
4t
τ
SB
τ
B
τ
SB
+ 13x
1 + 29x
2
. (255)
We then use the property of Big O notation: for any function f (t/τ
SB
, τ
B
/τ
SB
), if f(t/τ
SB
, τ
B
/τ
SB
) =
O
e
4t
τ
SB
τ
B
τ
SB
, then also f(t/τ
SB
, τ
B
/τ
SB
) = O(x
2
). We obtain:
kρ
true
t
− ρ
C
√
5
t
k
1
≤ x
2
12 + O
x
2
+ 13x
1 + 29x
2
. (256)
From this, the result of Eq. (28) immediately follows:
kρ
true
t
− ρ
C
√
5
t
k
1
≤ O(x) = O
e
6t
τ
SB
r
τ
B
τ
SB
. (257)
6.5 The Davies-Lindblad error
We write the CGME with an explicit dependence on T
a
:
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
] +
X
ω,ω
0
γ
T
a
ωω
0
(A
ω
ρ
C
A
ω
0
−
1
2
{ρ
C
, A
ω
0
A
ω
}) . (258)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 42
The error bound is known for arbitrary T
a
: see Eq. (242). We now wish to compare the CGME with finite
T
a
to the limit T
a
→ ∞ of the same equation:
˙ρ
C
(t) = −i[H + H
LS
, ρ
C
] +
X
ω,ω
0
γ
T
a
=∞
ωω
0
(A
ω
ρ
C
A
ω
0
−
1
2
{ρ
C
, A
ω
0
A
ω
}) + E
R
. (259)
In Appendix C, we show that this limit is just the Davies-Lindblad equation:
˙ρ
D
(t) = −i[H + H
LS
, ρ
D
] +
X
ω,ω
0
γ
T
a
=∞
ωω
0
(A
ω
ρ
D
A
ω
0
−
1
2
{ρ
D
, A
ω
0
A
ω
}) . (260)
We also show in Appendix C that the difference between the right-hand sides is bounded as:
kE
R
k
1
≤ 4
X
ωω
0
|γ
T
a
ωω
0
− γ
T
a
=∞
ωω
0
| ≤
4
n
τ
SB
O
1
T
a
δE
. (261)
Here the system Hilbert space dimension 2
n
appeared, and we note that 1/δE could contain an additional
exponential dependence on n. The accumulated error in the solution is bounded as:
kρ
true
(t) − ρ
D
(t)k
1
≤ (e
Λt
− 1)
kE
R
k
1
Λ
+ kρ
true
− ρ
C
k
1
. (262)
In the long time limit the leading terms are:
kρ
true
(t) − ρ
D
(t)k
1
≤ e
4t/τ
SB
4
n
O
1
T
a
δE
+
5T
a
τ
SB
e
4t/τ
SB
+ e
4t/τ
SB
τ
B
T
a
+ kρ
true
− ρ
BM
k
1
+ o(e
4t/τ
SB
) . (263)
Assuming that
4
n
δE
τ
B
, we find the optimal T
a
= 2
n
O(
p
τ
SB
/δE). Using that together with a bound for
t ≥ T
a
/2 we get:
kρ
true
(t) − ρ
D
(t)k
1
≤ 2
n
O
e
4t/τ
SB
√
τ
SB
δE
1 +
τ
B
τ
SB
!
+ kρ
true
− ρ
BM
k
1
≤ 2
n
O
e
12t/τ
SB
1
√
τ
SB
δE
+
τ
B
τ
SB
.
(264)
Assuming a constant system size n = O(1), this is the form presented in Eq. (27), our second main error
bound.
Then, assuming
4
n
δE
τ
B
, we can neglect the second term. Recall the discussion in Sec. 2.2: g =
1/
√
τ
B
τ
SB
, where g is the coupling strength g. Thus, we can express
p
1/τ
SB
δE = g
p
τ
B
/δE. If we set
t < cτ
SB
, we can write kρ
true
(t) − ρ
D
(t)k
1
= O(
p
1/τ
SB
δE) = O(g
p
τ
B
/δE), as stated in Eq. (1) in the
Introduction.
6.6 Redfield limits of integration error
There is one more error left to analyze: the error due to changing
R
t
→
R
∞
in the Redfield equation. We
already had the bound given in Eq. (47b):
kE
l
k
1
≤ 4c
BM
Z
∞
t
|C(t
0
)|dt
0
≤
4c
BM
t
Z
T
t
t
0
|C(t
0
)|dt
0
+ 4c
BM
Z
∞
T
|C(t
0
)|dt
0
≤
4c
BM
τ
B
tτ
SB
+
4c
BM
T
τ
SB
, (265)
where
T
= τ
SB
R
∞
T
|C(t
0
)|dt
0
[Eq. (48)]. We augment Eq. (265) by a trivial bound:
kE
l
k
1
≤ 4c
BM
Z
∞
t
|C(t
0
)|dt
0
≤
4c
BM
τ
SB
min(1, τ
B
/t +
T
) . (266)
Since the bound on kE
l
k
1
is time dependent, we cannot directly use Lemma 1 for Eq. (47a). Instead, we repeat
its proof and integrate the corresponding differential equations as was done in the proof of Lemma 2. We
define a function b
R
(t) ≥ kρ
BM
(t) − ρ
R
(t)k
1
(recall that ρ
BMl
= ρ
R
) that satisfies the following differential
equation, analogous to Eq. (144):
˙
b
R
=
4b
R
τ
SB
+
4c
BM
τ
SB
min(1, τ
B
/t +
T
) (267)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 43
Integrating it, we obtain:
kρ
BM
(t) − ρ
R
(t)k
1
≤ δ(t) = c
BM
e
4t/τ
SB
Z
4t/τ
SB
0
min(1,
4τ
B
xτ
SB
+
T
)e
−x
dx . (268)
As shown in Appendix I, this results in the error bound:
kρ
BM
(t) − ρ
R
(t)k
1
≤ c
BM
e
4t/τ
SB
4τ
B
τ
SB
1 − ln
1 − e
−
4τ
B
τ
SB
(1−
T
)
+
T
(269a)
≤ c
BM
e
4t/τ
SB
4τ
B
τ
SB
1 +
1
e
+ max
ln
τ
SB
(1 −
T
)
4τ
B
, 0
+
T
. (269b)
In other words, the O(τ
B
/τ
SB
) error scaling of the master equation just after Born-Markov approximation,
noted in Eq. (209) acquires a logarithmic correction.
We will now show how to simplify the total error of Redfield to obtain Eq.(26). We add the error of ρ
BM
from Eq. (211):
kρ
true
t
− ρ
R
t
k
1
≤
e
4t
τ
SB
− 1
e
8t
τ
SB
τ
B
12 + O
e
4t
τ
SB
τ
B
τ
SB
τ
SB
+ (270)
+ c
BM
e
4t/τ
SB
4τ
B
τ
SB
1 +
1
e
+ max
ln
τ
SB
(1 −
T
)
4τ
B
, 0
+
T
. (271)
We will set
T
= 0 in this estimate and use the following bound on c
BM
from Eq. (250):
c
BM
≤ 1 +
√
5X + 3X,
√
X ≤ 2
r
τ
B
τ
SB
e
4t
τ
SB
+1
. (272)
where we denoted x = e
12t
τ
SB
τ
B
τ
SB
. After a series of simplifications similar to the ones done after Eq. (255) we
arrive at:
kρ
true
t
− ρ
R
t
k
1
≤ O(x)(1 + max
ln
τ
SB
4τ
B
, 0
) . (273)
Since big O notation is defined for x ≤ for arbitrary small , it follows that τ
B
/τ
SB
≤ is also small, so the
logarithm is the leading term:
kρ
true
t
− ρ
R
t
k
1
≤ O(x)ln
τ
SB
τ
B
= O
e
12t
τ
SB
τ
B
τ
SB
ln
τ
SB
τ
B
. (274)
This is Eq. (26).
7 Summary and Conclusions
Our primary goal in this work was to derive a Markovian master equation with a complete positivity guarantee
that is in addition capable of incorporating arbitrary time-dependent driving, and that has a larger range
of applicability than the Davies-Lindblad Markovian master equation. The master equation that achieves
these desiderata is based on the coarse-graining approach introduced in Ref. [40], which only considered
the time-independent (no driving) case. We rederived this result by showing that it results from the non-
CP Redfield master equation after neglecting part of the time-integration domain, which restored complete
positivity. We then included arbitrary time-dependent driving and obtained a new CGME that is still CP,
and retains essentially the same form as its time-independent counterpart. This is our central new result,
given in Eq. (22).
This time-dependent CGME depends on one free parameter, the coarse graining time T
a
. We showed that
it is optimal to choose it as proportional to the geometric mean of two key timescales, the bath correlation
time τ
B
and the fastest system decoherence timescale τ
SB
, given in Eq. (11). These timescales are related to
the bath spectral density γ(ω) via the inequalities 2/τ
SB
≥ γ(ω) and 2τ
B
/τ
SB
≥ |γ
0
(ω)|. The 1-parameter
family of CGME equations involves a continuous integral, or equivalently, a sum over exponentially many
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 44
terms in the frequency domain. Alongside this 1-parameter family, we also derived an approximation in
terms of 3-parameter family (T
a
, ∆, k
∗
), given in Eq. (78). We showed that if [H, A] = O(1), where H is
the system Hamiltonian and A ⊗ B is the interaction with the bath [Eq. (3)], then there is no system-size
dependence in ∆, k
∗
, where k
∗
is the number of terms replacing the continuous integral.
Along the way we provided rigorous error bounds on three different master equations: Redfield, Davies-
Lindblad, and the CGME. These bounds [Eqs. (26)-(28)] are expressed in terms of the dimensionless ratio
τ
B
/τ
SB
, and in the Davies-Lindblad case also in terms of the quantity 1/
√
τ
SB
δE, where δE is the minimum
level spacing. These results are valid even when kBk diverges, so they apply in particular to, e.g., oscillator
baths. When the system-bath coupling strength g is finite we can write τ
B
/τ
SB
= (gτ
B
)
2
and 1/
√
τ
SB
δE =
g
p
τ
B
/δ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τ
B
1, but the Davies-
Lindblad master equation applies when in addition the minimum level spacing is large relative to g
2
τ
B
.
While these are all sufficient and not necessary conditions, if we informally view a large level spacing as a
requirement, then it is a very onerous one indeed, since in many cases of interest the spacing shrinks rapidly
(polynomially or even exponentially) with problem size, e.g., in adiabatic quantum computing [77]. Thus
the guaranteed range of applicability of the Davies-Lindblad master equation is severely limited, while the
CGME achieves complete positivity without imposing this cost.
To summarize the latter in the simplest possible terms, what we have shown is that the CGME and
Redfield master equations are valid for fast baths, in the sense of a very simple integral condition on the
bath correlation function. Ohmic and super-Ohmic baths satisfy this condition, provided the system-bath
coupling strength is weak but constant in system size (for a local system). For the Davies-Lindblad master
equation the coupling strength needs to be exponentially small in the system size. However, the Redfield
equation and the CGME are not on equal footing: the Redfield equation provides no guarantee of complete
positivity, while the CGME does.
We list a few potential applications of the CGME:
• The CGME opens the door to studies of open system dynamics with arbitrary time-dependent drives
H(t). In particular, the drive need not be adiabatic. In fact, the time-dependent CGME even correctly
describes a system driven by δ-function pulses, as demonstrated in Sec. 5.1 for dynamical decoupling.
Moreover, the time-dependent CGME overcomes an inherent problem in applying the Davies-Lindblad
ME to circuit model quantum computing, in particular in the context of fault-tolerant quantum error
correction. The latter ME is incompatible with the assumption of arbitrarily fast gates needed there [43],
but the time-dependent CGME can be consistently applied in this context. This opens the door to a
rigorous analysis of fault-tolerance in a first-principles Markovian setting.
• Earlier studies of the lifetime of the toric code at finite temperature took advantage of the fact that the
Davies generators are local for commuting models [78]. The CGME provides local generators of open
system finite temperature dynamics for non-commuting models (see Sec. 3.6). This opens the door to
studying the finite-temperature lifetime of, e.g., topological subsystem codes [79] in the thermodynamic
limit.
• Many-body localization effects in quantum annealing — previously studied only in the closed system
setting [80] — can now be extended to finite temperature open system quantum annealing.
We hope that this work, including the applications above, will stimulate further research into the uses of
the CGME, the only master equation we know of that provides a complete positivity guarantee, is locally
generated, and applies for arbitrarily driven open quantum systems.
Acknowledgements
The research is based upon work (partially) supported by the Office of the Director of National Intelligence
(ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office
contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and
should not be interpreted as necessarily representing the official policies or endorsements, either expressed or
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 45
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.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 46
A Proof that γ(ω) > 0
Trivial case. First consider V = A ⊗ B. The function γ(ω) is the Fourier transform of the correlation
function:
γ(ω) =
Z
∞
−∞
e
iωt
C(t)dt =
Z
∞
−∞
e
iωt
Tr[ρ
B
B(t)B]dt (275)
Our derivations in this paper are for a stationary bath ρ
B
(t) = ρ
B
, or [ρ
B
, H
B
] = 0. Note that [
√
ρ
B
, H
B
] = 0
follows.
γ(ω) =
X
kl
Z
∞
−∞
e
iωt
Tr[
√
ρ
B
e
iH
B
t
√
ρ
B
Be
−iH
B
t
B]dt (276)
Inserting the sum over bath eigenstates:
γ(ω) =
X
kl
Z
∞
−∞
e
i(ω+E
kl
)t
hk|
√
ρ
B
B|lihl|B
√
ρ
B
|kidt (277)
Integration results in delta-functions:
γ(ω) = 2π
X
kl
δ(ω −E
lk
)|hl|B
√
ρ
B
|ki|
2
≥ 0 (278)
This is enough for all the results in this paper. A more complicated interaction with environment, however,
requires a different proof.
General case. We consider the case V =
P
i
A
i
⊗ B
i
, so that γ(ω) becomes a matrix. ρ
B
will still be
assumed stationary: [ρ
B
, H
B
] = 0. Repeating the same steps, we arrive at
γ
ij
(ω) = 2π
X
kl
δ(ω −E
lk
)hk|
√
ρ
B
B
i
|lihl|B
j
√
ρ
B
|ki (279)
Thus we need to prove that the matrix
M
ij
= hk|
√
ρ
B
B
i
|lihl|B
j
√
ρ
B
|ki (280)
is positive. Let v
i
= hk|
√
ρ
B
B
i
|li, then
M
ij
= v
i
v
∗
j
. (281)
For positivity, we need to prove that for any vector ~w with elements w
i
it holds that
P
ij
w
i
M
ij
w
∗
j
≥ 0.
Indeed:
∀~w,
X
ij
w
i
M
ij
w
∗
j
=
X
ij
w
i
v
i
v
∗
j
w
j
=
X
i
v
i
w
i
2
≥ 0 . (282)
Thus
γ(ω) ≥ 0 (283)
as a matrix.
Bochner’s theorem Our proof above used the eigendecomposition of the bath, and γ(ω) was a sum of
Dirac δ functions. In practice, a smooth γ(ω) is used, corresponding to, e.g., a bath with a continuous
spectrum. In this sense it is desirable to avoid using bath eigenstates and δ-functions altogether. Here we
give another proof that uses Bochner’s theorem, as suggested, e.g., in the textbook [2], and which avoids the
issue mentioned above.
Lemma 3 (Bochner). If a function φ(t) is of positive type, i.e. ∀n ∈ Z, t
i
∈ R, ξ
i
∈ C
n
X
i,j=1
φ(t
i
− t
j
)ξ
i
ξ
∗
j
≥ 0, (284)
and also φ(0) is finite and φ(t) is continuous at t = 0, then its Fourier transform (if it exists) is non-negative:
Z
∞
−∞
φ(t)e
iωt
dt ≥ 0 (285)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 47
We again consider many interaction terms V =
P
i
A
i
⊗ B
i
(as in Sec. 3.7) and a stationary ρ
B
s.t.
[ρ
B
, H
B
] = 0. Since γ(ω) is Hermitian we can diagonalize it using a unitary transformation:
D ≡ UγU
†
⇒ D
αβ
=
X
i,j
U
αi
γ
ij
U
∗
βj
. (286)
D is diagonal so we only need to consider the diagonal elements (i.e., the eigenvalues of γ). Substituting
γ
ij
=
R
∞
−∞
e
iωs
C
ij
(s)ds, where C
ij
(s) = Tr[ρ
B
B
i
(s)B
j
], gives
D
α
=
Z
∞
−∞
e
iωs
X
i,j
U
αi
C
ij
(s)U
∗
αj
ds . (287)
We wish to show that D
α
is non-negative for each α. To do this we must consider the function in parenthesis.
D
α
is the Fourier transform of this function so if we can show that it is of positive type then D
α
must be
positive by Bochner’s theorem [81]. Define the following function with {t
i
} an arbitrary time partition:
f
α
mn
≡
X
i,j
U
αi
C
ij
(t
m
− t
n
)U
∗
αj
. (288)
Note that
B
αβ
(t, t − τ) ≡ hB
α
(t)B
β
(t − τ)i
B
= Tr
e
iH
B
t
B
α
e
−iH
B
t
e
iH
B
(t−τ)
B
β
e
−iH
B
(t−τ)
ρ
B
(289a)
= Tr
e
−iH
B
(t−τ)
e
iH
B
t
B
α
e
−iH
B
t
e
iH
B
(t−τ)
B
β
ρ
B
= Tr
e
iH
B
τ
B
α
e
−iH
B
τ
B
β
ρ
B
= hB
α
(τ)B
β
i
B
(289b)
= B
αβ
(τ, 0) ≡ B
αβ
(τ) , (289c)
Now use the property hB
α
(s)B
β
(0)i = hB
α
(t)B
β
(t − s)i [Eq. (289c)] to write f
α
mn
as
f
α
mn
=
X
i,j
U
αi
Tr [ρ
B
B
i
(t
m
)B
j
(t
n
)] U
∗
αj
= Tr
ρ
B
X
i
U
αi
B
i
(t
m
)
X
j
U
∗
αj
B
j
(t
n
)
. (290)
We need to show that f
α
is a positive matrix. For arbitrary |vi we have
hv|f
α
|vi =
X
m,n
v
∗
m
v
n
f
α
mn
= Tr
X
i,m
v
∗
m
U
αi
√
ρ
B
B
i
(t
m
)
X
j,n
v
n
U
∗
αj
B
j
(t
n
)
√
ρ
B
(291a)
= Tr
M
†
α
M
α
≥ 0 , (291b)
where M
α
≡
P
i,m
v
∗
m
U
αi
√
ρ
B
B
i
(t
m
).
We have established that hv|f
α
|vi ≥ 0 for any time partition {t
i
}. Therefore D
α
is non-negative by
Bochner’s theorem. Consequently, γ ≥ 0 as a matrix since all its eigenvalues are non-negative.
B Lamb shift simplification
The Lamb shift was given in Eq. (69) as H
LS
=
i
2T
a
R
T
a
/2
−T
a
/2
dt
1
R
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)A(t
2
)A(t
1
).
Recall also that A(t) =
P
ω
A
ω
e
−iωt
[Eq. (5)], so we can write:
H
LS
=
X
ω
1
ω
2
F
ω
1
ω
2
A
ω
2
A
ω
1
, (292)
where
F
ω
1
ω
2
=
i
2T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(t
1
− t
2
)C(t
2
− t
1
)e
−i(ω
2
t
2
+ω
1
t
1
)
(293a)
=
i
2T
a
Z
T
a
/2
−T
a
/2
dt
1
Z
T
a
/2
−T
a
/2
dt
2
sgn(θ
−
)C(−θ
−
)e
−i(ω
−
θ
−
+ω
+
θ
+
)
, (293b)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 48
where ω
±
=
1
2
(ω
1
±ω
2
) and θ
±
= t
1
±t
2
. This change of variables rotates the square by π/4 into a rhombus,
and the integration limits correspondingly change to:
θ
−
> 0 : −T
a
+ θ
−
≤ θ
+
≤ T
a
− θ
−
(294a)
θ
−
< 0 : −T
a
− θ
−
≤ θ
+
≤ T
a
+ θ
−
. (294b)
The Jacobian is 1/2. We can now perform the integral over θ
+
:
F
ω
1
ω
2
=
i
4T
a
−
Z
0
−T
a
dθ
−
Z
T
a
+θ
−
−T
a
−θ
−
dθ
+
+
Z
T
a
0
dθ
−
Z
T
a
−θ
−
−T
a
+θ
−
dθ
+
!
C(−θ
−
)e
−i(ω
−
θ
−
+ω
+
θ
+
)
(295a)
=
1
4T
a
ω
+
Z
0
−T
a
dθ
−
(e
−iω
+
(T
a
+θ
−
)
− e
iω
+
(T
a
+θ
−
)
) −
Z
T
a
0
dθ
−
(e
−iω
+
(T
a
−θ
−
)
− e
iω
+
(T
a
−θ
−
)
)
!
C(−θ
−
)e
−iω
−
θ
−
(295b)
=
1
4T
a
ω
+
Z
0
−T
a
dθ(e
−iT
a
ω
+
−iω
1
θ
− e
iT
a
ω
+
+iω
2
θ
)C(−θ) −
Z
T
a
0
dθ(e
−iT
a
ω
+
+iω
2
θ
− e
iT
a
ω
+
−iω
1
θ
)C(−θ)
!
(295c)
=
1
4T
a
ω
+
Z
T
a
0
dθ(e
−iT
a
ω
+
+iω
1
θ
− e
iT
a
ω
+
−iω
2
θ
)C(θ) −
Z
T
a
0
dθ(e
−iT
a
ω
+
+iω
2
θ
− e
iT
a
ω
+
−iω
1
θ
)C
∗
(θ)
!
(295d)
=
1
4T
a
ω
+
Z
T
a
0
dθ(e
−iT
a
ω
+
+iω
1
θ
− e
iT
a
ω
+
−iω
2
θ
)C(θ) + (e
−iT
a
ω
+
+iω
1
θ
− e
iT
a
ω
+
−iω
2
θ
)
∗
C
∗
(θ) (295e)
where we used the identity C(θ) = C
∗
(−θ) and a change of variables θ 7→ −θ to make all the integration
domains positive. When ω
+
= 0 the expression above should be understood as a lim
ω
+
→0
. Thus:
F
ω
1
ω
2
=
1
2T
a
ω
+
Re
Z
T
a
0
dθ
e
i(ω
1
θ−T
a
ω
+
)
− e
−i(ω
2
θ−T
a
ω
+
)
C(θ) , (296)
which is Eq. (21b) in the main text.
Recall that A
†
ω
= A
−ω
; for H
LS
=
P
ω
1
ω
2
F
ω
1
ω
2
A
ω
2
A
ω
1
to be Hermitian, it is thus sufficient if
F
ω
1
ω
2
= F
∗
ω
1
ω
2
= F
−ω
2
,−ω
1
. (297)
This is indeed satisfied by Eq. (296), thus ensuring the Hermiticity.
C The Davies-Lindblad equation is the infinite coarse-graining time limit of the CGME
The RWA used to derive the Davies-Lindblad equation leaves something to be desired. We simply dropped
terms with different Bohr frequencies, without a rigorous mathematical justification. Here we repeat, with
some extra clarifications, the derivation given in Ref. [40] which shows that the Davies-Lindblad equation is
the infinite coarse-graining timescale limit of the CGME.
First we rewrite Eq. (70) in the form
γ
ωω
0
(T
a
) =
1
T
a
Z
T
a
0
ds
Z
T
a
0
ds
0
e
i(ω
0
s−ωs
0
)
C(s −s
0
) , (298)
which we can do after changing variables ω
0
→ −ω
0
as remarked below Eq. (71). Our goal is to show that
lim
T
a
→∞
γ
ωω
0
(T
a
) = γ(ω)δ
ωω
0
.
Lemma 4. The following equivalent form holds for γ
ωω
0
(T
a
):
γ
ωω
0
(T
a
) =
1
T
a
e
i
ω
0
−ω
2
T
a
Z
T
a
0
dv cos
ω
0
− ω
2
(v − T
a
)
Z
v
−v
du e
i
ω+ω
0
2
u
C(u) . (299)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 49
Proof. Let us rewrite ω
0
s − ωs
0
in terms of a sum and difference of Bohr frequencies:
ω
0
s − ωs
0
=
1
2
(ω
0
− ω)v +
1
2
(ω
0
+ ω)u , (300)
where u = s − s
0
and v = s + s
0
. After this change of variables C(s − s
0
) = C(u), and since s = (v + u)/2
and s
0
= (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 ∈ [0, T
a
] and the lines u = 2T
a
−v
and v − 2T
a
for v ∈ [T
a
, 2T
a
]. Thus:
γ
ωω
0
(T
a
) =
1
2T
a
Z
T
a
0
dv e
i
ω
0
−ω
2
v
Z
v
−v
du e
i
ω+ω
0
2
u
C(u) +
1
2
Z
2T
a
T
a
dv e
i
ω
0
−ω
2
v
Z
2T
a
−v
−(2T
a
−v)
du e
i
ω+ω
0
2
u
C(u) . (301)
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:
γ
ωω
0
(T
a
) =
1
2T
a
Z
T
a
0
dv e
i
ω
0
−ω
2
[(v−T
a
)+T
a
]
Z
v
−v
du e
i
ω+ω
0
2
u
C(u) +
1
2
Z
T
a
0
dv e
−i
ω
0
−ω
2
[(v−T
a
)−T
a
]
Z
v
−v
du e
i
ω+ω
0
2
u
C(u)
(302a)
=
1
T
a
e
i
ω
0
−ω
2
T
a
Z
T
a
0
dv cos
ω
0
− ω
2
(v − T
a
)
Z
v
−v
du e
i
ω+ω
0
2
u
C(u) , (302b)
which is Eq. (298).
C.1 The ω = ω
0
case
For ω = ω
0
we now have:
γ
ωω
(t) =
1
T
a
Z
T
a
0
dv
Z
v
−v
du e
iωu
C(u) . (303)
Let U =
R
v
−v
du e
iωu
C(u). Then dU =
e
iωv
C(v) + e
−iωv
C(−v)
dv, and integrating by parts gives:
γ
ωω
(T
a
) =
1
T
a
v
Z
v
−v
du e
iωu
C(u)
T
a
0
−
1
T
a
Z
T
a
0
dv v
e
iωv
C(v) + e
−iωv
C(−v)
. (304)
Consider the second integral:
1
T
a
Z
T
a
0
dv ve
iωv
C(v)
≤
1
T
a
Z
T
a
0
dv v |C(v)| ≤
1
T
a
Z
∞
0
dv v |C(v)| (305a)
=
τ
B
τ
SB
T
a
T
a
→∞
−→ 0 , (305b)
where in the last step we used Eq. (11b). Since C(v) = C
∗
(−v), the third integral in Eq. (304) satisfies the
same bound and limit. We are thus left with
lim
T
a
→∞
γ
ωω
(T
a
) =
Z
∞
−∞
du e
iωu
C(u) = γ(ω) , (306)
with the last equality being due to Eq. (6b).
C.2 The ω 6= ω
0
case
For ω 6= ω
0
we also perform integration by parts of Eq. (299), but we shall see that this time the boundary
terms vanish. We write γ
ωω
0
(T
a
) =
1
T
a
e
i
ω
0
−ω
2
T
a
R
T
a
0
dV U(v), where now dV = cos
ω
0
−ω
2
(v − T
a
)
dv and
U(v) =
R
v
−v
du e
i
ω+ω
0
2
u
C(u). Then
V (v) =
2
ω
0
− ω
sin
ω
0
− ω
2
(v − T
a
)
(307a)
dU/dv = e
i
ω+ω
0
2
v
C(v) + e
−i
ω+ω
0
2
v
C(−v) (307b)
[U(v)V (v)]
T
a
0
= U(T
a
)V (T
a
) − U(0)V (0) = 0 . (307c)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 50
Therefore:
γ
ωω
0
(T
a
) = −
Z
T
a
0
V dU = −
2e
i
ω
0
−ω
2
T
a
(ω
0
− ω)T
a
Z
T
a
0
dv sin
(ω
0
− ω)
2
(v − T
a
)
e
i
ω+ω
0
2
v
C(v) + e
−i
ω+ω
0
2
v
C(−v)
.
(308)
Changing from v to −v in the second term we get
γ
ωω
0
(T
a
) = −
2e
i
ω
0
−ω
2
T
a
(ω
0
− ω)T
a
"
Z
T
a
0
dv sin
(ω
0
− ω)
2
(v − T
a
)
e
i
ω+ω
0
2
v
C(v) +
Z
0
−T
a
dv sin
(ω
0
− ω)
2
(−v − T
a
)
e
i
ω+ω
0
2
v
C(v)
#
(309a)
=
e
i
ω
0
−ω
2
T
a
(ω
0
− ω)T
a
Z
T
a
−T
a
dv
sin
ω
0
− ω
2
T
a
e
iωv
+ e
iω
0
v
+
sgn(v)
i
cos
ω
0
− ω
2
T
a
e
iωv
− e
iω
0
v
C(v) ,
(309b)
where we used the angle sum identity for the sine in the last equality. Thus:
lim
T
a
→∞
γ
ωω
0
(T
a
) = lim
T
a
→∞
e
i
ω
0
−ω
2
T
a
(ω
0
− ω)T
a
sin
ω
0
− ω
2
T
a
(γ(ω) + γ(ω
0
)) + 2 cos
ω
0
− ω
2
T
a
(S(ω) − S(ω
0
))
.
(310)
where S(ω) is the Lamb shift amplitude in Eq. (18b). Since nothing cancels with the overall T
−1
a
, we find
that the ω 6= ω
0
term vanishes.
A similar calculation could be done for the Lamb shift term (69), showing that Im(x
ωω
0
) → Im(x
ωω
)δ
ωω
0
.
Therefore, the Davies-Lindblad (RWA) equation can be understood as the T
a
→ ∞ limit of the CGME.
D Derivation of the discrete approximation Eq. (78)
The starting point is the time-independent CGME given in Eq. (19). We can obtain a discrete sum just by
discretizing
R
d. Below is the estimate of the number of terms (2k
∗
−1) in the sum for that case. The error
is:
E
k
=
Z
∞
−∞
d(A
ρ
C
A
†
−
1
2
n
ρ
C
, A
†
A
o
) − ∆
X
k,=∆k, |k|<k
∗
(A
ρ
C
A
†
−
1
2
n
ρ
C
, A
†
A
o
) , (311)
where A
is given in Eq. (20a). The discretization of an integral introduces errors as follows:
Z
df() − ∆
X
k,=∆k, |k|<k
∗
f()
1
≤
k
∗
− 0.5
2
max kf
0
k
1
∆
2
+
Z
−∆(k
∗
−0.5)
−∞
+
Z
∞
∆(k
∗
−0.5)
!
dkf()k
1
,
(312)
where f is an operator or operator-valued function.
Proof. By the triangle inequality:
Z
df() − ∆
X
k,=∆k, |k|<k
∗
f()
1
≤
Z
∆(k
∗
−0.5)
−∆(k
∗
−0.5)
df() − ∆
X
k,=∆k, |k|<k
∗
f()
1
(313a)
+
Z
−∆(k
∗
−0.5)
−∞
+
Z
∞
∆(k
∗
−0.5)
!
dkf()k
1
. (313b)
We apply the triangle inequality to pull the sum out:
Z
∆(k
∗
−0.5)
−∆(k
∗
−0.5)
df() − ∆
X
k,=∆k, |k|<k
∗
f()
1
≤
X
k,=∆k, |k|<k
∗
Z
+∆/2
−∆/2
dµf(µ) − ∆f ()
1
. (314)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 51
To bound this, we use O(t) = O(0) +
R
t
0
O
0
(θ)dθ, so that
R
1/2
−1/2
O(t)dt − O(0) =
R
1/2
−1/2
R
t
0
O
0
(θ)dθdt, which
implies
Z
1/2
−1/2
O(t)dt − O(0)
1
≤
1
4
max
θ∈[0,1]
kO
0
k
1
, (315)
which applied to the interval [, + ∆] reads:
Z
+∆/2
−∆/2
dµf(µ) − f ()
1
≤
∆
2
4
max
µ∈[,+∆]
kf
0
k
1
. (316)
Substituting this yields:
Z
∆(k
∗
−0.5)
−∆(k
∗
−0.5)
df() − ∆
X
k,=∆k, |k|<k
∗
f()
1
≤ (2k
∗
− 1)
∆
2
4
max
µ∈[−∆k
∗
,∆k
∗
]
kf
0
k
1
, (317)
which concludes the proof.
If we use the proven Eq. (315) for
f
= (A
ρ
C
A
†
−
1
2
n
ρ
C
, A
†
A
o
) , (318)
it will express the error E
k
as the sum of:
E
k
= E
k1
+ E
k2
, E
k1
=
k
∗
− 0.5
2
max kf
0
k
1
∆
2
, E
k1
=
Z
−∆(k
∗
−0.5)
−∞
+
Z
∞
∆(k
∗
−0.5)
!
dkf()k
1
. (319)
Let us start with the E
k2
term corresponding to truncating the integral at k
∗
. As becomes far removed
from all the frequencies ω in the system, this integral becomes small as 1/. To prove this, we integrate by
parts:
A
=
s
γ()
2πT
a
Z
T
a
/2
−T
a
/2
A(t)e
it
dt =
s
γ()
2πT
a
1
i
A(T
a
/2)e
iT
a
/2
− A(−T
a
/2)e
−iT
a
/2
+ i
Z
T
a
/2
−T
a
/2
[H, A(t)]e
it
dt
!
.
(320)
Taking the norm and recalling that kAk = 1, we get:
kA
k ≤
s
γ()
2πT
a
2 + k[H, A]kT
a
||
≤
s
1
πτ
SB
T
a
2 + k[H, A]kT
a
||
. (321)
Note that for a local n qubit Hamiltonian with O(1) local terms and a local operator A with kAk = 1, the
commutator k[H, A]k = O(1). Combining, we find:
kE
k2
k
1
≤
Z
−∆(k
∗
−0.5)
−∞
+
Z
∞
∆(k
∗
−0.5)
!
kA
k
2
d ≤ 2
Z
∞
∆(k
∗
−0.5)
1
πT
a
τ
SB
(2 + k[H, A]kT
a
)
2
||
2
d (322a)
= 2
1
πT
a
τ
SB
(2 + k[H, A]kT
a
)
2
∆(k
∗
− 0.5)
. (322b)
We enforce that the error is
p
τ
B
/τ
SB
1
τ
SB
, so as to match the other errors in our equation:
r
τ
B
τ
SB
1
τ
SB
= 2
1
πT
a
τ
SB
(2 + k[H, A]kT
a
)
2
∆(k
∗
− 0.5)
⇒ ∆(k
∗
− 0.5) = 2
1
πT
a
r
τ
SB
τ
B
(2 + k[H, A]kT
a
)
2
. (323)
Let us now focus on the term with the derivatives E
k1
:
kE
k1
k
1
≤ 2(k
∗
− 0.5) max kA
0
kkA
k∆
2
. (324)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 52
The necessary norm bounds are:
kA
k ≤
s
γ()T
a
2π
, kA
0
k ≤
1
2
s
|γ
0
()|T
a
2πγ()
+
T
a
4
s
γ()T
a
2π
. (325)
Using Eqs. (12) and (13) we have the bound
kA
kkA
0
k ≤
1
2
p
|γ
0
()|T
a
2π
+
T
a
4
γ()T
a
2π
≤
1
2
p
2τ
B
/τ
SB
T
a
2π
+
1
4
T
2
a
2πτ
SB
. (326)
Substituting T
a
=
p
τ
B
τ
SB
/5, we get:
kA
kkA
0
k ≤
p
2/5τ
B
4π
+
τ
B
40π
=
(10
p
2/5 + 1)τ
B
40π
, (327)
and the error bound is then:
kE
k1
k
1
≤ 2(k
∗
− 0.5)
(10
p
2/5 + 1)τ
B
10π
∆
2
. (328)
This leads to the step-size choice:
2(k
∗
− 0.5)
(10
p
2/5 + 1)τ
B
10π
∆
2
=
1
τ
SB
r
τ
B
τ
SB
, ∆ =
1
τ
SB
√
τ
B
τ
SB
(k
∗
− 0.5)∆
5π
(10
p
2/5 + 1)
. (329)
Combining this with Eq. (323), the solution is:
∆ =
1
τ
SB
r
τ
B
τ
SB
1
(2 + k[H, A]kT
a
)
2
√
5π
2
2(10
p
2/5 + 1)
, k
∗
= ceil
r
τ
SB
τ
B
3
(2 + k[H, A]kT
a
)
4
4(10
p
2/5 + 1)
π
3
+ 0.5
!
.
(330)
If [H, A] = O(1), there is no system-size dependence in ∆, k
∗
.
We remark that the bounds we used to derive Eq. (330) can be significantly tightened; our goal here was
to show that they are O(1). We do not actually recommend to use Eq. (330) for numerical calculations: one
can obtain accurate results with much larger ∆ and much smaller k
∗
.
E An error bound that is linear in t
Lemma 5. Let L be a linear superoperator and consider the pair of equations
˙x(t) = Lx(t) + E, ˙y(t) = Ly(t) . (331)
If any solution y(t) possesses the property that ∀t, y(0) such that ky(0)k
p
≤ 1 : ky(t)k
p
≤ c
y
, then for any
finite t:
kx(t) − y(t)k
p
≤ c
y
tkEk
p
. (332)
The result holds in the same form for the operator (p = ∞, omitted in our notation) and trace (p = 1)
norms.
Proof. The linear differential equation for y has a unique solution in a neighborhood of any initial condition.
A trivial corollary is that the evolution operator is reversible for any finite t:
e
Lt
e
−Lt
= . (333)
This is an equality between superoperators. Now we can express the difference as:
x(t) − y(t) = e
Lt
(e
−Lt
x(t) − y(0)) . (334)
The explicit form of the solution x(t) is:
x(t) = e
Lt
x(0) + e
Lt
Z
t
0
e
−Lτ
Edτ , (335)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 53
which yields, using Eq. (334):
x(t) − y(t) = e
Lt
Z
t
0
e
−Lτ
Edτ =
Z
t
0
e
L(t−τ)
Edτ . (336)
We now take the norm of both sides:
kx(t) − y(t)k
p
≤ tke
L(t−τ)
Ek
p
. (337)
By linearity of the time evolution and the norm:
ke
L(t−τ)
Ek
p
= kEk
p
ke
L(t−τ)
yk
p
, y =
E
kEk
p
(338)
Note that kyk
p
= 1. Since we can take y as the initial condition of the equation ˙y(t) = Ly(t), and since
e
L(t−τ)
y will then be the solution at time t −τ, we have, by assumption of the lemma:
ke
L(t−τ)
yk
p
≤ c
y
. (339)
Note that we had to make these assumptions for any y of bounded norm, not just density operators. It now
follows from Eq. (337) that:
kx(t) − y(t)k
p
≤ c
y
tkEk
p
. (340)
In particular, if L is a positive trace-preserving map then c
y
= 1.
The only error for which we can apply Lemma 5 in this work is the Markov error, which we present here:
kρ
BM
t
− ρ
B
t
k
1
≤ 16
τ
B
c
B
c
BM
t
τ
2
SB
, (341)
where c
y
= c
BM
and we used Eq. (203), which contains c
B
. We instead used a loose bound from Eq. (204)
since it results in a more compact expression. Apart from this we have no occasion to use Lemma 5 in
this work, since not all the equations appearing in our derivation are invertible. This is why we instead use
Lemma 1 or prove the necessary results independently.
F Proof of Lemma 2
Proof. The formal solutions allow us to write:
δρ(t) =
Z
t
0
L
θ
π(θ)dθ −
Z
t
0
L
θ
ρ(θ)dθ . (342)
Substituting π(t) = ρ(t) + δρ(t) gives:
δρ(t) =
Z
t
0
(L
θ
ρ(θ) − L
θ
ρ(θ))dθ +
Z
t
0
L
θ
δρ(θ)dθ . (343)
Making the time-averaging explicit:
δρ(t) =
1
T
a
Z
t
0
dθ
Z
T
a
/2
−T
a
/2
dτ[L
θ+τ
ρ(θ) − L
θ
ρ(θ)] +
Z
t
0
L
θ
δρ(θ)dθ . (344)
We would like to change the integration variables of the first term L
θ+τ
ρ(θ) in such a way that it becomes
L
θ
0
ρ(θ
0
−τ
0
). The change of variables is θ + τ = θ
0
, τ = τ
0
. In the second term −L
θ
ρ(θ) the change is trivial
θ = θ
0
, τ = τ
0
. 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 L
θ
0
ρ(θ
0
− τ
0
) to be the same as the second term. This is
illustrated in Fig. 11 and the corresponding integrals are:
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 54
Figure 11: Times contributing to δρ(t) for (a) t > T
a
/2 (b) t < T
a
/2.
Z
t
0
dθ
Z
T
a
/2
−T
a
/2
dτ =
Z
T
a
/2
−T
a
/2
dτ
0
Z
t+τ
0
τ
0
dθ
0
=
Z
T
a
/2
−T
a
/2
dτ
0
Z
t
0
dθ
0
+
Z
t+τ
0
t
dθ
0
−
Z
τ
0
0
dθ
0
!
(345)
After this change of variables and collecting all the contributions from the L
θ
ρ(θ) term in Eq. (344) into a
single term, we get:
δρ(t) =
1
T
a
Z
T
a
/2
−T
a
/2
dτ
0
Z
t
0
dθ
0
L
θ
0
(ρ(θ
0
−τ
0
)−ρ(θ
0
))+
1
T
a
Z
T
a
/2
−T
a
/2
dτ
0
Z
t+τ
0
t
dθ
0
−
Z
τ
0
0
dθ
0
!
L
θ
0
ρ(θ
0
− τ
0
)+
Z
t
0
L
θ
δρ(θ)dθ.
(346)
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 ≤ T
a
/2 the integrals in the second term cancel
over the area (T
a
/2 − t)
2
, as indicated by the yellow triangles in Fig. 11. In fact by splitting the integral
above differently we could get a slightly tighter (parametrically better for the transient) bound, but working
with rectangular regions is easier. Now we can take the norm:
kδρ(t)k
1
≤
Λ
T
a
Z
t
0
dθ
0
Z
T
a
/2
−T
a
/2
dτ
0
kρ(θ
0
− τ
0
) − ρ(θ
0
)k
1
+
c
ρ
Λ
2T
a
T
2
a
− (max (T
a
− 2t, 0))
2
+ Λ
Z
t
0
kδρ(θ)k
1
dθ ,
(347)
where Λ is a number such that
Λ ≥ sup
θ∈[0,t],kρk
1
=1
kL
θ
(ρ)k
1
. (348)
The difference between ρ’s at nearby time points is bounded trivially as:
˙ρ = Lρ ⇒ kρ(θ
0
− τ
0
) − ρ(θ
0
)k
1
≤ c
ρ
Λ|τ
0
| , (349)
which can be used to simplify:
kδρ(t)k
1
≤
c
ρ
Λ
T
a
Z
t
0
dθ
0
Z
T
a
/2
−T
a
/2
dτ
0
Λ|τ
0
| +
c
ρ
Λ min(T
a
, 4t − 4t
2
/T
a
)
2
+ Λ
Z
t
0
kδρ(θ)k
1
dθ (350a)
≤
c
ρ
Λ
2
tT
a
4
+
c
ρ
Λ min(T
a
, 4t − 4t
2
/T
a
)
2
+ Λ
Z
t
0
kδρ(θ)k
1
dθ . (350b)
We analyzed a similar inequality in the proof of Lemma 1, and as we saw there this amounts to solving
for b
2
(t) that saturates the inequality and upper bounds δρ(t). Note that despite the min function there
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 55
is no abrupt slope change on the r.h.s. and the derivative matches, so there is no discontinuity in
˙
b
2
. The
equations for b
2
such that kδρ(t)k
1
≤ b
2
(t) are:
˙
b
2
(t) =
c
ρ
Λ
2
T
a
4
+ c
ρ
Λ(2 − 4t/T
a
) + Λb
2
(t), b
2
(0) = 0, t ≤ T
a
/2 (351a)
˙
b
2
(t) =
c
ρ
Λ
2
T
a
4
+ Λb
2
(t), t ≥ T
a
/2 (351b)
Solving the first equation, we find:
b
2
(t) = 4c
ρ
t/T
a
+ c
ρ
4 − 2ΛT
a
− (ΛT
a
/2)
2
ΛT
a
(1 − e
Λt
) (352a)
b
2
(T
a
/2) = c
ρ
4 − (ΛT
a
/2)
2
ΛT
a
+ c
ρ
e
ΛT
a
/2
(ΛT
a
/2)
2
+ 2ΛT
a
− 4
ΛT
a
=
c
ρ
ΛT
a
2
+ O(Λ
2
T
2
a
) (352b)
where the last estimate is in the limit ΛT
a
1. Now we use that as the initial condition for the second
equation to find:
b
2
(t) = [b
2
(T
a
/2) + c
ρ
ΛT
a
/4]e
Λ(t−T
a
/2)
− c
ρ
ΛT
a
/4, t > T
a
/2 . (353)
Notice that
˙
b
2
> 0 always. We can express a bound kδρ(t)k
1
≤ b
2
(t) for all t as follows:
kδρ(t)k
1
≤ (b
2
(min(t, T
a
/2)) + c
ρ
ΛT
a
/4) max(e
Λ(t−T
a
/2)
, 1) − c
ρ
ΛT
a
/4 . (354)
The simplest observation is that for t ≤ c/Λ:
kδρ(t)k
1
≤ O(ΛT
a
) . (355)
G Derivation of the bound (245)
Here we discuss the derivation of the error introduced by time-averaging and restoring complete positivity,
starting from Eq. (239):
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ c
BM
e
Λ max(t,T
a
/2)
3T
a
τ
SB
−
T
a
τ
SB
+ e
Λt
2T
a
τ
SB
+
4τ
B
ΛT
a
τ
SB
max(e
Λ(t−T
a
/2)
− 1, 0) .
(356)
Similarly to Eq. (242), we loosen the bound by using e
−ΛT
a
/2
≤ 1. For t > T
a
/2 and c
Λ
= 4/Λτ
SB
≥ 1, we
obtain:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤ ((3c
BM
+ 2)e
Λt
− c
BM
)
T
a
τ
SB
+ (e
Λt
− 1)
c
Λ
τ
B
T
a
. (357)
At T
a
=
p
c
Λ
τ
B
τ
SB
/(3c
BM
+ 2) for t > T
a
/2 we have:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
Λt
−
1 +
c
BM
(3c
BM
+ 2)
. (358)
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:
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
Λ max(t,T
a
/2)
−
1 +
c
BM
(3c
BM
+ 2)
, T
a
=
r
c
Λ
τ
B
τ
SB
3c
BM
+ 2
,
(359)
and this holds for all t and all τ
B
/τ
SB
. We note that the more careful optimization in App. H allows us to
get a tighter bound, but the difference is most drastic for large τ
B
/τ
SB
when the bound is weaker than the
trivial kρ
BM,I
(t) − ρ
C,I
(t)k ≤ 1 + c
BM
.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 56
Note that for t ≤ T
a
/2 the expression in Eq. (359) grows exponentially with
p
τ
B
/τ
SB
. But we can have
a bound linear in
p
τ
B
/τ
SB
if we combine Eq. (359) with the trivial bound kρ
BM,I
(t) − ρ
C,I
(t)k ≤ 1 + c
BM
:
kρ
BM,I
(t) − ρ
C,I
(t)k ≤ min
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
Λ max(t,T
a
/2)
−
1 +
c
BM
(3c
BM
+ 2)
, 1 + c
BM
(360)
≤
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
Λt+1
− 6/5
. (361)
Proof. Here we prove Eq. (361). Note that the bound is monotonic in t, so the lowest value of τ
B
/τ
SB
=
∗2
for which the minimum is 2 for all t is given by equating the expressions at t = 0:
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
2T
a
/c
Λ
τ
SB
−
1 +
c
BM
(3c
BM
+ 2)
= 1 + c
BM
, (362a)
∗
q
c
Λ
(3c
BM
+ 2)
2e
2
∗
/
√
c
Λ
(3c
BM
+2)
−
1 +
c
BM
(3c
BM
+ 2)
= 1 + c
BM
(362b)
Relaxing this to an inequality, we find:
∗
q
c
Λ
(3c
BM
+ 2)
2 −
1 +
c
BM
(3c
BM
+ 2)
≤ 1 + c
BM
,
2
∗
√
c
Λ
√
3c
BM
+ 2
≤ 1, e
2
∗
/
√
c
Λ
(c
BM
+2)
≤ e
1/c
Λ
≤ e .
(363)
where we used c
Λ
≥ 1. Thus any bound on Eq. (359) for τ
B
/τ
SB
≤
∗2
that is monotonic in both τ
B
/τ
SB
and t will also be a bound for all τ
B
/τ
SB
. We choose to use the following in the exponent:
max(t, T
a
/2)
τ
SB
≤
t
τ
SB
+
1
2
s
c
Λ
τ
B
(3c
BM
+ 2)τ
SB
≤
t
τ
SB
+
∗
√
c
Λ
2
p
(3c
BM
+ 2)
, for
τ
B
τ
SB
≤
∗2
. (364)
Thus we arrive at
kρ
BM,I
(t) − ρ
C,I
(t)k
1
≤
s
c
Λ
(3c
BM
+ 2)τ
B
τ
SB
2e
2
∗
/
√
c
Λ
(3c
BM
+2)
e
Λt
−
1 +
c
BM
(3c
BM
+ 2)
,
which evaluates to the r.h.s of Eq. (361) using Eq. (363) and c
BM
≥ 1 ⇒ c
BM
/(3c
BM
+ 2) ≥ 1/5.
We will now redo this calculation for T
a
=
p
τ
B
τ
SB
/5 and c
Λ
= 1. Using this in Eq. (357):
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤
r
τ
B
5τ
SB
3e
Λt
− 1
c
BM
+
7e
Λt
− 5
. (365)
Now we find
∗
:
1 + c
BM
=
∗
√
5
3e
2
∗
/
√
5
− 1
c
BM
+
7e
2
∗
/
√
5
− 5
≥
∗
√
5
(2c
BM
+ 2), ⇒
2
∗
√
5
≤ 1 . (366)
This allows one to bound:
kρ
BM,I
(t) − ρ
C
√
5,I
(t)k
1
≤
r
τ
B
5τ
SB
3e
Λt+1
− 1
c
BM
+
7e
Λt+1
− 5
(367)
for all t.
H Optimal time
We begin with Eq. (357), and use it to write down the optimum t-dependent coarse-graining time:
T
a
(t) =
s
c
Λ
τ
B
τ
SB
(e
Λt
− 1)
(3c
BM
+ 2)e
Λt
− c
BM
. (368)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 57
In principle we could consider a family of equations with different T
a
(t) and solve each one of them just to
obtain ρ(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 kρ
BM,I
(t) − ρ
C,I
(t)k ≤ 1 + c
BM
:
1 + c
BM
(t
∗
) = ((3c
BM
(t
∗
) + 2)e
Λt
∗
− c
BM
(t
∗
))
T
a
τ
SB
+ (e
Λt
∗
− 1)
c
Λ
τ
B
T
a
. (369)
We would like to optimize T
a
at t = t
∗
, which would mean that t
∗
is maximized. The corresponding equation
for t
∗
is:
1 + c
BM
(t
∗
) = 2
r
c
Λ
τ
B
τ
SB
(e
Λt
∗
− 1)((3c
BM
(t
∗
) + 2)e
Λt
∗
− c
BM
(t
∗
)) , (370)
and solving the resulting quadratic equation yields:
e
Λt
∗
=
4c
BM
(t
∗
) + 2 + Y
2(3c
BM
(t
∗
) + 2)
, Y =
q
4(c
BM
(t
∗
) + 1)
2
+ 2(1 + c
BM
(t
∗
))
2
(3c
BM
(t
∗
) + 2)(τ
SB
/τ
B
c
Λ
) . (371)
Thus T
a
is:
T
a
=
s
c
Λ
τ
B
τ
SB
−2c
BM
(t
∗
) − 2 + Y
(3c
BM
(t
∗
) + 2)(2c
BM
(t
∗
) + 2 + Y )
. (372)
And the error is:
kρ
BM
(t) − ρ
C∗
(t)k
1
≤ ((3c
BM
(t
∗
) + 2)e
Λmax(t,T
a
/2)
− c
BM
(t
∗
))
s
c
Λ
τ
B
τ
SB
−2c
BM
(t
∗
) − 2 + Y
(3c
BM
(t
∗
) + 2)(2c
BM
(t
∗
) + 2 + Y )
(373a)
+ (e
Λmax(t,T
a
/2)
− 1)
s
c
Λ
τ
B
τ
SB
(3c
BM
(t
∗
) + 2)(2c
BM
(t
∗
) + 2 + Y )
−2c
BM
(t
∗
) − 2 + Y
(373b)
We should only use these expressions if t
∗
> T
a
/2. For τ
SB
τ
B
(when t
∗
is not likely to be > T
a
/2), this
changes the asymptotic form of T
a
to T
a
∼ τ
SB
. In the regime τ
SB
τ
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 ∼ τ
SB
.
I Proof of the bound (269b) for the Redfield limits of integration error
Starting from Eq. (268), we first send the integration limit to infinity:
Z
4t/τ
SB
0
min(1,
4τ
B
xτ
SB
+
T
)e
−x
dx ≤
Z
∞
0
min(1,
4τ
B
xτ
SB
+
T
)e
−x
dx =
T
+
Z
∞
0
min(1−
T
,
4τ
B
xτ
SB
)e
−x
dx . (374)
The min function evaluates to 1 −
T
when x ≤ x
∗
and to
4τ
B
xτ
SB
when x > x
∗
, where
x
∗
=
4τ
B
τ
SB
(1 −
T
)
. (375)
Therefore:
Z
∞
0
min(1 −
T
,
4τ
B
xτ
SB
)e
−x
dx =
Z
x
∗
0
(1 −
T
)e
−x
dx +
Z
∞
x
∗
4τ
B
xτ
SB
e
−x
dx (376a)
= (1 −
T
)(1 − e
−x
∗
) +
4τ
B
τ
SB
E
1
(x
∗
) , (376b)
where E
1
(x) denotes the exponential integral function: E
1
(x) ≡
R
∞
x
t
−1
e
−t
dt. A tight bound in terms of
elementary functions is given by [82][Thm. 2]:
E
1
(x) ≤ −ln(1 − e
−x
) (0 < x < ∞) . (377)
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 58
Thus
kρ
BM
(t) − ρ
R
(t)k
1
≤ c
BM
e
4t/τ
SB
T
+ (1 −
T
)(1 − e
−x
∗
) −
4τ
B
τ
SB
ln(1 −e
−x
∗
)
(378a)
≤ c
BM
e
4t/τ
SB
T
+
4τ
B
τ
SB
1 − ln
1 − e
−
4τ
B
τ
SB
(1−
T
)
, (378b)
where in the second line we used that ∀x: 1 − e
−x
< x.
Alternatively, and using more elementary methods, assuming x
∗
< 1, we obtain:
Z
∞
0
min(1,
4τ
B
xτ
SB
+
T
)e
−x
dx =
T
+ (1 −
T
)(1 − e
−x
∗
) +
Z
1
x
∗
4τ
B
xτ
SB
e
−x
dx +
Z
∞
1
4τ
B
xτ
SB
e
−x
dx (379a)
≤
T
+ (1 −
T
)x
∗
+
Z
1
x
∗
4τ
B
xτ
SB
dx +
Z
∞
1
4τ
B
τ
SB
e
−x
dx (379b)
=
T
+ (1 −
T
)x
∗
−
4τ
B
τ
SB
lnx
∗
+
4τ
B
eτ
SB
(379c)
=
T
+
4τ
B
τ
SB
+
4τ
B
τ
SB
ln
τ
SB
(1 −
T
)
4τ
B
+
4τ
B
eτ
SB
(x
∗
< 1) , (379d)
while when x
∗
≥ 1 the term with the logarithm is absent:
Z
∞
0
min(1,
4τ
B
xτ
SB
+
T
)e
−x
dx ≤
T
+
4τ
B
τ
SB
+
4τ
B
eτ
SB
(x
∗
≥ 1) . (380)
Combining these two cases directly yields Eq. (269b), which is an upper bound on the tighter bound presented
in Eq. (378b).
References
[1] R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications (Springer Science & Business
Media, 2007).
[2] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press,
Oxford, 2002).
[3] C.W. Gardiner and P. Zoller, Quantum Noise, Springer Series in Synergetics, Vol. 56 (Springer, Berlin,
2000).
[4] E. B. Davies, “Markovian master equations,” Communications in Mathematical Physics 39, 91–110
(1974).
[5] G. Lindblad, “On the generators of quantum dynamical semigroups,” Comm. Math. Phys. 48, 119–130
(1976).
[6] T. Albash, W. Vinci, A. Mishra, P. A. Warburton, and D. A. Lidar, “Consistency tests of classical and
quantum models for a quantum annealer,” Phys. Rev. A 91, 042314– (2015).
[7] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y.
Smirnov, M. Mohseni, and H. Neven, “Computational multiqubit tunnelling in programmable quantum
annealers,” Nat Commun 7 (2016).
[8] T. Albash, I. Hen, F. M. Spedalieri, and D. A. Lidar, “Reexamination of the evidence for entanglement
in a quantum annealer,” Physical Review A 92, 062328– (2015).
[9] A. Mishra, T. Albash, and D. A. Lidar, “Finite temperature quantum annealing solving exponentially
small gap problem with non-monotonic success probability,” Nature Communications 9, 2917 (2018).
[10] R. Feynman and F. Vernon, “The theory of a general quantum system interacting with a linear dissipative
system,” Annals of Physics 24, 118 – 173 (1963).
[11] A. Caldeira and A. Leggett, “Quantum tunnelling in a dissipative system,” Annals of Physics 149, 374
– 456 (1983).
[12] D. E. Makarov and N. Makri, “Path integrals for dissipative systems by tensor multiplication. condensed
phase quantum dynamics for arbitrarily long time,” Chemical Physics Letters 221, 482 – 491 (1994).
[13] E. Sim, “Quantum dynamics for a system coupled to slow baths: On-the-fly filtered propagator method,”
The Journal of Chemical Physics 115, 4450–4456 (2001).
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 59
[14] K. Kraus, States, Effects, and Operations (Springer, Berlin, 1983).
[15] J. M. Dominy and D. A. Lidar, “Beyond complete positivity,” Quant. Inf. Proc. 15, 1349 (2016).
[16] V. Gorini, A. Frigerio, M. Verri, A. Kossakowski, and E. C. G. Sudarshan, “Properties of quantum
Markovian master equations,” Reports on Mathematical Physics 13, 149–173 (1978).
[17] S. Nakajima, “On Quantum Theory of Transport Phenomena : Steady Diffusion,” Prog. Theor. Phys.
20, 948 (1958).
[18] R. Zwanzig, “Ensemble Method in the Theory of Irreversibility,” J. Chem. Phys. 33, 1338 (1960).
[19] A. Redfield, “The theory of relaxation processes,” in Advances in Magnetic Resonance, Advances in
Magnetic and Optical Resonance, Vol. 1, edited by J. S. Waugh (Academic Press, 1965) pp. 1 – 32.
[20] D. Bacon, D. A. Lidar, and K. B. Whaley, “Robustness of decoherence-free subspaces for quantum
computation,” Phys. Rev. A 60, 1944–1955 (1999).
[21] A. J. van Wonderen and K. Lendi, “Virtues and limitations of markovian master equations with a
time-dependent generator,” J. Stat. Phys. 100, 633–658 (2000).
[22] D. A. Lidar, Z. Bihary, and K. Whaley, “From completely positive maps to the quantum Markovian
semigroup master equation,” Chem. Phys. 268, 35 (2001).
[23] S. Daffer, K. Wodkiewicz, J.D. Cresser, J.K. McIver, “Depolarizing channel as a completely positive
map with memory,” Phys. Rev. A 70, 010304(R) (2004).
[24] A. Shabani and D. A. Lidar, “Completely positive post-markovian master equation via a measurement
approach,” Phys. Rev. A 71, 020101(R) (2005).
[25] S. Maniscalco and F. Petruccione, “Non-Markovian dynamics of a qubit,” Phys. Rev. A 73, 012111
(2006).
[26] J. Piilo, S. Maniscalco, K. Härkönen, and K.-A. Suominen, “Non-markovian quantum jumps,” Physical
Review Letters 100, 180402– (2008).
[27] H.-P. Breuer and B. Vacchini, “Quantum semi-markov processes,” Physical Review Letters 101, 140402–
(2008).
[28] R. S. Whitney, “Staying positive: going beyond lindblad with perturbative master equations,” Journal
of Physics A: Mathematical and Theoretical 41, 175304 (2008).
[29] L.-A. Wu, G. Kurizki, and P. Brumer, “Master equation and control of an open quantum system with
leakage,” Physical Review Letters 102, 080405– (2009).
[30] D. Chruściński and A. Kossakowski, “Non-markovian quantum dynamics: Local versus nonlocal,” Phys.
Rev. Lett. 104, 070406 (2010).
[31] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, “Quantum adiabatic Markovian master equations,”
New J. of Phys. 14, 123016 (2012).
[32] E. Mozgunov, “Local master equation for small temperatures,” arXiv:1611.04188 (2016).
[33] A. Y. Smirnov and M. H. Amin, “Theory of open quantum dynamics with hybrid noise,” New Journal
of Physics 20, 103037 (2018).
[34] R. Dann, A. Levy, and R. Kosloff, “Time-dependent markovian quantum master equation,” Phys. Rev.
A 98, 052129 (2018).
[35] L. C. Venuti and D. A. Lidar, “Error reduction in quantum annealing using boundary cancellation: Only
the end matters,” Phys. Rev. A 98, 022315 (2018).
[36] G. McCauley, B. Cruikshank, D. I. Bondar, and K. Jacobs, “Completely positive, accurate master
equation for weakly-damped quantum systems across all regimes,” arXiv:1906.08279 (2019).
[37] F. Benatti, R. Floreanini, and U. Marzolino, “Environment-induced entanglement in a refined weak-
coupling limit,” EPL (Europhysics Letters) 88, 20011 (2009).
[38] F. Benatti, R. Floreanini, and U. Marzolino, “Entangling two unequal atoms through a common bath,”
Phys. Rev. A 81, 012105 (2010).
[39] M. Merkli, “Quantum markovian master equations: Resonance theory shows validity for all time scales,”
Annals of Physics 412, 167996 (2020).
[40] C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, “Coarse graining can beat the rotating-wave
approximation in quantum markovian master equations,” Phys. Rev. A 88, 012103– (2013).
[41] T. S. Cubitt, A. Lucia, S. Michalakis, and D. Perez-Garcia, “Stability of local quantum dissipative
systems,” Communications in Mathematical Physics 337, 1275–1315 (2015).
[42] E. Knill, “Quantum computing with realistically noisy devices,” Nature 434, 39–44 (2005).
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 60
[43] R. Alicki, D. A. Lidar, and P. Zanardi, “Internal consistency of fault-tolerant quantum error correction
in light of rigorous derivations of the quantum markovian limit,” Phys. Rev. A 73, 052311 (2006).
[44] D. A. Lidar, “Lecture notes on the theory of open quantum systems,” arXiv preprint arXiv:1902.00967
(2019).
[45] T. Albash and D. A. Lidar, “Decoherence in adiabatic quantum computation,” Phys. Rev. A 91, 062320–
(2015).
[46] M. Žnidarič, “Dephasing-induced diffusive transport in the anisotropic heisenberg model,” New Journal
of Physics 12, 043001 (2010).
[47] M. V. Medvedyeva, T. Prosen, and M. Žnidarič, “Influence of dephasing on many-body localization,”
Phys. Rev. B 93, 094205 (2016).
[48] R. Bhatia, Matrix Analysis, Graduate Texts in Mathematics No. 169 (Springer-Verlag, New York, 1997).
[49] P. Gaspard and M. Nagaoka, “Slippage of initial conditions for the redfield master equation,” Journal of
Chemical Physics 111, 5668–5675 (1999).
[50] G. Vidal, “Efficient classical simulation of slightly entangled quantum computations,” Phys. Rev. Lett.
91, 147902 (2003).
[51] F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, “Matrix product density operators: Simulation of
finite-temperature and dissipative systems,” Phys. Rev. Lett. 93, 207204 (2004).
[52] E. H. Lieb and D.W. Robinson, “The finite group velocity of quantum spin systems,” Commun. Math.
Phys. 28, 251 (1972).
[53] J. Haah, M. Hastings, R. Kothari, and G. H. Low, “Quantum algorithm for simulating real time evolution
of lattice hamiltonians,” in 2018 IEEE 59th Annual Symposium on Foundations of Computer Science
(FOCS) (2018) pp. 350–360.
[54] H. Pichler, A. J. Daley, and P. Zoller, “Nonequilibrium dynamics of bosonic atoms in optical lattices:
Decoherence of many-body states due to spontaneous emission,” Phys. Rev. A 82, 063605 (2010).
[55] L.-M Duan and G.-C. Guo, “Reducing decoherence in quantum-computer memory with all quantum bits
coupling to the same environment,” Phys. Rev. A 57, 737 (1998).
[56] P. Zanardi and M. Rasetti, “Noiseless quantum codes,” Phys. Rev. Lett. 79, 3306–3309 (1997).
[57] D. A. Lidar, I. L. Chuang, and K. B. Whaley, “Decoherence-free subspaces for quantum computation,”
Phys. Rev. Lett. 81, 2594–2597 (1998).
[58] D. A. Lidar and K. B. Whaley, “Decoherence-free subspaces and subsystems,” in Irreversible Quantum
Dynamics, Lecture Notes in Physics, Vol. 622, edited by F. Benatti and R. Floreanini (Springer, Berlin,
2003) p. 83.
[59] P.G. Kwiat, A.J. Berglund, J.B. Altepeter, and A.G. White, “Experimental Verification of Decoherence-
Free Subspaces,” Science 290, 498 (2000).
[60] L. Viola, E. M. Fortunato, M. A. Pravia, E. Knill, R. Laflamme, and D. G. Cory, “Experimental
realization of noiseless subsystems for quantum information processing,” Science 293, 2059–2063 (2001).
[61] D. Kielpinski, V. Meyer, M. A. Rowe, C. A. Sackett, W. M. Itano, C. Monroe, and D. J. Wineland, “A
decoherence-free quantum memory using trapped ions,” Science 291, 1013–1015 (2001).
[62] J. E. Ollerenshaw, D. A. Lidar, and L. E. Kay, “Magnetic resonance realization of decoherence-free
quantum computation,” Phys. Rev. Lett. 91, 217904 (2003).
[63] L. Viola and S. Lloyd, “Dynamical suppression of decoherence in two-state quantum systems,” Phys.
Rev. A 58, 2733–2744 (1998).
[64] P. Zanardi, “Symmetrizing evolutions,” Physics Letters A 258, 77–82 (1999).
[65] D. Lidar and T. Brun, eds., Quantum Error Correction (Cambridge University Press, Cambridge, UK,
2013).
[66] D. Suter and G. A. Álvarez, “Colloquium: Protecting quantum information against environmental noise,”
Reviews of Modern Physics 88, 041001– (2016).
[67] H. K. Ng, D. A. Lidar, and J. Preskill, “Combining dynamical decoupling with fault-tolerant quantum
computation,” Phys. Rev. A 84, 012305– (2011).
[68] K. Szczygielski and R. Alicki, “Markovian theory of dynamical decoupling by periodic control,” Physical
Review A 92, 022349– (2015).
[69] J. E. Gough and H. I. Nurdin, “Can quantum markov evolutions ever be dynamically decoupled?” in
2017 IEEE 56th Annual Conference on Decision and Control (CDC) (2017) pp. 6155–6160.
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 61
[70] C. Addis, F. Ciccarello, M. Cascio, G. M. Palma, and S. Maniscalco, “Dynamical decoupling efficiency
versus quantum non-markovianity,” New Journal of Physics 17, 123004 (2015).
[71] C. Arenz, D. Burgarth, P. Facchi, and R. Hillier, “Dynamical decoupling of unbounded hamiltonians,”
Journal of Mathematical Physics, Journal of Mathematical Physics 59, 032203 (2018).
[72] L. Li, M. J. W. Hall, and H. M. Wiseman, “Concepts of quantum non-markovianity: A hierarchy,”
Concepts of quantum non-Markovianity: A hierarchy, Physics Reports 759, 1–51 (2018).
[73] I. de Vega, M. C. Bañuls, and A. Pérez, “Effects of dissipation on an adiabatic quantum search algo-
rithm,” New J. of Phys. 12, 123010 (2010).
[74] https://github.com/mvjenia/CGMEcode, code for the numerical section of the paper.
[75] L. Isserlis, “On certain probable errors and correlation coefficients of multiple frequency distributions
with skew regression,” Biometrika 11, 185 (1916).
[76] T. Albash, D. A. Lidar, M. Marvian, and P. Zanardi, “Fluctuation theorems for quantum processes,”
Phys. Rev. E 88, 032146– (2013).
[77] T. Albash and D. A. Lidar, “Adiabatic quantum computation,” Reviews of Modern Physics 90, 015002–
(2018).
[78] R. Alicki, M. Fannes, and M. Horodecki, “A statistical mechanics view on kitaev’s proposal for quantum
memories,” Journal of Physics A: Mathematical and Theoretical 40, 6451–6467 (2007).
[79] H. Bombin, “Topological subsystem codes,” Physical Review A 81, 032301– (2010).
[80] B. Altshuler, H. Krovi, and J. Roland, “Anderson localization makes adiabatic quantum optimization
fail,” Proceedings of the National Academy of Sciences 107, 12446–12450 (2010).
[81] M. Reed and B. Simon, Methods of Modern Mathematical Physics: Fourier analysis, self-adjointness,
Vol. 2 (Academic Press, 1975).
[82] H. Alzer, “On some inequalities for the incomplete gamma function,” Mathematics of Computation 66,
771 (1997).
Accepted in Quantum 2020-01-19, click title to verify. Published under CC-BY 4.0. 62