<EMAIL_ADDRESS>
# Hamiltonian Model for Fault Tolerant Singlet-Like Excitation: First
Principles Approach

Donghyun Jin  Georgia Institute of Technology, Atlanta, Georgia, United States
Grihith Manchanda  Georgia Institute of Technology, Atlanta, Georgia, United
States  Dragomir Davidović  Georgia Institute of Technology, Atlanta, Georgia,
United States  [

###### Abstract

Deriving quantum error correction and quantum control from the Schrödinger
equation for a unified qubit-environment Hamiltonian will give insights into
how microscopic degrees of freedom affect the capability to control and
correct quantum information beyond that of phenomenological theory. Here, we
investigate the asymptotic reduced state of two qubits coupled to each other
solely via a common heat bath of linear harmonic oscillators and search for
evidence of fault-tolerant excited qubit states. We vary the Hamiltonian
parameters, including the qubit-qubit and qubit-bath detuning, the bath
spectral density, and whether or not we use the Markov approximation in the
calculation of our dynamics. In proximity to special values of these
parameters, we identify these states as asymptotic reduced states that are
arbitrarily pure, excited, unique, and have high singlet fidelity. We
emphasize the central role of the Lamb-shift as an agent responsible for fault
tolerant excitations. To learn how these parameters relate to performance, we
discuss numerical studies on fidelity and error recovery time.

##  1 Introduction

Long-lasting coherence of quantum states is critical for quantum computation
and quantum information processing. [1] Alas, environmental coupling in
realistic quantum systems leads to rapid decoherence of quantum states.
Surmounting this decoherence and other quantum information errors like damping
remains the main challenge for the experimental realization of quantum
computers. Quantum error correction (QEC) [2, 3] opens a way to overcome
decoherence and damping. It works by encoding quantum states into quantum
codes, performing projective measurements of error syndromes, and recovering
the true state via fast gates. Other approaches to fight decoherence include
quantum feedback control for a continuously measured system [4, 5, 6] and
combination approaches. [7]

A measurement apparatus can always be regarded as part of a larger isolated
quantum system subject to unitary dynamics. At this point, we know that
overcoming decoherence can be derived from that unitary framework without
invoking any measurements. Indeed, coherent feedback control is an example of
that. It includes a Langevin equation approach to feedback mediated by
continuous measurements, which can be implemented without measurements. [4]
Unitary interactions alone between quantum systems can be utilized to
implement feedback control. [8, 9] Another example where explicit measurements
are not necessary comes from the principle of implicit measurement in quantum
computing and QEC, which states that at the end of computation, any qubits
which have not been measured may be assumed to have been measured. [1] The
result is that QEC can be performed without explicit measurement.

Much of the theory of QEC and quantum control is not based on first-principle
derivations. Rather, it is based upon phenomenological error models that begin
from a ”higher level, effective description, most notably that of Markovian
dynamics.” [10] As an illustration, let us consider two qubits immersed in a
heat bath of linear harmonic oscillators interacting solely via the bath,
sketched in Fig. 1 a). We want to determine if the reduced state of the qubits
in equilibrium can be all pure, excited, and fault-tolerant at the same time.
A black-box circuit model of the qubit-environment system that could lead to
such a state could be like the one sketched in Fig. 1 b). Fault-tolerance
would be a special interaction with the bath, where the qubit state is
teleported, transduced to the bath degrees of freedom, and swapped out by a
reverse-transduced and reverse-ported ancillary singlet. By the principle of
implicit measurement, quantum signal passed to the bath need not be measured
in order for this to work.

However, without the microscopic Hamiltonian of a system-bath universe, the
circuit model is only useful as a calculation tool, rather than a way to gain
insight into questions like: What is the recovery time from errors? How does
the recovery time affect state fidelity? What is the effect of the bath
spectral density on fidelity, especially at low frequencies? What are the
roles of the bath counterterms? The Hamiltonian model presented in this paper
opens a way to answering all these questions. Other approaches include quasi-
phenomenological descriptions of fault tolerant QEC [11, 12, 13, 14] and
continuous quantum feedback, [7, 15, 16] and fully microscopic descriptions of
dynamical decoupling. [17, 18]

Figure 1: a) Two qubits immersed in a common heath bath of linear harmonic
oscillators at zero temperature. b) Quantum circuit model for a heat bath as a
stabilizer of singlet. On the left, the inputs are an arbitrary two-qubit
state $|\Psi\rangle$ and the heat bath vacuum state $|0\rangle$ (red, lower-
left). On the right, in the qubit system, the output is the Bell or singlet
state $|\Psi^{-}\rangle=|S\rangle$, while in the bath the output is implicitly
measured (blue, lower-right). The large blue box is a unitary gate.

In most situations, the bath causes a damping in a system that it has contact
with, which leads to finite lifetimes of the system excited states. Our
central claim about the fault tolerant excited state (FTES) for a two-qubit
system is summarized as follows: by adjusting both the detuning of the qubits
and the asymmetry in their environmental coupling, an excited two-qubit state
can be decoupled from the dissipative effect of the bath. In the case of
equilibrium and zero temperature, the qubits are permanently in that pure and
excited state. Any fluctuation of the state is recovered by the action of the
bath, which means that the state is fault tolerant.

Additionally, the two lowest excited states of the two qubits are usually near
two-fold degenerate, making it easy to generate a highly entangled state by a
unitary transformation. This is provided here by the unitary environmental
effect. The net result is that the bath mediated interaction places the qubits
into a permanent, nearly maximally entangled, excited state.

We find the fault tolerant state not by carefully assembling the system-bath
to get the desired result, but by sweeping the microscopic parameters to see
if and when FTES emerges naturally on its own. The ab initio approach
distinguishes our work from previous works.

The advancements made in this paper were made possible by recent developments
of completely positive quantum master equations (MEs), namely the geometric-
arithmetic master equation (GAME), [19] similar to the universal Lindblad
equation (ULE) [20]. These two equations were developed independently and at
the same time. While they have identical dissipative dynamics, they differ in
the Lamb-shift. They both claim the same accuracy in the order for a system-
environment coupling, equivalent to that of the Redfield master equation (RE).
Indeed, a comparison of the solutions of ULE and GAME with the solutions of RE
leads to nearly identical errors with respect to RE solutions, [19] as well as
errors comparable to the estimated accuracy of RE. The Lamb-shift in ULE is
more self-consistent with the dissipator, but more difficult to calculate. The
mathematical form of the Lamb-shift is simpler in GAME, making implementation
faster. Rapid evaluation of the Lamb-shift will be paramount in extending FTES
to many-qubit systems.

Throughout the paper, we use the following pairs of words interchangeably:
heat bath and environment, and system and qubits. The outline of the paper is
as follows. We begin by introducing the system-environment Hamiltonian in Sec.
2. The dynamics of the reduced qubit system is derived from the Hamitonian in
Sec. 3, in which we discuss the four calculation methods we implemented. GAME
and its logic is discussed more succinctly relative to Ref. [19]. The Redfield
and time-convolutionless master equation (TCL4) are also reviewed. The section
is concluded by discussing an exact method for solving reduced dynamics, time
evolved matrix product operators (TEMPO), [21] which confirms and extends some
of the findings beyond the Markov approximation. Section 4 presents numerical
simulations that study the heat bath’s role as a stabilizer of qubit excited
states and the critical role of the Lamb-shift. Scaling properties of FTES and
the error recovery times are studied in Sec. 5. We conclude the paper by
discussion of what can be done next in Sec. 6.

##  2  Microscopic Model

Two qubits are immersed in a common heat-bath and are described by the spin-
boson Hamiltonian

$\displaystyle H_{T}=H_{S}+A\otimes
B+\sum_{\lambda}\omega_{\lambda}b_{\lambda}^{\dagger}b_{\lambda}$ (1)
$\displaystyle B=\sum_{\lambda}g_{\lambda}(b_{\lambda}+b_{\lambda}^{\dagger})$
(2) $\displaystyle
H_{S}=\frac{\Delta}{2}[\sigma_{z}^{1}+(1-\xi)\sigma_{z}^{2}]+H_{c}$ (3)
$\displaystyle
A=\frac{1}{2}[(1+\eta)\sigma_{x}^{1}+\sigma_{x}^{2}]-\frac{\zeta}{4}(\sigma_{y}^{1}\sigma_{z}^{2}-\sigma_{z}^{1}\sigma_{y}^{2}).$
(4)

Here, $b_{\lambda}$ and $b_{\lambda}^{\dagger}$ are boson annihilation and
creation operators, respectively, acting in the Fock-space of the heat bath,
$\sigma_{x,y,z}^{1,2}$ are the Pauli matrices with superscripts indicating the
qubit number, and $\Delta$ is the drive frequency. Parameters $\xi$ and
$\eta,\zeta$ are detunings in the qubit Hamiltonian and the couplings to the
bath, respectively.

We study integrable heat baths of linear harmonic oscillators with spectral
density

$\displaystyle J(\omega)$ $\displaystyle=$
$\displaystyle\pi\sum\limits_{\lambda}g_{\lambda}^{2}\delta(\omega-\omega_{\lambda})$
(5) $\displaystyle=$
$\displaystyle\frac{\pi}{2}\alpha\omega_{c}\left(\frac{\omega}{\omega_{c}}\right)^{s}\Theta(\omega)e^{-\frac{\omega}{\omega_{c}}},$

where $\alpha$ is the dimensionless coupling constant, $\Theta(\omega)$ is the
Heaviside step function and $\omega_{c}$ is the bath cut-off frequency.
Throughout the paper, we work at zero temperature ($T=0$K) and with the bath
frequency cutoff $\omega_{c}=10\Delta$. Here, $s=1$ and $s=3$ for an Ohmic and
super-Ohmic heat-bath, respectively.

$H_{c}$ in Eq. 3 is a contribution of the heat bath to the qubit Hamiltonian
known as the counterterm, [22]

$H_{c}=\frac{1}{2}\alpha\omega_{c}\Gamma(s)A^{2}.$ (6)

The counterterm is perhaps most easily understood in the simple case in which
both the system and environment are linear harmonic oscillators with
coordinate-coordinate coupling. The counterterm in that case arises naturally
after transformation of the bath Hamiltonian to the normal mode
representation; see for example the Rubin model in Ref. [23] More generally,
whether or not the counterterm should be included from the very beginning
depends on the system considered, as discussed in Ref [22].

Regardless of whether or not the counterterm is included, there is a second
unitary contribution to the system dynamics from the bath, known as the Lamb-
shift, [24] that is distinct from the counterterm. The Lamb-shift has a linear
dependence with the high frequency cut off $\omega_{c}$ of the bath. Without
the counterterm, the Lamb-shift would cause the reduced system dynamics to
significantly depend on $\omega_{c}$ for $\omega_{c}\gg\Delta$, which is
unphysical. Adding the counterterm cancels this linear dependence, [25] making
the reduced system dynamics independent of $\omega_{c}$. This cancellation is
studied in Sec. 4.1. In contrast to the counterterm, the Lamb-shift can have
higher order contributions in the power of $\alpha$, but those contributions
do not depend on $\omega_{c}$ for $\omega_{c}\gg\Delta$. This issue is
discussed in more detail in Appendix 8.

##  3  Reduced Quantum Dynamics

All equations that describe reduced qubit dynamics are derived from the
Liouville equation for the density matrix of the total system consisting of
the qubits and the heat bath,

$\frac{d\rho_{T}}{dt}=-i[H_{T},\rho_{T}].$ (7)

We determine the approximate reduced dynamics of two qubits using four
techniques, in increasing order of accuracy with respect to exact dynamics, as
described in the following four subsections.

###  3.1 Redfield Master Equation

After rotating to the interaction picture, e.g.,
$\varrho_{T}(t)=U_{0}(t)^{\dagger}\rho_{T}(t)U_{0}(t)$, where
$U_{0}=\exp[-i(H_{S}+\sum_{\lambda}\omega_{\lambda}b_{\lambda}^{\dagger}b_{\lambda})t]$
is the free propagator, Eq. 7 becomes

$\frac{d\varrho_{T}}{dt}=-i[A(t)\otimes B(t),\varrho_{T}].$ (8)

Throughout the paper, symbols $\rho$ and $\varrho$ will represent the density
matrix in the Schrödinger and interaction picture, respectively.

This equation is equivalent to integro-differential equation

$\displaystyle\frac{d\varrho_{T}}{dt}=-i[A(t)\otimes B(t),\varrho_{T}(0)]-$
$\displaystyle\int\limits_{0}^{t}d\tau\big{[}A(t)\otimes B(t),[A(\tau)\otimes
B(\tau),\varrho_{T}(\tau)]\big{]},$ (9)

where $\varrho_{T}(0)$ is the initial condition. We assume a factorized
initial state, e.g., $\varrho_{T}(0)=\rho(0)\otimes|0\rangle\langle 0|$, where
$\rho(0)$ is the initial qubit state and $|0\rangle\langle 0|$ is the bath
vacuum. Taking the partial trace over the bath, we arrive to the familiar
expression for exact quantum dynamics for the reduced qubit state,

