of 10
Supplemental Material for “Quasi-Lindblad pseudomode theory for open quantum
systems”
Gunhee Park (
박
ᄀ
ᅥ
ᆫ
희
),
1
Zhen Huang,
2
Yuanran Zhu,
3
Chao Yang,
3
Garnet Kin-Lic Chan,
4
and Lin Lin
2, 3
1
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California 91125, USA
2
Department of Mathematics, University of California, Berkeley, California 94720, USA
3
Applied Mathematics and Computational Research Division,
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
4
Division of Chemistry and Chemical Engineering,
California Institute of Technology, Pasadena, California 91125, USA
CONTENTS
SM-I. Proof of the equivalence condition for the bosonic Gaussian bath
1
SM-II. BCF expression for the spin-boson model
3
SM-III. BCF expression for the fermionic impurity model
4
SM-IV. Optimal gauge choice minimizing violation of CP condition
6
SM-V. Proofs of the stability condition in the noninteracting fermionic impurity model
6
SM-VI. Details of Fitting procedures
7
A. ESPRIT algorithm for complex exponential function fitting with complex weights
7
B. Complex exponential function fitting with positive weights
8
SM-VII. Details on DMRG calculations
8
References
9
SM-I. PROOF OF THE EQUIVALENCE CONDITION FOR THE BOSONIC GAUSSIAN BATH
In this section, we provide the derivation of the equivalence condition for the bosonic Gaussian bath based on
the bath correlation function (BCF). Our strategy here is different from that in [1], which is only applicable to the
Lindblad equation and relies on the dilation of the Lindblad system into an infinite dimensional unitary system.
We start from the following system-bath Liouvillian superoperator form introduced in the main text,
L
SA
=
−
i
X
j
S
j
F
j
+
i
X
j
e
S
j
e
F
j
.
(S1)
For example, in the case when there is only Hamiltonian coupling,
ˆ
H
SA
=
P
j
ˆ
S
j
ˆ
A
j
,
F
j
•
=
ˆ
A
j
•
and
e
F
j
•
=
•
ˆ
A
j
. In
the case of the quasi-Lindblad pseudomode with both Hamiltonian and Lindblad couplings with
D
SA
•
=
P
j
ˆ
L
′
j
•
ˆ
S
j
+
ˆ
S
j
•
ˆ
L
′†
j
−
1
2
{
ˆ
S
j
ˆ
L
′
j
+
ˆ
L
′†
j
ˆ
S
j
,
•}
,
F
j
•
=
ˆ
A
j
•
+
•
i
ˆ
L
′†
j
−
i
2
(
ˆ
L
′
j
+
ˆ
L
′†
j
)
•
,
e
F
j
•
=
•
ˆ
A
j
−
i
ˆ
L
′
j
•
+
•
i
2
(
ˆ
L
′
j
+
ˆ
L
′†
j
)
.
(S2)
For brevity, we will not explicitly distinguish between the non-tilde and tilde superoperators, and denote
L
SA
=
−
i
P
j
S
j
F
j
, treating
e
S
j
=
S
j
max
+
j
and
e
F
j
=
−F
j
max
+
j
. We will return to the explicit tilde superoperator notation
when the distinction is necessary.
2
In the interaction picture,
ˆ
e
ρ
SA
(
t
) =
e
−
(
L
S
+
L
A
)
t
ˆ
ρ
SA
(
t
),
S
j
(
t
) =
e
−L
S
t
S
j
e
L
S
t
,
F
j
(
t
) =
e
−L
A
t
F
j
e
L
A
t
, and
L
SA
(
t
) =
−
i
P
j
S
j
(
t
)
F
j
(
t
), the von Neumann equation of motion for
ˆ
e
ρ
SA
(
t
) is,
∂
t
ˆ
e
ρ
SA
(
t
) =
L
SA
(
t
)
ˆ
e
ρ
SA
(
t
)
.
(S3)
Note that the backward evolution operators such as
e
−L
A
t
are introduced mainly to simplify the notation, and these
operators do not appear explicitly after applying the time-ordering operations below. The formal solution of
ˆ
e
ρ
SA
(
t
)
can be obtained with the Dyson series expansion [2],
ˆ
e
ρ
SA
(
t
) =
∞
X
m
=0
1
m
!
T

Z
t
0
···
Z
t
0
L
SA
(
t
1
)
···L
SA
(
t
m
)
dt
1
···
dt
m

ˆ
ρ
SA
(0)
,
(S4)
where
T
refers to the time-ordering superoperator. We trace out the bath to obtain the system-reduced density
operator,
ˆ
e
ρ
S
(
t
) =
∞
X
m
=0
(
−
i
)
m
m
!
Z
t
0
···
Z
t
0
X
j
1
···
X
j
m
Tr
A
[
T
A
(
F
j
1
(
t
1
)
···F
j
m
(
t
m
)) ˆ
ρ
A
(0)]
T
S
(
S
j
1
(
t
1
)
···S
j
m
(
t
m
)) ˆ
ρ
S
(0)
.
(S5)
We separate the time-ordering superoperator into the system time-ordering superoperator
T
S
and the bath time-
ordering superoperator
T
A
to be applied within each Liouville operator space. We denote
C
j
1
,
···
,j
m
(
t
1
,
···
,t
m
) =
Tr
A
[
T
A
(
F
j
1
(
t
1
)
···F
j
m
(
t
m
)) ˆ
ρ
A
(0)]. From Wick’s theorem, if
m
is odd, it is zero, and if
m
is even (
m
= 2
n
),
C
j
1
,
···
,j
2
n
(
t
1
,
···
,t
2
n
) =
X
σ
∈
Π
2
n
n
Y
i
=1
C
j
σ
(2
i
−
1)
,j
σ
(2
i
)
(
t
σ
(2
i
−
1)
,t
σ
(2
i
)
)
.
(S6)
After inserting it into Eq. S5,
ˆ
e
ρ
S
(
t
) =
∞
X
n
=0
(
−
1)
n
(2
n
)!
Z
t
0
···
Z
t
0
X
j
1
···
X
j
n
X
σ
∈
Π
2
n
T
S
n
Y
i
=1
C
j
σ
(2
i
−
1)
,j
σ
(2
i
)
(
t
σ
(2
i
−
1)
,t
σ
(2
i
)
)
S
j
1
(
t
1
)
···S
j
n
(
t
n
)
!
ˆ
ρ
S
(0)
dt
1
···
dt
n
(S7)
Thanks to the system time ordering superoperator,
T
S
,
T
S
Q
n
i
=1
C
j
σ
(2
i
−
1)
,j
σ
(2
i
)
(
t
σ
(2
i
−
1)
,t
σ
(2
i
)
)
S
j
1
(
t
1
)
···S
j
n
(
t
n
)

