Quasi-Lindblad pseudomode theory for open quantum systems

PHYSICAL REVIEW B

110

, 195148 (2024)

Gunhee Park (

)

Zhen Huang

Yuanran Zhu,

Chao Yang,

Garnet Kin-Lic Chan

and Lin Lin

Division of Engineering and Applied Science,

California Institute of Technology

, Pasadena, California 91125, USA

Department of Mathematics,

University of California, Berkeley

, California 94720, USA

Applied Mathematics and Computational Research Division,

Lawrence Berkeley National Laboratory

, Berkeley, California 94720, USA

Division of Chemistry and Chemical Engineering,

California Institute of Technology

, Pasadena, California 91125, USA

(Received 29 August 2024; accepted 17 October 2024; published 25 November 2024)

We introduce a new framework to study the dynamics of open quantum systems with linearly coupled

Gaussian baths. Our approach replaces the continuous bath with an auxiliary discrete set of pseudomodes

with dissipative dynamics, but we further relax the complete positivity requirement in the Lindblad master

equation and formulate a quasi-Lindblad pseudomode theory. We show that this quasi-Lindblad pseudomode

formulation directly leads to a representation of the bath correlation function in terms of a complex weighted

sum of complex exponentials, an expansion that is known to be rapidly convergent in practice and thus leads to

a compact set of pseudomodes. The pseudomode representation is not unique and can differ by a gauge choice.

When the global dynamics can be simulated exactly, the system dynamics is unique and independent of the

specific pseudomode representation. However, the gauge choice may affect the stability of the global dynamics,

and we provide an analysis of why and when the global dynamics can retain stability despite losing positivity.

We showcase the performance of this formulation across various spectral densities in both bosonic and fermionic

problems, finding significant improvements over conventional pseudomode formulations.

DOI:

10.1103/PhysRevB.110.195148

I. INTRODUCTION

Open quantum systems (OQS) are central to many fields,

including quantum optics, condensed matter physics, chem-

ical physics, and quantum information science [

–

]. The

most widely studied case corresponds to a system coupled

to a bath with continuous (bosonic or fermionic) degrees of

freedom with the bath dynamics governed by a Hamiltonian

quadratic in the bath field operators (a Gaussian bath) and

with a system-bath coupling linear in the bath field operators.

The primary task in such setups is to compute the system’s

dynamical observables.

There are two main classes of approaches to simulate this

problem. The first approach integrates out the Gaussian bath

dynamics to yield a path integral for the reduced system

dynamics [

], which can then be evaluated using various

approximations [

–

]. The second approach, which is the

focus of this work, replaces the continuous physical bath with

a discrete mode representation [

–

], i.e., an auxiliary bath,

allowing the global state of the system and auxiliary bath to be

computationally propagated in time. In both approaches, the

influence of the bath on the reduced system dynamics is en-