$\frac{d\varrho}{dt}=-\text{Tr}_{B}\int\limits_{0}^{t}d\tau\big{[}A(t)\otimes
B(t),[A(\tau)\otimes B(\tau),\varrho_{T}(\tau)]\big{]}.$ (10)

The Redfield master equation is obtained by applying the Born-Markov
approximation, which approximates the reduced density matrix of the bath with
the initial state $|0\rangle\langle 0|$. In addition the integral neglects the
memory of the qubit reduced matrix. In other words, it replaces
$\varrho(\tau)$ with $\varrho(t)$. After the approximation, Eq. 10 becomes

$\displaystyle\frac{d\varrho}{dt}$ $\displaystyle=$ $\displaystyle
A(t)\varrho(t)\int\limits_{0}^{t}d\tau C(\tau-t)A(\tau)$ (11) $\displaystyle-$
$\displaystyle A(t)\int\limits_{0}^{t}d\tau C(t-\tau)A(\tau)\varrho(t)$
$\displaystyle-$ $\displaystyle\varrho(t)\int\limits_{0}^{t}d\tau
C(\tau-t)A(\tau)A(t)$ $\displaystyle+$ $\displaystyle\int\limits_{0}^{t}d\tau
C(t-\tau)A(\tau)\varrho(t)A(t).$

Here $C(t-\tau)=\langle 0|B(t)B(\tau)|0\rangle$ is the bath correlation
function (BCF), which satisfies $C(t-\tau)=C(\tau-t)^{\star}$. For the Ohmic
and super-Ohmic baths, we obtain the BCF as the inverse Fourier transform of
the SD given by Eq. 5,

$C(t)=\int\limits_{-\infty}^{\infty}d\omega\frac{J(\omega)}{\pi}e^{-i\omega
t}=\frac{\alpha\omega_{c}^{2}\Gamma(s+1)}{2(1+i\omega_{c}t)^{s+1}}.$ (12)

Equation 11 is next represented in the eigenbasis of $H_{S}$. After rotating
back to the Schrödinger picture, e.g.,
$\rho(t)=e^{-iH_{S}t}\varrho(t)e^{iH_{S}t}$, it becomes

$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$
$\displaystyle-i[H_{S},\rho]+A\rho\Lambda(t)^{\dagger}-\rho\Lambda(t)^{\dagger}A$
(13) $\displaystyle-$ $\displaystyle A\Lambda(t)\rho+\Lambda(t)\rho A.$

This is the time dependent Redfield equation. Here,

$\Lambda(t)=\int_{0}^{t}d\tau C(t-\tau)A(\tau-t)=A\circ\Gamma(t)^{T},$ (14)

where the superscript $T$ indicates transposition and $\circ$ is the Schur, or
Hadamard, product of matrices $A$ and $\Gamma^{T}$. Thus, we get a matrix with
elements $A_{ij}\Gamma(\omega_{ji},t)$. $E_{i}$ are the eigenenergies of
$H_{S}$ and $\omega_{ij}=E_{i}-E_{j}$ are the Bohr frequencies. The function

$\Gamma(\omega,t)=\int\limits_{0}^{t}d\tau C(\tau)e^{i\omega\tau}$ (15)

will be referred to here as the timed spectral density.

For the Ohmic and super-Ohmic bath, timed spectral densities can be determined
analytically, respectively as

$\displaystyle\Gamma(\omega,t)=-i\frac{\alpha\omega_{c}}{2}\Bigg{\\{}1-\frac{e^{i\omega
t}}{1+i\omega_{c}t}-\frac{\omega}{\omega_{c}}e^{-\frac{\omega}{\omega_{c}}}$
$\displaystyle\bigg{[}ei(\frac{\omega}{\omega_{c}})-ei(\frac{\omega}{\omega_{c}}+i\omega
t)-i\pi\Theta(-\frac{\omega}{\omega_{c}})\bigg{]}\Bigg{\\}}$ (16)

and

$\displaystyle\Gamma(\omega,t)=-i\frac{\alpha\omega_{c}}{12}\bigg{\\{}\frac{\omega^{2}}{\omega_{c}^{2}}\bigg{(}1-\frac{e^{i\omega
t}}{1+i\omega_{c}t}\bigg{)}$
$\displaystyle+\frac{\omega}{\omega_{c}}\bigg{(}1-\frac{e^{i\omega
t}}{(1+i\omega_{c}t)^{2}}\bigg{)}+2\bigg{(}1-\frac{e^{i\omega
t}}{(1+i\omega_{c}t)^{3}}\bigg{)}$
$\displaystyle-\frac{\omega^{3}}{\omega_{c}^{3}}e^{-\frac{\omega}{\omega_{c}}}\bigg{[}ei(\frac{\omega}{\omega_{c}})-ei(\frac{\omega}{\omega_{c}}+i\omega
t)-i\pi\Theta(-\frac{\omega}{\omega_{c}})\bigg{]}\bigg{\\}},$ (17)

where $ei(x)=\int_{-\infty}^{x}\frac{e^{t}}{t}dt$.

####  3.1.1 Markov Limit

Master equation 13 has time-dependent coefficients. However, these
coefficients become time independent in the limit $t\gg 1/\omega_{c}$. Then,
if the asymptotic state is unique, the initial time dependence of the
coefficients will have no effect on the asymptotic state. This can be obtained
from the master equation with asymptotic coefficients or the Markov limit of
the master equation. At time $t\gg 1/\omega_{c}$, the timed spectral density
can be replaced with

$\lim_{t\to\infty}\Gamma(\omega,t)=J(\omega)+iS(\omega).$ (18)

$J(\omega)$ is given by Eq. 5 and $S(\omega)$ is a real function we call
principle density, which can be obtained either directly from Eq. 18 or by
taking the Kramers-Kronig transform of $J(\omega)$.

For Ohmic bath, we find

$S(\omega)=-\frac{\alpha\omega_{c}}{2}\left[1-\frac{\omega}{\omega_{c}}e^{-\frac{\omega}{\omega_{c}}}ei\left(\frac{\omega}{\omega_{c}}\right)\right],$
(19)

while for the super-Ohmic one,

$S(\omega)=-\frac{\alpha\omega_{c}}{2}\Big{[}2+\frac{\omega}{\omega_{c}}+(\frac{\omega}{\omega_{c}})^{2}-(\frac{\omega}{\omega_{c}})^{3}e^{-\frac{\omega}{\omega_{c}}}ei(\frac{\omega}{\omega_{c}})\Big{]}.$
(20)

In the Markov limit, the Redfield equation becomes

$\frac{d\rho}{dt}=-i[H_{S},\rho]+A\rho\Lambda^{\dagger}-\rho\Lambda^{\dagger}A-A\Lambda\rho+\Lambda\rho
A,$ (21)

where $\Lambda=\lim_{t\to\infty}\Lambda(t)$. In a more symmetric form, the
equation can be written as

$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$
$\displaystyle-i[H_{S}+H_{L},\rho]$ (22) $\displaystyle+$ $\displaystyle
A\rho\Lambda^{\dagger}+\Lambda\rho
A-\frac{1}{2}\\{A\Lambda+\Lambda^{\dagger}A,\rho\\},$

where

$H_{L}=\frac{1}{2i}(A\Lambda-\Lambda^{\dagger}A)$ (23)

is the Lamb-shift, as we will see in the next section.

Matrix elements of Eq. 22 are

$\displaystyle\frac{d\rho_{nm}}{dt}$ $\displaystyle=$
$\displaystyle-i[H_{S}+H_{L},\rho]_{nm}+\sum\limits_{ij}\Big{[}\mathcal{G}_{ni,mj}$
(24)
$\displaystyle-\frac{1}{2}\sum\limits_{k}\big{(}\delta_{in}\mathcal{G}_{km,kj}+\mathcal{G}_{ki,kn}\delta_{jm}\big{)}\Big{]}\rho_{ij},$

where

$\mathcal{G}_{ni,mj}=A_{ni}A_{mj}^{\star}\\{\Gamma(\omega_{in})+[\Gamma(\omega_{jm})]^{\star}\\}.$
(25)

In the basis $I_{ij}=|i\rangle\langle j|$ on the Banach space of linear
operators acting in the qubit Hilbert space, we rewrite Eq. 24 as

$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$
$\displaystyle\mathcal{R}\rho=-i[H_{S}+H_{L},\rho]$ (26)
$\displaystyle+\frac{1}{2}\sum\limits_{ijnm}\mathcal{G}_{ni,mj}\Big{(}[I_{ni}\rho,I_{mj}^{\dagger}]+[I_{ni},\rho
I_{mj}^{\dagger}]\Big{)}.$

Here, $\mathcal{R}$ is the generator of the one-parameter semigroup on the
Banach space.

###  3.2 Geometric Arithmetic Master Equation

The generator $\mathcal{R}$ in Eq. 26 is easily compared to the standard form
of the generator of reduced quantum dynamics, given by Eq. 38 in Ref. [26].
The latter is the most general generator of a quantum dynamical contraction
semigroup derived by Gorini, Kossakowski, and Sudarshan, [27] and
independently by Lindblad (GKSL). [28]

It immediately follows that $H_{S}+H_{L}$ in the first line of Eq. 26 maps to
the effective qubit Hamiltonian, where $H_{L}$ is the unitary contribution of
the bath to reduced qubit dynamics, or the Lamb-shift. As for the second line
in Eq. 26, it will have the GKSL form if and only if $\mathcal{G}$ is positive
semidefinite. Unfortunately, a simple exercise in linear algebra shows that
for the Redfield equation, $\mathcal{G}$ will be positive semidefinite if and
only if the spectral density is flat. Thus, it is necessary to further
approximate $\mathcal{G}$ with a positive semidefinite matrix, to obtain the
quantum dynamical semigroup.

The GKSL equation was first rigorously derived from first principles by Davies
[29] in what is now known as the rotating wave approximation (RWA). Recently,
there have been many strides in deriving GKSL equations from first principles,
with a wider range of applicability than the RWA. [30] Those include coarse-
graining, [31, 32, 18], adiabatic master equations, [33] partial RWA, [34,
35], and the quasi-phenomenological approach. [36] In the year of 2020, two
papers were published claiming first principles derivation of the GKSL
equation with accuracy equivalent to that of the Redfield ME. They are the
Universal Lindblad Equation (ULE) [20] and the Geometric-Arithmetic Master
Equation (GAME) [19], and they were independently developed. While they have
identical dissipators, they differ in the Lamb-shift.

In the derivation of GAME, [19] we made a very crude positive semidefinite
approximation in Eq. 25,

$\mathcal{G}_{ni,mj}\approx
A_{ni}A_{mj}^{\star}2\sqrt{J(\omega_{in})J(\omega_{jm})}=M_{ni}M_{mj}^{\star},$
(27)

where $M_{ij}=A_{ij}\sqrt{2J(\omega_{ji})}$. The RHS is the positive definite
Kossakowski matrix for the GKSL master equation

$\frac{d\rho}{dt}=-i[H_{S}+H_{L},\rho]+M\rho
M^{\dagger}-\frac{1}{2}\\{M^{\dagger}M,\rho\\},$ (28)

where $M$ is

$M=A\circ\sqrt{2J^{T}},$ (29)

and $\sqrt{2J^{T}}$ is defined to be the matrix with elements
$[\sqrt{2J^{T}}]_{ij}=\sqrt{2J(\omega_{ji})}$. The explanation of the square-
root approximation of spectral densities (Eq. 27) is complicated, and is
relegated to appendix 7.

GAME and ULE have a pleasant feature. If we discard the Lamb-shift at $T=0$K,
then the asymptotic state will be the ground state of the qubit Hamiltonian
$H_{S}$, which is natural. Thus, if the asymptotic state is found to be
different from the ground state, this will be due to the unitary effect of the
environment expressed through the Lamb-shift. In this context, the Lamb-shift
becomes critical for the existence of the FTES. (By contrast, at $T>0$K, the
asymptotic state is no longer the Gibbs state, [37] and the above statement no
longer holds.)

###  3.3 Time-Convolutionless (TCL) Master Equation.

The TCL-master equation is an exact master equation obtained using the
Nakajima-Zwanzig projection-operator technique. [24] The main assumption is
that the time-evolution of the reduced density matrix of the system is
invertible. In that case, it is possible to eliminate the integral over the
history in the Nakajima-Zwanzig equation. It may be debatable if the
assumption is applicable for the evaluation of the asymptotic state. Namely,
since the asymptotic state results from an arbitrary initial state, it is
clearly not invertible. The asymptotic state is still worth exploring, since
the TCL is more accurate than the Redfield equation, which suffers from the
same problem.

The generator of the TCL-equation can be expressed as perturbation series in
$\alpha$. [24] In the lowest (linear) order of $\alpha$, the TCL equation is
equivalent to Eq. 13. In the literature, the names RE and TCL2 are used
interchangeably, but here we only use RE. In the second order of $\alpha$, TCL
is named TCL4. Implementation of TCL4 is complicated, and the discussion is
given in appendix 8.