be-
comes the same operator for all
σ
∈
Π
2
n
. The number of elements of Π
2
n
is (2
n
)!!, so Eq. S7 is reduced to,
ˆ
e
ρ
S
(
t
) =
∞
X
n
=0
(
−
1)
n
2
n
n
!
Z
t
0
···
Z
t
0
X
j
1
···
X
j
n
T
S
n
Y
i
=1
C
j
2
i
−
1
,j
2
i
(
t
2
i
−
1
,t
2
i
)
S
j
1
(
t
1
)
···S
j
n
(
t
n
)
!
ˆ
ρ
S
(0)
dt
1
···
dt
n
(S8)
=
T
S
exp
−
1
2
Z
t
0
Z
t
0
X
j
1
,j
2
C
j
1
,j
2
(
t
1
,t
2
)
S
j
1
(
t
1
)
S
j
2
(
t
2
)
dt
1
dt
2
ˆ
ρ
S
(0)
(S9)
=
T
S
exp
−
Z
t
0
Z
t
1
0
X
j
1
,j
2
C
j
1
,j
2
(
t
1
,t
2
)
S
j
1
(
t
1
)
S
j
2
(
t
2
)
dt
1
dt
2
ˆ
ρ
S
(0)
.
(S10)
From the second to the third line, we fix the ordering between
t
1
and
t
2
as
t
1
> t
2
. After returning to the original
frame and parameterizing the integrand in the time non-local exponential superoperator,
ˆ
ρ
S
(
t
) =
e
L
S
t
T
S
exp

Z
t
0
Z
t
1
0
W
(
t
1
,t
2
)
dt
1
dt
2

ˆ
ρ
S
(0)
,
(S11)
which is termed the influence superoperator [3–5]. The expression for
W
(
t
1
,t
2
) after reviving the tilde superoperators
becomes,
W
(
t
1
,t
2
) =
−
X
jj
′

⟨F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
A
S
j
(
t
1
)
S
j
′
(
t
2
)
−⟨F
j
(
t
1
)
e
F
j
′
(
t
2
)
⟩
A
S
j
(
t
1
)
e
S
j
′
(
t
2
)
−⟨
e
F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
A
e
S
j
(
t
1
)
S
j
′
(
t
2
) +
⟨
e
F
j
(
t
1
)
e
F
j
′
(
t
2
)
⟩
A
e
S
j
(
t
1
)
e
S
j
′
(
t
2
)

,
(S12)
3
where
⟨•⟩
A
= Tr
A
[
•
ˆ
ρ
A
(0)]. This expression can be further simplified by using the following property, Tr
A
[
F
j
•
] =
Tr
A
[
e
F
j
•
]. For the Hamiltonian coupling part, it is clear to see Tr
A
[
ˆ
A
j
•
] = Tr
A
[
•
ˆ
A
j
]. For the Lindblad coupling part
in Eq. S2, it can be seen from direct calculation,
Tr
A
[
i
•
ˆ
L
′†
j
−
i
2
(
ˆ
L
′
j
+
ˆ
L
′†
j
)
•
] = Tr
A
[
i
2
(
ˆ
L
′†
j
−
ˆ
L
′
j
)
•
] = Tr
A
[
−
i
ˆ
L
′
j
•
+
i
2
•
(
ˆ
L
′
j
+
ˆ
L
′†
j
)]
.
(S13)
Combining this property with the trace-preserving property, Tr
A
[
L
A
•
] = 0 or Tr
A
[
e
−L
A
t
•
] = Tr
A
[
•
], we have,
⟨F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
A
=
⟨
e
F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
A
,
⟨F
j
(
t
1
)
e
F
j
′
(
t
2
)
⟩
A
=
⟨
e
F
j
(
t
1
)
e
F
j
′
(
t
2
)
⟩
A
.
(S14)
Furthermore, by using the fact that (
F
j
•
)
†
=
e
F
j
(
•
)
†
, we have,
⟨
e
F
j
(
t
1
)
e
F
j
′
(
t
2
)
⟩
A
=
⟨F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
∗
A
.
(S15)
Therefore, by defining the BCF as,
C
jj
′
(
t
1
,t
2
) =
⟨F
j
(
t
1
)
F
j
′
(
t
2
)
⟩
A
= Tr
A
[
F
j
e
L
A
(
t
1
−
t
2
)
F
j
′
e
L
A
(
t
2
)
ˆ
ρ
A
(0)]
,
(S16)
the superoperator
W
(
t
1
,t
2
) is expressed as,
W
(
t
1
,t
2
) =
−
X
jj
′

S
j
(
t
1
)
−
e
S
j
(
t
1
)

×

C
jj
′
(
t
1
,t
2
)
S
j
′
(
t
2
)
−
C
∗
jj
′
(
t
1
,t
2
)
e
S
j
′
(
t
2
)

