Completely Positive, Simple, and Possibly Highly Accurate
Approximation of the Redfield Equation
Dragomir Davidović
Georgia Institute of Technology, Atlanta, Georgia, United States
Here we present a Lindblad master equation that approximates the Redfield equation, a well
known master equation derived from first principles, without significantly compromising the range
of applicability of the Redfield equation. Instead of full-scale coarse-graining, this approximation
only truncates terms in the Redfield equation that average out over a time-scale typical of the
quantum system. The first step in this approximation is to properly renormalize the system
Hamiltonian, to symmetrize the gains and losses of the state due to the environmental coupling.
In the second step, we swap out an arithmetic mean of the spectral density with a geometric
one, in these gains and losses, thereby restoring complete positivity. This completely positive
approximation, GAME (geometric-arithmetic master equation), is adaptable between its time-
independent, time-dependent, and Floquet form. In the exactly solvable, three-level, Jaynes-
Cummings model, we find that the error of the approximate state is almost an order of magnitude
lower than that obtained by solving the coarse-grained stochastic master equation. As a test-bed,
we use a ferromagnetic Heisenberg spin-chain with long-range dipole-dipole coupling between up
to 25-spins, and study the differences between various master equations. We find that GAME
has the highest accuracy per computational resource.
1 Introduction
Quantum correlations between particles in quantum systems are counterintuitive, giving us exotic algorithms
that can greatly accelerate some calculations using quantum computing as opposed to conventional comput-
ing. It remains an unsettled issue, however, to describe correlated quantum dynamics in large many-body
quantum systems, with the precision required for quantum computing, if the systems are coupled to an envi-
ronment. In the Born-Markov approximation, the time dependence of the reduced state of the system can be
approximated by a first-order differential equation, (e.g., the master equation), which greatly simplifies this
description. [1, 2] The Born-Markov approximation is justified if the time in which the environment changes
the system state in the interaction picture, τ
r
, is much longer than the time at which the correlations in
the environment decay, τ
c
. [2] We will use the terms bath and environment interchangeably throughout this
paper.
The master equation does not necessarily preserve the positivity of the reduced quantum state of the
system. There is only one form of the master equation that is a completely positive map, the deeply respected
Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation. [24] The GKSL equation is a consequence of
the Dirac–von Neumann axioms, while the physical scope of the equation rests in the coarse-graining of the
reduced system dynamics over a timescale much longer than τ
c
. [5]
It has been very difficult to derive the GKSL equation from first principle calculations using the total
Hamiltonian of the system and environment as a starting point. [6] The operator-projection technique has
been the standard way to derive quantum master equations [2, 79] initially pioneered by Redfield. [10, 11]
The Redfield equation, however, suffers from a lack of positivity preservation leading to negative probabilities
of observables [6] which violates the mathematical definition of probability. Negative states also pose problems
in terms of how to quantify the entanglement of the system, which is important for quantum computing.
One measure of entanglement is the negativity of the state’s partial transpose, [12] which raises the question:
Dragomir Davidović: dragomir.davidovic@physics.gatech.edu
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 1
arXiv:2003.09063v4 [quant-ph] 17 Sep 2020
Can we trust negativity of the partial transpose as a genuine measure of entanglement if the approximate
state itself is negative?
Although it has this problem of producing negative probabilities, the Redfield equation has major ad-
vantages as well. It is a highly accurate master equation that is simple to apply and not axiomatic. Thus,
the Redfield equation and its modifications have found significant use in quantum chemistry, condensed
matter physics, and quantum optics. [6, 1320] The issue of negative states and their relevance remains de-
bated. [21, 22] The Redfield equation with time dependent coefficients (TDC) has reduced negativity. The
nonpositivity of that equation can even be seen as a benefit of a red-flag for the breakdown of the master
equation concept. [22] This benefit does not extend well to the Redfield equation with asymptotic coefficients
as we will discuss in this paper. The latter is the widely used master equation that we are approximating
into a GKSL form in this paper.
Our main focus here is to find a simple, completely positive approximation of the Redfield equation
without significantly compromising its accuracy. Note that in the same weak coupling regime where the
Redfield equation is valid, quantum trajectory description is available which guarantees positivity. [23, 24]
There are also quite a few non-perturbative approaches such as the exact quasi adiabatic path integral
method (QUAPI), [2527] the hierarchical equations of motion (HEOM), [2832], the multiconfiguration
time-dependent Hartree (MCTDH) [33, 34], and the multilayer formulation ML-MCTDH. [35, 36] These are
powerful approaches but can be very challenging. Exact methods aiming at improved efficiency include quan-
tum trajectory based hierarchy of stochastic pure states (HOPS) [3739] and tensor network methods. [40, 41]
Other strong coupling approaches involve unitary transformations that result in a non-perturbative bath-
renormalization of the system Hamiltonian, followed by a quantum master equation in the weak coupling
limit after the renormalization. [42, 43]
There have been several attempts to cure the issue of negative states by modifying the Redfield master
equation into a GKSL form. This was done for the first time by Davies, [44] who gave us the rotating-wave
approximation (RWA), also known as the secular approximation. The RWA amounts to coarse-graining of
the reduced state on a timescale comparable to the Heisenberg time. Since this time scales exponentially
with system size, the RWA does not work well for large many-body quantum systems that generally have
exponentially suppressed level spacings with increasing size of the system. More recently, coarse-grained
GKSL master equations have been derived from first principles [4547] capable of capturing correlation effects
that the RWA cannot. Another notable completely positive approximation is the partial-RWA (PRWA). [48,
49] An alternative approach to restoring complete positivity includes the dynamically coarse-grained (DCG)
master equation. [50, 51] Phenomenological models have also worked well to briskly set up a completely
positive master equation. [5255]
The approximate, completely positive GAME works in two steps. First, it identifies an implicit renormal-
ization of the Hamiltonian (e.g., the Lamb-shift) in the Redfield equation. After separating the Lamb-shift
by adding it to the system Hamiltonian, the second step separates the dissipator into two parts: a completely
positive dissipator and the remainder. The latter has the property that it averages-out on a short time-scale
compared to the system relaxation time, and is therefore dropped, while the renormalized Hamiltonian and
the completely positive dissipator remain intact.
Practically, in terms of the spectral density (SD), we swap out an arithmetic mean of the SD at two
system frequencies, with a geometric mean of the SD. The difference between the two means has the sought
after property that it averages out on a typical time scale of the system. We find that the state error added
by the approximation is linear in system-bath coupling and generally comparable to the accuracy of the
Redfield equation (when applicable).
GAME belongs to a relatively new group of approximations that we refer to here as the
SD-approximations.
[55, 56] In this paper we compare many of the aforementioned GKSL-equations with GAME, by applying
them on a spin-1/2 Heisenberg chain, and have not found a single one that outperforms GAME in terms of
the closeness of its solutions to those of the Redfield equation.
The renormalized Hamiltonian plays the key role and is the primary reason that GAME is capable of
accounting for the correlation effects that the bath introduces into the system dynamics. The effects of the
renormalization are studied in detail in the case of the 3-level Jaynes-Cummings model, where the eigenstates
of the renormalized Hamiltonian are identified as the relaxation modes of the system coupled to the heat
bath, those that decay with a single relaxation time.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 2
By construction, GAME retains the simplicity of the Redfield equation, which is one of its assets. Even
more important is its combination of both simplicity and accuracy.
The paper is organized as follows. In Sec. 2 we introduce the notation relevant for this work. In Sec. 3, 3.1,
and 3.2 we review the Redfield, Coarse-Grained Redfield, and the Davies-Lindblad master equations, respec-
tively, as a prelude to GAME, which is presented in Sec. 4. In Sec. 5 we implement GAME on an exactly
solvable 3-level Jaynes–Cummings model. We follow by the detailed investigation of a 25-body spin-1/2 fer-
romagnetic spin-chain in Sec. 6. We end the analysis by presenting a detailed comparison with other GKSL
equations in Sec. 7, followed by discussion and conclusion.
2 Introducing the problem
Perturbative master equations have been derived many times in the literature. For a tutorial and some recent
examples we refer the reader to Refs. [21, 22, 47, 57].
Here we consider a quantum system coupled to a heat-bath, represented by the total Hamiltonian
H
tot
= H
0
1
B
+ 1
S
H
B
+ A B. (1)
H
0
and A are the respective Hamiltonian and arbitrary operators performing in the Hilbert space of the
quantum system, H
B
and B are the corresponding operators for the heat bath, and A B is the interaction
Hamiltonian between the system and the heat bath. All the operators are assumed to be Hermitian. For
notational clarity we only consider one coupling term in the interaction Hamiltonian. The results can be
straightforwardly extended to a sum of interaction terms assuming uncorrelated heat baths.
In the Hilbert space of the quantum system, we use the eigenbasis of H
0
, H
0
|ni = E
n
|ni, where |ni
are the eigenvectors and E
n
are the eigenenergies. Then, in the interaction picture A(t) = e
iH
0
t
Ae
iH
0
t
is
decomposed as
A(t) =
X
n,m
A
nm
e
mn
t
|nihm|, (2)
where A
nm
= hn|A|mi, ω
mn
= E
m
E
n
, t is time, and we use ~ = 1. Operator A(t) can be represented in
the matrix form as
A(t) = A Ω(t), (3)
where is the Hadamard product (the element-wise product) and [Ω(t)]
nm
= e
mn
t
.
The bath correlation function is defined as C(t) = T r[ρ
B
B(t)B(0)], where B(t) = e
iH
B
t
Be
iH
B
t
and
ρ
B
is the reduced density matrix of the heat bath. We work under the premise that the excitations in the
heat-bath decay on a time-scale much smaller than the time scale of the system state in the interaction
picture. In that case, the change in the reduced density matrix of the bath due to the coupling to the system
can be neglected, in the leading order of the approximation, and the density matrix of the environment can
be approximated by the thermal state where C(t) = C
?
(t). [2] In our notation we reserve the symbols ρ
and %, for the density matrix in the Schrödinger and interaction picture, respectively.
3 Redfield equation
The equation is a widely used fundamental master equation describing fluctuations and dissipation in quantal
systems. [9] It is derived from first principles, by applying the Born-Markov approximation and projection
operator technique. We begin here from Eq. 16 from a recent re-derivation, [47]
dt
= i[H
0
, ρ] AA
f
ρ ρA
f
A + A
f
ρA + AρA
f
, (4)
where
A
f
=
Z
0
C(τ)A(τ) A Γ
?
(5)
is termed the "filtered" operator, and adopt the notation from that reference throughout this paper. The
star symbol indicates complex-conjugate. The matrix elements of Γ can be expressed as
Γ
nm
=
Z
0
C(τ)e
nm
τ
=
1
2
γ
nm
+ iS
nm
. (6)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 3
C(t) γ(ω) S(ω)
1
2π
R
−∞
e
it
γ(Ω)d
R
−∞
e
t
C(t)dt
1
2π
P
R
−∞
γ(Ω)
ω
d
gω
2
c
(1+
c
t)
2
2πgωe
ω
ω
c
Θ(
ω
ω
c
) gω
c
h
1
ω
ω
c
e
ω
ω
c
Ei(
ω
ω
c
)
i
gω
2
c
[ln |ω
c
t| + γ
+i
π
2
sgn(t)] + O[t
2
ln(t)]
2πgω
ω
2
c
ω
2
c
+ω
2
Θ(
ω
ω
c
) gω
c
h
π/2(ω/ω
c
) ln(|ω|
c
)
1+(ω
c
)
2
i
6gω
2
c
(1+
c
t)
4
2πg
ω
3
ω
2
c
e
ω
ω
c
Θ(
ω
ω
c
)
gω
c
h
2 +
ω
ω
c
+ (
ω
ω
c
)
2
(
ω
ω
c
)
3
e
ω
ω
c
Ei(
ω
ω
c
)
i
.
Table 1: Integral transforms between bath correlation function [C(t)], and spectral functions γ(ω) (spectral density) and
S(ω) (principal density) at zero temperature. P indicates the principal value of the integral. The second and third row
display these functions in case of an Ohmic bath with exponential and Drude-Lorentz cutoff at frequency ω
c
. Here, Θ(x)
is the Heaviside step function, Ei(x) is the exponential integral
R
x
−∞
e
t
t
dt, g is the dimensionless system-bath coupling
constant, and γ is the Euler–Mascheroni constant. The forth row exhibits the super-Ohmic spectral density.
In this work we also use functions Γ(x) =
1
2
γ(x) + iS(x), so that Γ
nm
= Γ(ω
nm
), γ
nm
= γ(ω
nm
), and
S
nm
= S(ω
nm
). Here γ(x) is the spectral density (SD), equal to the Fourier transform of the bath correlation
function, and S(x) is the principal density (PD), equal to the imaginary part of the half-range Fourier
transform of the bath correlation function. The symbols Γ, γ, and S mean matrices with respective matrix
elements Γ
nm
, γ
nm
and S
nm
. Table 1 summarizes key formulae for the bath correlation functions and its
integral transforms, and three examples of heat baths at zero temperature, that we apply in this paper.
Eq. 4 can be recast in a different form,
dt
= i[H
0
, ρ]
1
2
[AA
f
A
f
A, ρ]
1
2
{AA
f
+ A
f
A, ρ}+ A
f
ρA + AρA
f
. (7)
With the renormalization of the system Hamiltonian,
H = H
0
i
2
(AA
f
A
f
A) H
0
+ H
L
, (8)
the Redfield equation becomes
dt
+ i[H, ρ] =
1
2
{AA
f
+ A
f
A, ρ}+ AρA
f
+ A
f
ρA, (9)
or in terms of matrix elements,
nm
dt
+ i[H, ρ]
nm
=
X
ij
A
?
in
ρ
ij
A
jm
γ
in
+ γ
jm
2
i(S
jm
S
in
)
1
2
X
i,j
ρ
ni
A
?
mj
A
ij
γ
ij
+ γ
mj
2
i(S
ij
S
mj
)
1
2
X
i,j
A
?
ji
A
ni
ρ
jm
γ
ni
+ γ
ji
2
i(S
ni
S
ji
)
. (10)
The matrix elements of the renormalized Hamiltonian 8 are
H
nm
= E
n
δ
nm
+
X
i
H(ω
ni
, ω
mi
)A
ni
A
im
, (11)
where we introduce the kernel for the unitary component of the reduced system dynamics (from now on the
unitary kernel):
H(ω, ω
0
) =
1
2
S(ω) + S(ω
0
) + i
γ(ω) γ(ω
0
)
2
. (12)
The role of the renormalization is to rewrite the loss terms in the Redfield equation so that all three
summands in Eq. 10 involve this one and only function of system frequencies,
G(ω, ω
0
) =
1
2
γ(ω) + γ(ω
0
)
i
S(ω
0
) S(ω)
, (13)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 4
that we refer to here as the dissipative kernel. That is to say, we reformulate the Redfield equation as
nm
dt
+ i[H, ρ]
nm
=
X
ij
A
?
in
ρ
ij
A
jm
G
in,jm
1
2
X
i,j
ρ
ni
A
ij
A
?
mj
G
mj,ij
1
2
X
i,j
A
ni
A
?
ji
ρ
jm
G
ji,ni
, (14)
with one simple kernel
G
in,jm
= G(ω
in
, ω
jm
) = Γ
in
+ Γ
?
jm
. (15)
The importance of this form of the Redfield equation is that we can approximate the kernel only. The
approximation will then carry over to the state gains and losses unbiasedly, maintaining the balance between
them. Furthermore, if the kernel in Eq. 15 were positive-semidefinite, it would lead to a GKSL-equation.
See for example Ref. [50]. So restoring positivity of the Redfield equation now amounts to approximating its
kernel with a positive-semidefinite one.
As written above, the Redfield equation has no dependence on the initial time and the history of the state
is fully encoded in the present state. In the derivation of the Redfield equation, before the last approximation
is applied, there is a master equation with time-dependent coefficients (TDCs), that depend on the choice of
the initial time. The filtered operator in that case has time dependence
A
f
(t) =
Z
t
0
C(τ)A(τ) A Γ
?
(t), (16)
which leads to the time-dependent SDs and PDs vis-à-vis
Γ
t
(ω
nm
) =
1
2
γ
t
(ω
nm
) + iS
t
(ω
nm
), (17)
while the master equation is
dt
= i[H
0
, ρ] AA
f
(t)ρ ρA
f
(t)A + A
f
(t)ρA + AρA
f
(t). (18)
This equation depends on the initial time, but, as t becomes much longer than τ
c
, the coefficients Γ
nm
(t)
approach the asymptotic value, and the equation no longer depends on the history of the system. [2] We also
call this the
R
t
R
approximation.
In this paper we refer to equations 18 and 4 as the TDC-Redfield equation and the Redfield equation,
respectively. The Redfield equation does not resolve the system dynamics over a time scale of the bath
correlation time. [2]
A time-dependent Redfield equation in form equivalent to 9 has been studied recently in context of
embedding of the quantum system into an expanded system that displays completely positive quantum
dynamics by Breuer. [58] The density matrix of the subsystem, after the embedding, is presented by the
off-diagonal block of the expanded state, and such embedding does not resolve negativity of the reduced
state. The renormalization of the Hamiltonian is mentioned, but not written down explicitly.
The range of applicability of the Redfield equation is given by the condition of the validity of the Born-
Markov approximation. The failure of the Born-Markov approximation can be identified by comparing the
second and fourth order terms in the Dyson expansion of the system state, [59, 60] which leads to the condition
τ
r
τ
c
. The Born-Markov approximation is justified if τ
r
τ
c
.
For any given initial state of the system, the error bound of the Redfield equation can be expressed in
terms of the properties of the heat bath only, [47]
1/2||ρ
RED
(t) ρ
exact
(t)||
1
O
τ
b
τ
sb
e
12t
τ
sb
ln
τ
sb
τ
b
. (19)
Here ρ
exact
(t) is the unapproximated quantum state, τ
b
and τ
sb
are the bath correlation time and the smallest
system relaxation time, respectively, and ||...||
1
is the trace-distance. τ
sb
is roughly τ
b
/g. Ignoring the
logarithmic correction, the error bound is linear in system-bath coupling constant g. Recent simulations on
a two qubit-system find linear scaling of the actual error (as opposed to the error bound) with g. [22]
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 5
In the numerical simulations in this paper, the state relaxation time will be significantly longer than
τ
sb
, as it is in many cases in quantum optics and atomic physics. [61, 62] In that case the error bound 19 is
exponentially large at the system relaxation time. Only if the system relaxation time is comparable to τ
sb
the
bound will be tight enough to present a good error estimate. So here we adopt a system-dependent criterion
to determine the range of applicability of the Redfield master equation. If, for example, we find that the
Redfield states become too negative, we consider it as a sign that we are outside the range of applicability. A
similar point of view has been expressed in Ref. [22]. In addition, if the Redfield equation has an instability,
which it usually does at large system-bath coupling, then it will clearly not be applicable.
3.1 Coarse-grained Redfield Equation
Here we coarse-grain Eq. 14, as follows. First, we transform to the interaction picture [%
nm
= exp(
nm
t)ρ
nm
],
where the Redfield equation reads
d%
nm
dt
= i[e
iH
0
t
H
L
e
iH
0
t
, %]
nm
+
X
ij
A
?
in
%
ij
A
jm
e
i(ω
jm
ω
in
)t
G(ω
in
, ω
jm
)
1
2
X
i,j
h
%
ni
A
ij
A
?
mj
e
im
t
G(ω
mj
, ω
ij
) + A
ni
A
?
ji
%
jm
e
nj
t
G(ω
ji
, ω
ni
)
i
. (20)
Next, we change the time symbol from t to τ, and time average on a rectangular window. That is, we apply
the integral operator to the equation,
ˆ
G =
1
T
0
Z
t+T
0
/2
tT
0
/2
(21)
while keeping % at time t. We will refer to T
0
as the coarse-graining time. This leads to
d%
nm
dt
+ i[
˜
H
L
(t), %]
nm
=
X
ij
A
?
in
%
ij
A
jm
e
i(ω
jm
ω
in
)t
˜
G(ω
in
, ω
jm
)
1
2
X
i,j
h
%
ni
A
ij
A
?
mj
e
im
t
˜
G(ω
mj
, ω
ij
) + A
ni
A
?
ji
%
jm
e
nj
t
˜
G(ω
ji
, ω
ni
)
i
, (22)
where
˜
G(ω, ω
0
) =
1
2
γ(ω) + γ(ω
0
)
+ i
S(ω) S(ω
0
)
sinc
(ω ω
0
)T
0
2
. (23)
Here sinc(x) = [sin(x)]/x, and
˜
H
L
(t) =
ˆ
G e
iH
0
τ
H
L
e
iH
0
τ
. (24)
In the third step, we transform out of the interaction picture, [ρ
nm
= exp(
nm
t)%
nm
], and obtain the
coarse-grained Redfield equation,
nm
dt
+ i[
˜
H, ρ]
nm
=
X
ij
A
?
in
ρ
ij
A
jm
˜
G(ω
in
, ω
jm
)
1
2
X
i,j
h
ρ
ni
A
ij
A
?
mj
˜
G(ω
mj
, ω
ij
) + A
ni
A
?
ji
ρ
jm
˜
G(ω
ji
, ω
ni
)
i
. (25)
Similarly, the matrix elements of the coarse-grained renormalized Hamiltonian are
˜
H
nm
E
n
δ
nm
+
˜
H
L,nm
= E
n
δ
nm
+
X
i
˜
H(ω
ni
, ω
mi
)A
ni
A
im
, (26)
where
˜
H(ω, ω
0
) = H(ω, ω
0
)sinc
(ω ω
0
)T
0
2
. (27)
Coarse-graining increases the error of the approximate state. Error bound exists for the state error added
by coarse-graining, and increases linearly with the coarsegraining time. [47]
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 6
3.2 Rotating wave and CGSE approximations
The RWA follows from Eqs. 23, 25 and 26 in the limit T
0
1/|ω
nm
|, ω
nm
6= 0. In that limit the sinc-function
becomes the Kronecker-delta thereby truncating the coarse-grained Redfield equation. The truncation causes
cancellations of the principle densities S on the RHS of Eq. 10, leading to the historically important Davies-
Lindblad master equation:
nm
dt
+ i[
˜
˜
H, ρ]
nm
=
X
i,j
δ
ω
in
jm
γ
in
A
?
in
ρ
ij
A
jm
1
2
X
ij
γ
ij
(ρ
ni
A
ij
A
jm
δ
E
i
,E
m
+ A
nj
A
ji
δ
E
n
,E
i
ρ
im
) (28)
where
˜
˜
H is the system Hamiltonian renormalized due to the environmental coupling,
˜
˜
H
nm
= E
n
δ
n,m
+
X
i
˜
˜
H(ω
ni
, ω
mi
)A
ni
A
im
, (29)
which is also known as the Lamb-shift. The unitary kernel in the RWA is
˜
˜
H(ω, ω
0
) =
δ
ω
0
2
S(ω) + S(ω
0
)
. (30)
In the interaction picture, the Davies-Lindblad master equation has time independent coefficients. The
equation and the Lamb-shift simplify further, respectively as
nm
dt
= i
nm
ρ
nm
+
X
i,j
δ
ω
in
jm
γ
in
A
?
in
ρ
ij
A
jm
1
2
ρ
nm
X
i
|A
ni
|
2
γ
ni
+ |A
mi
|
2
γ
mi
, (31)
and
n
= E
n
+
X
i
|A
n,i
|
2
S
n,i
, (32)
if the energy levels of the system are non-degenerate.
RWA is valid when the coarse-graining time is longer than the Heisenberg time τ
h
= 1E, where δE is
the smallest level spacing in the system. In that respect, the system relaxation time also needs to be of that
order or longer, so that the coarse-grained equation can be of value. Therefore, the main disadvantage to
the RWA is in its highly limited range of applicability. Note that at zero temperature, the dissipative part
of the system dynamics under the RWA only involves positive frequencies that correspond to the emission of
quanta from the system into the bath (since the SD at negative frequency is zero). [61]
More recently, a coarse-grained stochastic equation was derived [45] with better error and range of appli-
cability than that of Davies’,
1/2||ρ
cgse
(t) ρ
exact
(t)||
1
O
r
τ
c
τ
sb
, (33)
but not as low as the Redfield error bound given by Eq. 19.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 7
4 Geometric Arithmetic Master Equation (GAME)
Since coarse-graining increases the approximate state error, as discussed at the end of Sec. 3.1, here we will
attempt to restore positivity by picking and choosing which part of the Redfield equation to coarse-grain and
by how much. First, we rewrite the Redfield equation (Eq. 10) in terms of the geometric mean SD, instead
of the arithmetic one. For example, we write (γ
in
+ γ
jm
)/2 =
γ
in
γ
jm
+ (
γ
in
γ
jn
)
2
/2. The Redfield
equation 10 becomes
nm
dt
= i[H, ρ]
nm
+
X
i,j
ρ
ij
A
?
in
A
jm
γ
in
γ
jm
1
2
ρ
ni
A
ij
A
?
mj
γ
ij
γ
mj
1
2
A
?
ji
A
ni
ρ
jm
γ
ji
γ
ni
+
X
ij
ρ
ij
A
?
in
A
jm
f (ω
in
, ω
jm
)
1
2
ρ
ni
A
ij
A
?
mj
f(ω
ij
, ω
mj
)
1
2
A
?
ji
A
ni
ρ
jm
f(ω
ji
, ω
ni
)
, (34)
where we introduce the "detuning" function,
f(ω, ω
0
) =
1
2
q
γ(ω)
q
γ(ω
0
)
2
+ i
S(ω) S(ω
0
)
. (35)
Note that f(ω, ω) = 0, which is why we call it the detuning function. We can write ω = ω
a
+ (ω ω
0
)/2 and
ω
0
= ω
a
(ω ω
0
)/2, where ω
a
= (ω + ω
0
)/2 is the center frequency and ω ω
0
is the detuning. For small
detuning, defined as |ω ω
0
| < |ω
a
|, we apply the Tyler expansion and find
f(ω, ω
0
) = i(ω ω
0
)
S(ω
a
)
ω
a
+ O
h
(ω ω
0
)
2
i
. (36)
For Ohmic bath at zero temperature and exponential cutoff, the images in Fig. 1(a) and (d) display
the geometric mean g(ω, ω
0
) =
p
γ(ω)γ(ω
0
) and the detuning function magnitude versus system frequencies,
respectively. The entities f and g are independent of the system operator A. Without any coarse-graining,
f and g are overall comparable in magnitude, which is unfortunate, because, if the detuning function were
negligibly small compared to the geometric mean, we could neglect the second line in Eq. 34, and end up
with a GKSL master equation.
In spite of the significant value of the detuning function, we can make a compromise following the steps
discussed next. We recast the coarse-grained Eq. 25 as
nm
dt
= i[
˜
H, ρ]
nm
+
X
i,j
ρ
ij
A
?
in
A
jm
˜g(ω
in
, ω
jm
)
1
2
ρ
ni
A
ij
A
?
mj
˜g(ω
ij
, ω
mj
)
1
2
A
?
ji
A
ni
ρ
jm
˜g(ω
ji
, ω
ni
)
+
X
ij
ρ
ij
A
?
in
A
jm
˜
f (ω
in
, ω
jm
)
1
2
ρ
ni
A
ij
A
?
mj
˜
f(ω
ij
, ω
mj
)
1
2
A
?
ji
A
ni
ρ
jm
˜
f(ω
ji
, ω
ni
)
, (37)
where ˜g and
˜
f are the coarse-grained geometric mean and detuning function,
˜g(ω, ω
0
) =
q
γ(ω)γ(ω
0
) sinc
(ω ω
0
)T
0
2
, (38)
˜
f(ω, ω
0
) =
n
1
2
q
γ(ω)
q
γ(ω
0
)
2
+ i
S(ω) S(ω
0
)
o
sinc
(ω ω
0
)T
0
2
. (39)
The characteristic detuning frequency of these functions is governed by the sinc function. Irrespective of
the center frequency, if |ω ω
0
| 1/T
0
, then the sinc function will suppress
˜
f in Eq. 39 by prefactor of 1/T
0
.
Lets us next consider the other frequency range, |ω ω
0
| < 1/T
0
. If in this range |ω
a
| 1/T
0
, then the
condition |ω ω
0
| |ω
a
| will hold. Eq. 36 is then applicable and leads to
˜
f(ω, ω
0
) i(ω ω
0
)
×
S(ω
a
)
ω
a
sinc
(ωω
0
)T
0
2
2i
T
0
S(ω
a
)
ω
a
sinc
(ωω
0
)T
0
2
.
(40)
So the corresponding terms in Eq. 37 are also suppressed by the inverse coarse-grain time. Vice-versa,
if |ω
a
| < 1/T
0
, then Eq. 36 may not hold, and the corresponding terms in Eq. 37 are not suppressed.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 8
Figure 1: a) The effect of coarse-graining. a)-c) Coarse-grained geometric mean
SD, at the respective coarse-graining times corresponding to ω
c
T
0
= 0, 6.8, and
25. Yellow-Hot scale range is the same between the panels. d)-f) Coarse-grained
detuning function, at ω
c
T
0
= 0, 6.8, and 25, respectively. Yellow-Hot scale range
is the same between the panels, but reduced by factor of two relative to that in
panels a)-c). g = 1.
Putting it all together, |
˜
f(ω, ω
0
)|
is not suppressed in the frequency
range |ω ± ω
0
| < 1/T
0
, which
corresponds to a region bordered
by a rotated square in the fre-
quency space [Fig. 1(e)]. As the
coarsegraining time increases, the
square shrinks towards the origin
and the entire function
˜
f(ω, ω
0
) is
suppressed.
Let us visually inspect these
coarse-graining effects at times
T
0
= 6.8
c
and 25
c
. The
coarse-grained geometric mean
SD, displayed in Figs. 1(b,c), is
significantly suppressed in the fre-
quency range |ω ω
0
| > 1/T
0
,
but remains mostly unchanged at
frequencies near the diagonal ω =
ω
0
. The coarse-grained detuning
function is similarly suppressed in
the frequency range |ω ω
0
| >
1/T
0
, as shown in Figs. 1(e,f). In
contrast to the geometric mean,
the coarse-grained detuning func-
tion is strongly suppressed in its
entirety. The white-dashed ro-
tated square in Fig. 1(e) is for il-
lustrative purpose only, display-
ing the region inside which the de-
tuning function is not suppressed
by coarse-graining. The square
shrinks as T
0
increases. An im-
portant point to make here is that
the characteristic coarse-graining
time that suppresses the detun-
ing function relative to the geo-
metric mean is independent of the
strength of the coupling of the sys-
tem to the environment and of the
details of the system operator A.
For the Ohmic bath with ex-
ponential frequency cutoff, we quantify the effect of coarse-graining in terms of the trace-norm of functions
˜g and
˜
f. Fig. 2 displays the ratio of the norms of the coarse-grained detuning function and geometric
mean. While at small T
0
, the norm ratio is constant, above a characteristic time it fits well to inverse time
dependence as shown by the dashed line in Fig. 2. There is a crossover time of order bath-correlation time
1
c
at which the time dependence changes from constant to inverse.
This analysis conveys a picture of dissipative quantum dynamics. The effects of the detuning function
average out relative to the effects of the geometric mean density, on timescale governed by τ
c
, and possibly
a system time-scale encoded in the operator A. Namely, the second line of the master equation 37 has a
crossover coarse-graining time T
C
above which its contribution starts decreasing as inverse coarse-graining
time. Then if we consider system dynamics on timescale longer than T
c
, we can neglect the detuning function
altogether, but retain all other contributions in their primitive (e.g., not coarse-grained) form.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 9
Figure 2: a) Ratio of the tracenorms of two kernels versus coarse-graining time. Black-full curve is the ratio between the
detuning-function and the geometric mean SD. Red-dashed line is the best inverse fit at large T
0
. Diagram: Four master
equations divided into quadrants based on complete positivity (CP) and time dependence of the coefficients. GAME is
initially obtained by approximating the Redfield equation. Its accuracy is assessed by comparison with the TDC-Redfield
equation, eider directly or via an intermediate TDC-GAME equation.
Next we introduce the generator matrix,
L = A
γ, (41)
(where
γ is defined as the matrix with elements
γ
ij
), which leaves us with a GKSL master equation that
we call GAME,
dt
= i[H, ρ]
1
2
{LL
, ρ}+ L
ρL, (42)
that is more general than the Davies-Lindblad equation. Without the renormalized Hamiltonian H, this
equation would be identical to the phenomenological PERLind master equation. [55] It is interesting that
the coarse-graining time is absent in Eq. 42, just as it is absent in the Davies-Lindblad equation.
In this paper we study by how much the solutions of different master equations can diverge from each
other for the same initial condition. In the diagram in Fig. 2, the Redfield and GAME are shown in the
upper row. They both have asymptotic coefficients, but differ in complete positivity. As described above we
restore complete positivity of the Redfield equation directly, by dropping the detuning function, which is the
approximation step indicated by the arrow in the upper row of the diagram.
We also include the TDC-Redfield equation into the analysis, e.g., Eq. 18, because it is a higher level of
approximation relative to the Redfield equation, in terms of accuracy. [22, 47] So we can directly compare
GAME and TDC-Redfield equation, along the diagonal arrow in the diagram in Fig. 2.
Alternatively, we first restore complete positivity of the TDC-Redfield equation by an arithmetic geometric
mean SD approximation. This leads to an intermediate master equation with TDCs that we refer to as TDC-
GAME. Specifically, we introduce a time-dependent renormalization of the system Hamiltonian,
H(t) = H
0
i
2
[AA
f
(t) A
f
(t)A] H
0
+ H
L
(t). (43)
The matrix elements of H(t) are expressed in analogy with the time-independent case,
H
nm
(t) = E
n
δ
nm
+
X
i
H
t
(ω
ni
, ω
mi
)A
ni
A
im
, (44)
where
H
t
(ω, ω
0
) =
1
2
S
t
(ω) + S
t
(ω
0
) + i
γ
t
(ω) γ
t
(ω
0
)
2
. (45)
Now we rewrite Eq. 18 analogous to how we rewrite the Redfield equation in Sec. 3. Then we separate the
geometric mean of the time-dependent SDs, and drop the time-dependent detuning function. The resulting
TDC-GAME equation in the Schrödinger picture is
dt
= i[H(t), ρ]
1
2
{L(t)L(t)
, ρ}+ L(t)
ρL(t), (46)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 10
where in terms of matrix elements γ
t
(ω
nm
) defined in Eq. 17, we have
L(t) = A
γ
t
. (47)
Notice that the matrix elements γ
t
(ω
nm
) can be negative and therefore
γ
t
can be imaginary. Nevertheless
Eq. 46 retains the Lindblad form. GAME is now derived from TDC-GAME, by applying the
R
t
R
approximation, that is, by replacing the TDCs with their asymptotic values. We will study the relation
between these four equations in Sec. 7.4.
5 3-level Jaynes-Cummings model
Quantum correlations mediated by the heat bath can be exploited to generate nonlocal unitary transforma-
tions including entanglement between qubits. [22, 45, 50, 51, 63, 64] GAME can well describe such correlations.
This is primarily due to the unitary component of the system evolution, governed by the renormalized Hamil-
tonian given by Eq. 8. The correlations become particularly strong near the crossings of system frequencies,
which turn to avoided crossings when the frequency difference becomes lower than 1
r
as we study in this
section.
We apply GAME on an exactly solvable Jaynes-Cummings type of model, complex enough to display
these quantum correlations. Specifically, we utilize a V-type 3-level system studied in Ref. [45].
The unperturbed system eigenstates are |0i, |1i, and |2i, with the respective sorted eigenenergies E
0
= 0,
E
1
, and E
2
. The total Hamiltonian of the system and the bath of noninteracting linear harmonic oscillators,
is
H
tot
=
2
X
i=0
E
i
|iihi| +
X
k
ω
k
b
k
b
k
+ H
sb
, (48)
where H
sb
is the system-bath coupling Hamiltonian
H
sb
= C B + C
B
. (49)
Here
C = |1ih0| + |2ih0| =
0 0 0
1 0 0
1 0 0
, (50)
B =
X
k
g
k
b
k
, (51)
b
k
are boson annihilation operators, and ω
k
0. The appropriate bath-correlation function for this system
is C(t) = T r[ρ
b
B(t)B
(0)] =
P
k
|g
k
|
2
e
k
t
. [45]
The Redfield master equation is
dt
= i[H
0
, ρ] CC
f
ρ ρC
f
C
+ C
f
ρC + C
ρC
f
, (52)
where C
f
= C (γ/2iS). Note that the Hermitian conjugate of C is now explicit, in contrast to that in Eq. 4
where it would be redundant. Similarly, the renormalized Hamiltonian features the Hermitian conjugate of
C
H = H
0
i
2
(CC
f
C
f
C
), (53)
or explicitly
H =
0 0 0
0 E
1
+ S(E
1
) S
i
4
γ
0 S +
i
4
γ E
2
+ S(E
2
)
, (54)
where S =
1
2
[S(E
1
) + S(E
2
)] and γ = γ(E
2
) γ(E
1
). Lastly, the generator in equation 42 has matrix
elements L
nm
= C
nm
γ
nm
.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 11
Figure 3: Time dependence of the populations, for initial condition
ρ(0) = |1ih1|. Blue and red lines are from the exact solution. Black
lines are from the approximate solution of Eq. 42. a) and b): Case A:
E
1
= 0.095, and E
2
= 0.105. c) and d): Case B: E
1
= 0.09975, and
E
2
= 0.10025. g = 0.001.
For the initial condition, we suppose
that the system is in state |1i. Repeating
from Ref. [45], the time dependent state of
the total system in the interaction picture
is
|ψ(t)i = [c
1
(t)|1i + c
2
(t)|2i] |0i
b
+ |0i
X
k
d
k
(t)b
k
|0i
b
, (55)
where |0i
b
is the bath vacuum. The ex-
act solution follows by numerically solving
a system of two Volterra integro-differential
equations,
˙c
1
= f
1
? c
1
e
12
t
f
2
? c
2
, (56)
˙c
2
= f
2
? c
2
e
12
t
f
1
? c
1
, (57)
where f
1
(t) = exp(
10
t)C(t), f
2
(t) =
exp(
20
t)C(t), f ? g =
R
t
0
f(τ)g(t τ).
For the system parameters we use g =
0.001, and E
1
= 0.095ω
c
, and E
2
= 0.105ω
c
in case A, and E
1
= 0.09975ω
c
and E
2
=
0.10025ω
c
in case B. These are chosen to be
the same parameters as in Ref. [45], so that
we can directly compare the accuracy of the
master equations. Figure 3 displays numerically obtained exact populations of levels 1 and 2 versus time, in
parameter cases A and B, respectively. The blue and red curves reproduce previous work. [45]
5.1 Accuracy of the approximate solution
Next we investigate the error of the approximate state. We solve Eq. 42 as explained in appendix 9. The
approximate populations of levels 1 and 2 versus time are displayed by the black lines in Fig. 3. The difference
between the approximate and the exact solution is barely discernible. Not shown is the coherence ρ
12
(t),
which also agrees well with the exact solution.
The trace-distance between the approximate and exact solutions is displayed by the black lines in Fig. 4.
Also shown by the blue line in Fig. 4(a), is the error of the optimal CGSE, obtained by digitizing the blue
line of data in Fig. 3 of Ref. [45]. GAME is almost an order of magnitude more accurate than the CGSE,
and almost completely captures the population oscillations as seen by the diffuse error with respect to time.
We also determine the trace-distance between the solutions of the Redfield and GAME equations, and
find that the distance is ten times lower than the black lines in Fig. 4. Additionally, we solve Eq. 42 without
renormalizing the system Hamiltonian (Eq. 53), which is equivalent to PERLind approach. [55] In that case
the state error is much larger than the black lines in Fig. 4. The error increases by a factor of several hundred
and demonstrates the necessity to renormalize the Hamiltonian to properly account for quantum correlation.
The dissipator can induce correlations as well, [55] but it is insufficient for high accuracy. We will return to
the subject of accuracy of PERLind in Sec. 7.3.
Even though the interaction Hamiltonian (Eq. 49) is of a RWA form, the Redfield equation in the inter-
action picture has oscillations and therefore it is not in the RWA. Although the large oscillation frequencies
ω
10
and ω
20
are absent, the oscillations at detuning frequency ω
21
remain.
The accuracy of the RWA in this system was found to be poor for the system parameters we apply. [45]
The error is exacerbated in the regime where the relaxation rate is larger than detuning ω
21
. This can be
traced to a known property of the RWA [22, 63, 65] that there is an artificial discontinuity of the equation
at zero detuning (level degeneracy). As an example, Eqs. 29 and 32 do not merge at the crossings of levels n
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 12
Figure 4: Trace-distance between the exact quantum state and the approximate states obtained by solving the master
equations. Black lines: Error of the solutions of the Lindblad equation (Eq. 42). Blue line: Error of the optimal CGSE
obtained by scanning and digitizing data in Ref. [45]. g = 0.001. a) E
1
= 0.095ω
c
, E
2
= 0.105ω
c
. b) E
1
= 0.09975ω
c
,
E
2
= 0.10025ω
c
.
and m. In the three level system in particular, at finite detuning, the RWA of the Hamiltonian 54 will be its
diagonal part. At zero detuning, however, the RWA groups the frequency duplicates under same generator,
leading to the renormalized Hamiltonian 54.
The artifact has been remedied in case of two qubit systems coupled to a heat bath, by applying a coarse-
grained approach, [63, 65] which leads to a continuous description of entanglement creation as function of the
frequencies of the qubits. The entanglement between two-level atoms can also be created in the dissipative
process of spontaneous emission. [66] By similar mechanism, GAME replaces the discontinuity of the RWA
with a continuous avoided crossing, as will be shown next, thereby creating superpositions between levels 1
and 2 at finite detuning.
5.2 Results and discussion
Interestingly, Figs. 3 c) and d) show a few initial oscillations that decay rapidly followed by a much slower
relaxation process without any oscillations. In case A, as can be seen in Figs. 3 a) and b), the oscillations
and populations decay on a similar scale. At first sight, this property is puzzling since the Fermi-golden rule
rates are approximately the same in cases A and B.
Thinking further, however, the disparity between the relaxation rates is not utterly surprising. Namely,
if levels E
1
and E
2
were degenerate, one could find a basis such that one state has a zero matrix element
in the coupling Hamiltonian. For example, if in one basis the coupling Hamiltonian is represented by the
matrix in Eq. 50, we can change the basis to
h
|0i, (|1i + |2i)/
2, (|1i |2i)/
2
i
. (|1i|2i)/
2 is the dark
state. [45] If the levels are not exactly, but nearly degenerate, the bath will renormalize the states so that
one state approaches the dark state resulting in the suppressed relaxation rate.
Here we present a simple procedure that determines the system oscillation frequencies and relaxation
rates based entirely on the properties of the renormalized Hamiltonian given by Eq. 54. The renormalized
energies are defined as the nonzero eigenvalues of the matrix in Eq. 54,
E
0
1,2
= E + S
s
S
2
+
E + S
2
2
+
γ
4
2
, (58)
where E = (E
1
+ E
2
)/2, S = S(E
2
) S(E
1
), and the other terms are identical to those in Eq. 54.
Let us introduce a parameter dependence of the system Hamiltonian,
H
0
= E
1
|1ih1| + [E
2
λ(E
2
E
1
)] |2ih2|. (59)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 13
Figure 5: Renormalized eigenenergies and relaxation rates in units of ω
c
, versus parameter λ of the unperturbed Hamiltonian.
The higher level (red) has suppressed effective width and becomes a dark state with zero relaxation rate at λ = 1. The
energy gap imposed by the heat bath is 2S(E
1
), where S is the PD. At λ = 0, E
1
= 0.1ω
c
, E
2
= 0.105ω
c
. g = 0.001.
Figure 6: Populations of the renormalized states versus time, in case B. E
1
= 0.09975ω
c
, E
2
= 0.10025ω
c
, and g = 0.001.
The initial conditions are ρ(0) = |1
0
ih1
0
| and ρ(0) = |2
0
ih2
0
| in a) and b), respectively. Dashed black line is the best fit to
exp(Γ
0
t). Γ
0
is the effective (renormalized) width.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 14
level 1 level 2
bare widths Γ
1
= γ(E
1
) Γ
2
= γ(E
2
)
renormalized energies E
0
1
= E 2S κ
(∆E)
2
S
E
0
2
= E + κ
(∆E)
2
S
renormalized widths Γ
0
1
= Γ
1
+ Γ
2
Γ
0
2
Γ
0
2
=
Γ
1
2
4
E
2S(E)
2
Table 2: Bare and renormalized parameters in the strong coupling case.
As a function of λ, the energy levels of H
0
cross at λ = 1, while the renormalized eigenenergies (Eq. 58) have
avoided crossing as shown in Fig. 5(a) with the energy gap 2S(E
1
). The avoided crossing regime is found
in the region of energy where the widths γ(E
1,2
) are comparable to the spacing E
2
E
1
. In the regime of
strong anticrossing, e.g., |E
2
E
1
| γ(E
1,2
), Eq. 58 reduces to
E
0
1
=
E 2S κ
(∆E)
2
S
, E
0
2
= E + κ
(∆E)
2
S
, (60)
where
κ =
"
1
16
γ
E
2
+
1
4
1 +
S
E
2
#
E=0
.
At the center of the crossing, the renormalized states are
|1
0
i =
|1i + |2i
2
, (61)
|2
0
i =
|1i |2i
2
. (62)
The last is the dark state.
Here, the avoided crossing is caused by the coupling between the system and the environment of linear
harmonic oscillators. We find that not only the states, but also the relaxation rates are strongly affected
near such avoided crossings. To obtain the relaxation rates, we consider the initial states of the system
ρ(0) = |1
0
ih1
0
| and |2
0
ih2
0
| and determine the populations of the renormalized states versus time again by
solving Eq. 42. The results are shown in Fig. 6 for case B. In contrast to the time dependence of the
unperturbed states shown in Fig. 3, there are no oscillations in Fig. 6. The goodness of fit to e
Γ
0
t
in Fig. 6
indicates that the populations of the renormalized states decay with single effective or renormalized widths.
The width Γ
0
versus λ is displayed in Fig. 5(b). The higher energy state has a strongly suppressed width
and becomes dark at λ = 1, while the lower one has enhanced width and increases in brightness by factor of
two at λ = 1.
Table 2 displays the effective level widths of the reduced 3-level system, in the regime near the avoided
crossings. The accuracy of the prefactor in the relaxation rate of the dark state is approximately 2-3%. Note
the "sum rule" where the sum of the widths Γ
0
2
+ Γ
0
3
of the renormalized states is independent of the coupling
strength.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 15
6 Heisenberg Ferromagnetic Spin-Chain
As a testbed for various master equations, here we consider a spin-chain of up to n = 25 spin-1/2 particles.
Each particle couples to its nearest neighbor via the ferromagnetic exchange interaction and to each other
through the r
3
magnetic dipole-dipole interaction. This system has a rapidly increasing density of states with
energy and is rich in many-body phenomena, including, the energy gap (ferromagnetic resonance), collective
excitations (standing spin-waves), and topological excitations (domain walls). There is a corresponding
plethora of relaxation processes, including Gilbert damping, emission of spin-waves, and domain wall motion,
spanning wide range of relaxation times, that can be exponential in system size. The focus here will be on
relaxation of the density matrix in response to a macroscopic displacement of the magnetization, which is
initially set perpendicular to the chain (easy) axis.
6.1 Magnetic Hamiltonian
The system is modeled by the Heisenberg model with dipole-dipole coupling,
H
0
= J
n1
X
i=1
~
S
i
~
S
i+1
d
n1
X
i=1
n
X
j=i+1
3S
z,i
S
z,j
~
S
i
~
S
j
(j i)
3
. (63)
The Hamiltonian commutes with S
z
=
P
i
S
iz
, but does not commute with
~
S
2
, where
~
S =
P
i
~
S
i
. Here
n = 25 and the eigenvalues display Kramers degeneracy. In addition, the Hamiltonian is invariant with
inversion of spins about i = (n + 1)/2. We use the Kramers, S
z
, and the inversion symmetry to reduce the
Hamiltonian, and find the eigenvalues and eigenstates using the Lanzos algorithm on sparse matrices. Most
simulations were done on a 16-core Intel Xeon processor and 192GB RAM, while the heaviest simulations
in Sec. 7 were done on a 32-core AMD EPYC and 384GB RAM. The eigenstates have quantum numbers
S
z
= S, S + 1, ..., S and are even or odd functions of spin. The states with S
z
< 0 are obtained by
time-reversal operation on the states with S
z
> 0.
For parameters we choose J = 400 and
d
= 6. With no anisotropy, (e.g.,
d
= 0), the spacing between
the lowest lying spin-multiplets is δ = 2.5. Each multiplet represents a different standing spin wave, while
the spin-wave number varies within the corresponding multiplet. In the case
d
= 6, the ground state is
that of a uniaxial ferromagnet, |S, S
z
= ±Si, (chain axis is along z-direction), the ground state energy is
E
g
= 2485.3, and there is an anisotropy gap ∆ = 20.1 in the excitation spectrum, as can be seen in Fig. 7.
The gap is the ferromagnetic resonance. We work in the regime where > δ, which is the regime common
in devices.
Fig. 7 displays the lowest 1728 energy levels versus S
z
. There are approximately 3 ×10
6
Bohr frequencies
spanning the range 158.5 < E
ij
< 158.5. In the macrospin approximation, also known as the Stoner-
Wohlfarth model, [67] H
0,sf
= KS
2
z
/S. The uniaxial anisotropy K can be estimated from the energy
gap as K = /2 = 10. The blue parabola represents the energy of the magnet vs. S
z
in the macrospin
approximation. In that approximation, the Néel energy barrier for this uniaxial 1D magnet is S/2 = 125.
We are interested in the low energy collective magnetic dynamics. It suffices to find the eigenvalues and
eigenstates, spanning the energy range between the ground state to not far over the Néel barrier, as we will
see shortly. In particular, let us consider the initial state with the magnetization fully oriented along the
x-direction. To obtain this initial state, we truncate S
x
to the vector space on the lowest 1728 eigenstates and
find the eigenstate of the truncated S
x
with the maximum eigenvalue. In that state, we obtain S
x
12.41,
which is lower than 12.5 because of the truncation error. Since this value is still very close to 12.5, it shows
that the vector space is sufficiently large to describe the uniformly magnetized state, and the states accessible
by relaxation from that state. The expectation value of the energy in the perpendicularly magnetized state
is approximately the same as the energy at the top of the blue parabola in Fig. 7(a), while the energy
uncertainty is rms(H
0
˙) = 11.6.
Note that we have applied a very weak magnetic field along the z-direction to break the S
z
degeneracy.
The difference in energy between |S, Si and |S, Si is 0.000125∆.
The initial uniformly magnetized state is decomposed into the energy eigenstates as shown in Fig. 7(b). It
shows a wide spread of the perpendicularly magnetized state over many eigenstates, with the spacing between
different components much lower than the gap. Often these spacings will be smaller than the coupling to
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 16
Figure 7: a) Low lying energy levels in a 25 site spin-1/2 chain, in the Heisenberg model with short range ferromagnetic
exchange and 1/r
3
dipole-dipole interactions. J = 400 and
d
= 6, respectively. Blue parabola: Energy barrier between
the ground states in the macrospin approximation. b) Decomposition of the uniformly magnetized state (|Ψ
x
i) with
magnetization oriented perpendicular to the easy axis of the spin-chain, in the energy eigenbasis. The weight on the y-axis
is defined as log |hE|Ψ
x
i|
2
. The black double-arrow line indicates the energy range centered at the top of the Néel barrier,
at system-bath coupling g
tot
= 1. n = 25 and N = 1728.
the heat bath, which can introduce avoided crossings as we discussed in the example of the 3-level system.
These properties make our system well suited for testing completely positive master equations in many-body
regime.
The energy levels under the parabola in Fig. 7(a) are noteworthy. They are grouped into branches, which
have weak dependence of energy with S
z
. On these branches, energies have two-fold near degeneracy [not
visible in Fig. 7(a)] between even and odd spin states. At small |S
z
|, these states are similar to domain walls
located at distance S
z
from the middle of the chain. The energy of the domain wall doesn’t change much if
it moves near the middle of the chain, explaining the flatness of the branches. At this time, it would be a
distraction to delve into these states but they may be interesting for future reasearch.
6.2 Coupling to the heat bath
Each spin is coupled to a noisy magnetic field, analogous to the environment in the celebrated Brown
model. [68] In our case, there are 75 independent baths of linear harmonic oscillators, 3 for each spin, one
for each direction of the field. The coupling Hamiltonian is
H
sb
=
n
X
i=1
X
α=x,y,z
S
B
. (64)
Here B
i,α
=
P
k
g
k
(b
k,iα
+ b
k,iα
) are Hermitian bath operators, and we assume slow bath at zero temperature
with Ohmic SD and cutoff frequency ω
c
= 6∆. All baths have the same spectral functions.
The Redfield equation for the spin-chain is
dt
= i[H
0
, ρ] +
X
i,α
S
S
iα,f
ρ ρS
iα,f
S
+ S
iα,f
ρS
+ S
ρS
iα,f
, (65)
where S
iα,f
= S
Γ, and Γ is given by Eq. 6.
The corresponding equation in GAME is
dt
= i[H
0
, ρ] +
X
i,α
L
ρL
1
2
{L
L
, ρ}
, (66)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 17
Figure 8: Almost unitary magnetization dynamics of the spin-chain in the regime of weak coupling to the bath. a) Rapid
decay of the initially perpendicular magnetization due to the effective phase randomization of the many-body states of the
spin-chain. b) Energy of the spin-chain versus time. n = 25, N = 1728, g
tot
= 0.01, T
fm
= 2π/, ∆ = 20.1.
with the renormalized Hamiltonian,
H
0
= H
0
i
2
X
(S
S
iα,f
S
iα,f
S
). (67)
The Lindblad generators are obtained from Eq. 41, e.g. L
= S
γ. We solve these equations by numeric
iterations in the truncated Hilbert space, as explained in appendix 9.
6.3 Nearly Unitary Magnetization Dynamics
Let us first familiarize with the system by considering what we call very weak coupling to the heat bath, e.g.,
we assume g = 0.000133 for the SD in each bath. Since there are 75 independent baths, the total coupling
is g
tot
= 75g = 0.01.
Fig. 8(a) displays a rapid decay of the perpendicular magnetization. The bottom axis is time in units
of T
fm
, where T
fm
= 2π/ is the period of the ferromagnetic resonance. We consider the dynamics as
mostly unitary, in a sense that T
2
changes from one at t = 0 to 0.744 after thirty FMR-periods. The
decomposition of the initially uniformly magnetized state into eigenstates, displayed in Fig. 7(b), results in
effective loss of magnetic coherence over the period of the FMR. Over that time, T r(ρ
2
) drops by only 0.01,
so the system retains quantum coherence even when the magnetization drops close to zero.
The energy of the system hH
0
i is shown in Fig. 8(b). It decays slowly due to the dissipative coupling to
the heat bath. The small drop of hH
0
i, relative to the initial energy, indicates that the figure displays an
early stage of the relaxation cascade.
6.4 Dissipation Time Scales
Now we examine the time-scale on which the environment changes the system properties in the interaction
picture. First we evaluate the density matrix relaxation rate 1
r
for two initial conditions: a uniform
perpendicular magnetization and the first excited state.
To determine τ
r
, we rotate into the interaction picture % = exp(iH
0
t)ρ exp(iH
0
t), and determine τ
r
as
this,
||%
RED
||
1
τ
r
=
d%
RED
dt
1
, (68)
where the subscript RED indicates the state is obtained by solving the Redfield equation. From now, we
approximate this as
1
τ
r
=
d%
RED
dt
1
(69)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 18
Figure 9: a) State relaxation rate versus time, for g
tot
= 0.005 (blue-solid) and 0.01 (red-dotted). b) Energy relaxation of
the FMR at g
tot
= 1. Dashed thick-black line is the best fit to exp(Γ
0
t). Γ
0
= 0.220T
1
fm
= 0.035∆ is the FMR width.
c) Trace-distance between the solutions of the Redfield equation and the approximate states in the RWA (red-dotted) and
GAME (blue-solid), g
tot
= 0.01. In all panels, n = 25, N = 1728, ω
c
= 6∆, T
fm
= 2π/, ∆ = 20.1.
since ||%
RED
||
1
is only slightly larger than one. (The trace norm of the state can be larger than one if there
are negative eigenvalues, but the effect is usually small in the range of applicability of the Redfield equation.)
Fig. 9(a) displays 1
r
versus time, for the initially perpendicular magnetization, for g
tot
= 0.005 and
0.01. The state relaxation rate 1
r
scales with g
tot
as expected. We find the rate 1
r
by time-averaging,
which calibrates the rate according to 1
r
= 0.293g
tot
. We similarly find that for the system initially
prepared in the first excited state, the state relaxation rate is 1
E
= 0.08g
tot
.
In many cases in micromagnetics, the parameter of interest is the energy relaxation rate, rather than the
density-matrix relaxation rate. To give a reference point to where our spin-chain system lies in that context,
we find the energy relaxation rate of the first excited state of the system, e.g., the width of the ferromagnetic
resonance. We calculate the system energy versus time, assuming the initial condition to be the first excited
state. We fit energy versus time to exp(Γt), where Γ = 1
E
is the width we seek. The one parameter fit is
shown in Fig. 9(b), leading to the FMR width 1
E
= 0.0348g
tot
. The FMR is well resolved at our highest
value of g
tot
= 1. In fact, the FMR width in that case corresponds to the Gilbert-damping parameter of
1/(4τ
E
∆) = 0.0087, which is not uncommon in transition metal ferromagnets. [69]
In Fig. 9(c), we plot the trace-distance between the Redfield state and those of GAME and RWA, at
g
tot
= 0.01. GAME clearly produces a state much closer to the Redfield state relative to that of the
RWA. In this weak coupling regime 1
r
= 0.00293∆, and the Born-Markov condition is well satisfied,
1
c
τ
r
5 × 10
4
, assuring high accuracy of the Redfield equation. Still, the RWA has poor accuracy since
the smallest level spacing is much smaller than 1
r
.
We note in passing on the ability of trace-distance to distinguish quantum states. Fig. 8 would naively
suggest that the RWA accounts for the dynamical properties of the system spectacularly well, yet Fig. 9(c)
tells a very different story, that the trace-distance between the RWA and the Redfield state in the same
situation increases with time and is far from saturation.
To study the scaling of the trace-distance with system-bath coupling, in Fig. 10(a) we display the trace-
distance
1
2
||%
RED
(t) %
GAME
(t)||
1
versus time and the coupling strength g
tot
. As the coupling decreases,
we see that the trace distance versus time saturates on a time scale independent of the coupling g
tot
. This
introduces a new time scale of the system, the saturation time T
s
. It is indicated by the arrow in Fig. 10(a),
showing that it is comparable to T
fm
. In appendix 10 we show that increasing the bath cutoff frequency
above 6∆ has a weak effect on T
s
.
In Fig. 10(a) and at low g
tot
, the trace-distance versus time saturates because of the transient nature of
the detuning term that we drop from the Redfield equation. That is, the integral of %
RED
(t) %
GAME
(t)
will saturate on time-scale over which the detuning superoperator averages out to zero. Thus restoration of
complete positivity can be viewed as a state slip on time scale T
c
, with no further trace-distance increase
after that time. A different state-slip at the correlation time scale of the bath is a known phenomenon in
stochastic quantum dynamics, [7074] but here the slippage to regain complete positivity takes a longer,
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 19
Figure 10: a) and b) display the trace-distance between the solutions of the Redfield and GAME equations versus time, in
the units of the period of ferromagnetic resonance. g
tot
is the total dimensionless coupling constant between the spin-chain
and the heat bath. b) Maximum trace-distance between the solutions of the two equations, versus rate 1
r
in units of the
cutoff frequency ω
c
. The blue line is the linear scaling discussed in the text. n = 25, N = 1728, ω
c
= 6∆, T
fm
= 2π/,
∆ = 20.1.
system dependent time. The trace-distance accumulated over time T
c
is T
c
r
= gT
c
c
. It is linear with g,
since neither T
c
nor τ
c
depend on g.
In Fig. 10(b) we show the maximum trace-distance between the solutions of the Redfield and GAME
equations versus state relaxation rate 1
r
. The figure displays asymptotic linear scaling between the trace-
distance and the state relaxation rate at low 1
r
. The best linear fit on the lowest four points is
1
2
||%
RED
%
GAME
||
1,max
6.41
1
ω
c
τ
r
=
1.07
τ
r
, (70)
and is shown by the blue line.
7 Comparisons Between Various Master Equations
.
There have been several recently derived CP master equations with the claim that they perform well in
the limit of small level spacings of the quantum system. In addition, there are master equations with TDCs,
that we discussed in Secs. 3 and 4 that have fewer problems than the Redfield equation with asymptotic
coefficients. In this section we evaluate how much the solutions of different master equations can diverge
from each other for the same initial condition.
For the initial condition we use the same uniformly polarized state perpendicular to chain axis. In Sec. 7.4
we also study the decay of the FMR. For the coupling to the Ohmic bath, we use g = 0.0133 per bath (g
tot
= 1
for the 25-site chain), the highest value in Fig. 10, where the break-down of the master equation is hinting.
Also we study the case g = 0.000133 per bath (g
tot
= 0.01 for the 25-site chain), which we refer to in Sec. 6.3
as "nearly unitary" dynamics. While the higher value of g
tot
is more realistic in terms of damping in magnetic
materials, the Redfield approach is perhaps on shaky grounds. For the smaller g
tot
, the Redfield equation is
valid.
In all examples we use the same exponential cutoff with frequency of ω
c
= 6∆ = 120. In this case
τ
c
r
= 0.0493 for the larger g
tot
at n = 25, so we are comfortably but perhaps not too comfortably within
the limit of the validity of the Bohr-Markov approximation. Since different master equation have different
complexities, we adjust the code for fastest performance, depending on the approximation we use. The details
will be provided as we go through the approximations.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 20
7.1 Partial Rotating Wave Approximation.
The partial RWA (PRWA) extends the range of applicability of the RWA. [19, 22, 48, 49] It begins by
observing that if the SD is flat, then the Redfield equation will be in Lindblad form. The approximation
works as follows. First, the system Bohr energies or frequencies ω
nm
are sorted and divided into bins.
In the next step, the Redfield equation (Eq. 10) is truncated, but less draconically so than in the RWA.
In the first line of Eq. 10, for example, all terms except those where ω
jm
and ω
in
belong to the same bin
are discarded. At the same time, the SD within the bin is approximated by an average (and therefore flat)
SD, one for each bin. As a result, each bin contributes to a Lindblad generator specific to that bin, and
different bins are decoupled. In our case we find the average frequency within each bin, (we call it the coarse-
grained frequency), and calculate the SD at that frequency. We checked that two other frequency averaging
techniques give approximately the same results. This binning is also performed on the loss terms in Eq. 10.
The bin width (δb) is an adjustable parameter, but the least compromising value seems to be set by the
relaxation rate τ
1
r
. If δb 1
r
, then the PRWA will have the same problem as the RWA: the oscillations
at the coarse-grained frequency differences will not average out on time scale τ
r
, resulting in poor accuracy.
If on the other hand δb 1
r
, the assumption of the flat SD becomes less valid.
The number of Lindblad generators increases with decreasing bin width. The RWA is the limit of the
PRWA when there is one frequency per bin. Depending on the bin width, it may take significant resources to
numerically perform the binning and obtain the generators. For example, it takes approximately four hours
to count the frequency duplicates and find the generators in the RWA, for the 25-site Heisenberg chain with
N = 1728. Fortunately, the sparseness of the generators also increases with the decreasing bin width, so the
issue of time and memory is managed using sparse matrices.
Note that for the n = 25 spin-chain, the relaxation rate of the state in the interaction picture 1
r
= 5.91
at g
tot
= 1, while the frequency span is 316. So we expect that the optimum number of bins of 316/5.91 53
for g
tot
= 1. In that case the average number of system frequencies per bin is approximately 6 × 10
4
, and it
fluctuates widely between the bins. In the case of g
tot
= 0.01, the optimum bin width should be much lower.
But we do not solve the PRWA equation over time scale τ
r
, because it would take unreasonable simulation
time.
The renormalization of the system Hamiltonian given by Eq. 11 in a sense completes the PRWA, as fol-
lows. Eq. 10 shows that flat spectral functions result in the cancellation of the PDs on the RHSs. Had we not
renormalized the Hamiltonian according to Eq. 11, the principal density would not cancel in the loss term.
The renormalization therefore completes the PRWA by accounting the PDs via the renormalized Hamilto-
nian, thereby insuring correctness. Previous approaches involved either discarding the Lamb shift, [48] or
coarsegraining the RWA Lamb-shift to the representative bin frequency. [22] The renormalized Hamiltonian
in Eq. 8 does not need any binning so we apply it here in that primitive form.
In Fig. 11 we vary the number of bins between 11 and 201 and find the corresponding solutions in the
PRWA. In Fig. 11(a) numerical calculations show that the overall magnetic dynamics is well accounted for
by the PRWA, especially when compared to the RWA.
The traces distance between the solutions of the RWA, PRWA, GAME and the Redfield equation is shown
in Fig. 11(b,c). Initially all trace-distances are zero because the states are the same. As time increases, the
states under different master equations evolve differently and the trace distances increase. At long times,
however, as the states relax towards the ground state the trace-distances become suppressed. Thus, a trace-
distance versus time exhibits a maximum.
The optimum bin-width corresponds to the trace-distance with the smallest maximum. For g
tot
= 1,
the optimum bin width is in the expected range. At weak coupling, (g
tot
= 0.01), the trace-distances to
the PRWA solutions versus time do not saturate, while the trace-distance between GAME and the Redfield
equation saturates at time similar to the FMR-period.
In summary, we find that GAME is a significantly more accurate approximation of the Redfield equation
than the PRWA. Since both GAME and PRWA are derived by approximating the same Redfield equation,
the former is most likely also a significantly better approximation of true reduced quantum dynamics.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 21
Figure 11: Partial Rotating Wave Approximation. a) Magnetization of the 25-site spin chain versus time, in units of the
ferromagnetic resonance period. The Redfield and GAME solutions are virtually the same. g
tot
= 1. b) Trace norm distance
to the solution of the Redfield equation versus time, from the states solved using the RWA, GAME, and PRWA. g
tot
= 1.
c) The same as b), but with g
tot
= 0.01. The trace-distance in c) is far from saturation in all approximations except GAME.
N = 1728, ∆ = 20, ω
c
= 6∆. Ohmic bath with exponential frequency cutoff.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 22
7.2 Dynamical Coarse-Graining Approximation
As an example of completely positive coarse-grained master equations, here we study the properties of the
dynamically coarse-grained (DCG) master equation. [50, 51] The equation is derived by coarse-graining of the
Liouville equation. The derivation applies the Born, but not Markov approximation. In the vectorized from
of the density matrix, the DCG equation utilizes the Liouville superoperator, which is provided explicitly by
Schaller and Brandes. [50] After translating to the notation we use in this article, we find that the DCG-
master equation reads the same as the Redfield equation in the form of Eq. 14, except that the dissipative
kernel G(ω, ω
0
) in Eq. 15 is replaced with
G
dc
τ
(ω, ω
0
) =
τ
2π
e
i
(ωω
0
)τ
2
Z
−∞
γ(Ω) sinc
(Ω ω)τ
2
sinc
(Ω ω
0
)τ
2
d, (71)
where τ is the coarse-graining time. Note that our double-indices in Eq. 15 are different from Ref. [50], i.e.,
G
dc,ab
(τ) = γ
cd,ba
(τ), where γ
cd,ba
(τ) is the dissipative kernel given by Eq. 23 in that reference.
As alluded to by Schaller and Brandes, [50] this expression can be decomposed in terms of single-variable
functions only,
G
dc
τ
(ω, ω
0
) = e
i
(ωω
0
)τ
2
(
1
2
[γ
τ
(ω
0
) + γ
τ
(ω)]sinc
(ω
0
ω)τ
2
S
τ
(ω
0
) S
τ
(ω)
(ω
0
ω)
τ
2
cos
(ω
0
ω)τ
2
)
(72)
which is critical to our numerical calculation. Namely, to determine the N
4
values of G
dc
τ
(ω, ω
0
), it suffices
to calculate the functions γ
τ
(ω), S
τ
(ω), and S
τ
(ω)/∂ω, thereby reducing the computational complexity
from O(N
4
) to O(N
2
). We show in appendix 12 that γ
τ
(ω) are S
τ
(ω) are the time-dependent spectral and
principal densities, respectively, as defined by Eq. 17. The kernel G
dc
τ
(ω, ω
0
) is a positive semidefinite N
2
×N
2
matrix, which assures that the DCG equation is in a Lindblad form. [50]
In addition, the renormalized Hamiltonian H in Eq. 11 is replaced with
H
dc
τ,nm
= E
n
δ
nm
+
X
i
H
dc
τ
(ω
ni
, ω
mi
)A
ni
A
im
, (73)
with the unitary kernel
H
dc
τ
(ω, ω
0
) =
τ
2π
e
i(ωω
0
)τ/2
Z
−∞
S(Ω) sinc
(Ω ω)τ
2
sinc
(Ω ω
0
)τ
2
d. (74)
We decompose the unitary kernel similar to the dissipative one, as described in more detail in the appendix,
which yields the explicit formula for the kernel that governs the unitary contribution of the heat bath to
system quantum dynamics:
H
dc
τ
(ω, ω
0
) = e
i
(ωω
0
)τ
2
(
1
2
[S
τ
(ω
0
) + S
τ
(ω)]sinc
(ω
0
ω)τ
2
+
1
4
γ
τ
(ω
0
) γ
τ
(ω)
(ω
0
ω)
τ
2
cos
(ω
0
ω)τ
2
)
. (75)
Here we utilize a 23-site spin-chain and solve the DCG-master equation at various coarse-graining times,
assuming a perpendicularly magnetized state as the initial condition. In Fig 12, we display the trace-
distances from the solution of the Redfield equation, for Ohmic SD with exponential frequency cutoff. The
trace-distances versus time exhibit maxima as discussed in the previous section. The curve with the smallest
maximum defines the optimum coarse graining time. For the weak environmental coupling [Fig 12(b)], it
takes unreasonable amount of time to calculate how the curves merge at t and the maximum may
not be reached within the time range of the simulation. In that case we apply the maximum trace-distance
within the range of the figure. As τ decreases from to τ
opt
, the maximum trace-distance of the DCG-state
decreases. However, below τ
opt
the maximum trace-distance increases with decreasing coarse-graining time.
[At τ = 0, the trace-distance varies as 1/2 ||%
RED
(t) ρ
0
||
1
]. At g
tot
= 0.92 and 0.0092, τ
opt
= 1.3T
fm
and
16T
fm
, respectively, in good agreement with square root dependence on g, which makes it consistent with
the coarse-graining approximation. [45, 47] The corresponding values of
τ
c
τ
r
(the optimum coarse-graining
time in the CGSE, see Sec. 7.4) are approximately factor of 10 smaller than the values we find.
It is striking how in Fig 12(b) the trace-distance of GAME saturates at a time scale τ
s
much smaller
than in any other approximation. This results in improved accuracy of GAME, assuming that the figure
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 23
Figure 12: Trace-distances between the solutions of various master equations and that of the Redfield equation, versus time.
a) Regime of moderate environmental coupling g
tot
= 69g = 0.92. Thick lines: RWA (dashed-blue), DCG
opt
(dashed-black),
and GAME (red-full). The approximate optimum coarsegraining time is τ
opt
= 1.3. Thin-full lines: coarse-graining times
τ = 64, 32, 16, 6.4, 3.2, and 1.6, follow the maxima between the RWA and DCG
opt
from top to bottom, respectively. Thin-
dashed lines: τ = 0.96, 0.64, 0.32, and 0.16, follow the maxima from DCG
opt
, bottom to top, respectively. b) Analogous to
a) in the weak environmental coupling regime g
tot
= 0.0092. τ
opt
= 16. Thin-full lines, top-to-bottom: τ = 160, 64, 32,
and 22.4. Thin-dashed lines, bottom-to-top: τ = 11.2, 6.4, 3.2, 1.6, and 0.96. All time scales are in units of T
fm
. Ohmic
SD with exponential cutoff, n = 23, N = 694, ∆ = 20.1, and ω
c
= 6∆.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 24
of merit is the trace distance to the Redfield state. The rapid saturation of the trace-distance (also seen
in Fig. 11) is a specific trait of GAME. It is due to the transient nature of the term that we dropped
from the Redfield equation in order to restore CP. The saturation is less pronounced for moderately weak
environmental coupling [Fig 12(a)], due to the more rapid relaxation which prevents us from reaching the
flat part of the trace-distance curve.
We also perform numerical calculations assuming Ohmic SD with the Drude-Lorentz frequency cutoff. In
the Drude-Lorentz case we obtain the coefficients of the DCG master equation using a different numerical
method, yet we find that results are essentially the same, e.g., see Fig. 17 in appendix 12.
To summarize, we again find that GAME is a significantly more accurate approximation of the Redfield
equation. However, in contrast to the previous section, now the compared approximations operate under
different premises, because they are differently derived from first principles. Here the improved approximation
of the Redfield state does not necessarily imply that GAME is an improved approximation of the true reduced
quantum dynamics. Rather, it opens a question of why there is such a difference between the Redfield and
DCG approximations in the first place. Since the focus of this paper is to find a simple and accurate CP-
approximation of the Redfield equation, we do not speculate on possible answers to this question at this time.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 25
Figure 13: Comparison of four approximations: PERLind, PERLind with the added RWA Lamb-shift, ULE, and GAME.
Trace norm distances are taken with respect to the Redfield equation solutions. a) g
tot
= 1 and b) g
tot
= 0.01. n = 25,
N = 1728, ω
c
= 6∆, T
fm
= 2π/, ∆ = 20.1.
7.3
SD-Approximations
The Lamb-shifts being the only difference between various
SD-approximations, let us investigate them in
detail. The PERLind approach does not utilize a Lamb-shift. [55] Out of curiosity, we add to it the Lamb-shift
from the RWA to see if it improves the accuracy. ULE, on the other hand, does provide the Lamb-shift [56],
that translates to our notation as
H
nm
= E
n
δ
nm
+
X
i
H
ule
(ω
ni
, ω
mi
)A
ni
A
im
, (76)
with the kernel given by
H
ule
(ω, ω
0
) =
1
2π
P
Z
−∞
d
q
γ(Ω + ω)γ(Ω + ω
0
). (77)
The diagonals of the RWA, GAME, and ULE Lamb-shifts are identical. (In ULE the diagonal calculation
involves a Kramers-Kronig transform given in table 1.) So the differences between various Lamb-shifted
SD-approximations are all in the off-diagonal matrix elements of the Lamb-shift.
In terms of numerics, the ULE Lamb-shift requires N
3
numerical integrations, as we could not find a
way to decompose the kernel in Eq. 77 into single-variable functions of frequency, while GAME needs only
N
2
analytical calculations of the spectral functions. After optimizing the codes for ULE and GAME, in the
example n = 25 and N = 1728 it takes approximately 61 hours to compute the Lamb-shift of ULE, compared
to approximately 30 seconds of that in GAME. Similarly, in terms of memory, ULE requires storage of a
rank 3 tensor (Eq. 76), while GAME only stores 2D matrices of SDs and PDs at Bohr frequencies.
In Fig. 13 we present the trace-distance to the solution of the Redfield equation, from the solutions of the
PERLind, PERLind plus LS from the RWA, ULE, and GAME. In Fig. 13(a) the coupling to the heat bath
g
tot
= 1 is the highest one from Fig. 10, while in Fig. 13(b) the coupling is g
tot
= 0.01
GAME and ULE are clearly the closest to the Redfield solutions. The Lamb-shifts of GAME and ULE
are weakly off-diagonal in this example. Nevertheless, the small off-diagonal elements matter significantly,
since PERLind with the added RWA Lamb-shift (which is diagonal) has significantly higher error than both
ULE and GAME.
GAME is slightly closer to the Redfield state than ULE, as can be clearly seen in Fig. 13(b). This is not
surprising because GAME uses the implicit Lamb-shift of the Redfield equation. It is nevertheless surprising
how close the two approximations are, even though the Lamb-shift calculations are quite different.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 26
7.4 Master Equations with Time Dependent Coefficients.
Figure 14: Decay of the perpendicularly magnetized state. a) and b):
Trace-distances between the solutions of equations X and the TDC-Redfield
equation, at g
tot
= 1 and 0.01, respectively. There is a separation in time
scales for divergence of the trace-distance in the Redfield versus TDC-
GAME. The former diverges on a time scale of the bath, while the latter
diverges on a time scale comparable to T
fm
. c) and d): Sum of the
negative eigenvalues of the TDC-Redfield and Redfield solutions versus
time, for g
tot
= 1 and 0.01, respectively. n = 25, N = 1728, ω
c
= 6∆,
T
fm
= 2π/, ∆ = 20.1.
Here we compare the solutions of the
Redfield, GAME, and TDC-GAME
equations to that of the TDC-Redfield
equation (recall diagram in Fig. 2).
Since the TDC-Redfield equation is the
least approximated among the four, it
is likely the most accurate.
In Fig. 14 we study the decay
of the initially perpendicularly mag-
netized state. In particular, (a) and
(b) show numerically obtained trace-
distances versus time at g
tot
= 1
and g
tot
= 0.01, respectively. At
early times, defined here as t < T
fm
,
the trace-distance between the TDC-
Redfield and the Redfield state in-
creases rapidly with a characteristic
time scale much smaller than T
fm
,
independent of g. In comparison,
the trace-distance between the TDC-
Redfield and TDC-GAME increases
with a delay, on a time scale com-
parable to the characteristic system
time (FMR-period), also independent
of g. This suggests there is a sepa-
ration of time scales that govern the
two approximations: one that restores
CP and the other that replaces TDC-
coefficients with their asymptotic val-
ues. We interpret the rapid initial in-
crease in the trace-distance on time
scale T
fm
in terms of the initial state
slip we mentioned earlier. [7074]
Note that now is not the first time that a CP-restoration is governed by a system-specific time-scale,
that is independent of the system-bath coupling strength. The RWA is a well known example, where that
time-scale is the Heisenberg time. In GAME CP restoration, the time scale is similar to the characteristic
system time, in this case the FMR-period, which is much smaller than the Heisenberg time. This suggests
that GAME may be most suitable for applications on quantum systems with a dominant system frequency,
which is the case in ferromagnets or superconductors, or any other gapped system.
Other coarsegraining approximations such as CGSE and PRWA work quite differently. They usually
involve some optimization process that depends on g. For example, the CGSE equation has optimum coarse-
graining time of order
τ
r
τ
c
but it seems that there is no system dependence of the optimum time. Thus
CGSE may be applicable even if 1
r
is larger than the dominant system frequency, if ω
c
is large enough.
We note that in such regime the Redfield equation is unlikely to be stable. The question of which master
equation works best in what regimes is a complicated matter and will remain the subject of future reasearch.
Fig. 14(c,d) presents the sum of the negative eigenvalues of the solutions of the Redfield equations,
corresponding to g
tot
= 1 and 0.01, respectively. In (c), both Redfield equations develop significant negativity.
This would suggest that there is a breakdown of the perturbative approach, according to the analysis by
Hartmann and Strunz. [22] It suggests that quantum theory of magnetics may need to extend beyond the
master-equation platform and involve regime of strong coupling to the heat bath. On the other hand, in
Fig. 14(d), the TDC-Redfield equation has strongly suppressed negativity, which implies the validity of the
weak coupling regime.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 27
Figure 15: Decay of the ferromagnetic resonance. a) and b): Trace-
distances between the solutions of equations X and the TDC-Redfield
equation, at g
tot
= 1 and 0.01, respectively. c) and d): Sum of the
negative eigenvalues of the TDC-Redfield and Redfield solutions versus
time, for g
tot
= 1 and 0.01, respectively. n = 25, N = 1728, ω
c
= 6∆,
T
fm
= 2π/, ∆ = 20.1.
We have investigated the properties
displayed in Fig. 14 for different ini-
tial conditions. In the specific exam-
ple of the lowest excited state (FMR),
the TDC-GAME is always closest to
the TDC-Redfield state as shown in
Fig. 15(a,b), in contrast to Fig. 14(a,b)
where TDC-GAME and Redfield trace-
distances cross much sooner.
The asymptotic GAME is also
closer to the TDC-Redfield state than
the Redfield state at g
tot
= 1. However,
at g
tot
= 0.01, GAME and the Red-
field states have essentially the same
distance to the TDC-Redfield state. At
this time it is too early to comment on
this, but it suggests that the Redfield
and GAME approximations are on sim-
ilar, if not the same, level of approxima-
tion. Negativity of the TDC-Redfield
state is strongly suppressed at g
tot
=
0.01, analogous to that in Fig. 14(d).
More exotic initial states including
the Dicke and the Greenberger Horne
Zeilinger state will be the subject of fu-
ture reasearch, to establish which equa-
tion is most suitable to describe quan-
tum correlations and their decay. We would like to note at this time that the perpendicularly magnetized
state, which is initially a product state, evolves into a nearly maximally entangled state at approximately one
period of the FMR. The study of multipartite entanglement determined by these equations is well outside
the scope of the present paper, but a hint is given in Appendix 11.
8 Conclusion
The key to an effective master equation is in the combination of high accuracy and low complexity, and GAME
accomplishes both with a complete positivity guarantee. The accuracy is the result of the combination of
two steps involved in approximating the Redfield equation into a complete positive form. The first step is
the renormalization of the system Hamiltonian in the Redfield equation, due to the coupling to the heat
bath. This renormalization is done before any approximation, to ensure balance between dissipative gains
and losses of the state. The second step is to approximate the gain and loss terms to restore complete
positivity. This is accomplished by swapping out an arithmetic mean of the SD with a geometric mean of
the SD. The difference between the two means amounts to a transient term in the master equation. The
error accumulated due to its neglect saturates on time scale of the transient, which is on the order of typical
system frequency.
We test GAME on a 25-site ferromagnetic Heisenberg spin-chain with dipole-dipole magnetic anisotropy,
by numerically simulating the relaxation of an initially perpendicularly magnetized state. The density matrix
is solved directly, by iterating the master equation. We compare GAME to several other master equations,
and find that GAME produces the closest state to the solution of the Redfield equation.
Over the past several years, we have studied GAME in time dependent driven open quantum systems.
In fact, we derived GAME initially for time-dependent Hamiltonians, but reduced the scope of this paper
to time-independent ones. The extension of GAME to time-dependent Hamiltonians is both straightforward
and retains the simplicity of the (now) Floquet-Redfield equation. In appendix 11 we sketch the extension
by presenting an example of the ferromagnetic Heisenberg spin-chain driven by a transverse time dependent
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 28
magnetic field. We study how the periodic field drives this uniaxial magnet into a long-range quantum
coherent state at zero temperature.
Possible future applications of GAME include the following:
As an extension of GAME, we plan to investigate quantum many-body systems with moderately strong
coupling to the heat-bath. In such regime, τ
c
τ
r
1/E
max
, with E
max
the largest Bohr frequency
of the system. These studies aim to understand overdamped ground states and quantum dynamics in
solid-state systems like magnets and superconducting quantum circuits. Although we expect that the
Redfield equation in that regime develops states with significant negative eigenvalues and a possible
instability, there is no fundamental reason not to investigate such systems within the master-equation
paradigm, because the Born-Markov approximation is still valid. A specific question is what type of
CP-restoration of the Redfield theory produces a steady state that is a reasonable approximation of
the reduced ground state, broadened by significant vacuum fluctuations, as well as the dynamics above
that state. Possibly, this investigation could later on be extended to the true strong-coupling regime,
by including the Hamiltonian renormalization and CP restoration into a larger unitary transformation
such as the polaron transformation [75].
Practically, GAME opens possibilities to model solid state devices that utilize long range entanglement
as resource. Complete positivity cures the headache of negative eigenvalues before taking the partial
transpose, while low complexity of the master equation enables studies of large number of entangled
particles. One challenge will be to determine the local-local and local-global spin correlation functions
on density matrices such as those in appendix 11. Quantum regression theorem is well established
and can be applied in that task. [1, 2] These computations will be relevant for quantum sensing of
spins without an immediate spin-charge conversion, opening a way to possible applications in neutron
scattering spectroscopy, high-energy physics, and gravitational astrophysics.
Another direction will be to apply well-established density matrix renormalization group (DMRG)
methods to the study of finite and infinite spin systems with long range 1/(r
α
) interactions. The
infinite DMRG has been around for 15 or so years, and such applications have been well developed. [76]
Tensor methods will apply well-established infinite DMRG methods to the study of this system out of
equilibrium. [77]
Beyond quantum-information science, it is our wish that GAME will find appeal in broader range of
fields, including quantum chemistry, energy transfer dynamics, and fuel efficiency. GAME is the only quantum
master equation that assures complete positivity and retains the simplicity and accuracy of the widely applied
63-year old Redfield equation.
We thank Jason Dark, Elyana Crowder, Brian Kennedy, and Glen Evenbly for discussions, and Ugo
Marzolino for a comment about the discontinuity of the RWA pertinent to Sec 5. This research was supported
by DOE contract DE-FG02-06ER46281, including development of the numerical simulations of quantum
dynamics in magnets. Additional support from the Georgia Tech Quantum Alliance (GTQA), a center
funded by the Georgia Tech Institute of Electronics and Nanotechnology was used to develop exact numerical
methods on entangled low-dimensional magnetic systems.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 29
9 Appendix: Numerical Iteration of Master Equations
We discuss how we solve a master equation
dt
= L(ρ), (78)
within an arbitrary error () in norm distance, given the initial condition ρ(0). First we consider the case
where L is a time independent superoperator. The algorithm we use simulates the power series expansion of
e
Lt
, e.g., as follows. Given a state ρ(t) at time t, and time increment dt, let ρ
d
= ρ
0
= ρ(t). Start the loop
m = 1, 2, ... ,. For each m, let ρ
d
L{ρ
d
}, followed by ρ
0
ρ
0
+ (dt
m
)ρ
d
/m!. Terminate the loop when
||(dt
m
)ρ
d
/m!||
2
< . For long simulations we use the threshold = 10
7
, which leads to a typical error per
step in the 10
9
range. For the ferromagnetic chain we use the time increment of dt = T
fm
/64.
In case of TDC master equations, we loose the above property. So we apply Runge-Kutta RK4 method
in the interaction picture, and use the same time increment dt. (That is, a time step at time t calculates the
RK4-coefficients at times t, t + dt/2, and t + dt.)
10 Appendix: Saturation Time versus Bath Cutoff Frequency.
Trace distance between the
solutions of the Redfield
and GAME master equa-
tions versus time, at three
values of ω
c
. n = 23,
N = 694, = 20.1,
T
fm
= 2π/. g
tot
=
0.0092, 0.0055, and 0.0033,
at ω
c
= 6∆, 60∆, and 600∆,
respectively. The satura-
tion time (T
s
) of the trace-
distance is weakly depen-
dent on ω
c
in this range.
Fits to a single exponen-
tial (not shown) lead to
T
s
= 1.4/1.2, at ω
c
=
6∆/60∆.
11 Appendix: Extension Of GAME to time-dependent Hamiltonian
GAME can be straightforwardly extended to study the dynamics of open quantum systems with time depen-
dent system Hamiltonian. In that case, H
0
in Eq. 1 is replaced with time dependent periodic Hamiltonian
H
0
(t) with period T . The periodicity does not necessarily narrow the scope. For example, if we want to
study the effect of pulses, we could make the period of H
0
(t) much longer than other relevant time scales
and rely on fast Fourier transform algorithm to rapidly calculate the spectral functions as shown below. We
generalize GAME by rotating into the Floquet basis.
Here consider a 19-site spin-chain with the same exchange energy and anisotropy as in the main text. In
this example, the dimension of the truncated Hilbert space is N = 448, while H
0
(T ) = H
0
z
S
x
h(t), where
h(t) is a shape function bound by one, and
z
is the Zeeman splitting at h(t) = 1.
We apply a time dependent square wave of magnetic field perpendicular to the chain axis with period
T = 10T
fm
. In experimental setting the frequency would be roughly 1GHz, for typical magnets used in
magneto-electronics. Fig. 16 displays one simulated magnetic field period, with amplitude in units of the
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 30
Figure 16: Dynamics of a driven ferromagnetic spin-chain switches between mostly dissipative and mostly unitary during
one period of the drive. a) A nearly square field wave is applied perpendicular to the chain with period 10 × T
fm
. b)
X-component of the magnetization versus time. c) and d): The magnitude of the density matrix at half-period is displayed
in the S
x
and S
y
basis and demonstrates coherence between oppositely magnetized states.
FMR gap . The magnetic field switches from high to zero smoothly, is zero for approximately half the
period, and then it switches back to the large value. In the high field, the ground state magnetization in the
adiabatic approximation is close to aligned with the applied field.
We use a super-Ohmic spectral density given in Table 1, which is appropriate for relaxation mediated by
the phonon bath, with the same cutoff frequency as for the Ohmic bath (6∆) at zero temperature, and the
coupling constant g 0.2660 per heat bath. At this g, the super-Ohmic spectral density leads to a similar
energy relaxation time at zero field, as in the example studied in Secs. 7.1 and 7.3 at g
tot
= 1.
The Floquet theorem expresses the states of the system as
|φ(t)i = e
i
φ
t
|u
φ
(t)i, (79)
where
φ
are the quasi-energies and |u
φ
(t)i = |u
φ
(t + T )i. The quasi-energies are determined by finding
the eigenvalues e
i
φ
T
of the unitary time evolution operator (U) for one period of the drive. We find U
numerically using RK4 with a time increment of 3.8147 × 10
6
T and verify that ||UU
1||
1
< 10
11
.
Next, the density matrix and operator A are represented in the Floquet basis, hu
φ
|ρ|u
ξ
i and hu
φ
|A|u
ξ
i,
respectively. As such, A becomes time dependent and periodic, so we write it down as A(t). The Floquet
Redfield equation [59, 78, 79] generalizes the Redfield Eq. 4:
dt
= i[diag(
φ
), ρ] A(t)A
f
(t)ρ ρA
f
(t)A(t) + A
f
(t)ρA(t) + A(t)ρA
f
(t). (80)
A
f
(t) is the time dependent filtered operator, which is obtained as
A
f
(t) = IFFT
Γ
m
FFT
A(t)