tirely determined by the bath correlation function (BCF) [

As long as the auxiliary BCF matches the original BCF, the

reduced system dynamics can be shown to be identical [

Therefore the primary consideration in this second approach

is to begin with a compact auxiliary bath representation of the

BCF.

The conventional discrete mode representation is to di-

rectly discretize the continuous Hamiltonian [

–

] without

introducing any dissipation. However, this discretization re-

quires a large number of modes in the long-time limit, which

scales linearly with the propagated time [

]. In most physical

cases, the continuous bath provides dissipation to the system.

It is thus natural to introduce an ansatz where the auxiliary

bath dynamics is dissipative, in which case, the bath modes are

referred to as pseudomodes [

–

]. The Liouvillian

for the system and auxiliary bath can then be chosen to satisfy

the constraint of physical dynamics, taking a Lindblad form

and generating a completely positive trace-preserving (CPTP)

dynamics. The corresponding BCF (assuming independent

pseudomodes) is then a sum of

complex

exponential functions

weighted by

positive

numbers. This BCF exactly represents

the continuous spectral density as a positive weighted sum

of Lorentzian functions. Unfortunately, many studies have

found that this natural Lorentzian pseudomode expansion is

inefficient, achieving only a slow rate of convergence as the

number of pseudomodes increases [

]. Improving this

convergence rate remains an active area of research, with

approaches including the introduction of additional intra-bath

Hamiltonian terms [

–

] or non-Hermitian pseudomode

formulations [

–

In the current work, we demonstrate that by relaxing the

complete positivity (CP) of the global dynamics of the system

and pseudomodes, we gain the flexibility to use a broader

set of functions to approximate the BCF, significantly en-

hancing the convergence rate of the pseudomode expansion.

The generator of the dynamics is a slight modification of

the Lindblad form, which we refer to as a quasi-Lindblad

pseudomode representation. The BCF is expressed as a sum

complex

exponentials weighted by

complex

numbers. This

ansatz is well-known in the signal processing and applied

mathematics communities for yielding an exponentially fast

convergence rate for many functions [

–

], a result we also

2469-9950/2024/110(19)/195148(12)

195148-1

PARK, HUANG, ZHU, YANG, CHAN, AND LIN

PHYSICAL REVIEW B

110

, 195148 (2024)

observe in our numerical tests. Additionally, this ansatz for

the BCF has been adopted in non-pseudo-mode approaches

to OQS dynamics, such as hierarchical equations of motion

(HEOM) [

–

] and tensor network influence functionals

[

], providing a connection among these various methods.

It is worth noting that the BCF alone does not uniquely

determine the pseudomode formulation. Meanwhile, the sys-

tem dynamics is uniquely determined by the BCF and is

independent of the pseudomode formulation as long as the

global dynamics can be simulated exactly (i.e., without further

numerical approximations such as truncating the number of

bosonic modes or even finite precision arithmetic operations).

However, certain pseudomode formulations can be unstable

if the corresponding Liouvillian has eigenvalues with posi-

tive real parts, leading to exponentially growing modes. This

phenomenon has also been observed in other contexts, such

as HEOM [

–

]. Therefore, in practice, the stability of

the pseudomode formulation can indeed impact the quality

of system dynamics. Nevertheless, with a careful choice of

the pseudomode formulation that minimally violates the CP

condition, we observe that the global dynamics can be stable

across a range of physically relevant OQS in our simulations.

We attribute this stability to the mixing of the divergent and

dissipative parts of the global dynamics via the coherent terms

in the Liouvillian.

II. QUASI-LINDBLAD PSEUDOMODE CONSTRUCTION

A. Problem setup: Original unitary OQS

We consider a quantum system

linearly coupled to a con-

tinuous Gaussian bath

, either bosonic or fermionic, evolving

under the following Hamiltonian,

phy

ωω

†

(1)

where

denotes the system-only Hamiltonian and ̄

1. ˆ

†

denotes the creation operator for the bath mode with energy

In the system-bath coupling term,

and

are operators

within the system and bath, respectively, and can be assumed

to be Hermitian operators, without loss of generality [

The bath coupling operator is given by

(

)ˆ

∗

(

)ˆ

†

. We assume the initial state to take a factorized form,

(0)

⊗

(0), where ˆ

(0) is a Gaussian state. We

will also assume that ˆ

(0) is a number-conserving operator

(and also for all initial bath reduced density operators hence-

forth) for simplicity.

Time-evolution of the density operator is described by

the von Neumann equation,

∂

=−

[

phy

], and the

system-reduced density operator is obtained by tracing out the

bath, ˆ

(

)

[ˆ

(

)]. The influence of the Gaussian bath

on the system-reduced density operator is determined by the

bath correlation function (BCF),

(

)

[

(

)

(

)ˆ

(0)]

(2)

for

0 with

(

)

−

. When ˆ

(0) is sta-

tionary, i.e.,

−

(0)

(0), the BCF only depends

on the time difference

−

, i.e.,

(

)

(

−

B. Conventional pseudomode

In many cases of interest, the above unitary dynamics using

a continuous bath leads to a dissipative bath correlation func-

tion, i.e., a decaying

(

−

) with lim

−

→∞

(

−

)

→

0. To reproduce the associated dynamics of ˆ

(

) with a dis-

crete bath, we now introduce a pseudomode description. First,

we review the conventional pseudomode formulation [

which introduces an auxiliary bosonic bath

with a Lindblad

dissipator

applied within the bath

. The density operator

across the system and auxiliary bath, ˆ

, evolves under the

following Lindblad master equation,

∂

=−

[

aux

]

(3)

where the Hamiltonian

aux

is composed

of the system (bath)-only Hamiltonian

(

) and coupling

Hamiltonian

has a Lindblad form

•=

•

†

−

{

†

•}

(4)

which guarantees the CPTP property of the global dynamics

consisting of the system and the pseudomodes. Here,

can

be an arbitrary operator within the bath

. The Lindblad

dissipator can be expressed in terms of a pseudomode basis,

•=

•

†

−

{

†

•}

(5)

where

is an operator within the

th pseudomode, and

is a positive semidefinite (PSD) matrix from the CPTP condi-

tion.

and

are quadratic in the creation and annihilation

operators of the pseudomode, and

and

are linear, to

maintain the Gaussian properties of the bath.

has negative

eigenvalues, and thus, we can interpret it as adding dissipation

to the reduced system dynamics.

The BCF for the pseudomode can be defined using

the bath-only Liouvillian

•=−

[

•

]

•

and initial

pseudomode reduced density operator ˆ

(0),

(

)

(

−

)

(0)

(6)

Reference [

] established the equivalence condition between

the unitary bath dynamics and the pseudomode dynamics in

terms of the BCF: if

(

)

(

)for

∀

then Tr

[ˆ

(

)]

[ˆ

(

)] [

]. As in the unitary case in

Sec.

II A

, when ˆ

(0) is stationary, i.e.,

(0)

(0),

the BCF depends solely on

−

C. Quasi-Lindblad pseudomode

As described in the introduction, the conventional pseudo-

mode description can reproduce a dissipative bath correlation

function within a finite discretization, but it typically requires

a large number of pseudomodes for an accurate approximation

to the BCF. To remedy this, we now provide a more general

pseudomode description that introduces an additional system-

bath coupling

, expressed in terms of Lindblad dissipators,

•=

•

†

−

{

†

•}

(7)

195148-2

QUASI-LINDBLAD PSEUDOMODE THEORY FOR OPEN ...

PHYSICAL REVIEW B

110

, 195148 (2024)

FIG. 1. Pseudomode theory defines a discrete auxiliary dissi-

pative bath

to reproduce a system-reduced density operator.

(a) The conventional pseudomode assumes bath-only dissipators,

. (b) Our quasi-Lindblad pseudomode relaxes the complete pos-

itivity condition by introducing an additional system-bath Lindblad

coupling with dissipators,

where

. We call this coupling a Lindblad

coupling, as illustrated in Fig.

1(b)

. The total dissipator

can be compactly written as

•=

•

†

−

{

†

•}

(8)

where the index

and

refer to both the coupling index

and the pseudomode bath index

, i.e.,

∈{

}∪{

}

, and if

. The coefficient

†

(9)

where the bold font denotes matrices. We note that the matrix

is not PSD in general, unlike the PSD matrix

, and

hence, this pseudomode equation of motion violates the CP

condition in the Lindblad master equation. Therefore we call

it a quasi-Lindblad pseudomode representation.

To take into account both the Hamiltonian and Lindblad

coupling in the BCF, we define the system-bath Liouvillian

•=−

[

•

]

•

and further decompose it in the

following form:

=−

(10)

where the system superoperators,

and

are defined as

•=

•

and

•=•

. The superoperators,

and

,are

obtained as

•=

•+•

†

−

(

†

)

•

•=•

−

•+•

(

†

)

(11)

The BCF is defined with the superoperators

(

)

(

−

)

(0)

(12)

We remark that

reduces to

when

0, and then

the BCF in Eq. (

) also reduces to the BCF in Eq. (

Similar to the pseudomode description in Sec.

II B

, the above

BCF leads to the equivalence of the reduced system dynamics

to that generated by the original unitary formulation under

the equivalence of the BCFs. We provide the derivation of

the equivalence condition in the Supplemental Material (SM)

[

In this framework, violation of the CP condition leads the

dissipator Liouvillian

to have eigenvalues with positive real

parts. Thus, in the absence of any coherent term (i.e., the com-

mutator term of the form

−

[

aux

•

] for some Hamiltonian

aux

), this would lead to unstable global dynamics with modes

growing exponentially in time. Note that when simulations are

performed exactly, the reduced system dynamics, after tracing

out the bath, can still be perfectly reproduced from the BCF

equivalence condition. However, this cannot be guaranteed

in practical simulations involving finite precision arithmetic

operations and physical approximations, such as truncating

the number of bosonic modes. Therefore the stability of the

Liouvillian is an important issue to study.

Interestingly, we find that the Liouvillian, including both

the coherent and dissipator terms,

can

generate stable dynam-

ics despite violating the CP condition. We will analyze this

effect in more detail in Sec.

D. Ansatz for quasi-Lindblad pseudomode

We next introduce an ansatz for the quasi-Lindblad pseu-

domode that yields a simple form for the BCF, facilitating

numerical fitting. We describe this ansatz for both bosonic

and fermionic baths. We mainly target the reproduction of the

BCF of the unitary dynamics in Sec.

II A

with the initial bath

being a stationary state, such as a thermal state ˆ

(0)

∝

−

where

represents the inverse temperature.

We start with a bosonic bath, setting the initial reduced

density operator to a vacuum state, ˆ

(0)

. The dissi-

pator

is set as

, where

is the bosonic annihilation

operator of the

th pseudomode, and the Hamiltonian

†

. This ansatz satisfies the stationary condition

(0)

0. Then,

(0)

(0) and

(

) is only

dependent on the time difference

(

)

(

−

)

(

), where

−

. The system-bath coupling is set

∗

†

We obtain the BCF based on the above ansatz:

(

)

(

−

)

−

(

)

†

(13)

where

. This BCF expression has a gauge redun-

dancy:

→

GZG

−

→

(

−

)

−

→

(

)

†

, for arbitrary invertible matrix

. Therefore we

focus on a diagonal matrix

diag(

{

}

), which can be

transformed to a more general diagonalizable

through a

unitary transformation [

]. In this case, we define

†

and

, and

. Even with

a diagonal matrix

and

still have a gauge redundancy

associated with a diagonal

, since

remains invariant under

the transformation

GZG

−

. It is important to note that

while different choices of gauge yield the same BCF and

thus the same reduced system dynamics (in the absence of

simulation errors), the total system-bath Liouvillian differs,

leading to different global system-bath dynamics.

The case of a fermionic bath can be introduced simi-

larly. In particular, we focus on number-conserving fermionic

impurity models where the coupling Hamiltonian is given

†

(

)ˆ

(ˆ

†

is the

195148-3

PARK, HUANG, ZHU, YANG, CHAN, AND LIN

PHYSICAL REVIEW B

110

, 195148 (2024)

creation operator of the impurity with index

, which can

involve both spin and orbital degrees of freedom). The

number-conserving property of the fermionic impurity model

reduces the BCFs to two different BCFs, namely, the greater

and lesser BCFs:

and

. For any fermionic Gaussian

state ˆ

(0)

and

are defined as

(

)

[

(

)

†

(

)ˆ

(0)]

(

)

[

†

(

)

(

)ˆ

(0)]

(14)

Corresponding to these two BCFs, we use two different aux-

iliary baths with an initial reduced density operator ˆ

(0)

⊗

(0)

|⊗|

. The bath-only Hamiltonian

and dissipator are set to be diagonal without loss of gen-

erality, as discussed in the bosonic case. The Hamiltonian

†

where

and

refer

to the pseudomode index in each bath, respectively, and the

dissipator is set to be

•=

•

†

−

†

•

•=

†

•

−

†

•

(15)

which satisfies

(0)

0 for both

2. The system-

bath Hamiltonian coupling is

†

(

)

(

)

and the system-bath Lind-

blad coupling is set to be

•=

•

†

•

†

−

{

†

•}

•=

†

•

†

•

−

{

†

•}

(16)

where

)

and

)

The BCF with this ansatz becomes

(

)

(

−

)

−

(

)

†

(

)

(

−

)

∗

−

(

)

(17)

where

and

are diagonal matrices with elements

and

−

. We note that the baths

and

only contribute to the BCF

and

, respectively. We

remark that the BCF from the fermionic bath is also expressed

as a complex weighted sum of complex exponential functions

and has the same gauge redundancy as the bosonic bath.

It is instructive to consider a model with a single coupling

term (with only

1). In this case, the BCF simplifies to a

scalar function without indices:

(

)

w

−

(18)

where

w

(

−

)(

∗

−

∗

). This expression is a sum

complex

exponential functions with

complex

weights. In

contrast, the conventional pseudomode description corre-

sponds to

0, resulting in

positive

weights

w

Physically, this exactly represents the BCF from a positive

weighted sum of the Lorentzian spectrum, and thus, we also

refer to the conventional pseudomode as a Lorentzian pseudo-

mode. We see, however, that the quasi-Lindblad pseudomode

provides for a more flexible and general representation, whose

influence on the compactness of the pseudomode description

we will examine in the numerics below.

In addition to the enhanced flexibility, the use of complex

exponential expressions also simplifies the numerical fitting

of the BCF. In this work, we utilize the ESPRIT algorithm

[

], which has previously been reported as one of the

most competitive methods for BCF fitting [

], to determine

the complex exponents,

. We then perform a least-squares

fitting of the complex weights

w

. We note that this simple un-

constrained least-squares approach is possible because we do

not restrict

w

0 as would be required in the conventional

pseudomode approach. Once

w

is determined, we choose the

gauge for

and

. In this limit, the gauge choice is param-

eterized as

−

√

w

∗−

w

∗

where

∈

. Here, we follow a simple heuristic of minimizing the

magnitude of

, as this term determines the extent of CP

condition violation. We describe the result for the bosonic bath

expression in Eq. (

)for

and

, and its adaptation to the

fermionic BCF expression is straightforward. The solution of

and

with minimum

is then given by

−

√

w

with

∈

1) which we provide its proof in

Ref. [

]. We will examine other choices of

and

their effect on the stability of the system-bath dynamics in

Sec.

There are several approaches to relax the CPTP condi-

tion in pseudomode representations. One approach involves

a non-Hermitian pseudomode formulation [

–

], where the

system-bath coupling is given by a non-Hermitian Hamilto-

nian (without Lindblad coupling). This formulation results in

real weights,

w

∈

[

], or complex weights from a pair

of pseudomodes [

]. Another approach, presented in [

is a Lindblad-like pseudomode formulation derived from the

HEOM with a similarity transformation. While this approach

allows for complex weights, it is important to emphasize

that the Lindblad-like equation in [

] represents one gauge

choice for

[

] within our quasi-Lindblad pseudomode

framework.

For general couplings with multiple

, a similar

strategy can be adopted. In this case, we first find a sum

of exponential functions with a matrix-valued weight

(

)

≈

−

. We determine

by fitting a

scalar quantity, such as the sum of the BCF entries,

(

)

≈

(

)

−

, or the trace

(

)

≈

(

)

−

. The matrix-valued

complex weights

are then obtained through the

least-squares fitting.

Finding

and

from

involves additional matrix de-

compositions. However, the representation from Eq. (

) has a

rank-1 matrix factorization form, (

)

(

−

)

(

)

∗

, and

from the least-squares fitting is generally not

a rank-1 matrix, but of rank

, after denoting the range of

. To address this, we represent the rank-

indexing the pseudomode with

(

), where the index

is an additional index with a range of

. Then, we assign

the diagonal elements of

(

)

. This gives the

195148-4

QUASI-LINDBLAD PSEUDOMODE THEORY FOR OPEN ...

PHYSICAL REVIEW B

110

, 195148 (2024)

FIG. 2. Exponential function fitting of BCF for [(a) and (d)] subohmic, [(b) and (e)] semicircular, and [(c) and (f)] two-site impurity spectral

densities. The real part of the fitted BCF is shown in (a)–(c) from the complex-weighted fit to complex exponentials (red solid, quasi-Lindblad

pseudomode formulation in this work), positive-weighted fit to complex exponentials (blue dash-dotted, conventional Lorentzian pseudomode

formulation), and positive-weighted fit to real exponentials (green dotted, discrete Hamiltonian) compared to the exact BCF (black dashed).

The inset shows the zoomed-in data at

1 (only the results from complex exponential fits). The number of exponential functions,

exp

,used

for the fitting was

exp

4 and 6 for (a) and (b), (c), respectively. [(d)–(f)] Fitting errors,

fit

,asafunctionof

exp

from the fitting with complex

weights (red dots), positive weights (blue stars), and discrete Hamiltonian (green crosses).

following expression:

(

)

(

)

(

−

)

(

)

∗

−

(19)

where (

)

(

−

)

(

)

∗

. Thus, given

exp

exponents

and

coupling operators,

exp

pseu-

domodes are used in the BCF fitting. With this setting, we

decompose

with a singular value decomposition (SVD),

†

, where

is a diagonal matrix. One possi-

ble choice for

and

is (

−

)

(

)

(

)

(

)

, which reduces to the gauge choice

above for the single coupling limit [

III. NUMERICAL RESULTS

A. Fitting the bath correlation function

The quasi-Lindblad pseudomode formulation provides the

BCF as a complex weighted sum of complex exponential

functions. This type of BCF representation is widely used,

for example, in HEOM [

–

] and tensor network influ-

ence functionals [

]. It is known in many applications

to show a fast exponential convergence rate with the number

of pseudomodes. In this section, we focus on comparing the

performance of BCF fitting using the quasi-Lindblad ansatz to

that achieved with the conventional Lorentzian pseudomode

ansatz and a discrete Hamiltonian ansatz obtained from direct

discretization.

Given a spectral density

(

)

(

)

∗

(

), we will fit

the associated zero-temperature BCF [

(

)

(

)

−

ω,

(20)

using the three different ansatz. We examine the performance

of these ansatz for three spectral densities, as illustrated in

Fig.

. The first case is a subohmic spectral density,

(

)

−

ω/ω

,ω

∈

∞

)

(21)

the second, a semicircular spectral density,

(

)

−

,ω

∈

]

(22)

and the third corresponds to a two-site spectral density,

(

)

(

)

∗

(

)

(23)

where

(

)

exp(

−

ω/

) and

(

) is the semicircular

spectral density specified above. For the subohmic case, we

use

5 and

1, and for the semicircular case,

we set

1 and

Figures

2(a)

–

2(c)

compare the fitted BCF to the exact

expressions shown by the black dashed lines. Three different

fits are compared: (1) complex exponential function fits with

complex weights (red solid), (2) with positive (or PSD for the

multisite case) weights (blue dash-dotted), and (3) a discrete

Hamiltonian (green dotted). Fitting scheme (2) represents the

Lorentzian pseudomode. Here, the parameters are optimized

195148-5

PARK, HUANG, ZHU, YANG, CHAN, AND LIN

PHYSICAL REVIEW B

110

, 195148 (2024)

FIG. 3. Time evolution of coherence

from the quench dy-

namics of the spin-boson model. The analytical result from the

continuous bath with

5 subohmic spectral density is compared

to results from a two-mode discretized bath using quasi-Lindblad

pseudomodes, Lorentzian pseudomodes, and a discrete Hamiltonian.

(Inset) Errors of coherence compared to the analytical result.

through numerical minimization of the fitting error. For fit-

ting (3), we used a discretization scheme based on Gaussian

quadrature using Legendre polynomials [

]. Details of these

fitting schemes are in Ref. [

]. The fits to the quasi-Lindblad

ansatz with

exp

6) for Fig.

2(a)

[Figs.

2(b)

and

2(c)

], respectively, show high accuracy already with a small

number of exponentials. In Figs.

2(d)

–

2(f)

, the maximum fit-

ting error as a function of

exp

is illustrated. This comparison

shows the clear advantage of the quasi-Lindblad pseudomode

compared to the Lorentzian pseudomode and direct discretiza-

tion schemes and illustrates an exponential decrease of error

with respect to

exp

in both cases.

B. Quench dynamics simulation

We then test the accuracy of the quasi-Lindblad pseudo-

mode formulation in simulating the quench dynamics of a

spin-boson model and a fermionic impurity model. First, we

demonstrate the effectiveness of this approach in the spin-

boson model. We choose the Hamiltonian

2 with

4 and

and the initial system-reduced density

operator ˆ

(0)

=|+ +|

, with

|+ =

√

). While

the analytical solution for this setup is known [

], obtain-

ing accurate numerical results remains nontrivial [

particularly due to the requirement for a precise description of

the BCF.

We choose the

5 subohmic spectral density as above

(

1) and consider a zero temperature bath initial

condition. The BCF is fitted using two bosonic modes, and we

compare the quasi-Lindblad pseudomode, Lorentzian pseudo-

mode, and a discrete Hamiltonian from Gaussian quadrature.

Figure

shows the time evolution of the coherence

(

)

using these three different bath discretizations. We see that the

time evolution from the discrete Hamiltonian shows unphys-

ical oscillations in the long-time limit, while the Lorentzian

FIG. 4. Quench dynamics of spinless (a) single-site and (b) two-

site impurity models with

with

=−

, (c) spinful

SIAM and (d) spinless two-site impurity model with

with

=−

. [(a) and (b)] Time-dependence of

impurity occupations,

(

)

, computed from Gaussian calculations

from our quasi-Lindblad pseudomode, Lorentzian pseudomode, dis-

crete Hamiltonian, and exact unitary dynamics with

exp

4. (Inset)

A maximum error of

(

)

from pseudomode dynamics as a func-

tion of

exp

for three different fitting schemes in Fig.

. [(c) and

(d)]

(

)

of our pseudomode compared to the original unitary

dynamics with large bath discretization, computed from tdDMRG.

pseudomode qualitatively captures the dissipative dynamics

but at a much lower accuracy than the quasi-Lindblad pseudo-

mode representation.

We next study the dynamics of the spinless noninteract-

ing single-site and two-site models with

(ˆ

†

). Here, we can carry out efficient Gaussian calculations,

which we do for both a continuous unitary bath and for the

Lindblad pseudomode dynamics [

]. We choose the semi-

circular spectral density in Eqs. (

) and (

) with

but with

∈

[

−

]

,ε

=−

, and consider the initial bath

as a fully occupied (empty) state for

ω<

ω>

0). We

fit the two BCFs,

(

) and

(

), with

exp

exponential

functions each and represent the bath with in total 2

exp

) pseudomodes for the single-site (two-site) impurity

model, respectively. The initial impurity state is assumed to

be an empty state, ˆ

(0)

, for the single-site model,

and ˆ

(0)

for the two-site model. Figures

4(a)

and

4(b)

show the time-dependence of the impurity occupation,

(

)

, with

exp

4 for the exact dynamics, the quasi-

Lindblad pseudomode, Lorentzian pseudomode, and discrete

Hamiltonian dynamics. These results clearly show the im-

proved accuracy of the quasi-Lindblad pseudomode relative to

the Lorentzian pseudomode, which reaches the wrong steady

state. The inset in Figs.

4(a)

and

4(b)

shows the maximum

error of the impurity occupation as a function of

exp

, which

also shows an exponential convergence in

exp

195148-6