.
(S17)
We see, therefore, that the influence of the bath on the system dynamics, as expressed through the influence super-
operator, is entirely parametrized by the BCF. Thus, any two baths with the same BCF give rise to the same system
dynamics, which is the equivalence condition.
SM-II. BCF EXPRESSION FOR THE SPIN-BOSON MODEL
In this section, we derive the BCF expression for the spin-boson model,
C
A
(∆
t
) = (
V
−
i
M
)
e
−
Z
∆
t
(
V
+
i
M
)
†
.
(S18)
This can be shown from a direct calculation. First,
e
L
A
t
′
ˆ
ρ
A
(0) =
e
L
A
t
′
|
0
⟩⟨
0
|
=
|
0
⟩⟨
0
|
from
L
A
|
0
⟩⟨
0
|
= 0. We have,
F
j
′
|
0
⟩⟨
0
|
=
X
k
(
V
∗
j
′
k
−
iM
∗
j
′
k
)
ˆ
d
†
k
|
0
⟩⟨
0
|
,
(S19)
from
ˆ
d
k
|
0
⟩⟨
0
|
=
|
0
⟩⟨
0
|
ˆ
d
†
k
= 0. Then, from
L
A
[
ˆ
d
†
k
|
0
⟩⟨
0
|
] =
P
l
(
−
iH
A
lk
−
Γ
lk
)
ˆ
d
†
l
|
0
⟩⟨
0
|
, we get
e
L
A
t
h
P
k
(
V
∗
j
′
k
−
iM
∗
j
′
k
)
ˆ
d
†
k
|
0
⟩⟨
0
|
i
=
P
lk
(
e
−
Z
t
)
lk
(
V
∗
j
′
k
−
iM
∗
j
′
k
)
ˆ
d
†
l
|
0
⟩⟨
0
|
where
Z
=
Γ
+
i
H
A
. Using the property in Eq. S13,
Tr
A
[
F
j
•
] = Tr
A
"
(
X
k
(
V
jk
−
iM
jk
)
ˆ
d
k
+ (
V
∗
jk
+
iM
∗
jk
)
ˆ
d
†
k
)
•
#
.
(S20)
Therefore, this leads to,
C
A
jj
′
(
t
+
t
′
,t
′
) = Tr
A
h
F
j
e
L
A
t
F
j
′
e
L
A
t
′
ˆ
ρ
A
(0)
i
=
X
kl
(
V
jl
−
iM
jl
)(
e
−
Z
t
)
lk
(
V
∗
j
′
k
−
iM
∗
j
′
k
)
,
(S21)
or in the matrix representation,
C
A
(
t
) = (
V
−
i
M
)
e
−
Z
t
(
V
+
i
M
)
†
.
(S22)
4
SM-III. BCF EXPRESSION FOR THE FERMIONIC IMPURITY MODEL
In this section, we discuss the quasi-Lindblad pseudomode formulation for the fermionic impurity model. The
pseudomode model contains two different auxiliary baths,
A
1
and
A
2
with the initial state, ˆ
ρ
SA
(0) = ˆ
ρ
S
(0)
⊗
ˆ
ρ
A
(0)
with ˆ
ρ
A
(0) = ˆ
ρ
A
1
(0)
⊗
ˆ
ρ
A
2
(0) =
|
0
⟩⟨
0
|⊗|
1
⟩⟨
1
|
. We define index sets,
I
1
and
I
2
, to denote the auxiliary mode index
k
1
∈
I
1
(
k
2
∈
I
2
) for the fermions in the auxiliary bath
A
1
(
A
2
), respectively. The density operator evolves under the
following Lindblad equation,
∂
t
ˆ
ρ
SA
=
−
i
[
ˆ
H
aux
,
ˆ
ρ
SA
] + (
D
A
1
+
D
A
2
+
D
SA
1
+
D
SA
2
) ˆ
ρ
SA
.
(S23)
The auxiliary Hamiltonian is given by,
ˆ
H
aux
=
ˆ
H
S
+
X
k
∈
I
1
∪
I
2
E
k
ˆ
d
†
k
ˆ
d
k
+
X
j
ˆ
a
†
j
ˆ
A
j
+
ˆ
A
†
j
ˆ
a
j
,
ˆ
A
j
=
X
k
∈
I
1
(
V
1
)
jk
1
ˆ
d
k
1
+
X
k
∈
I
2
(
V
2
)
jk
2
ˆ
d
k
2
.
(S24)
The dissipators are given by,
D
A
1
•
=
X
k
1
∈
I
1
2
γ
k
1

ˆ
d
k
1
•
ˆ
d
†
k
1
−
1
2
{
ˆ
d
†
k
1
ˆ
d
k
1
,
•}

,
D
A
2
•
=
X
k
2
∈
I
2
2
γ
k
2

ˆ
d
†
k
2
•
ˆ
d
k
2
−
1
2
{
ˆ
d
k
2
ˆ
d
†
k
2
,
•}