The TCL4 master equation reads as Eq. 13 with an extra term,

$\displaystyle\frac{d\rho_{nm}}{dt}$ $\displaystyle=$
$\displaystyle\\{-i[H_{S},\rho]+A\rho\Lambda(t)^{\dagger}-\rho\Lambda(t)^{\dagger}A$
(30) $\displaystyle-$ $\displaystyle A\Lambda(t)\rho+\Lambda(t)\rho A\\}_{nm}$
$\displaystyle+$
$\displaystyle\sum\limits_{i,j}\delta\mathcal{G}_{ni,mj}(t)\rho_{ij}.$ (31)

Tensor $\delta\mathcal{G}_{ni,mj}(t)$ is quadratic with $\alpha$ and time
dependent. It is given by triple integrals of various four-point bath
correlation functions, which can be found in the appendix 8.

Since our focus is to study the asymptotic state, we consider Eq. 31 in the
Markov limit $t\to\infty$. We find $\delta\mathcal{G}_{ni,mj}(t)$ converges
well for both Ohmic and Super-Ohmic SDs. Calculating the tensor
$\delta\mathcal{G}_{ni,mj}(\infty)$ is the hard part, after which the
numerical calculation of system dynamics and the asymptotic state can be
determined as easily as with the RE and GAME master equations.

###  3.4 Time Evolving Matrix Product Operators (TEMPO).

The Born-Markov approximation used in the previous sections is usually valid
if the coupling between the system and the heat bath is weak. If the coupling
increases, at some point the correlations in the heat bath imparted by the
system dynamics persist long enough to feed back into the system, altering the
system dynamics. This feedback is nonlocal in time and can dramatically change
the asymptotic state. In an extreme example, one qubit coupled to Ohmic bath
displays the spin-boson quantum phase transition at the critical coupling to
the bath. [38, 23, 39] Even if the coupling is far below critical, the
entanglement of the asymptotic state of the system, calculated using exact
quantum dynamics, can differ significantly from that calculated by Markov
master equations. [25]

Since the phenomenon of FTES involves feedback from the bath, there is an
immediate concern if the phenomenon is well represented by the Markov master
equations. Thus, we will investigate if the phenomenon persists in non-
Markovian dynamics. In the path-integral formulation of quantum mechanics, the
effect of the environment of linear harmonic oscillators can be exactly
accounted for by the Feynman-Vernon influence functional [40]. Its
implementation includes the exact quasi adiabatic path integral method
(QUAPI), [41, 42, 43, 44] where the central quantity is the augmented density
tensor (ADT). This carries the probability distribution of possible histories
and auto-correlations of the system. The density matrix is obtained by tracing
over those histories. Another path integral method commonly used is Monte
Carlo [45]. More recently, tensor network methods have been utilized [46, 47,
21, 48, 49, 50, 51, 52] to simplify and extend the QUAPI method.

Other methods of exact dynamics include the hierarchical equations of motion,
[53, 54, 55, 56, 57], the multiconfiguration time-dependent Hartree (MCTDH)
algorithm [58, 59], and the multilayer formulation of the previous algorithm.
[60, 61] The Nakajima-Zwanzig equation is also exact, [62, 63] but it is as
difficult to solve as the Liouville equation. [24] The efficiency of the
method has recently been improved using quantum trajectory based hierarchy of
stochastic pure states, [64, 65, 66] including investigation of entanglement
between two spins mediated by the bath. [25]

In this paper, we apply the time evolving matrix product operators (TEMPO)
[21] technique. The main idea is to compress the ADT in the form of a matrix
product state [67, 68] using singular-value-decomposition (SVD), which can
significantly extend the history time of the approximate path integral. The
introduction to the method and the python code we use in most of this paper
are available in Ref. [21]. The TEMPO algorithm can also be modified to
calculate the process tensor [69]. This alternative formulation has been used
in both optimizing quantum control procedures [70] and in calculating exact
bath dynamics. [71] An open source python package [72] improves the
performance and includes new approaches.

##  4 Heat-Bath as State Stabilizer

It is theoretically established that two qubits interacting solely via a
common heat bath can exhibit permanent entanglement, depending on the system-
bath coupling and bath parameters. [73, 74, 75, 76, 77, 25] The entanglement
usually results from the unitary transformation that the heat bath applies on
the qubit system. [31, 78, 32, 79, 80, 19] However, asymptotic entanglement of
qubits is typically weak, as dissipative dynamics drives the system close to
the ground state, which is not easily entangled.

Our claim is that the unitary effect of the bath can, in addition, create
fault tolerant excitations. In the context of Eq. 28, an excited state
$|FT\rangle$ of the Hamiltonian $H_{S}+H_{L}$ will be fault tolerant if two
conditions are met: 1) the state is annihilated by the generator, e.g.,
$M|FT\rangle=0$; 2) the asymptotic state of the master equation
$\rho_{\infty}$ must be unique.

Under the first condition, $|FT\rangle\langle FT|$ is a solution of Eq. 28,
while under the second condition, the state can recover from errors. Let us
assume that an error corrupts $|FT\rangle\langle FT|$, which can be due to an
environment not accounted for by the model. Since $|FT\rangle\langle FT|$ is
both asymptotic and unique, the corrupted state will spontaneously recover
under Eq. 28, assuming no additional error will take place during the recovery
time.

In general, the asymptotic state need not be unique. For example, noiseless
quantum codes [81] or decoherence free subspaces (DFS) [82] have degenerate
asymptotic states. A state in a DFS is annihilated by the generator. However,
the purity of the asymptotic state is not guaranteed, because a convex sum of
different asymptotic states is a mixed asymptotic state. Noiseless quantum
codes and DFS occur if there is a symmetry of the Hamiltonian with respect to
permutations between qubits. All such symmetries are broken in our model,
guaranteeing the uniqueness of the asymptotic state.

###  4.1 Bath Induced Renormalization

The unitary effect of the heat-bath transforms the two lowest excited states
of two weakly detuned qubits into singlet-like and triplet-like
superpositions. The singlet-like state exhibits suppressed damping because of
the destructive interference caused by the singlet superposition. The
interference is at the heart of suppressing damping to create a stable excited
state. Thus, we begin by finding out how the bath transforms quantum states of
the qubit system.

Figure 2: Renormalizations of the qubit energy levels by the Ohmic heat bath.
a) Counterterm only. b) Lamb-shift only. c) Both the counterterm and the Lamb-
shift. Full lines show 2-qubit system energy levels versus $\alpha$, on the
left axes. Dashed-lines display singlet fidelity of the appropriate (e.g.,
singlet-like) excited state versus $\alpha$, on the right axis. Insets: zoom-
in at avoided level-crossings. $\Delta=1$, $\xi=0.05$, $\eta=0.179$,
$\zeta=0$, $\omega_{c}=10$, $s=1$.

In our recent work, [19] we studied bath induced renormalization, or avoided
crossing, between two weakly detuned excited levels in a three-level Jaynes-
Cummings model. Excluding the counterterm, the avoided crossing energy is
$2S(\Delta)$, where $\Delta$ is the drive frequency of the model. The pair of
weakly detuned excited levels is transformed into a singlet-like state
$|S\rangle^{\prime}\approx|S\rangle=(|1,0\rangle-|0,1\rangle)/\sqrt{2}$ at
higher energy, and a triplet like state
$|T0\rangle^{\prime}\approx|T0\rangle=(|1,0\rangle+|0,1\rangle)/\sqrt{2}$ at
lower energy. The avoided crossing energy can be quite large for large
$\omega_{c}$, suggesting that the counterterms should be considered.

Here we investigate two qubit system Hamiltonian renormalizations. If we take
the limit $\omega/\omega_{c}\to 0$ in Eqs. 19-20, then $A_{f}\approx-iS(0)A$,
and the Lamb-shift in Eq. 23 becomes
$S(0)A^{2}=-\alpha\omega_{c}\Gamma(s)A^{2}/2$, equal and opposite to the
counterterm given by Eq. 6. At finite $\omega/\omega_{c}$, there is incomplete
cancellation of the terms, the effects of which have recently been studied in
Ref. [25]

Figure 3: Renormalizations of the qubit energy levels by an Ohmic and a super-
Ohmic heat bath. a) $\Delta=1$, $\xi=0.05$, $\eta=0.179$, $\zeta=0$, $s=1$. b)
$\Delta=1$, $\xi=0.05$, $\eta=0.37$, $\zeta=0$, $s=3$. c) Singlet fidelity
versus scaled damping $\alpha/\xi$, at four values of $\xi$. $s=1$,
$\eta=0.179$, $\zeta=0$, $s=1$. In all figures, $\omega_{c}=10$.

Note that in our model, the renormalized system Hamiltonian is always block-
diagonal. In the computational basis, it is

$\left(\begin{array}[]{cccc}a&0&0&b\\\ 0&c&d&0\\\ 0&d^{\star}&e&0\\\
b^{\star}&0&0&f\\\ \end{array}\right).$ (32)

We focus first on the subspace spanned by the basis $|01\rangle$ and
$|10\rangle$. Without any detuning, the Hamiltonian in the subspace is
explicitly given as

$H_{S}+H_{L}=\bigg{[}q\Delta+\frac{S(\mathcal{E}_{-})\mathcal{E}_{+}-S(\mathcal{E}_{+})\mathcal{E}_{-}}{\Delta\sqrt{q^{2}+4}}\bigg{]}|T0\rangle\langle
T0|$ (33)

where $q=\Gamma(s)\alpha\omega_{c}/2\Delta$ and

$\mathcal{E}_{\pm}=\frac{\Delta}{2}\left(q\pm\sqrt{q^{2}+4}\right)$ (34)

are the extremal eigenenergies of $H_{c}$. In the singlet-triplet basis, the
renormalized Hamiltonian is diagonal, with the singlet-triplet splitting
depending on the delicate balance between the counterterm and Lamb-shift.

As an example, for $\alpha=0.2$ and $\omega_{c}=10\Delta$, we find a singlet-
triplet splitting of $-0.067\Delta$ [$0.0108\Delta$] for the Ohmic [super-
Ohmic] bath. In comparison, the norm of the counterterm for the Ohmic bath is
half that of the super-Ohmic. The larger counterterm with a smaller level
repulsion implies that the super-Ohmic bath has a much better cancellation of
the counterterm and the Lamb-shift, which is consistent with the cancellation
trends studied recently in Ref. [25]

The regime at finite system and system-bath detunings is numerically examined.
We determine the eigenvalues and eigenstates of Hamiltonians $H_{S}$,
$H_{S}-H_{c}+H_{L}$, and $H_{S}+H_{L}$. In Fig. 2, we plot the renormalized
qubit levels versus damping $\alpha$. In a) and b), the renormalizations are
only due to the counterterm and Lamb-shift, respectively. The levels shift
asymmetrically with increasing $\alpha$, with approximately equal and opposite
centers of gravity. The insets on the bottom left show repulsion between
weakly detuned-levels with an avoided crossing at $\alpha\approx 0.005$. At
larger $\alpha$, the singlet-like level does not shift with $\alpha$. The
singlet fidelity, indicated by the dashed line (right axis), ranges from one-
half at small $\alpha$ to one at large $\alpha$.

In comparison, Fig. 2 c) displays the qubit levels renormalized by both the
counterterm and the Lamb-shift as a function of $\alpha$. The center of
gravity of the spectrum remains close to zero, showing an overall cancellation
of the two terms. Despite this, some weak shifts remain.

Figs. 3 a) and b) display avoided crossings between weakly detuned levels
versus damping for the Ohmic and super-Ohmic bath, respectively, with both
terms included. The minimum spacing between levels 2 and 3 is located at
$\alpha=0.0673$ and $0.1908$ for a) and b), respectively. The avoided crossing
energy is approximately three times larger for the Ohmic bath, even though in
this case the super-Ohmic bath has a 6.7 times larger counterterm norm. This
again shows that the counterterm and the Lamb-shift cancel more effectively
with large $s$. Jumping ahead, stronger level repulsion for the Ohmic bath,
relative to the super-Ohmic bath, will make it easier to prepare FTES and
gives an advantage to slower baths.

If we plot the singlet fidelity of the singlet-like excited state versus
scaled parameter $\alpha/\xi$, then the curves will coalesce at very small
$\alpha$, as shown in Fig. 3 c). Thus, in the weak coupling regime, singlet
fidelity versus $\alpha$ and $\xi$ is a function of only one variable,
$\xi/\alpha$. We will encounter this scaling again and again.

###  4.2 Suppression of Singlet Damping by System-Bath Detuning

Throughout the paper, the initial condition will be the singlet state. Master
equations are solved by decomposing the initial vectorized density matrix into
a linear combination of the eigenvectors of the superoperator, after which the
solutions can be found by exponentiation in a fraction of a second. TEMPO
calculations are much slower and can take anywhere between several days to a
month per curve, as detailed in appendix 9.