. (81)
FFT and IFFT are the forward and inverse discrete Fourier transforms in MATLAB, respectively. The
frequencies we use are (m 1)ω
p
, m = 1, 2, .., 2048, where ω
p
= 2π/T . Γ
m
is a matrix, related to the
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 31
half ranged integral Fourier transform of the bath correlation function given by Eq. 5. That is, Γ
mφξ
=
Γ(m
0
ω
p
+
φ
ξ
). FFT in MATLAB does not center frequency around zero. The easiest way to fix this issue
is to shift m as this: m
0
= m, if m 1025, and m
0
= m 2049 if m 1026. Virtually the entire spectrum
of A(t) is accounted for within our double precision calculation, because the norm of FFT[A(t)] at the half
frequency is approximately 10
12
of that at zero frequency (not shown).
To establish the time-dependent version of GAME, the first step is to extract the Lamb-shift without
changing the Floquet-Redfield equation. We find a time dependent renormalized Hamiltonian, by generalizing
Eq. 8,
H(t) = diag(
Φ
)
i
2
[A(t)A
f
(t) A
f
(t)A(t)] diag(
Φ
) + H
L
(t), (82)
in the Floquet basis. The next step is to approximate the kernel (implicit in Eq. 80), and we do this by
generalizing the derivation of Eq. 41. That is, we calculate a time dependent Lindblad generator
L(t) = IFFT
γ
m
FFT
A(t)