,
D
SA
1
•
=
X
j
ˆ
a
j
•
ˆ
L
†
1
j
+
ˆ
L
1
j
•
ˆ
a
†
j
−
1
2
{
ˆ
L
†
1
j
ˆ
a
j
+ ˆ
a
†
j
ˆ
L
1
j
,
•}
,
D
SA
2
•
=
X
j
ˆ
a
†
j
•
ˆ
L
2
j
+
ˆ
L
†
2
j
•
ˆ
a
j
−
1
2
{
ˆ
L
2
j
ˆ
a
†
j
+ ˆ
a
j
ˆ
L
†
2
j
,
•}
,
(S25)
with
ˆ
L
1
j
= 2
P
k
1
∈
I
1
(
M
1
)
jk
1
ˆ
d
k
1
and
ˆ
L
2
j
= 2
P
k
2
∈
I
2
(
M
2
)
jk
2
ˆ
d
k
2
.
We adopt a super-fermion representation [6–9] of Lindblad dynamics to formulate the superoperator formalism
for the fermionic system. This representation maps a density operator to a state vector, ˆ
ρ
7→ |
ρ
⟩⟩
, in the so-called
‘super-Fock’ space with twice the number of original orbitals. This is done with a left-vacuum vector,
|
I
⟩⟩
=
Y
k
(
ˆ
d
†
k
+
ˆ
̃
d
†
k
)
|
0
⟩⊗|
0
⟩
,
(S26)
where
̃
d
†
k
is a fermionic creation operator in the ‘tilde’ space. The density operator is then represented as
|
ρ
⟩⟩
=
( ˆ
ρ
⊗
I
)
|
I
⟩⟩
. For example, the initial bath density operator, ˆ
ρ
A
(0) =
|
0
⟩⟨
0
|⊗|
1
⟩⟨
1
|
becomes,
ˆ
ρ
A
(0)
7→|
ρ
A
(0)
⟩⟩
=
|
01
⟩⟩
⊗
I
1
⊗|
10
⟩⟩
⊗
I
2
.
(S27)
The trace operation can also be expressed with the left-vacuum vector,
Tr[
ˆ
O
ˆ
ρ
] =
⟨⟨
I
|
ˆ
O
⊗
I
|
ρ
⟩⟩
.
(S28)
The following conjugation rules for the left-vacuum vector will become useful,
ˆ
d
k
|
I
⟩⟩
=
ˆ
̃
d
k
|
I
⟩⟩
,
ˆ
d
†
k
=
−
ˆ
̃
d
†
k
|
I
⟩⟩
.
(S29)
We will assume that the initial density operator is an even-parity density operator for simplicity. Any superoperator
L
can be cast to the operator in the super-Fock space,
ˆ
L
, or
L
ˆ
ρ
7→
ˆ
L|
ρ
⟩⟩
. For example, the following dissipator maps
to,
D•
= 2
γ
k

ˆ
d
k
•
ˆ
d
†
k
−
1
2
{
ˆ
d
†
k
ˆ
d
k
,
•}

7→
ˆ
D
= 2
γ
k
ˆ
̃
d
†
k
ˆ
d
k
−
γ
k
ˆ
d
†
k
ˆ
d
k
+
γ
k
ˆ
̃
d
†
k
ˆ
̃
d
k
−
γ
k
.
(S30)
This relation can be shown from the direct calculation of (
D
ˆ
ρ
)
|
I
⟩⟩
=
ˆ
D|
ρ
⟩⟩
with the conjugation rule in Eq. S29. The
bath Liouvillian superoperator is written as,
ˆ
L
A
=
X
k
1
∈
I
1
(
−
iE
k
1
−
γ
k
1
)
ˆ
d
†
k
1
ˆ
d
k
1
+ (
−
iE
k
1
+
γ
k
1
)
ˆ
̃
d
†
k
1
ˆ
̃
d
k
1
+ 2
γ
k
1
ˆ
̃
d
†
k
1
ˆ
d
k
1
+
iE
k
1
−
γ
k
1
+
X
k
2
∈
I
2
(
−
iE
k
2
+
γ
k
2
)
ˆ
d
†
k
2
ˆ
d
k
2
+ (
−
iE
k
2
−
γ
k
2
)
ˆ
̃
d
†
k
2
ˆ
̃
d
k
2
+ 2
γ
k
2
ˆ
d
†
k
2
ˆ
̃
d
k
2
+
iE
k
2
−
γ
k
2
.
(S31)
5
Note that
ˆ
L
A
|
ρ
A
(0)
⟩⟩
= 0. With this representation, we decompose the coupling superoperators in the following form,
ˆ
L
SA
=
−
i
X
j

ˆ
a
†
j
ˆ
F
−
j
+
ˆ
F
+
j
ˆ
a
j
+
ˆ
̃
a
†
j
ˆ
e
F
−
j
+
ˆ
e
F
+
j
ˆ
̃
a
j