Figure 4 [5] displays singlet fidelity $\langle S|\rho|S\rangle$ versus time,
with [without] the counterterm, at fixed qubit detuning $\xi$ and varied
qubit-bath detuning $\eta$. Overall, different methods lead to qualitatively
similar results. They all exhibit a strong suppression of singlet damping as a
function of a bath detuning $\eta$, which is part of the main result of the
paper.

Singlet fidelity decays incoherently in Fig. 4. As $\eta$ increases in the
range $[0,0.35]$, the renormalized energy spacing between singlet and triplet-
like states changes, corresponding to the oscillation period in the interval
$[60/\Delta,78/\Delta]$. The period is too large compared to the damping and
decoherence time, explaining the incoherent singlet decay. In the absence of
the counterterm, the splitting is much larger [compare the spacing between
levels 2 and 3 in Fig. 2 b) versus 2 c)], leading to a much shorter
oscillation period, in quantitative agreement with the oscillations in Fig. 5.

Figure 4: Singlet fidelity versus time for Ohmic heat bath including the
counterterm effect. (a-d) Numerical results obtained, respectively, using
GAME, RE, TCL4, and TEMPO. Full red lines, bottom-to-top:
$\eta=0,0.1,0.15,0.179$, respectively. Dashed black lines, top-to-bottom:
$\eta=0.19,0.22,0.25,0.3,0.35$, respectively. Other parameters: $\Delta=1$,
$\xi=0.05$, $\zeta=0.0437$, $\alpha=0.2$, $\omega_{c}=10$, and
$\rho(0)=|S\rangle\langle S|$. TEMPO parameters: time step $\tau_{s}=0.02$,
number of steps stored in the ADT $K_{max}=150$, $\epsilon_{SVD}=10^{-6}$.

Depending on the qubit-bath detuning $\eta$, Figs. 4 b-c) and 5 b-c) can
display negative rates. That is, the singlet probability can increase with
time and exceed one. We checked that this unphysical behavior is due to
density matrix positivity violation, which is not that uncommon in the case of
the RE and TCL4 master equations. The asymptotic state is not bound in the
range of $\eta$ with negative rates. In comparison, the TEMPO curves in Fig. 4
have significant negative curvature when the singlet fidelity increases versus
time, suggesting a bound asymptotic state. We also verified density matrix
positivity violations in our TEMPO calculations.

TEMPO represents exact quantum dynamics, [70] which must be completely
positive. We lack the computational resources needed to maintain positivity of
the density matrix in TEMPO. We vary the simulation parameters, such as memory
cutoff time, precision of the SVD compression, the time step, and the boundary
condition at memory cut-off, and find that positivity violations are very
stubborn within reasonable calculation time-frames.

Figure 5: The same as Fig. 4, but without the counterterm. (a-d) Full red
lines, bottom-to-top: $\eta=0.02,0,-0.02,-0.04,-0.06,-0.075,-0.09$,
respectively. Dashed black lines, top-to-bottom:
$\eta=-0.13,-0.16,-0.18,-0.2,-0.25,-0.3$, respectively. Other parameters:
$\Delta=1$, $\xi=0.05$, $\zeta=0$, $\alpha=0.2$, $\omega_{c}=10$, and
$\rho(0)=|S\rangle\langle S|$. TEMPO parameters: time step $\tau_{s}=0.02$,
number of steps stored in the ADT $K_{max}=150$, compression precision
$\epsilon_{SVD}=10^{-6}$.

One likely explanation of this problem is that the Ohmic bath has a slow decay
of bath correlations. In fact, we find that the problem goes away if we
consider the super-Ohmic bath. With a smaller memory cut-off for the ADT, we
still get a positive semidefinite density matrix, thereby making the
simulation more reliable. We can even explore singlet decay in the strong
coupling regime, where master equations fail and the dynamics is not
Markovian, e.g., as follows.

Fig. 6 a-d) displays singlet fidelity versus scaled time $\alpha t$ for the
super-Ohmic bath. Between a) and d), $\alpha$ increases by four orders of
magnitude, while $\xi/\alpha=0.2$, and $\zeta$, $\Delta$, and $\omega_{c}$ are
constant. In each sub-figure, $\eta$ changes within the same set of values.

The curves in Fig. 6 a-c) are obtained by solving the Redfield master
equation. For $\alpha=0.001$, Fig. 6 a), we also find the solutions of the
GAME master equation, and obtain curves indistinguishable from those in the
figure. Since the RE and GAME master equations have equivalent accuracy that
scales linearly with the coupling parameter $\alpha$, [19] this
indistinguishability confirms correctness of both master equations. In
addition, we determine the solutions of the RE master equation at
$\alpha=10^{-4}$, and also find curves indistinguishable from those in a).
This implies scaling of qubit dynamics in the weak coupling limit, similar to
the invariance in the singlet fidelity in Fig. 3 c).

Figure 6: Singlet decay for the super-Ohmic bath. (a-c) show the RE results,
for $\alpha=10^{-3},10^{-2},10^{-1}$, respectively. d) shows the exact
solution by TEMPO. Full red lines, bottom-to-top:
$\eta=0,0.20,0.28,0.33,0.36$, respectively. Dashed black lines, top-to-bottom:
$\eta=0.37,0.38,0.39,0.42$, respectively. Other parameters: $\alpha=1$.
$\Delta=1$, $\xi/\alpha=0.2$, $\zeta=0$, $\omega_{c}=10$, and
$\rho(0)=|S\rangle\langle S|$. TEMPO: $\tau_{s}=0.02$, $K_{max}=75$,
$\epsilon_{SVD}=10^{-6}$.

Comparing Figs. 6 a) and b), we see that the curves are slightly different.
They no longer scale, which is attributed to the non-negligible coupling to
the bath in b). With another order of magnitude increase in damping, Fig. 6 c)
displays unphysical singlet fidelity, larger than one or less than zero. In
this regime, the weak-coupling approximation to reduced quantum dynamics is
failing. With further increase in orders of magnitude of $\alpha$ not shown,
the RE master equation is very unstable.

On the other hand, the TEMPO simulation is stable at $\alpha=1$ and exhibits a
positive definite density matrix and singlet damping as shown in Fig. 6 d).
The damping curves in the weak and strong coupling regimes are overall very
similar. For example, the singlet damping rate is strongly suppressed near the
same value of $\eta$ in a) and d) ($\eta\approx 0.37$). At $\eta=0$, the
oscillation periods measured in the damping curves are $6.11/\xi$ and
$8.08/\xi$ in a) and d), respectively. We compare those periods with the
detuning of the renormalized qubit Hamiltonian, and obtain the respective
numbers of $6.12/\xi$ and $7.32/\xi$. The discrepancy between the frequencies
in the strong coupling regime will be the subject of future work. Regardless,
the main conclusion of this section is that the strong suppression versus bath
detuning in singlet damping remains valid in the strong coupling regime.

We conclude that the suppression of the singlet damping rate is overall weakly
dependent on the coupling strength $\alpha$ if $\xi/\alpha$ is constant. It
behooves that the suppression of damping of the singlet-like excited state
persists in the strong-coupling regime. Phenomenological theory of quantum
error correction is based on the operator-sum representation of quantum
operations, [1] which is correct regardless of the strength of the coupling.

###  4.3 Cancelation of ”Singlet” Damping

Consistency of the solutions between the four different methods is
encouraging. Out of these methods, only the geometric-arithmetic master
equation has an easily calculable and positive semidefinite asymptotic density
matrix. We thus focus on the GAME method.

The asymptotic state is determined as the eigenstate of the superoperator with
eigenvalue zero. Fig. 7 displays the singlet fidelity and purity of the
asymptotic state as a function of bath detuning $\eta$ at qubit detuning
$\xi=0.001$ and damping $\alpha=0.004$. Dashed lines are calculated using
$\zeta=0$. The black curves with peaks at negative $\eta$ correspond to the
system Hamiltonian without the counterterm, and vice versa for the red and
blue curves with peaks at positive $\eta$. The insets zoom-in the peaks with
the counterterm.

Figure 7: a) and b): Singlet fidelity $\langle S|\rho_{\infty}|S\rangle$ and
purity $Tr\,\rho_{\infty}^{2}$ of the asymptotic state, respectively, for
$\alpha=0.004$. Full and barely visible dashed black line: $\zeta=7.5\times
10^{-5}$ and $\zeta=0$, respectively, with no counterterm. Full red and dashed
blue line: $\zeta=1.368\times 10^{-3}$ and $\zeta=0$, respectively, with
counterterm included. Insets: Zoom-in peak fidelity and purity, with
counterterm included. $\xi=0.001$, $s=1$, $\omega_{c}=10$.

One observes from the insets in Fig. 7 that if $\zeta=0$, the maximum singlet
fidelity and state purity are only $0.66$ and $0.55$, respectively. So, we
optimize both $\eta$ and $\zeta$ to maximize the asymptotic state purity. This
can be done efficiently using standard optimization methods. The full black
($\zeta=7.5\times 10^{-5}$) and red ($\zeta=-1.36\times 10^{-3}$) lines show
singlet fidelity and state purity versus $\eta$ at the optimum $\zeta$. The
optimization significantly increases the maximum singlet fidelity and purity
if the counterterm is included. The asymptotic state purity reaches one,
within the precision of at least $10^{-8}$, limited by the numerical accuracy
of the optimizer. We thus consider the asymptotic state to be pure. We verify
that the asymptotic state is an excited eigenstate of the qubit Hamiltonian.

Figure 8: a) and b): Singlet fidelity $\langle S|\rho_{\infty}|S\rangle$ and
purity $Tr\,\rho_{\infty}^{2}$ of the asymptotic state, respectively, for
$\alpha=0.2$. Full and barely visible dashed black line: $\zeta=3.6\times
10^{-3}$ and $\zeta=0$, respectively, with no counterterm. Full red and dashed
blue line: $\zeta=0.0437$ and $\zeta=0$, respectively, with counterterm
included. $\xi=0.05$, $s=1$, $\omega_{c}=10$.

For a weak system-bath coupling, the FTES occurs within a narrow range of
system-environment parameters. This range increases with increasing damping
and detuning, as shown by the broader peaks in Fig. 8 a) and b). At the
maximum, the purity is again larger than $1-10^{-8}$, and the asymptotic state
is the singlet-like excited eigenstate of the qubit Hamiltonian.

###  4.4 FTES Recovery

In QEC, the faster an unknown quantum state recovers from an error transient,
the higher the state fidelity. Similarly, one figure of merit of FTES is the
recovery rate, which can be obtained from the smallest magnitude nonzero
eigenvalue of the superoperator of the quantum dynamical semigroup.

Figure 9 a) displays the smallest relaxation rate ($\Gamma_{r}$) versus bath
detuning $\eta$ at $\alpha=0.004$. Looking at the larger range of $\eta$,
there is virtually no difference between GAME and RE curves. However, the
inset zooms in at the minimum and shows differences. Critically, the
relaxation rates of RE and TCL4 become negative, and the asymptotic state
diverges in the neighborhood of the FTES identified by GAME. In the latter
approximation, $\Gamma_{r}$ is positive, guaranteed by the contracting
semigroup. [26] Since RE and GAME have comparable accuracy with respect to
exact quantum dynamics, we take the recovery rate of the FTES to be only an
estimate, with an error comparable to the estimate.

In the case of stronger damping, Figure 9 b) includes the recovery rate
obtained from TEMPO simulations, where $\Gamma_{r}$ is obtained by fitting
singlet fidelity versus time to the exponential, as discussed in appendix 9.
The magnitudes of the minima are comparable between the four methods,
confirming the validity of the GAME recovery rate within the factor of order
one. The ability to estimate the recovery rate easily makes GAME the preferred
method to search for and study FTES.

Figure 9: FTES recovery rate $\Gamma_{r}$. a) $\alpha=0.004$. Full-red line:
GAME-values. Thick dashed black line: RE values. Inset: Zoom in at minimum
rate. Blue circles and line: TCL4 values and best parabolic fit.
$\zeta=1.36\times 10^{-3}$, $\xi=0.001$. b) $\alpha=0.2$. Full-red line: GAME-
values. Thick dashed black line: RE values. Blue circles and line: TCL4 values
and best parabolic fit. Pink squares and line: TEMPO-values and best parabolic
fit. $\zeta=0.0437$, $\xi=0.05$. All $\Gamma_{r}$ are in units of drive
frequency $\Delta$. $s=1$, $\omega_{c}=10$, counterterm included.

###  4.5 Wrapping it up

We proceed to qualitatively explain the dynamics of FTES. The sketch in Fig.
10 a) displays four energy levels of two weakly detuned qubits, without taking
into account any unitary effect induced by the heat bath. We assume that the
Fermi-golden rule relaxation rates (equal to $2J$ and indicated by the blue
arrows) are much larger than the detuning $\xi\Delta$ between states
$|0,1\rangle$ and $|1,0\rangle$, but much smaller than $\Delta$. In terms of
the GKSL-equation 28, all the information about these rates is in the
generator $M$.