. (83)
γ
m
is the SD matrix with elements γ
mφξ
= γ(m
0
ω
p
+
φ
ξ
). In this case the transformation of the
time dependent master equation into a Lindblad form is assisted by the FFT-convolution theorem. This
approximation sequence leads to a time dependent GKSL-equation,
dt
= i[H(t), ρ] +
L(t)
ρL(t)
1
2
{L(t)L(t)
, ρ}
. (84)
.
The magnetization x-component versus time is shown in Fig. 16(b). When the magnetic field is high, the
level spacing of the system is high due to Zeeman splitting. The chain rapidly relaxes towards the quantum
ground state, facilitated by the enhanced super-Ohmic spectral density. We could say that the high field
part of the wave constitutes the state preparation stage.
At zero magnetic field, the level spacing is suppressed and set by the anisotropy , leading to a much
longer τ
r
. Analogous enhancement of τ
r
at low field is well established in spin-1/2 quantum dots. [80, 81]
Now we arrive to the key observation. The magnetic dynamics at zero magnetic field quickly creates a highly
correlated many-body state. At a half period, Figs. 16(c,d) display the magnitude of the density matrix
in the Sx and S
y
eigenbasis, respectively, within the truncated Hilbert space. The off-diagonal elements
connect states with opposite magnetization and are spread out over the entire range of Sx and Sy. Similar
coherences are found when the density matrix is represented in the Sz basis (not shown). It will be interesting
to calculate various correlation functions, to test if they may violate a Bell inequality in this macroscopic
setting.
12 Appendix: Time-Dependent Spectral Functions
The decomposition of the kernel in Eq. 71 into single-variable functions proceeds as discussed in appendix
F of Ref. [50]. After the decomposition, the dissipative kernel 72 has time-dependent coefficients γ
t
(ω) and
S
t
(ω) defined as
γ
t
(ω) =
1
π
Z
−∞
γ
ω +
x
t
sinc(x)dx (85)
and
S
t
(ω) =
1
π
Z
−∞
γ
ω +
x
t
(sin
x
2
)
2
x
dx. (86)
First we show that these are the time-dependent SD and PD introduced in Eq. 17.
Proof. The time-dependent spectral functions γ
t
(ω) and S
t
(ω) are defined as
Γ
t
(ω) =
1
2
γ
t
(ω) + iS
t
(ω) =
Z
t
0
C(τ)e
τ
. (87)
Expressing the bath-correlation function as the inverse Fourier transform of the SD, as shown in Table 1,
and changing the order of integrals, we find
Γ
t
(ω) =
1
2π
Z
−∞
γ(Ω)d
Z
t
0
e
i[(ωΩ)τ]
=
t
2π
Z
−∞
γ(Ω)e
i
[
(ωΩ)
t
2
]
sinc
(ω Ω)
t
2
d. (88)
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 32
After substituting Ω = ω + x/τ, this equation leads to
1
2
γ
t
(ω) + iS
t
(ω) =
1
2π
Z
−∞
γ
ω +
x
t
e
i
x
2
sinc
x
2
dx. (89)
By taking the real and imaginary parts, we arrive at Eqs. 85 and 86, respectively. QED.
Now we demonstrate the equivalency of Eqs. 74 and Eq. 75.
Proof: Introduce the real functions R
t
(ω) and W
t
(ω),
R
t
(ω) =
1
π
Z
−∞
S
ω +
x
t
sinc(x)dx (90)
and
W
t
(ω) =
1
π
Z
−∞
S
ω +
x
t
(sin
x
2
)
2
x
dx. (91)
Notice that Eqs. 71 and 74 map to each other, if we swap out the SD with PD, and change the sign of the
exponential in the prefactor. We apply this mapping to Eq. 72, by swapping [γ
t
(ω), S
t
(ω)] with [R
t
(ω), W
t
(ω)]
and changing the sign in the exponential. This leads to the unitary kernel in terms of the single variable
functions of frequency R
t
(ω) and W
t
(ω):
H
dc
t
(ω, ω
0
) = e
i
(ωω
0
)τ
2
(
1
2
[R
τ
(ω
0
) + R
τ
(ω)]sinc
(ω
0
ω)τ
2
W
τ
(ω
0
) W
τ
(ω)
(ω
0
ω)
τ
2
cos
(ω
0
ω)τ
2
)
(92)
To evaluate R
t
(ω) and W
t
(ω), express the PD in terms of the bath correlation function (Table 1):
R
t
(ω) =
1
π
Im
Z
−∞
sinc(x)dx
Z
0
C(τ)e
i(ωτ+
τ
t
x)
, (93)
and
W
t
(ω) =
1
π
Im
Z
−∞
(sin
x
2
)
2
x
dx
Z
0
C(τ)e
i(ωτ+
τ
t
x)
. (94)
After some algebra, changing the order of integration, and applying the identity
Z
−∞
sin(ax) sin
2
(x/2)/xdx =
1
2
Z
−∞
cos(ax)sinc(x)dx = Θ(1 |a|)sign(a)π/2, (95)
we arrive at:
R
t
(ω) = S
t
(ω), (96)
W
t
(ω) =
1
4
γ
t
(ω). (97)
Substituting these into Eq. 92, we obtain Eq. 75. QED.
For the Ohmic bath with exponential frequency cutoff, the time dependent spectral functions can be
determined explicitly,
Γ
t
(ω) = igω
c
(
1
e
t
1 +
c
t
ω
ω
c
e
ω
ω
c
Ei(
ω
ω
c
) Ei(
ω
ω
c
+ t) Θ(
ω
ω
c
)
)
. (98)
However, we find it is faster to evaluate the spectral functions by numerical integration of the correlation
function, using the function "integrate" in MATLAB. Hence, after verifying that the numerical and analytical
results agree, we numerically obtain all the results related to the DCG-approximation.
To check the validity of these transforms, we have solved the Redfield, GAME, and DCG master equations
for Ohmic spectral density with the Drude-Lorentz cutoff, presented in table 1. In this calculation, we find
the time-dependent SD and PD directly from the asymptotic SD according to Eqs. 85 and 86. That way
we bypass the difficult correlation function of the Drude-Lorentz SD. Similarly, we calculate the functions
Z
t
(ω) and W
t
(ω) by numerical integration of Eqs. 90 and Eqs. 91. Since the Drude PD is a relatively simple
function of frequency as shown in table 1, and contains no special functions, these integrals are relatively fast
to evaluate. That is, we do not take advantage of the explicit unitary kernel given by Eq. 75, but integrate
the asymptotic coefficients over frequency. The results of the simulation of the system state versus time are
displayed in Fig. 17, showing a very reasonable agreement with Fig. 12.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 33
Figure 17: The same as Fig. 12 but for Ohmic SD with the Drude-Lorentz cutoff. Trace-distances between the solutions
of various master equations and that of the Redfield equation, versus time. a) Regime of moderate environmental coupling
g
tot
= 69g = 0.92. Thick lines: RWA (dashed-blue), DCG
opt
(dashed-black), and GAME (red-full). The approximate
optimum coarsegraining time is τ
opt
= 1.3. Thin-full lines: coarse-graining times τ = 64, 32, 16, 6.4, 3.2, and 1.6, follow
the maxima between the RWA and DCG
opt
from top to bottom, respectively. Thin-dashed lines: τ = 0.96, 0.64, 0.32, and
0.16, follow the maxima from DCG
opt
bottom to top, respectively. b) Analogous to a) in the weak environmental coupling
regime g
tot
= 0.0092. τ
opt
= 16. Thin-full lines, top-to-bottom: τ = 160, 64, and 32. Thin-dashed lines, bottom-to-top:
τ = 6.4, 3.2, 1.6, and 0.96. All time scales are in units of T
fm
. n = 23, N = 694, ∆ = 20.1, and ω
c
= 6∆.
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 34
References
[1] C. W.Gardiner, Quantum noise,” (1991).
[2] H.-P. Breuer and F. Petruccione, The theory of open quantum systems,” (2007).
[3] G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
[4] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976),
https://aip.scitation.org/doi/pdf/10.1063/1.522979 .
[5] D. A. Lidar, Z. Bihary, and K. Whaley, Chemical Physics 268, 35 (2001).
[6] D. Kohen, C. C. Marston, and D. J. Tannor, The Journal of Chemical Physics 107, 5236 (1997),
https://doi.org/10.1063/1.474887 .
[7] R. Zwanzig, The Journal of Chemical Physics 33, 1338 (1960), https://doi.org/10.1063/1.1731409 .
[8] S. Nakajima, Progress of Theoretical Physics 20, 948 (1958), https://academic.oup.com/ptp/article-
pdf/20/6/948/5440766/20-6-948.pdf .
[9] K. Ryogo, T. Morikazu, and H. Natsuki, Statistical physics ii,” (1998).
[10] A. G. Redfield, IBM Journal of Research and Development 1, 19 (1957).
[11] A. REDFIELD, in Advances in Magnetic Resonance, Advances in Magnetic and Optical Resonance,
Vol. 1, edited by J. S. Waugh (Academic Press, 1965) pp. 1 32.
[12] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
[13] W. T. Pollard, A. K. Felts, and R. A. Friesner, “The redfield equation in condensed-phase quan-
tum dynamics,” in Advances in Chemical Physics (John Wiley & Sons, Ltd, 2007) pp. 77–134,
https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470141526.ch3 .
[14] I. Kondov, U. Kleinekathöfer, and M. Schreiber, The Journal of Chemical Physics 114, 1497 (2001),
https://doi.org/10.1063/1.1335656 .
[15] D. Egorova, M. Thoss, W. Domcke, and H. Wang, The Journal of Chemical Physics 119, 2761 (2003),
https://doi.org/10.1063/1.1587121 .
[16] M. Schröder, M. Schreiber, and U. Kleinekathöfer, Journal of Luminescence 125, 126 (2007), festschrift
in Honor of Academician Alexander A. Kaplyanskii.
[17] A. Montoya-Castillo, T. C. Berkelbach, and D. R. Reichman, The Journal of Chemical Physics 143,
194108 (2015), https://doi.org/10.1063/1.4935443 .
[18] C. Timm, Phys. Rev. B 77, 195416 (2008).
[19] J. Jeske, D. J. Ing, M. B. Plenio, S. F. Huelga, and J. H. Cole, The Journal of Chemical Physics 142,
064104 (2015), https://doi.org/10.1063/1.4907370 .
[20] W. P. Bricker, J. L. Banal, M. B. Stone, and M. Bathe, The Journal of Chemical Physics 149, 024905
(2018), https://doi.org/10.1063/1.5036656 .
[21] R. S. Whitney, Journal of Physics A: Mathematical and Theoretical 41, 175304 (2008).
[22] R. Hartmann and W. T. Strunz, Phys. Rev. A 101, 012103 (2020).
[23] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Phys. Rev. A 60, 91 (1999).
[24] I. de Vega, D. Alonso, P. Gaspard, and W. T. Strunz, The Journal of Chemical Physics 122, 124106
(2005), https://doi.org/10.1063/1.1867377 .
[25] N. Makri and D. E. Makarov, The Journal of Chemical Physics 102, 4600 (1995),
https://doi.org/10.1063/1.469508 .
[26] M. Thorwart, E. Paladino, and M. Grifoni, Chemical Physics 296, 333 (2004), the Spin-Boson Problem:
From Electron Transfer to Quantum Computing ... to the 60th Birthday of Professor Ulrich Weiss.
[27] P. Nalbach and M. Thorwart, Phys. Rev. B 81, 054308 (2010).
[28] Y. Tanimura and R. Kubo, Journal of the Physical Society of Japan 58, 1199 (1989),
https://doi.org/10.1143/JPSJ.58.1199 .
[29] Y. Tanimura, Journal of the Physical Society of Japan 75, 082001 (2006),
https://doi.org/10.1143/JPSJ.75.082001 .
[30] Y. Tanimura, The Journal of Chemical Physics 141, 044114 (2014), https://doi.org/10.1063/1.4890441
.
[31] Z. Li, N. Tong, X. Zheng, D. Hou, J. Wei, J. Hu, and Y. Yan, Phys. Rev. Lett. 109, 266403 (2012).
[32] Y. Cheng, W. Hou, Y. Wang, Z. Li, J. Wei, and Y. Yan, New Journal of Physics 17, 033009 (2015).
[33] H.-D. Meyer, U. Manthe, and L. Cederbaum, Chemical Physics Letters 165, 73 (1990).
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 35
[34] M. Beck, A. Jäckle, G. Worth, and H.-D. Meyer, Physics Reports 324, 1 (2000).
[35] H. Wang and M. Thoss, The Journal of Chemical Physics 119, 1289 (2003),
https://doi.org/10.1063/1.1580111 .
[36] J. Zheng, Y. Xie, S. Jiang, and Z. Lan, The Journal of Physical Chemistry C 120, 1375 (2016),
https://doi.org/10.1021/acs.jpcc.5b09921 .
[37] D. Suess, A. Eisfeld, and W. T. Strunz, Phys. Rev. Lett. 113, 150403 (2014).
[38] P.-P. Zhang and A. Eisfeld, The Journal of Physical Chemistry Letters 7, 4488 (2016), pMID: 27775345,
https://doi.org/10.1021/acs.jpclett.6b02111 .
[39] R. Hartmann and W. T. Strunz, Journal of Chemical Theory and Computation 13, 5834 (2017), pMID:
29016126, https://doi.org/10.1021/acs.jctc.7b00751 .
[40] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Nature Communications 9, 3322
(2018).
[41] F. A. Y. N. Schröder, D. H. P. Turban, A. J. Musser, N. D. M. Hine, and A. W. Chin, Nature
Communications 10, 1062 (2019).
[42] S. Jang, The Journal of Chemical Physics 131, 164101 (2009), https://doi.org/10.1063/1.3247899 .
[43] D. P. S. McCutcheon, N. S. Dattani, E. M. Gauger, B. W. Lovett, and A. Nazir, Phys. Rev. B 84,
081305 (2011).
[44] E. B. Davies, Comm. Math. Phys. 39, 91 (1974).
[45] C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, Phys. Rev. A 88, 012103 (2013).
[46] D. Farina and V. Giovannetti, Phys. Rev. A 100, 012107 (2019).
[47] E. Mozgunov and D. Lidar, 4, 227 (2020), 1908.01095 [Quantum] .
[48] N. Vogt, J. Jeske, and J. H. Cole, Phys. Rev. B 88, 174514 (2013).
[49] T. V. Tscherbul and P. Brumer, The Journal of Chemical Physics 142, 104107 (2015),
https://doi.org/10.1063/1.4908130 .
[50] G. Schaller and T. Brandes, Phys. Rev. A 78, 022106 (2008).
[51] F. Benatti, R. Floreanini, and U. Marzolino, Phys. Rev. A 81, 012105 (2010).
[52] W. J. Munro and C. W. Gardiner, Phys. Rev. A 53, 2633 (1996).
[53] J. Wilkie, Phys. Rev. E 62, 8808 (2000).
[54] B. Palmieri, D. Abramavicius, and S. Mukamel, The Journal of Chemical Physics 130, 204512 (2009),
https://doi.org/10.1063/1.3142485 .
[55] G. Kiršanskas, M. Franckié, and A. Wacker, Phys. Rev. B 97, 035432 (2018).
[56] F. Nathan and M. S. Rudner, Phys. Rev. B 102, 115109 (2020).
[57] S. Kryszewski and J. Czechowska-Kryszk, Master equation - tutorial approach,” (2008),
arXiv:0801.1757 [quant-ph] .
[58] H.-P. Breuer, Phys. Rev. A 70, 012106 (2004).
[59] D. W. Hone, R. Ketzmerick, and W. Kohn, Phys. Rev. E 79, 051129 (2009).
[60] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, New Journal of Physics 14, 123016 (2012).
[61] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82,
1155 (2010).
[62] A. Daley, Advances in Physics 63 (2014), 10.1080/00018732.2014.933502.
[63] F. Benatti, R. Floreanini, and U. Marzolino, EPL (Europhysics Letters) 88, 20011 (2009).
[64] A. Rivas, Phys. Rev. A 95, 042104 (2017).
[65] F. Benatti, R. Floreanini, and M. Piani, Phys. Rev. Lett. 91, 070402 (2003).
[66] R. Tana and Z. Ficek, Journal of Optics B: Quantum and Semiclassical Optics 6, S90 (2004).
[67] S. E. Clifton and W. E. P., Philosophical Transactions of the Royal Society of London 240, 599 (1948).
[68] W. F. Brown, Phys. Rev. 130, 1677 (1963).
[69] K. Gilmore, Y. U. Idzerda, and M. D. Stiles, Phys. Rev. Lett. 99, 027204 (2007).
[70] F. Haake and M. Lewenstein, Phys. Rev. A 28, 3606 (1983).
[71] P. Gaspard and M. Nagaoka, The Journal of Chemical Physics 111, 5668 (1999),
https://doi.org/10.1063/1.479867 .
[72] Y. C. Cheng and R. J. Silbey, The Journal of Physical Chemistry B 109, 21399 (2005), pMID: 16853776,
https://doi.org/10.1021/jp051303o .
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 36
[73] A. Suarez, R. Silbey, and I. Oppenheim, The Journal of Chemical Physics 97, 5101 (1992),
https://doi.org/10.1063/1.463831 .
[74] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Physics Letters A 265, 331 (2000).
[75] R. Silbey and R. A. Harris, The Journal of Chemical Physics 80, 2615 (1984),
https://doi.org/10.1063/1.447055 .
[76] Z.-X. Gong, M. F. Maghrebi, A. Hu, M. Foss-Feig, P. Richerme, C. Monroe, and A. V. Gorshkov, Phys.
Rev. B 93, 205115 (2016).
[77] G. Evenbly and G. Vidal, Phys. Rev. Lett. 115, 180405 (2015).
[78] M. Grifoni and P. Hänggi, Physics Reports 304, 229 (1998).
[79] T. Shirai, J. Thingna, T. Mori, S. Denisov, P. Hänggi, and S. Miyashita, New Journal of Physics 18,
053008 (2016).
[80] J. Elzerman, R. Hanson, and L. W. van Beveren et al., Nature 430, 431–435 (2004).
[81] A. Morello, J. Pla, and F. Z. et al., Nature 467, 687 (2010).
Accepted in Quantum 2020-08-20, click title to verify. Published under CC-BY 4.0. 37