.
(S32)
where
ˆ
F
−
j
=
ˆ
A
j
−
i
2
ˆ
L
1
j
+
i
(
ˆ
e
L
2
j
+
1
2
ˆ
L
2
j
)
,
ˆ
F
+
j
=
ˆ
A
†
j
+
i
(
ˆ
e
L
†
1
j
−
1
2
ˆ
L
†
1
j
) +
i
2
ˆ
L
†
2
j
,
ˆ
e
F
−
j
=
ˆ
e
A
j
+
i
(
ˆ
L
1
j
+
1
2
ˆ
e
L
1
j
)
−
i
2
ˆ
e
L
2
j
,
ˆ
e
F
+
j
=
ˆ
e
A
†
j
+
i
2
ˆ
e
L
†
1
j
+
i
(
ˆ
L
†
2
j
−
1
2
ˆ
e
L
†
2
j
)
,
(S33)
with
ˆ
e
A
j
=
P
k
∈
I
1
∪
I
2
V
jk
ˆ
̃
d
k
,
ˆ
̃
L
1
j
= 2
P
k
1
∈
I
1
(
M
1
)
jk
1
ˆ
̃
d
k
1
and
ˆ
L
2
j
= 2
P
k
2
∈
I
2
(
M
2
)
jk
2
ˆ
̃
d
k
2
. Note that the
ˆ
F
operators
with superscripts
−
(+) are composed of a linear combination of
ˆ
d
k
and
ˆ
̃
d
k
(
ˆ
d
†
k
and
ˆ
̃
d
†
k
), respectively. The associated
superoperator BCFs are defined by,
C
>A
jj
′
(
t
1
,t
2
) =
⟨⟨
I
|
ˆ
F
−
j
e
ˆ
L
A
(
t
1
−
t
2
)
ˆ
F
+
j
′
e
ˆ
L
A
t
2
|
ρ
A
(0)
⟩⟩
,
(S34)
C
<A
jj
′
(
t
1
,t
2
) =
⟨⟨
I
|
ˆ
F
+
j
e
ˆ
L
A
(
t
1
−
t
2
)
ˆ
F
−
j
′
e
ˆ
L
A
t
2
|
ρ
A
(0)
⟩⟩
.
(S35)
Explicit calculation leads to the BCF expression in terms of a sum of complex exponential functions,
C
>A
(
t
+
t
′
,t
′
) = (
V
1
−
i
M
1
)
e
−
Z
1
t
(
V
1
+
i
M
1
)
†
,
C
<A
(
t
+
t
′
,t
′
) = (
V
2
−
i
M
2
)
∗
e
−
Z
2
t
(
V
2
+
i
M
2
)
T
,
(S36)
where
Z
1
and
Z
2
are diagonal matrices with elements
z
k
1
=
γ
k
1
+
iE
k
1
and
z
k
2
=
γ
k
2
−
iE
k
2
.
From
ˆ
L
A
|
ρ
A
(0)
⟩⟩
= 0,
e
ˆ
L
A
t
′
|
ρ
A
(0)
⟩⟩
=
|
ρ
A
(0)
⟩⟩
. We have,
ˆ
F
+
j
|
ρ
A
(0)
⟩⟩
= (
ˆ
A
†
j
−
i
2
ˆ
L
†
1
j
)
|
ρ
A
(0)
⟩⟩
=
X
k
1
∈
I
1
((
V
1
)
∗
jk
1
−
i
(
M
1
)
∗
jk
1
)
ˆ
d
†
k
1
|
ρ
A
(0)
⟩⟩
,
ˆ
F
−
j
|
ρ
A
(0)
⟩⟩
= (
ˆ
A
j
+
i
2
ˆ
L
2
j
)
|
ρ
A
(0)
⟩⟩
=
X
k
2
∈
I
2
((
V
2
)
jk
2
+
i
(
M
2
)
jk
2
)
ˆ
d
k
2
|
ρ
A
(0)
⟩⟩
.
(S37)
For
k
1
∈
I
1
,
ˆ
L
A
ˆ
d
†
k
1
|
ρ
A
(0)
⟩⟩
= (
−
iE
k
1
−
γ
k
1
)
ˆ
d
†
k
1
|
ρ
A
(0)
⟩⟩
, and for
k
2
∈
I
2
,
ˆ
L
A
ˆ
d
k
2
|
ρ
A
(0)
⟩⟩
= (
iE
k
2
−
γ
k
2
)
ˆ
d
k
2
|
ρ
A
(0)
⟩⟩
.
Hence,
e
ˆ
L
A
t
ˆ
F
+
j
|
ρ
A
(0)
⟩⟩
=
X
k
1
∈
I
1
e
(
−
iE
k
1
−
γ
k
1
)
t
((
V
1
)
∗
jk
1
−
i
(
M
1
)
∗
jk
1
)
ˆ
d
†
k
1
|
ρ
A
(0)
⟩⟩
,
e
ˆ
L
A
t
ˆ
F
−
j
|
ρ
A
(0)
⟩⟩
=
X
k
2
∈
I
2
e
(
iE
k
2
−
γ
k
2
)
t
((
V
2
)
jk
2
+
i
(
M
2
)
jk
2
)
ˆ
d
k
2
|
ρ
A
(0)
⟩⟩
.
(S38)
From the conjugation rule in Eq. S29,
⟨⟨
I
|
ˆ
d
†
k
=
⟨⟨
I
|
ˆ
̃
d
†
k
and
⟨⟨
I
|
ˆ
d
k
=
−⟨⟨
I
|
ˆ
̃
d
k
. Therefore,
⟨⟨
I
|
ˆ
F
−
j
=
⟨⟨
I
|

ˆ
A
j
−
i
2
ˆ
L
1
j
−
i
2
ˆ
L
2
j

=
X
k
∈
I
1
∪
I
2
⟨⟨
I
|
(
V
jk
−
iM
jk
)
ˆ
d
k
,
⟨⟨
I
|
ˆ
F
+
j
=
⟨⟨
I
|

ˆ
A
†
j
+
i
2
ˆ
L
†
1
j
+
i
2
ˆ
L
†
2
j