In the presence of the unitary contribution from the bath, levels
$|0,1\rangle$ and $|1,0\rangle$ are strongly admixed, since their relaxation
rate was assumed to be much larger than the detuning. On the other hand, the
matrix 32 admixes levels $|0,0\rangle$ and $|1,1\rangle$ weakly, because their
spacing is much larger than the relaxation rates.

Figure 10: Dynamics of the FTES $|S\rangle^{\prime}$. a) Energy levels of a
two qubit system in the absence of a unitary transformation by the heat-bath.
The arrows represent the relaxation process due to a dissipative coupling to
the bath at zero temperature. The rates are proportional to the thickness of
the arrows and given by the Fermi-golden rule, where only energy-decreasing
transitions are allowed. b) Energy levels after the unitary transformation.
The top-most and the lowest levels are weakly renormalized. The weakly detuned
levels, $|0,1\rangle$ and $|1,0\rangle$, are strongly renormalized into
singlet and triplet like states $|S\rangle^{\prime}$ and
$|T0\rangle^{\prime}$. The triplet relaxation rate is enhanced, while that of
the singlet is suppressed. Energy increasing processes are weakly allowed, due
to the admixing.

We consider the GKSL equation in the basis of the renormalized states, which
are the eigenstates of $H_{S}+H_{L}$. The effect of the unitary transformation
on the generator is

$M=\left(\begin{array}[]{cccc}0&0&0&0\\\ \bullet&0&0&0\\\ \bullet&0&0&0\\\
0&\bullet&\bullet&0\\\
\end{array}\right)\to\left(\begin{array}[]{cccc}0&\bullet&\bullet&0\\\
\bullet&0&0&\bullet\\\ \bullet&0&0&\bullet\\\ 0&\bullet&\bullet&0\\\
\end{array}\right),$ (35)

while that on the damping term is

$M^{\dagger}M=\left(\begin{array}[]{cccc}\bullet&0&0&0\\\
0&\bullet&\bullet&0\\\ 0&\bullet&\bullet&0\\\ 0&0&0&0\\\
\end{array}\right)\to\left(\begin{array}[]{cccc}\bullet&0&0&\bullet\\\
0&\bullet&\bullet&0\\\ 0&\bullet&\bullet&0\\\ \bullet&0&0&\bullet\\\
\end{array}\right),$ (36)

where full circles indicate nonzero matrix elements. The arrows in Fig. 10 b)
sketch the relaxation rates in the renormalized basis: the triplet-like
relaxation rate is approximately doubled, while the singlet-like one is
strongly suppressed relative to the Fermi-golden rule, due to constructive and
destructive interference, respectively.

Additionally, the rates of the reverse processes (yellow arrows) are nonzero
(but weak). This can be seen, for example, by the nonzero matrix element
$(4,4)$ on the RHS of Eq. 36, which indicates a nonzero loss rate of the
renormalized ground state $|1,1\rangle^{\prime}$.

If we tweak the system-bath interaction Hamiltonian through parameters
($\eta,\zeta$) we can make the singlet relaxation rate be the bottleneck rate.
As a result, the asymptotic state of the qubits will not be far from the
singlet-like excitation. In the regime of fault tolerance, the bottleneck rate
is precisely zero, and the asymptotic state is the pure ”singlet” excitation.
In that regime, the recovery rate is given by a nontrivial combination of
various rates indicated in picture 10 b), which we obtain by finding the
smallest magnitude nonzero eigenvalue of the superoperator of the GAME master
equation.

##  5 Scaling Properties of the FTES

Here, we vary the system-environment parameters $\xi,\eta,\zeta,\alpha$ and
$s$, and identify the parameter range that exhibits FTES, e.g., where the
asymptotic state $\rho_{\infty}$ has a fidelity of one to an excited
eigenstate of the renormalized qubit Hamiltonian. We also explore the
parameter range where the FTES has high singlet-fidelity. All results in this
section are obtained by finding the asymptotic solutions of the GAME master
equation and maximizing its purity versus system parameters, so that
$\text{Tr}\rho_{\infty}^{2}>1-10^{-8}$.

###  5.1 The effects of qubit detuning on FTES.

Here, we vary the detuning $\xi$ at fixed values of $\alpha$, $s$, and
$\omega_{c}$. For each $\xi$, we optimize $\eta$ and $\zeta$ to find the FTES,
which leads to the unique value of $\eta$ and $\zeta$. The optimum $\eta$ and
$\zeta$ versus $\xi$ are displayed in Figs. 11 a) and b) for the Ohmic and
super-Ohmic baths, respectively. Note that here $\eta$ can be larger than $1$,
which means there’s strong asymmetry between the couplings of the qubits to
the environment, where the damping of one qubit is enhanced by $(1+\eta)$.

Let us first examine the regime of weak detuning, which we define as the range
where $\xi<0.2$, $\eta<0.2$, and $\zeta<0.2$. At small $\xi$, $\eta$ and
$\zeta$ increase linearly with $\xi$, as shown in Figs. 11 a) and b). The
slope $\eta$ versus $\xi$ is proportional to $1/\alpha$. The property
$\eta\approx\text{const}\,\xi/\alpha$ is consistent with the analysis of the
suppression of singlet damping in Fig. 6.

Figure 11: System-bath coupling parameters $\eta$ and $\zeta$ versus system
detuning $\xi$, at the FTES, for Ohmic a) and super-Ohmic b) heat baths. The
top three (black) curves show $\eta$. The bottom three (red) curves show
$\zeta$. The dashes correspond to different dampings $\alpha$, as displayed in
the legend. $\omega_{c}=10$. Figure 12: Maximum and minimum relaxation rate
versus system detuning for Ohmic (a) and super-Ohmic (b) heat baths. The top
three (black) curves show the maximum relaxation rate. The bottom three (red)
curves show the recovery (e.g., minimum relaxation) rate. The dashes
correspond to different dampings $\alpha$, as displayed in the legends.
$\omega_{c}=10$ and all frequencies are in units of the drive frequency
$\Delta$.

The recovery rate and the maximum qubit relaxation rate are obtained by
diagonalizing the superoperator, finding the eigenvalues with the largest and
smallest real part, taking the real part, and multiplying by $-1$. The rates
are displayed in Fig. 12. The recovery rate for the Ohmic bath scales as
$\Gamma_{r}\simeq\alpha\xi^{2}$ over the entire range in Fig. 12 a). By
contrast, $\Gamma_{max}\simeq\alpha$ and is independent of $\xi$ at small
$\xi$, with a crossover to a more complicated behavior at high $\xi$. The
complicated behavior will be discussed in the next section. The quadratic
scaling of $\Gamma_{r}$ with $\xi$ is the general property of decoherence-free
proximity spaces. [82] Namely, if the permutational symmetry of the qubits is
broken by a weak inhomogeneity (e.g., the detuning), the relaxation rate out
of the decoherence free subspace scales quadratically with the inhomogeneity.

Note that the master equation is likely valid in most of the parameter range
in Fig. 12 since the maximum relaxation rate is at worst somewhat smaller than
the drive $\Delta$. The highest recovery rate is found at the highest damping.
For the Ohmic bath, $\Gamma_{r}$ is in the $10^{-2}\Delta$ range, which is a
rapid error correction rate. Namely, if there is an extrinsic noise source
that drives the qubit-state out of its asymptotic state at the rate
$\Gamma_{error}$, then the state fidelity will be approximately
$1-\exp(-\Gamma_{r}/\Gamma_{error})$. For example, if $\Gamma_{r}=0.01\Delta$
and $\Gamma_{error}=10^{-3}\Delta$, then the FTES fidelity will be $>0.9999$.

###  5.2 The effects of damping on FTES.

Figure 13: System-bath coupling parameters $\eta$ and $\zeta$ at fault
tolerance versus damping $\alpha$ at fixed $\xi$, for Ohmic (a) and super-
Ohmic (b) heat baths. The top three (black-dashed) curves show $\eta$. The
bottom three (red-full) curves show $\zeta$. $\omega_{c}=10$. Figure 14:
Maximum and minimum relaxation rate versus damping, for Ohmic (a) and super-
Ohmic (b) heat baths. The top three (black-dashed) curves show maximum
relaxation rate. The bottom three (red-full) curves show the recovery (e.g.,
the minimum relaxation) rate. $\omega_{c}=10$ and all frequencies in unit of
drive frequency $\Delta$.

We obtain further insights by following the fault-tolerant excited state as a
function of damping at a fixed system detuning $\xi$. The bath coupling
parameters $\eta$ and $\zeta$ are obtained as explained in the previous
section, and follow the curves in Fig. 13. The bath detuning $\eta$ exceeds
one at small $\alpha$. As discussed in the previous section, FTES follows the
constraint $\alpha\eta\simeq\xi$. So, for fixed $\xi$, $\eta$ increases as
$1/\alpha$ and can become larger than 1.

The recovery and the maximum relaxation rate versus $\alpha$ are shown in Fig.
14. As in the previous section, $\Gamma_{r}\simeq\alpha\xi^{2}$, while
$\Gamma_{max}$ has a crossover from $\simeq\alpha$ at large $\alpha$ to
$\simeq$const at small $\alpha$. The flattening of $\Gamma_{max}$ versus
$\alpha$ at small $\alpha$ is also related to the constraint
$\eta\simeq\xi/\alpha$. If $\alpha$ is small, $\eta$ can be larger than one,
and the damping rate of one qubit will be enhanced to
$(1+\eta)\alpha\approx\eta\alpha\propto\xi$. This is constant in this figure,
thereby explaining why the rate $\Gamma_{max}$ becomes flat at small $\alpha$.
In comparison, the recovery rate is not affected by large $\eta$. The
following table summarizes the various scalings of the FTES.

$\eta<1$ $\eta>1$ $\eta$ $\xi/\alpha$ $\xi/\alpha$ $\Gamma_{max}$ $\alpha$
const $\Gamma_{r}$ $\alpha\xi^{2}$ $\alpha\xi^{2}$

###  5.3 Singlet fidelity of the FTES.

High singlet fidelity would be important if the goal were to prepare and
preserve singlet states as resource for quantum information processing. Fig.
15 displays singlet fidelity of the FTES versus system detuning. The arrows in
the figure indicate the points where the bath-coupling detuning $\eta$ is one.
That is, $\eta>1$ to the right of the arrows. We conclude that the condition
for significant singlet fidelity of the FTES is that the detuning of the
qubit-coupling to the heat bath is significantly less than 1, e.g., that the
couplings are close to symmetric. In that regime, the microscopic details of
the Hamiltonian and spectral density have only a weak effect on the FTES.

Figure 15: Singlet fidelity of the FTES versus system detuning $\xi$. The
higher the damping, the higher the characteristic system detuning, above which
the FTES becomes non-universal and not singlet like.

##  6 Discussion

In this paper we have presented a theoretical argument for the fault-tolerant
excited state (FTES) under the umbrella of a unified system-environment
Hamiltonian. The system is two qubits coupled to a common bath of linear
harmonic oscillators. The FTES is an excited state of qubits that behaves like
the ground state. It is excited because its energy is higher than the lowest
eigenenergy of the qubit Hamiltonian. It behaves like the ground state in that
the coupling to the environment at zero temperature produces equilibrium,
where the reduced qubit state is pure.

The FTES is not the Gibbs state at an effective temperature, and exemplifies
break-down of the eigenstate thermalization hypothesis. [83, 84] There is a
universality aspect to the FTES, in that over a wide range of system-
environment parameters in the regime of strong repulsion between weakly
detuned qubit levels caused by the unitary effect of the heat bath, the FTES
has high singlet fidelity and therefore weakly depends on those parameters.
This is analogous to the universality of the Boltzmann distribution, which is
also independent of system-environment parameters over a significant range.

The eigenstate thermalization hypothesis can break down if the system is
integrable. [83] Since the bath of linear harmonic oscillators is integrable,
this may be one reason that the FTES is in violation of the hypothesis. A
possible next step of this research could be to examine if the FTES can occur
with an anharmonic heat bath, which is not integrable.

We have shown that the FTES is the consequence of the GKSL master equation if
two conditions are met: first, an excited eigenstate of the qubit Hamiltonian
is annihilated by the dissipative generator of the equation; second, the
asymptotic state of the equation is unique. A counterargument for the
existence of the FTES would thus be the counterargument for the validity of
the GKSL master equation, which would be challenging to find.

A peculiar feature of our model in the weak system-bath coupling regime is
that the slower the bath, the less sensitive the FTES with respect to changes
in the model parameters. This will make it easier to experimentally find the
parameter range for a FTES if the bath is slower. The reason is that slower
baths have a weaker compensation by the counterterm, leading to stronger net
unitary effect on the reduced qubit system.

An interesting question is what kind of entangled qubit states can be created
if several qubits are coupled to the same environment. We have found fault-
tolerant manifolds of excited states in systems with four and six qubits,
which we will publish in future. Thus, fault tolerant excitations can be
extended to multiple qubits. This leads to a topic worth exploring if FTES can
be used to prepare entangled states as initial states for quantum computation.
This capability would bypass the need for the large number of unitary gates
necessary for state preparation. However, experimental tuning of the
environmental couplings for the FTES will be very challenging as well.

The question that remains to be answered is if the FTES can occur in non-
Markovian quantum dynamics. Our simulations based on TEMPO are consistent with
the master equations, but only on a short time scale compared to the recovery
time of the FTES. We do not have the computational resources to determine if
non-Markovian dynamics can produce a FTES. With better computational resources
and improvements in the algorithm, a study of this question could be doable.
Alternatively, note that FTES shares some properties with decoherence-free-
subspaces, as discussed in the text. The latter are valid for exact quantum
dynamics. It will be worth exploring if FTES can be understood in terms of
some decoherence free proximity space, which in turn would imply the exactness
of FTES for arbitrary strengths of system-bath coupling from the exactness of
decoherence-free-subspace.

FTES is not QEC, at least not yet. QEC protects unknown quantum states, while
FTES is a known quantum state. Expanding the approach presented in this paper
could lead to a Hamiltonian model of QEC.

The authors would like to thank Gerald Fux for discussion and adding a feature
to the open source package on our request, and Brayden Ware for bringing to
our attention the eigenstate thermalization hypothesis. This research was
supported by DOE contract DE-FG02-06ER46281. 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 system.

##  7 Appendix 1. Approximation in Eq. 27

To approximate a master equation well, it is neither necessary nor sufficient
that the coefficients of the before/after equations be very close. Such is the
case with the approximation in Eq. 27. The accuracy of GAME stems from the
coarse-graining properties of master equations, rather than any similarity in
coefficients.

To see how coarse-graining works, we first vectorize the master equation. A
matrix like $\rho$ is replaced by a column vector $|\rho\rangle$ made up of
the columns of the matrix appended one after the other. In this format,
$\mathcal{G}$ becomes a Hermitian matrix,
$\mathcal{G}_{ni,mj}=\mathcal{G}_{mj,ni}^{\star}$, while the ME is

$\frac{d|\rho\rangle}{dt}=\mathcal{R}|\rho\rangle.$ (37)

Note that here the generator $\mathcal{R}$ is more general than in Eq. 26. For
example, it may be the generator of GAME. The generator is often referred to
as a superoperator, and can be obtained from the master equation as described,
for example, in Ref. [85]. Unlike $\mathcal{G}$, $\mathcal{R}$ need not be
Hermitian.

A coarse-grained master equation is obtained by time averaging the
coefficients in the interaction picture. In particular, the coarse-grained
superoperator is

$\tilde{\mathcal{R}}=\frac{1}{T_{c}}\int\limits_{-T_{c}/2}^{T_{c}/2}e^{-\mathcal{R}_{0}\tau}\mathcal{R}e^{\mathcal{R}_{0}\tau}d\tau,$
(38)

where $\mathcal{R}_{0}=-i[H_{S},\bullet]$ is the qubit free Liouvillian, and
$\mathcal{R}_{int}=e^{-\mathcal{R}_{0}\tau}\mathcal{R}e^{\mathcal{R}_{0}\tau}$
is the superoperator in the interaction picture. In terms of matrix elements,
coarse-graining amounts to
$\tilde{\mathcal{R}}_{nm,ij}=\mathcal{R}_{nm,ij}\text{sinc}[(\omega_{nm}-\omega_{ij})T_{c}/2]$,
where $\text{sinc}(x)=\sin(x)/x$.

Figure 16: Conjoining of superoperators of two coarse-grained master equations
versus coarse-graining time $T_{c}$. Lines, top-to bottom: $\alpha=2\times
10^{-k}$, $k=1,2,3,4$. a) $s=1$, $\xi=\alpha/4$, $\eta=0.1788$,
$\zeta=-0.04349$. b) $s=3$, $\xi=\alpha/5$, $\eta=0.37$, $\zeta=0$.
$\tilde{\mathcal{F}}$ is the Frobenius distance between the superoperators of
the RE and GAME master equations. Schrödinger picture. $\Delta=1$,
$\omega_{c}=10$.

The Frobenius distance between coarse-grained superoperators of RE and GAME
decreases inversely with coarse-graining time, as shown in Fig. 16. In two
qubit system, we see that the onset coarse-grain time for the inverse
dependence is of order inverse drive frequency $1/\Delta$.

On a time scale much shorter than the relaxation time, the dynamics of the
density matrix in the interaction picture, under Markov ME, is approximately
the same as the effect of coarse-graining. To see this, note that the dynamics
is

$\varrho(t/2)=\varrho(-t/2)+\int\limits_{-t/2}^{t/2}d\tau\mathcal{R}_{int}(\tau)\varrho(\tau).$
(39)

If $\tau$ is much smaller than the relaxation time, we can make the
approximation $\varrho(\tau)\approx\varrho(0)$, and Eq. 39 becomes

$\varrho(t/2)=\varrho(-t/2)+t\tilde{\mathcal{R}}\varrho(0).$ (40)

Just as the distance between coarse-grained superoperators decreases above the
onset time in Fig. 16, the dynamics under the RE and GAME expressed through
Eq. 40 become closer versus time, validating the approximation given by Eq.
27. Further details can be found in Ref. [19] Our numerical studies find
accuracy comparable to that of the RE Master equation. [19] Equivalent
accuracy is analytically established for ULE. [20]

##  8 Appendix 2. TCL4 Master Equation

We closely follow the derivation of the TCL4 master equation given by Breuer,
Kappler and Petruccione in Ref. [86]. The master equation is obtained as a
perturbative expansion of the Nakajima-Zwanzig equation in $\alpha$. Note that
$\alpha$ in this article is not the same as in that reference. The $\alpha$
that we use is quadratic versus $\alpha$ from Ref. [86]. In the lowest order,
the perturbative expansion gives the time dependent Redfield master equation
13. As discussed in Sec. 3.1, the Markov limit of the equation emerges by
replacing the timed spectral density with its limit at $t\to\infty$.

The next term in the perturbative expansion is second order in $\alpha$. In
the interaction picture, it is given by Eq. 29 in Ref. [86],

$\displaystyle\left(\frac{d\varrho}{dt}\right)_{2}=\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}$
(41) $\displaystyle\times\Big{\\{}\langle 02\rangle\langle
13\rangle\left[\hat{0},\left[\hat{1},\hat{2}\right]\hat{3}\varrho\right]-\langle
02\rangle\langle
31\rangle\left[\hat{0},\left[\hat{1},\hat{2}\right]\varrho\hat{3}\right]-\langle
20\rangle\langle
13\rangle\left[\hat{0},\hat{3}\varrho\left[\hat{1},\hat{2}\right]\right]+\langle
20\rangle\langle
31\rangle\left[\hat{0},\varrho\hat{3}\left[\hat{1},\hat{2}\right]\right]$
$\displaystyle+\langle 03\rangle\langle
12\rangle\left(\left[\hat{0},\left[\hat{3},\hat{2}\right]\varrho\hat{1}\right]+\left[\hat{0},\left[\hat{1}\hat{2},\hat{3}\right]\varrho\right]\right)+\langle
30\rangle\langle
21\rangle\left(\left[\hat{0},\hat{1}\varrho\left[\hat{3},\hat{2}\right]\right]+\left[\hat{0},\varrho\left[\hat{2}\hat{1},\hat{3}\right]\right]\right)$
$\displaystyle-\langle 03\rangle\langle
21\rangle\left[\hat{0},\left[\hat{1},\hat{3}\right]\varrho\hat{2}\right]-\langle
30\rangle\langle
12\rangle\left[\hat{0},\hat{2}\varrho\left[\hat{1},\hat{3}\right]\right]\Big{\\}}.$

Here the short-hand notion is introduced,

$\displaystyle\hat{i}$ $\displaystyle=$ $\displaystyle A(t_{i}),\,\,i\doteq
0,1,2,3$ $\displaystyle\langle ij\rangle$ $\displaystyle=$ $\displaystyle
C(t_{i}-t_{j}),\,\,i,j\doteq 0,1,2,3.$

Evaluating the expression 41 is immensely complicated. For completeness, we
present the final result for $\delta\mathcal{G}_{in,jm}(t)$ in Eq. 31:

$\displaystyle\delta\mathcal{G}_{ni,mj}(t)$ $\displaystyle=$
$\displaystyle\sum\limits_{a,b}A_{ni}A_{ja}A_{ab}A_{bm}\int\limits_{0}^{t}d\tau
e^{i\omega_{mj}t+i\omega_{jb}\tau}H_{bm,ja}(t,\tau)-e^{i\omega_{mj}t+i(\omega_{ja}+\omega_{bm})\tau})H_{ab,ja}(t,\tau)$
$\displaystyle-$
$\displaystyle\sum\limits_{a,b,k}\delta_{ni}A_{ja}A_{ab}A_{bk}A_{km}\int\limits_{0}^{t}d\tau[e^{i\omega_{kj}t+i\omega_{jb}\tau}H_{bk,ja}(t,\tau)-e^{i\omega_{kj}t+i(\omega_{ja}+\omega_{bk})\tau}H_{ab,ja}(t,\tau)]$
$\displaystyle-$
$\displaystyle\sum\limits_{a,b}A_{na}A_{ab}A_{bi}A_{jm}\int\limits_{0}^{t}d\tau
e^{i(\omega_{ij}+\omega_{ma})t+i\omega_{jm}\tau}[e^{i\omega_{ab}\tau}F_{bi,jm}(t,\tau)-e^{i\omega_{bi}\tau}F_{ab,jm}(t,\tau)]$
$\displaystyle+$
$\displaystyle\sum\limits_{a,b}A_{nb}A_{bi}A_{ja}A_{am}\int\limits_{0}^{t}\tau
e^{i(\omega_{ij}+\omega_{an})t+i\omega_{ja}\tau}[e^{i\omega_{nb}\tau}F_{bi,ja}(t,\tau)-e^{i\omega_{bi}\tau}F_{nb,ja}(t,\tau)]$
$\displaystyle+$
$\displaystyle\sum\limits_{a,b}A_{na}A_{ai}A_{jb}A_{bm}\int\limits_{0}^{t}d\tau[e^{i(\omega_{ma}+\omega_{ib})t+i\omega_{ai}\tau}I_{bm,jb}(t,\tau)-e^{i(\omega_{ba}+\omega_{ij})t+i\omega_{ai}\tau}I_{jb,bm}(t,\tau)]$
$\displaystyle-$
$\displaystyle\sum\limits_{a,b}A_{ni}A_{jb}A_{ba}A_{am}\int\limits_{0}^{t}d\tau[e^{i(\omega_{an}+\omega_{ib})t+i\omega_{ni}\tau}I_{ba,jb}(t,\tau)-e^{i(\omega_{bn}+\omega_{ij})t+i\omega_{ni}\tau}I_{jb,ba}(t,\tau)]$
$\displaystyle+$
$\displaystyle\sum\limits_{a,b}A_{ni}A_{jb}A_{ba}A_{am}\int\limits_{0}^{t}d\tau[e^{i\omega_{aj}t+i\omega_{ba}\tau}I_{jb,am}(t,\tau)-e^{i\omega_{aj}t+i\omega_{jb}\tau}I_{ba,am}(t,\tau)]$
$\displaystyle-$
$\displaystyle\sum\limits_{a,b,k}\delta_{ni}A_{jk}A_{ka}A_{ab}A_{bm}\int\limits_{0}^{t}d\tau[e^{i\omega_{aj}t+i\omega_{ka}\tau}I_{jk,ab}(t,\tau)-e^{i\omega_{aj}t+i\omega_{jk}\tau}I_{ka,ab}(t,\tau)]$
$\displaystyle-$
$\displaystyle\sum\limits_{a,b}A_{na}A_{ai}A_{jb}A_{bm}\int\limits_{0}^{t}d\tau\bigg{[}e^{i(\omega_{ba}+\omega_{ij})t+i\omega_{jb}\tau}Y_{ai,bm}(t,\tau)-e^{i(\omega_{ma}+\omega_{ib})t+i\omega_{bm}\tau}Y_{ai,jb}(t,\tau)\bigg{]}$
$\displaystyle+$
$\displaystyle\sum\limits_{a,b}A_{ni}A_{jb}A_{ba}A_{am}\int\limits_{0}^{t}d\tau\bigg{[}e^{i(\omega_{bn}+\omega_{ij})t+i\omega_{jb}\tau}Y_{ni,ba}(t,\tau)-e^{i(\omega_{an}+\omega_{ib})t+i\omega_{ba}\tau}Y_{ni,jb}(t,\tau)\bigg{]}.$

Here we introduce the bath correlation functions