=
X
k
∈
I
1
∪
I
2
⟨⟨
I
|
(
V
jk
−
iM
jk
)
∗
ˆ
d
†
k
.
(S39)
Combining Eq. S38 and Eq. S39, the BCFs are expressed as,
C
>A
(
t
+
t
′
,t
′
) = (
V
1
−
i
M
1
)
e
−
Z
1
t
(
V
1
+
i
M
1
)
†
,
C
<A
(
t
+
t
′
,t
′
) = (
V
2
−
i
M
2
)
∗
e
−
Z
2
t
(
V
2
+
i
M
2
)
T
,
(S40)
The equivalence condition with the BCF in the fermionic bath can be similarly constructed with the Dyson series
expansion as in the bosonic bath. For a rigorous construction of the fermionic influence superoperator, we refer to
Ref. [5].
6
SM-IV. OPTIMAL GAUGE CHOICE MINIMIZING VIOLATION OF CP CONDITION
In the single coupling limit with a diagonal
Z
, the relationship between the weight
w
and the coupling parameters
V
and
M
is established as (
V
−
iM
)(
V
∗
−
iM
∗
) =
w
. Here,
w,V,M
∈
C
. Once
w
is given, the coupling parameters
V
and
M
have a gauge choice parameterized as
V
−
iM
=
κ
√
w,V
+
iM
=
κ
∗−
1
√
w
∗
where
κ
∈
C
. We show
that the solution of
V
and
M
with minimum
|
M
|
is then given by
V
−
iM
=
√
w
with
V,M
∈
R
, (
κ
= 1), or
V
= Re(
√
w
)
, M
=
−
Im(
√
w
).
Proposition 1.
Given
w
∈
C
, for
∀
V,M
∈
C
satisfying the constraint, (
V
−
iM
)(
V
∗
−
iM
∗
) =
w
, and its real
solution,
V
r
,M
r
∈
R
,
V
r
= Re(
√
w
)
,M
r
=
−
Im(
√
w
), then
|
M
|≥|
M
r
|
.
Proof.
With a parameterization,
V
=
|
V
|
e
iθ
1
,
M
=
|
M
|
e
iθ
2
, the constraint becomes,
|
V
|
2
−|
M
|
2
−
2
i
|
V
||
M
|
cos(
θ
1
−
θ
2
) =
w
or
|
V
|
2
−|
M
|
2
= Re(
w
) and
|
V
||
M
|
cos(
θ
1
−
θ
2
) =
−
Im(
w
)
/
2. It leads to
|
M
|
p
Re(
w
) +
|
M
|
2
|
cos(
θ
1
−
θ
2
)
|
=
|
Im(
w
)
|
/
2
.
Since
|
M
|
p
Re(
w
) +
|
M
|
2
is an increasing function in
|
M
|
, we have minimum
|
M
|
when
|
cos(
θ
1
−
θ
2
)
|
= 1, which
corresponds to
θ
1
=
θ
2
or
θ
1
=
θ
2
±
π
. The real solution
V
r
,M
r
satisfies this condition on
θ
1
and
θ
2
, so
|
M
|≥|
M
r
|
.
SM-V. PROOFS OF THE STABILITY CONDITION IN THE NONINTERACTING FERMIONIC
IMPURITY MODEL
Recall that for a noninteracting fermionic impurity model, its one-particle reduced density matrix (1-RDM)
P
pq
(
t
) =
Tr[ˆ
c
†
q
ˆ
c
p
ˆ
ρ
SA
(
t
)] satisfies a continuous differential Lyapunov equation [10, 11]:
∂
t
P
=
XP
+
P X
†
+
Y
,
(S41)
where
X
=
−
i
H
S
−
i
V
1
−
M
1
−
i
V
2
−
M
2
−
i
V
†
1
−
M
†
1
−
Z
1
0
−
i
V
†
2
−
M
†
2
0
−
Z
2
,
Y
=
0
0
2
M
2
0
0
0
2
M
†
2
0
2
Γ
2
.
(S42)
Proposition 2.
Eq. S41 is asymptotically stable if and only if the real part of all eigenvalues of
X
is strictly negative.
Proof.
Let Vec(
P
) denote the vectorization of
P
. Then Eq. S41 can be rewritten as:
∂
t
Vec(
P
) = (
X
⊗
I
+
I
⊗
X
∗
) Vec(
P
) + Vec(
Y
)
.
(S43)
Therefore,
P
(
t
) is asymptotically stable if and only if all eigenvalues of
X
⊗
I
+
I
⊗
X
∗
have strictly negative real
parts. Note that (
λ,
u
) is an eigen-pair of
X
⊗
I
if and only if
u
=
u
1
⊗
u
∗
2
and
X
u
1
=
λu
1
, and similarly for
I
⊗
X
∗
.
Also, note that
X
⊗
I
and
I
⊗
X
∗
commutes. Therefore
e
λ
ij
=
λ
i
+
λ
∗
j
is the eigenvalue of
X
⊗
I
+
I
⊗
X
∗
where
λ
i
,λ
j
are the eigenvalues of
X
. Therefore, Eq. S41 is asymptotically stable if and only if all eigenvalues of
X
have
strictly negative real parts.
For simplicity, we now assume a single-site system, i.e.,
ˆ
H
S
=
h
ˆ
a
†
ˆ
a
where
h
∈
R
, and
X
is
X
=

−
ih
−
i
V
−
M
−
i
V
†
−
M
†
−
Z

,
(S44)
where we do not separate baths 1 and 2 since there is no difference in the expression for
X
. The matrices
V
and
M
are now row vectors of size 1
×
N
b
, where
N
b
is the total number of bath modes. When needed, we will interchangeably
refer to these matrices as vectors.
Proposition 3.
When
∥
M
∥
<
1
2
min
γ
k
for
∀
k
, there exists
m >
0, such that if
∥
V
∥
> m
, the dynamics Eq. S41
with
X
being Eq. S44 is asymptotically stable.
7
Proof.
With Proposition 2, Eq. S41 being asymptotically stable is equivalent to
∂
t
u
=
X
u
being asymptotically stable
for any initial vector
u
(0). To study the stability of
∂
t
u
=
X
u
, let us define a Lyapunov function
V
(
u
) =
1
2
u
†
Θ
u
, for
which the quadratic form
Θ
will be specified later. Then we have
̇
V
(
u
(
t
)) =
1
2
u
†
(
Θ
X
+
X
†
Θ
)
u,
Then, by the Lyapunov stability theorem, if
Θ
X
+
X
†
Θ
is negative definite, then
∂
t
u
=
X
u
is asymptotically stable.
Let us consider
Θ
=
1
−
i
ε
V
1+
|
V
|
2
i
ε
V
†
1+
|
V
|
2
I
!
.
(S45)
Here
ε
are parameters that will be determined later. We define
Q
(
V
) =
−
(
Θ
X
+
X
†
Θ
). Let
V
=
c
e
V
where
e
V
is a
unit vector, then we have
lim
c
→∞
Q
(
c
e
V
) =
−
lim
c
→∞
(
Θ
X
+
X
†
Θ
) =
2
ε
2
M
2
M
†
2
Γ
−
2
ε
e
V
†
e
V
!
:=
Q
∞
.
If we let
ε
=
1
2
min
γ
k
,
Γ
≥
2
ε
I
, and
Q
∞
≥
2
ε
2
M
2
M
†
4
ε
I
−
2
ε
e
V
†
e
V
!
≥

2
ε
2
M
2
M
†
2
ε
I