$\displaystyle F_{ab,jm}(t,\tau)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\tau}duC(t-u)e^{i\omega_{ab}u}\big{[}\Gamma_{jm}(\tau)-\Gamma_{jm}(\tau-u)\big{]}^{\star}$
(42) $\displaystyle H_{ab,jm}(t,\tau)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\tau}duC(u-t)e^{i\omega_{ab}u}\big{[}\Gamma_{jm}(\tau)-\Gamma_{jm}(\tau-u)\big{]}^{\star}$
(43) $\displaystyle I_{ab,jm}(t,\tau)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\tau}duC(u-\tau)e^{i\omega_{ab}u}\big{[}\Gamma_{jm}(t)-\Gamma_{jm}(t-u)\big{]}^{\star}$
(44) $\displaystyle Y_{ab,jm}(t,\tau)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\tau}duC(\tau-u)e^{i\omega_{ab}u}\big{[}\Gamma_{jm}(t)-\Gamma_{jm}(t-u)\big{]}^{\star},$
(45)

where $\Gamma_{jm}(t)=\Gamma(\omega_{jm},t)$ is the timed spectral density
given by Eqs. 16 or 17. We must also add the Hermitian-conjugate, e.g.,
substitute
$\delta\mathcal{G}_{ni,mj}(t)\to\delta\mathcal{G}_{ni,mj}(t)+\delta\mathcal{G}_{mj,ni}(t)^{\star}$
before inserting in Eq. 31.

The most time consuming step is to find the correlation functions 42-45. $I$
and $Y$ decay more slowly with time than $F$ and $H$ and require longer $t$ to
reach the Markov limit ($t\to\infty$). For the Ohmic heat bath, $I$ and $Y$
decay inversely with $t-t_{1}$, and calculating the integrals is time
consuming.

First we determine the timed spectral densities analytically using Eqs. 16 or
17, and apply Simpson’s rule with time increment $dt=0.04$ and $n_{0}=2^{17}$
time-steps to find $\delta G_{in,jm}(t)$, where $t=n_{0}dt\approx
52420/\omega_{c}$. This leads to the accuracy of approximately five
significant digits. For $F$ and $H$, $n_{0}=2^{14}$ time steps are sufficient
to reach that accuracy. For super-Ohmic baths, the integrals converge much
faster, and we use $n_{0}=2^{14}$ time steps with $dt=0.02$. To verify the
accuracy of our TCL4 implementation, we reproduced the low temperature results
from Fig. 10.15 in Ref. [24].

An interesting question is should there be a second or higher order
counterterms in order to cancel any remaining $\omega_{c}$ dependence of
$\delta\mathcal{G}$ in the Markov limit? Such higher order counterterms have
been implied by the dots in the sum given in Eq. 4.6 in Ref. [87], but no
analysis of the higher order terms were given. The derivation of the most
general coupling between the system and bath of linear harmonic oscillators
leads to the counterterm being linear with $\alpha$. [22].

So let us examine the second order in $\alpha$ terms of the TCL4 master
equation, and find out if they depend on $\omega_{c}$ for
$\omega_{c}\gg\Delta$. Eqs. 42-45 show that we integrate over difference
between two timed spectral densities at two different times. In the Markov
limit $t\to\infty$, the terms in Eqs. 16 and 17 linear with $\omega_{c}$
effectively cancel after taking the integrals over both $\tau$ and $u$.
Indeed, we calculate the Frobenius distance between $\delta\mathcal{G}$ at
cutoff frequencies $10^{k}\omega_{c}$ and $10^{k+1}\omega_{c}$, for
$k=0,1,2,3$, and find that the distance decreases inversely with the cutoff
frequency. Thus, in the Markov limit, $\delta\mathcal{G}$ is independent of
$\omega_{c}$ at $\omega_{c}\gg\Delta$. For the same reason, we expect that all
higher order terms in the ME expansion of ordered cumulants [24] are also
independent of $\omega_{c}$ under those conditions, which is consistent with
the absence of higher order counterterms.

The Lamb-shift portion independent of $\omega_{c}$ for $\omega_{c}\gg\Delta$
has higher order terms with $\alpha$. In particular, we can identify the
unitary contribution, or the Lamb-shift, in Eq. 8 analogously to how we
identified it in the RE, which leads to

$H_{L}^{(2)}=\frac{\mathcal{F}^{\dagger}-\mathcal{F}}{2i},$ (46)

where

$\mathcal{F}=\langle 20\rangle\langle
31\rangle\hat{3}[\hat{1},\hat{2}]\hat{0}+\langle 30\rangle\langle
21\rangle[\hat{2}\hat{1},\hat{3}]\hat{0},$ (47)

with the matrix elements

$\displaystyle\mathcal{F}_{jm}$ $\displaystyle=$
$\displaystyle\sum\limits_{a,b,k}A_{ja}A_{ab}A_{bk}A_{km}\int\limits_{0}^{t}dt_{1}\Big{[}e^{i\omega_{kj}t+i\omega_{jb}t_{1}}H_{bk,ja}(t,t_{1})-e^{i\omega_{kj}t+i(\omega_{ja}+\omega_{bk})t_{1}}H_{ab,ja}(t,t_{1})$
$\displaystyle+$ $\displaystyle
e^{i\omega_{bj}t+i\omega_{ab}t_{1}}I_{ja,bk}(t,t_{1})-e^{i\omega_{bj}t+i\omega_{ja}t_{1}}I_{ab,bk}(t,t_{1})\Big{]}.$

##  9 Appendix 3. TEMPO Simulation Parameters

Throughout this paper we use the TEMPO software package downloaded from Ref.
[21]. We made some modifications, such as updating the system-bath coupling
Hamiltonian according to Eq. 4 and including the counterterm in Eq. 6. There
are three sources of errors: 1) The Trotter splitting error, which is due to
the finite time step $\tau_{s}$. In all simulations presented here, the time
step is $0.02/\Delta$; 2) The finite memory approximation, originally
introduced by Makri and Makarov, [41, 42] limits the time over which the
histories of the system are stored. This error is measured here by the finite
number of time steps $K_{max}$ in the ADT. That is, the various histories of
the system are stored over the last $K_{max}$ time steps. In our simulations,
$K_{max}=150$ for the Ohmic bath and $K_{max}=75$ for the super-Ohmic bath; 3)
Precision parameter of the SVD-compression. In all our simulations, the
precision is $\epsilon_{SVD}=10^{-6}$, which means that in the tensor network
representation of the ADT, we neglect the singular vectors with singular
values $<10^{-6}$.

We have varied $\tau_{s}$, $K_{max}$, and $\epsilon$ until the numerically
obtained quantum dynamics is roughly converged. Here, that means that changing
the parameters in the direction of higher accuracy, by factors of two, will
change a curve in Figs. 5-6 by approximately five percent. However, the array
of curves in Figs. 5-6 remain very similar as the accuracy is increased,
unambiguously confirming the suppression of singlet damping. The positivity
violations in the density matrix can lead to singlet fidelity exceeding one
near fault tolerance conditions due to these errors.

Fig. 17 shows singlet fidelity versus time at $\eta=0.22$ from Fig. 4 d). The
data between $70/\Delta$ and $100/\Delta$ is fit to an exponential with offset
and a linear function. The differences between best fits and data are shown in
the inset. The exponential fit is much better than the linear, showing that we
can reliably infer the curvature from which we estimate the recovery time. The
latter is $T_{r}=2286/\Delta$, an unfeasible long time for TEMPO, while the
asymptotic state singlet fidelity is $1.29$. Thus, while the solution appears
to be bound, the fidelity exceeding one confirms a positivity violation, which
we have also confirmed by the negative eigenvalues of the density matrix.

Figure 17: Tempo simulation: Red line: Singlet fidelity versus time, for
$\eta=0.22$ in Fig. 4 d). Black line:best fit to exponential, in the interval
$t>70/\Delta$. Inset: Difference between best fit and TEMPO simulation, for
linear (red) and exponential fit (black).

## References

  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. (Cambridge University Press, USA, 2011).
  * Shor [1995] P. W. Shor, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 52,“ “bibinfo pages R2493 (“bibinfo year 1995).
  * Steane [1996] A. M. Steane, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 77,“ “bibinfo pages 793 (“bibinfo year 1996).
  * Wiseman [1994] H. M. Wiseman, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 49,“ “bibinfo pages 2133 (“bibinfo year 1994).
  * Goetsch  _et al._ [1996] P. Goetsch, P. Tombesi, and D. Vitali, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 54,“ “bibinfo pages 4519 (“bibinfo year 1996).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn, “emph “bibinfo title Quantum Measurement and Control (Cambridge University Press, 2009).
  * Ahn  _et al._ [2002] C. Ahn, A. C. Doherty, and A. J. Landahl, “bibfield journal “bibinfo journal Physical Review A“ “textbf “bibinfo volume 65 (“bibinfo year 2002),“ 10.1103/physreva.65.042301.
  * Lloyd [2000] S. Lloyd, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 62,“ “bibinfo pages 022108 (“bibinfo year 2000).
  * Nelson  _et al._ [2000] R. J. Nelson, Y. Weinstein, D. Cory, and S. Lloyd, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 85,“ “bibinfo pages 3045 (“bibinfo year 2000).
  * Alicki  _et al._ [2006] R. Alicki, D. A. Lidar, and P. Zanardi, “bibfield journal “bibinfo journal Physical Review A“ “textbf “bibinfo volume 73 (“bibinfo year 2006),“ 10.1103/physreva.73.052311.
  * Mohseni and Lidar [2005] M. Mohseni and D. A. Lidar, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 94,“ “bibinfo pages 040507 (“bibinfo year 2005).
  * Terhal and Burkard [2005] B. M. Terhal and G. Burkard, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 71,“ “bibinfo pages 012336 (“bibinfo year 2005).
  * Aliferis  _et al._ [2006] P. Aliferis, D. Gottesman, and J. Preskill, Quantum Info. Comput. 6, 97–165 (2006).
  * Aharonov  _et al._ [2006] D. Aharonov, A. Kitaev, and J. Preskill, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 96,“ “bibinfo pages 050504 (“bibinfo year 2006).
  * Sarovar  _et al._ [2004] M. Sarovar, C. Ahn, K. Jacobs, and G. J. Milburn, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 69,“ “bibinfo pages 052324 (“bibinfo year 2004).
  * Sarovar and Milburn [2005] M. Sarovar and G. J. Milburn, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 72,“ “bibinfo pages 012306 (“bibinfo year 2005).
  * Khodjasteh and Lidar [2005] K. Khodjasteh and D. A. Lidar, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 95,“ “bibinfo pages 180501 (“bibinfo year 2005).
  * Mozgunov and Lidar [2020] E. Mozgunov and D. Lidar, “ “textbf “bibinfo volume 4,“ “bibinfo pages 227 (“bibinfo year 2020), 1908.01095 [Quantum] .
  * Davidović [2020] D. Davidović, “bibfield journal “bibinfo journal Quantum“ “textbf “bibinfo volume 4,“ “bibinfo pages 326 (“bibinfo year 2020).
  * Nathan and Rudner [2020] F. Nathan and M. S. Rudner, “bibfield journal “bibinfo journal Phys. Rev. B“ “textbf “bibinfo volume 102,“ “bibinfo pages 115109 (“bibinfo year 2020).
  * Strathearn  _et al._ [2018] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, “bibfield journal “bibinfo journal Nat. Commun.“ “textbf “bibinfo volume 9,“ “bibinfo pages 3322 (“bibinfo year 2018).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, “bibfield journal “bibinfo journal Annals of Physics“ “textbf “bibinfo volume 149,“ “bibinfo pages 374 (“bibinfo year 1983).
  * Weiss [2012] U. Weiss, (2012).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione, “The theory of open quantum systems,” (2007).
  * Hartmann and Strunz [2020a] R. Hartmann and W. T. Strunz, “bibfield journal “bibinfo journal Quantum“ “textbf “bibinfo volume 4,“ “bibinfo pages 347 (“bibinfo year 2020“natexlaba).
  * Alicki and Lendi [2014] R. Alicki and K. Lendi,  _Quantum Dynamical Semigroups and Applications_ (Springer Publishing Company, Incorporated, 2014).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, “bibfield journal “bibinfo journal Journal of Mathematical Physics“ “textbf “bibinfo volume 17,“ “bibinfo pages 821 (“bibinfo year 1976), https://aip.scitation.org/doi/pdf/10.1063/1.522979 .
  * Lindblad [1976] G. Lindblad, “bibfield journal “bibinfo journal Comm. Math. Phys.“ “textbf “bibinfo volume 48,“ “bibinfo pages 119 (“bibinfo year 1976).
  * Davies [1974] E. B. Davies, “bibfield journal “bibinfo journal Comm. Math. Phys.“ “textbf “bibinfo volume 39,“ “bibinfo pages 91 (“bibinfo year 1974).
  * Chen and Lidar [2020] H. Chen and D. A. Lidar, “Hoqst: Hamiltonian open quantum system toolkit,” (2020), arXiv:2011.14046 [quant-ph] .
  * Majenz  _et al._ [2013] C. Majenz, T. Albash, H.-P. Breuer, and D. A. Lidar, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 88,“ “bibinfo pages 012103 (“bibinfo year 2013).
  * Schaller and Brandes [2008] G. Schaller and T. Brandes, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 78,“ “bibinfo pages 022106 (“bibinfo year 2008).
  * Albash  _et al._ [2012] T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, “bibfield journal “bibinfo journal New Journal of Physics“ “textbf “bibinfo volume 14,“ “bibinfo pages 123016 (“bibinfo year 2012).
  * Vogt  _et al._ [2013] N. Vogt, J. Jeske, and J. H. Cole, “bibfield journal “bibinfo journal Phys. Rev. B“ “textbf “bibinfo volume 88,“ “bibinfo pages 174514 (“bibinfo year 2013).
  * Tscherbul and Brumer [2015] T. V. Tscherbul and P. Brumer, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 142,“ “bibinfo pages 104107 (“bibinfo year 2015), https://doi.org/10.1063/1.4908130 .
  * Kiršanskas  _et al._ [2018] G. Kiršanskas, M. Franckié, and A. Wacker, “bibfield journal “bibinfo journal Phys. Rev. B“ “textbf “bibinfo volume 97,“ “bibinfo pages 035432 (“bibinfo year 2018).
  * Lee and Yeo [2020] J. S. Lee and J. Yeo, “Comment on ”universal lindblad equation for open quantum systems”,” (2020), arXiv:2011.00735 [quant-ph] .
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, “bibfield journal “bibinfo journal Rev. Mod. Phys.“ “textbf “bibinfo volume 59,“ “bibinfo pages 1 (“bibinfo year 1987).
  * Hur [2008] K. L. Hur, “bibfield journal “bibinfo journal Annals of Physics“ “textbf “bibinfo volume 323,“ “bibinfo pages 2208 (“bibinfo year 2008).
  * Feynman and Vernon [1963] R. Feynman and F. Vernon, “bibfield journal “bibinfo journal Annals of Physics“ “textbf “bibinfo volume 24,“ “bibinfo pages 118 (“bibinfo year 1963).
  * Makri and Makarov [1995a] N. Makri and D. E. Makarov, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 102,“ “bibinfo pages 4600 (“bibinfo year 1995“natexlaba), https://doi.org/10.1063/1.469508 .
  * Makri and Makarov [1995b] N. Makri and D. E. Makarov, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 102,“ “bibinfo pages 4611 (“bibinfo year 1995“natexlabb), https://doi.org/10.1063/1.469509 .
  * Thorwart  _et al._ [2004] M. Thorwart, E. Paladino, and M. Grifoni, “bibfield journal “bibinfo journal Chemical Physics“ “textbf “bibinfo volume 296,“ “bibinfo pages 333 (“bibinfo year 2004), the Spin-Boson Problem: From Electron Transfer to Quantum Computing … to the 60th Birthday of Professor Ulrich Weiss.
  * Nalbach and Thorwart [2010] P. Nalbach and M. Thorwart, “bibfield journal “bibinfo journal Phys. Rev. B“ “textbf “bibinfo volume 81,“ “bibinfo pages 054308 (“bibinfo year 2010).
  * Mühlbacher  _et al._ [2005] L. Mühlbacher, J. Ankerhold, and A. Komnik, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 95,“ “bibinfo pages 220404 (“bibinfo year 2005).
  * Chin  _et al._ [2010] A. W. Chin, n. Rivas, S. F. Huelga, and M. B. Plenio, “bibfield journal “bibinfo journal Journal of Mathematical Physics“ “textbf “bibinfo volume 51,“ “bibinfo pages 092109 (“bibinfo year 2010), https://doi.org/10.1063/1.3490188 .
  * Prior  _et al._ [2010] J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 105,“ “bibinfo pages 050404 (“bibinfo year 2010).
  * Strathearn [2020] A. Strathearn, “emph “bibinfo title Modelling Non-Markovian Quantum Systems Using Tensor Networks, Springer Theses (Springer International Publishing, Cham, 2020).
  * Schröder  _et al._ [2019] F. A. Y. N. Schröder, D. H. P. Turban, A. J. Musser, N. D. M. Hine, and A. W. Chin, “bibfield journal “bibinfo journal Nature Communications“ “textbf “bibinfo volume 10,“ “bibinfo pages 1062 (“bibinfo year 2019).
  * Wall  _et al._ [2016] M. L. Wall, A. Safavi-Naini, and A. M. Rey, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 94,“ “bibinfo pages 053637 (“bibinfo year 2016).
  * Pollock  _et al._ [2018] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 97,“ “bibinfo pages 012127 (“bibinfo year 2018).
  * Luchnikov  _et al._ [2019] I. A. Luchnikov, S. V. Vintskevich, H. Ouerdane, and S. N. Filippov, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 122,“ “bibinfo pages 160401 (“bibinfo year 2019).
  * Tanimura and Kubo [1989] Y. Tanimura and R. Kubo, “bibfield journal “bibinfo journal Journal of the Physical Society of Japan“ “textbf “bibinfo volume 58,“ “bibinfo pages 1199 (“bibinfo year 1989), https://doi.org/10.1143/JPSJ.58.1199 .
  * Tanimura [2006] Y. Tanimura, “bibfield journal “bibinfo journal Journal of the Physical Society of Japan“ “textbf “bibinfo volume 75,“ “bibinfo pages 082001 (“bibinfo year 2006), https://doi.org/10.1143/JPSJ.75.082001 .
  * Tanimura [2014] Y. Tanimura, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 141,“ “bibinfo pages 044114 (“bibinfo year 2014), https://doi.org/10.1063/1.4890441 .
  * Li  _et al._ [2012] Z. Li, N. Tong, X. Zheng, D. Hou, J. Wei, J. Hu, and Y. Yan, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 109,“ “bibinfo pages 266403 (“bibinfo year 2012).
  * Cheng  _et al._ [2015] Y. Cheng, W. Hou, Y. Wang, Z. Li, J. Wei, and Y. Yan, “bibfield journal “bibinfo journal New Journal of Physics“ “textbf “bibinfo volume 17,“ “bibinfo pages 033009 (“bibinfo year 2015).
  * Meyer  _et al._ [1990] H.-D. Meyer, U. Manthe, and L. Cederbaum, “bibfield journal “bibinfo journal Chemical Physics Letters“ “textbf “bibinfo volume 165,“ “bibinfo pages 73 (“bibinfo year 1990).
  * Beck  _et al._ [2000] M. Beck, A. Jäckle, G. Worth, and H.-D. Meyer, “bibfield journal “bibinfo journal Physics Reports“ “textbf “bibinfo volume 324,“ “bibinfo pages 1 (“bibinfo year 2000).
  * Wang and Thoss [2003] H. Wang and M. Thoss, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 119,“ “bibinfo pages 1289 (“bibinfo year 2003), https://doi.org/10.1063/1.1580111 .
  * Zheng  _et al._ [2016] J. Zheng, Y. Xie, S. Jiang, and Z. Lan, “bibfield journal “bibinfo journal The Journal of Physical Chemistry C“ “textbf “bibinfo volume 120,“ “bibinfo pages 1375 (“bibinfo year 2016), https://doi.org/10.1021/acs.jpcc.5b09921 .
  * Nakajima [1958] S. Nakajima, “bibfield journal “bibinfo journal Progress of Theoretical Physics“ “textbf “bibinfo volume 20,“ “bibinfo pages 948 (“bibinfo year 1958), https://academic.oup.com/ptp/article-pdf/20/6/948/5440766/20-6-948.pdf .
  * Zwanzig [1960] R. Zwanzig, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 33,“ “bibinfo pages 1338 (“bibinfo year 1960), https://doi.org/10.1063/1.1731409 .
  * Suess  _et al._ [2014] D. Suess, A. Eisfeld, and W. T. Strunz, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 113,“ “bibinfo pages 150403 (“bibinfo year 2014).
  * Zhang and Eisfeld [2016] P.-P. Zhang and A. Eisfeld, “bibfield journal “bibinfo journal The Journal of Physical Chemistry Letters“ “textbf “bibinfo volume 7,“ “bibinfo pages 4488 (“bibinfo year 2016), pMID: 27775345, https://doi.org/10.1021/acs.jpclett.6b02111 .
  * Hartmann and Strunz [2017] R. Hartmann and W. T. Strunz, “bibfield journal “bibinfo journal Journal of Chemical Theory and Computation“ “textbf “bibinfo volume 13,“ “bibinfo pages 5834 (“bibinfo year 2017), pMID: 29016126, https://doi.org/10.1021/acs.jctc.7b00751 .
  * Orús [2014] R. Orús, “bibfield journal “bibinfo journal Annals of Physics“ “textbf “bibinfo volume 349,“ “bibinfo pages 117 (“bibinfo year 2014).
  * Cirac  _et al._ [2020] I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, “Matrix product states and projected entangled pair states: Concepts, symmetries, and theorems,” (2020), arXiv:2011.12127 [quant-ph] .
  * Jørgensen and Pollock [2019] M. R. Jørgensen and F. A. Pollock, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 123,“ “bibinfo pages 240602 (“bibinfo year 2019).
  * Fux  _et al._ [2021] G. E. Fux, E. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling, “ (“bibinfo year 2021), arXiv:2101.03071 .
  * [71] D. Gribben, G. E. Fux, A. Strathearn, P. Kirton, and B. W. Lovett, .
  * The TEMPO collaboration [2020] The TEMPO collaboration, “enquote “bibinfo title TimeEvolvingMPO: A Python 3 package to efficiently compute non-Markovian open quantum systems.“ (2020).
  * Braun [2002] D. Braun, “bibfield journal “bibinfo journal Phys. Rev. Lett.“ “textbf “bibinfo volume 89,“ “bibinfo pages 277901 (“bibinfo year 2002).
  * McCutcheon  _et al._ [2009] D. P. S. McCutcheon, A. Nazir, S. Bose, and A. J. Fisher, “bibfield journal “bibinfo journal Physical Review A“ “textbf “bibinfo volume 80 (“bibinfo year 2009),“ 10.1103/physreva.80.022337.
  * Benatti  _et al._ [2009] F. Benatti, R. Floreanini, and U. Marzolino, “bibfield journal “bibinfo journal EPL (Europhysics Letters)“ “textbf “bibinfo volume 88,“ “bibinfo pages 20011 (“bibinfo year 2009).
  * Kast and Ankerhold [2014] D. Kast and J. Ankerhold, “bibfield journal “bibinfo journal Phys. Rev. B“ “textbf “bibinfo volume 90,“ “bibinfo pages 100301 (“bibinfo year 2014).
  * Sahrapour and Makri [2013] M. M. Sahrapour and N. Makri, “bibfield journal “bibinfo journal The Journal of Chemical Physics“ “textbf “bibinfo volume 138,“ “bibinfo pages 114109 (“bibinfo year 2013), https://doi.org/10.1063/1.4795159 .
  * Benatti  _et al._ [2010] F. Benatti, R. Floreanini, and U. Marzolino, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 81,“ “bibinfo pages 012105 (“bibinfo year 2010).
  * Rivas [2017] A. Rivas, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 95,“ “bibinfo pages 042104 (“bibinfo year 2017).
  * Hartmann and Strunz [2020b] R. Hartmann and W. T. Strunz, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 101,“ “bibinfo pages 012103 (“bibinfo year 2020“natexlabb).
  * Zanardi and Rasetti [1997] P. Zanardi and M. Rasetti, “bibfield journal “bibinfo journal Physical Review Letters“ “textbf “bibinfo volume 79,“ “bibinfo pages 3306–3309 (“bibinfo year 1997).
  * Lidar  _et al._ [1998] D. A. Lidar, I. L. Chuang, and K. B. Whaley, “bibfield journal “bibinfo journal Physical Review Letters“ “textbf “bibinfo volume 81,“ “bibinfo pages 2594–2597 (“bibinfo year 1998).
  * Deutsch [2018] J. M. Deutsch, “bibfield journal “bibinfo journal Reports on Progress in Physics“ “textbf “bibinfo volume 81,“ “bibinfo pages 082001 (“bibinfo year 2018).
  * Garrison and Grover [2018] J. R. Garrison and T. Grover, “bibfield journal “bibinfo journal Phys. Rev. X“ “textbf “bibinfo volume 8,“ “bibinfo pages 021026 (“bibinfo year 2018).
  * Havel [2003] T. F. Havel, “bibfield journal “bibinfo journal Journal of Mathematical Physics“ “textbf “bibinfo volume 44,“ “bibinfo pages 534 (“bibinfo year 2003).
  * Breuer  _et al._ [1999] H.-P. Breuer, B. Kappler, and F. Petruccione, “bibfield journal “bibinfo journal Physical Review A“ “textbf “bibinfo volume 59,“ “bibinfo pages 1633–1643 (“bibinfo year 1999).
  * Alicki [1989] R. Alicki, “bibfield journal “bibinfo journal Phys. Rev. A“ “textbf “bibinfo volume 40,“ “bibinfo pages 4077 (“bibinfo year 1989).