. Since
∥
M
∥
< ε
=
1
2
min
γ
k
,
then
Q
∞
>
0. Therefore, for any unit vector
e
V
0
, there exists a constant
c
0
>
0 depending on
e
V
0
such that if
c > c
0
,
then
Q
=
Q
(
c
e
V
0
)
>
0. Since both
Q
and
c
0
depend continuously on
e
V
, and the unit sphere is compact, there exists
a positive constant
m
such that when
c
=
∥
V
∥
> m
, we have
Q >
0 and the dynamics is asymptotically stable.
Proposition 3 implies that the dynamics, though violating the positivity condition, is stable with sufficiently large
system-bath Hamiltonian coupling
ˆ
H
SA
and small Lindblad coupling
D
SA
.
SM-VI. DETAILS OF FITTING PROCEDURES
A. ESPRIT algorithm for complex exponential function fitting with complex weights
In this section, we briefly describe the Estimation of Signal Parameters via Rotational Invariant Techniques (ES-
PRIT) algorithm [12] used for fitting the BCF as a sum of complex exponential functions with complex weights,
C
(
t
)
≈
P
k
w
k
e
−
z
k
t
, where
w
k
and
z
k
are complex-valued. The ESPRIT algorithm has been widely adopted in a
signal processing community [13–17], and has recently been applied to quantum computing applications [18, 19]. In
this work, we take the ESPRIT algorithm based on the following Hankel matrix [16]:
H
=
C
(
t
0
)
C
(
t
1
)
C
(
t
2
)
···
C
(
t
n
−
1
)
C
(
t
1
)
C
(
t
2
)
C
(
t
3
)
···
C
(
t
n
)
C
(
t
2
)
C
(
t
3
)
C
(
t
4
)
···
C
(
t
n
+1
)
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
C
(
t
n
−
1
)
C
(
t
n
)
C
(
t
n
+1
)
···
C
(
t
2
n
−
2
)
(S46)
or
H
ab
=
C
(
t
a
+
b
), where
t
a
=
aδt
after discretizing time with a timestep
δt
. The ESPRIT algorithm is summarized
in the following steps to obtain the complex exponent
z
k
.
•
Step 1: Perform the singular-value decomposition (SVD) of the matrix
H
:
H
=
U
Σ
V
, where the singular values
(diagonal entries of the diagonal matrix Σ) are ordered non-increasingly.
•
Step 2: let
U
1
=
U
[1 :
n
−
1
,
1 :
N
exp
],
U
2
=
U
[2 :
n,
1 :
N
exp
]. Calculate
W
=
U
+
1
U
2
, where
U
+
1
= (
U
†
1
U
1
)
−
1
U
†
1
is the Moore-Penrose pseudo inverse of
U
1
.
•
Step 3: Calculate the eigenvalue
λ
1
,
···
,λ
N
exp
of the matrix
W
. Let
z
k
=
−
log(
λ
k
)
/
∆
t
, where we select the
branch of log(
·
) such that Im(log(
λ
k
))
∈
[
−
π,π
).
8
After we get the complex exponent
z
k
, the complex weights
w
k
are obtained via the following least-squares fitting:
min
{
w
k
}
X
a
C
(
t
a
)
−
X
k
w
k
e
−
z
k
t
a
2
.
(S47)
In the multi-site case, the BCFs are given by a matrix-valued function,
C
jj
′
(
t
)
≈
P
k
(
W
k
)
jj
′
e
−
z
k
t
. We construct a
scalar quantity
C
(
t
) by summing over all the entries in
C
jj
′
(
t
), i.e.,
C
(
t
) =
P
jj
′
C
jj
′
(
t
), to estimate
z
k
as above.
Then, we obtain the matrix-valued complex weights, (
W
k
)
jj
′
, with the elementwise least-squares fitting with Eq. S47.
B. Complex exponential function fitting with positive weights
As a comparison to the complex exponential function fitting with complex weights, we performed the complex
exponential function fitting with positive weights (i.e.,
w
k
∈
R
+
), which is a BCF decomposition for the Lorentzian
pseudomode model, with numerical fitting using the L-BFGS-B algorithm. We minimized the following fitting error
function,
E
fit
=
δt
v
u
u
t
X
a
C
(
t
a
)
−
X
k
w
k
e
−
z
k
t
a
2
.
(S48)
During the optimization steps, we only used the gradient with respect to
z
k
and determined
w
k
from the non-negative
least-squares fitting following [20]. As an initial guess for the optimization problem, we followed the deterministic
procedure from Ref. [10, 21–23], which determines the imaginary part of
z
k
from the direct discretization of the
unitary bath and sets the real part of
z
k
to some finite values. We chose the imaginary part of
z
k
from a set of
Gauss-Legendre quadrature nodes on the real axis, and the real part of
z
k
as Re[
z
k
] = 0
.
1
×
W/N
exp
where
W
is a
width of the spectrum.
In the multi-site case, the fitting error function is modified as,
E
fit
=
δt
v
u
u
t
X
jj
′
X
a
C
jj
′
(
t
a
)
−
X
k
(
W
k
)
jj
′
e
−
z
k
t
a
2
,
(S49)
with the constraint that
W
k
≥
0.
SM-VII. DETAILS ON DMRG CALCULATIONS
We describe some details of the time-dependent density matrix renormalization group (tdDMRG) calculations,
focusing on the fermionic impurity models. We first express the fermionic Lindblad master equation in the super-
fermion representation. This allows us to represent the density operator as a state vector with Schr ̈odinger equation-
like propagation so that the two-site time-dependent variational principle, implemented in the Block2 [24, 25] package,
can be directly used for the state vector propagation. Furthermore, the Liouville operator has a number-conserving
form in the super-fermion representation, so we can take advantage of
U
(1) symmetry adaptation in the tdDMRG
calculations.
Ordering of orbitals is important in DMRG calculations. We employed two different ordering schemes illustrated
in Fig. S1. The illustration is based on the spinless single-site impurity model with two bath orbitals in each bath
A
1
and
A
2
. We have ˆ
ρ
SA
(0) = ˆ
ρ
S
(0)
⊗
ˆ
ρ
A
1
(0)
⊗
ˆ
ρ
A
2
(0), where ˆ
ρ
A
1
(0) =
|
00
⟩⟨
00
|
and ˆ
ρ
A
2
(0) =
|
11
⟩⟨
11
|
. In the super-
fermion representation,
|
ρ
SA
(0)
⟩⟩
=
|
ρ
S
(0)
⟩⟩⊗|
ρ
A
1
(0)
⟩⟩⊗|
ρ
A
2
(0)
⟩⟩
where
|
ρ
A
1
(0)
⟩⟩
=
|
0101
⟩⟩
and
|
ρ
A
2
(0)
⟩⟩
=
|
1010
⟩⟩
.
The orbitals from the ‘tilde’ space are drawn with hatching lines in Fig. S1. The two orbital ordering schemes are
distinguished by whether the bath
A
1
and
A
2
are on the same side (Fig. S1a) or different sides (Fig. S1b) of the
impurity. We chose the first scheme for the spinful impurity model to put two different spins on different sides of the
impurity [26]. We chose the second scheme for the spinless impurity model. The remaining ordering of the orbitals
within each bath is determined from the imaginary part of
z
k
, or equivalently, the energy part,
E
k
, following the
physical insights in [27, 28].
9
impurity
Bath
A
1
Bath
A
2
|0
|1
|1
|0
|0
|1
|1
|0
(a)
impurity
Bath
A
2
Bath
A
1
|0
|1
|1
|0
|0
|1
|1
|0
(b)
FIG. S1. Two different orbital ordering schemes for tdDMRG calculations.
[S1] D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Nonperturbative treatment of non-markovian dynamics of
open quantum systems, Phys. Rev. Lett.
120
, 030402 (2018).
[S2] A. Rivas and S. F. Huelga,
Open quantum systems
(Springer, 2012).
[S3] H.-P. Breuer and F. Petruccione,
The Theory of Open Quantum Systems
(Oxford University Press, 2007).
[S4] E. Aurell, R. Kawai, and K. Goyal, An operator derivation of the feynman-vernon theory, with applications to the
generating function of bath energy changes and to an-harmonic baths, Journal of Physics A: Mathematical and Theoretical
53
, 275303 (2020).
[S5] M. Cirio, P.-C. Kuo, Y.-N. Chen, F. Nori, and N. Lambert, Canonical derivation of the fermionic influence superoperator,
Phys. Rev. B
105
, 035121 (2022).
[S6] M. Schmutz, Real-time green’s functions in many body problems, Zeitschrift f ̈ur Physik B Condensed Matter
30
, 97
(1978).
[S7] A. A. Dzhioev and D. S. Kosov, Super-fermion representation of quantum kinetic equations for the electron transport
problem, The Journal of Chemical Physics
134
, 10.1063/1.3548065 (2011), 044121.
[S8] A. Dorda, M. Nuss, W. von der Linden, and E. Arrigoni, Auxiliary master equation approach to nonequilibrium correlated
impurities, Phys. Rev. B
89
, 165105 (2014).
[S9] G. Park, N. Ng, D. R. Reichman, and G. K.-L. Chan, Tensor network influence functionals in the continuous-time limit:
Connections to quantum embedding, bath discretization, and higher-order time propagation, Phys. Rev. B
110
, 045104
(2024).
[S10] M. Lotem, A. Weichselbaum, J. von Delft, and M. Goldstein, Renormalized lindblad driving: A numerically exact
nonequilibrium quantum impurity solver, Phys. Rev. Res.
2
, 043052 (2020).
[S11] T. Barthel and Y. Zhang, Solving quasi-free and quadratic lindblad master equations for open fermionic and bosonic
systems, Journal of Statistical Mechanics: Theory and Experiment
2022
, 113101 (2022).
[S12] A. Paulraj, R. Roy, and T. Kailath, A subspace rotation approach to signal parameter estimation, Proceedings of the
IEEE
74
, 1044 (1986).
[S13] R. Roy and T. Kailath, Esprit-estimation of signal parameters via rotational invariance techniques, IEEE Transactions
on Acoustics, Speech, and Signal Processing
37
, 984 (1989).
[S14] S. Rouquette and M. Najim, Estimation of frequencies and damping factors by two-dimensional esprit type methods,
IEEE Transactions on Signal Processing
49
, 237 (2001).
[S15] S. Shahbazpanahi, S. Valaee, and M. Bastani, Distributed source localization using esprit algorithm, IEEE Transactions
on Signal Processing
49
, 2169 (2001).
[S16] W. Li, W. Liao, and A. Fannjiang, Super-resolution limit of the esprit algorithm, IEEE Transactions on Information
Theory
66
, 4593 (2020).
[S17] Z. Ding, E. N. Epperly, L. Lin, and R. Zhang, The ESPRIT algorithm under high noise: Optimal error scaling and noisy
super-resolution, in
2024 Symposium on Foundations of Computer Science (FOCS 2024)
(2024) arXiv:2404.03885.
[S18] M. E. Stroeks, J. Helsen, and B. M. Terhal, Spectral estimation for hamiltonians: a comparison between classical
imaginary-time evolution and quantum real-time evolution, New Journal of Physics
24
, 103024 (2022).
[S19] H. Li, H. Ni, and L. Ying, Adaptive low-depth quantum algorithms for robust multiple-phase estimation, Phys. Rev. A
108
, 062408 (2023).
[S20] Z. Huang, E. Gull, and L. Lin, Robust analytic continuation of green’s functions via projection, pole estimation, and
semidefinite relaxation, Phys. Rev. B
107
, 075151 (2023).
[S21] M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T. Mitchison, S. R. Clark, and J. Goold, Tensor-network method