Quasi-Lindblad pseudomode theory for open quantum systems

We consider a quantum system $S$ linearly coupled to a continuous Gaussian bath $B$, either bosonic or fermionic, evolving under the following Hamiltonian,

	$\displaystyle\hat{H}_{\text{phy}}$	$\displaystyle=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}$
		$\displaystyle=\hat{H}_{S}+\int\!d\omega\ \omega\ \hat{c}_{\omega}^{\dagger}\hat{c}_{\omega}+\sum_{j}\hat{S}_{j}\hat{B}_{j},$		(1)

where $\hat{H}_{S}$ denotes the system-only Hamiltonian and $\hbar=1$. $\hat{c}_{\omega}^{\dagger}$ denotes the creation operator for the bath mode with energy $\omega$. In the system-bath coupling term, $\hat{H}_{SB}$, $\hat{S}_{j}$ and $\hat{B}_{j}$ are operators within the system and bath, respectively, and can be assumed to be Hermitian operators, without loss of generality [1]. The bath coupling operator is given by $\hat{B}_{j}=\int d\omega\ g_{j}(\omega)\hat{c}_{\omega}+g_{j}^{*}(\omega)\hat{c}_{\omega}^{\dagger}$. We assume the initial state to take a factorized form, $\hat{\rho}_{SB}(0)=\hat{\rho}_{S}(0)\otimes\hat{\rho}_{B}(0)$, where $\hat{\rho}_{B}(0)$ is a Gaussian state. We will also assume that $\hat{\rho}_{B}(0)$ is a number-conserving operator (and also for all initial bath reduced density operators henceforth) for simplicity.

Time-evolution of the density operator is described by the von Neumann equation, $\partial_{t}\hat{\rho}_{SB}=-i[\hat{H}_{\text{phy}},\hat{\rho}_{SB}]$, and the system-reduced density operator is obtained by tracing out the bath, $\hat{\rho}_{S}(t)=\operatorname{Tr}_{B}[\hat{\rho}_{SB}(t)]$. The influence of the Gaussian bath on the system-reduced density operator is determined by the bath correlation function (BCF),

$$C_{jj^{\prime}}^{B}(t,t^{\prime})=\operatorname{Tr}_{B}\left[\hat{B}_{j}(t)\hat{B}_{j^{\prime}}(t^{\prime})\hat{\rho}_{B}(0)\right],$$

(2)

for $t\geq t^{\prime}\geq 0$ with $\hat{B}_{j}(t)=e^{i\hat{H}_{B}t}\hat{B}_{j}e^{-i\hat{H}_{B}t}$. When $\hat{\rho}_{B}(0)$ is stationary, i.e., $e^{-i\hat{H}_{B}t}\hat{\rho}_{B}(0)e^{i\hat{H}_{B}t}=\hat{\rho}_{B}(0)$, the BCF only depends on the time difference $t-t^{\prime}$, i.e., $C_{jj}^{B}(t,t^{\prime})=C_{jj}^{B}(t-t^{\prime})$.

In many cases of interest, the above unitary dynamics using a continuous bath leads to a dissipative bath correlation function, i.e., a decaying $C_{jj}^{B}(t-t^{\prime})$ with $\lim_{t-t^{\prime}\to\infty}C_{jj}^{B}(t-t^{\prime})\to 0$. To reproduce the associated dynamics of $\hat{\rho}_{S}(t)$ with a discrete bath, we now introduce a pseudomode description. First, we review the conventional pseudomode formulation [20], which introduces an auxiliary bosonic bath $A$ with a Lindblad dissipator $\mathcal{D}_{A}$ applied within the bath $A$. The density operator across the system and auxiliary bath, $\hat{\rho}_{SA}$, evolves under the following Lindblad master equation,

$$\partial_{t}\hat{\rho}_{SA}=-i\left[\hat{H}_{\text{aux}},\hat{\rho}_{SA}\right]+\mathcal{D}_{A}\hat{\rho}_{SA},$$

(3)

where the Hamiltonian $\hat{H}_{\text{aux}}=\hat{H}_{S}+\hat{H}_{A}+\hat{H}_{SA}$ is composed of the system (bath)-only Hamiltonian $\hat{H}_{S}$ ($\hat{H}_{A}$) and coupling Hamiltonian $\hat{H}_{SA}=\sum_{j}\hat{S}_{j}\hat{A}_{j}$. $\mathcal{D}_{A}$ has a Lindblad form

$$\mathcal{D}_{A}\bullet=\sum_{\mu}\hat{L}_{\mu}\bullet\hat{L}_{\mu}^{\dagger}-\frac{1}{2}\{\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu},\bullet\},$$

(4)

which guarantees the CPTP property of the global dynamics consisting of the system and the pseudomodes. Here, $\hat{L}_{\mu}$ can be an arbitrary operator within the bath $A$. The Lindblad dissipator can be expressed in terms of a pseudomode basis,

$$\mathcal{D}_{A}\bullet=2\sum_{kk^{\prime}}\Gamma_{kk^{\prime}}\left(\hat{F}_{k^{\prime}}\bullet\hat{F}_{k}^{\dagger}-\frac{1}{2}\{\hat{F}_{k}^{\dagger}\hat{F}_{k^{\prime}},\bullet\}\right),$$

(5)

where $\hat{F}_{k}$ is an operator within the $k$th pseudomode, and $\Gamma_{kk^{\prime}}$ is a positive semidefinite (PSD) matrix from the CPTP condition. $\hat{H}_{A}$ and $\mathcal{D}_{A}$ are quadratic in the creation and annihilation operators of the pseudomode, and $\hat{A_{j}}$ and $\hat{F}_{k}$ are linear, to maintain the Gaussian properties of the bath. $\mathcal{D}_{A}$ 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 $\mathcal{L}_{A}\bullet=-i[\hat{H}_{A},\bullet]+\mathcal{D}_{A}\bullet$ and initial pseudomode reduced density operator $\hat{\rho}_{A}(0)$,

$$C_{jj^{\prime}}^{A}(t,t^{\prime})=\operatorname{Tr}_{A}\left[\hat{A}_{j}e^{\mathcal{L}_{A}(t-t^{\prime})}\hat{A}_{j^{\prime}}e^{\mathcal{L}_{A}t^{\prime}}\hat{\rho}_{A}(0)\right].$$

(6)

Ref. [20] established the equivalence condition between the unitary bath dynamics and the pseudomode dynamics in terms of the BCF: if $C^{B}_{jj^{\prime}}(t,t^{\prime})=C^{A}_{jj^{\prime}}(t,t^{\prime})$ for $\forall t\geq t^{\prime}\geq 0$, then $\operatorname{Tr}_{B}[\hat{\rho}_{SB}(t)]=\operatorname{Tr}_{A}[\hat{\rho}_{SA}(t)]$ ¹¹1In the original formulation in [20], the additional conditions on the one-point bath expectation values $\operatorname{Tr}[\hat{B}_{j}(t)\hat{\rho}_{B}(0)]$ were imposed. In the main text, we further assumed that the initial bath reduced density operators, $\hat{\rho}_{B}(0)$ and $\hat{\rho}_{A}(0)$, are number-conserving operators. It leads to the one-point bath expectation values becoming zero for both the original unitary and pseudomode formulations.. As in the unitary case in Sec. II.1, when $\hat{\rho}_{A}(0)$ is stationary, i.e., $e^{\mathcal{L}_{A}t^{\prime}}\hat{\rho}_{A}(0)=\hat{\rho}_{A}(0)$, the BCF depends solely on $t-t^{\prime}$.

As described in the introduction, the conventional pseudomode 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 $\mathcal{D}_{SA}$, expressed in terms of Lindblad dissipators,

$$\mathcal{D}_{SA}\bullet=\sum_{j}\hat{L}^{\prime}_{j}\bullet\hat{S}_{j}+\hat{S}_{j}\bullet\hat{L}^{\prime{\dagger}}_{j}-\frac{1}{2}\{\hat{S}_{j}\hat{L}^{\prime}_{j}+\hat{L}^{\prime{\dagger}}_{j}\hat{S}_{j},\bullet\},$$

(7)

where $\hat{L}^{\prime}_{j}=\sum_{k}2M_{jk}\hat{F}_{k}$. We call this coupling a Lindblad coupling, as illustrated in Fig. 1b. The total dissipator $\mathcal{D}=\mathcal{D}_{A}+\mathcal{D}_{SA}$ can be compactly written as

$$\mathcal{D}\bullet=2\sum_{pq}\widetilde{\Gamma}_{pq}\left(\hat{F}_{q}\bullet\hat{F}^{\dagger}_{p}-\frac{1}{2}\{\hat{F}^{\dagger}_{p}\hat{F}_{q},\bullet\}\right),$$

(8)

where the index $p$ and $q$ refer to both the coupling index $j$ and the pseudomode bath index $k$, i.e., $p\in\{j\}\cup\{k\}$, and if $p=j$, $\hat{F}_{j}=\hat{S}_{j}$. The coefficient $\widetilde{\Gamma}_{pq}$ is

$$\widetilde{\bm{\Gamma}}=\left(\begin{array}[]{cc}\bm{0}&\bm{M}\\ \bm{M}^{\dagger}&\bm{\Gamma}\end{array}\right),$$

(9)

where the bold font denotes matrices. We note that the matrix $\widetilde{\Gamma}_{pq}$ is not PSD in general, unlike the PSD matrix $\Gamma_{pq}$, 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 $\mathcal{L}_{SA}\bullet=-i[\hat{H}_{SA},\bullet]+\mathcal{D}_{SA}\bullet$ and further decompose it in the following form,

$$\mathcal{L}_{SA}=-i\sum_{j}\mathcal{S}_{j}\mathcal{F}_{j}+i\sum_{j}\widetilde{\mathcal{S}}_{j}\widetilde{\mathcal{F}}_{j},$$

(10)

where the system superoperators, $\mathcal{S}_{j}$ and $\widetilde{\mathcal{S}}_{j}$ are defined as $\mathcal{S}_{j}\bullet=\hat{S}_{j}\bullet$ and $\widetilde{\mathcal{S}}_{j}\bullet=\bullet\hat{S}_{j}$. The superoperators, $\mathcal{F}_{j}$ and $\widetilde{\mathcal{F}}_{j}$, are obtained as,

$$\begin{gathered}\mathcal{F}_{j}\bullet=\hat{A}_{j}\bullet+\bullet i\hat{L}^{\prime{\dagger}}_{j}-\frac{i}{2}(\hat{L}^{\prime}_{j}+\hat{L}^{\prime{\dagger}}_{j})\bullet,\\ \widetilde{\mathcal{F}}_{j}\bullet=\bullet\hat{A}_{j}-i\hat{L}^{\prime}_{j}\bullet+\bullet\frac{i}{2}(\hat{L}^{\prime}_{j}+\hat{L}^{\prime{\dagger}}_{j}).\end{gathered}$$

(11)

The BCF is defined with the superoperators $\mathcal{F}_{j}$,

$$C_{jj^{\prime}}^{A}(t,t^{\prime})=\operatorname{Tr}_{A}\left[\mathcal{F}_{j}e^{\mathcal{L}_{A}(t-t^{\prime})}\mathcal{F}_{j^{\prime}}e^{\mathcal{L}_{A}t^{\prime}}\hat{\rho}_{A}(0)\right].$$

(12)

We remark that $\mathcal{F}_{j}$ reduces to $\hat{A}_{j}$ when $\hat{L}^{\prime}_{j}=0$, and then the BCF in Eq. 12 also reduces to the BCF in Eq. 6. Similar to the pseudomode description in Sec. II.2, 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) [52].

In this framework, violation of the CP condition leads the dissipator Liouvillian $\mathcal{D}$ to have eigenvalues with positive real parts. Thus, in the absence of any coherent term (i.e., the commutator term of the form $-i[\hat{H}_{\rm aux},\bullet]$ for some Hamiltonian $\hat{H}_{\rm 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 dynamics despite violating the CP condition. We will analyze this effect in more detail in Sec. IV.

We next introduce an ansatz for the quasi-Lindblad pseudomode 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.1 with the initial bath being a stationary state, such as a thermal state $\hat{\rho}_{B}(0)\propto e^{-\beta\hat{H}_{B}}$, where $\beta$ represents the inverse temperature.

We start with a bosonic bath, setting the initial reduced density operator to a vacuum state, $\hat{\rho}_{A}(0)=\mathinner{|{\bm{0}}\rangle}\!\!\mathinner{\langle{\bm{0}}|}$. The dissipator $\mathcal{D}_{A}$ is set as $\hat{F}_{k}=\hat{d}_{k}$, where $\hat{d}_{k}$ is the bosonic annihilation operator of the $k$th pseudomode, and the Hamiltonian $\hat{H}_{A}=\sum_{kk^{\prime}}H^{A}_{kk^{\prime}}\hat{d}_{k}^{\dagger}\hat{d}_{k^{\prime}}$. This ansatz satisfies the stationary condition $\mathcal{L}_{A}\hat{\rho}_{A}(0)=0$. Then, $e^{\mathcal{L}_{A}t^{\prime}}\hat{\rho}_{A}(0)=\hat{\rho}_{A}(0)$ and $C^{A}_{jj^{\prime}}(t,t^{\prime})$ is only dependent on the time difference $C^{A}_{jj^{\prime}}(t,t^{\prime})=C^{A}_{jj^{\prime}}(t-t^{\prime})=C^{A}_{jj^{\prime}}(\Delta t)$, where $\Delta t=t-t^{\prime}$. The system-bath coupling is set as $\hat{A}_{j}\ =\sum_{k}V_{jk}\hat{d}_{k}+V_{jk}^{*}\hat{d}_{k}^{\dagger}$.

We obtain the BCF based on the above ansatz:

$$\bm{C}^{A}(\Delta t)=(\bm{V}-i\bm{M})e^{-\bm{Z}\Delta t}(\bm{V}+i\bm{M})^{\dagger},$$

(13)

where $\bm{Z}=\bm{\Gamma}+i\bm{H}^{A}$. This BCF expression has a gauge redundancy: $\bm{Z}\rightarrow\bm{G}\bm{Z}\bm{G}^{-1}$, $\bm{V}-i\bm{M}\rightarrow(\bm{V}-i\bm{M})\bm{G}^{-1}$, $\bm{V}+i\bm{M}\rightarrow(\bm{V}+i\bm{M})\bm{G}^{\dagger}$, for arbitrary invertible matrix $\bm{G}$. Therefore, we focus on a diagonal matrix $\bm{Z}=\mathrm{diag}(\{z_{k}\})$, which can be transformed to a more general diagonalizable $\bm{Z}$ through a unitary transformation ²²2In principle, we have a general expression by considering non-diagonalizable $\bm{Z}$. However, numerically, fitting with arbitrary $\bm{Z}$ yields the same optimal solution as with diagonalizable $\bm{Z}$.. In this case, we define $\hat{H}_{A}=\sum_{k}E_{k}\hat{d}^{\dagger}_{k}\hat{d}_{k}$ and $\Gamma_{kk^{\prime}}=\gamma_{k}\delta_{kk^{\prime}}$, and $z_{k}=\gamma_{k}+iE_{k}$. Even with a diagonal matrix $\bm{Z}$, $\bm{V}$ and $\bm{M}$ still have a gauge redundancy associated with a diagonal $\bm{G}$, since $\bm{Z}$ remains invariant under the transformation $\bm{G}\bm{Z}\bm{G}^{-1}=\bm{Z}$. 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 similarly. In particular, we focus on number-conserving fermionic impurity models where the coupling Hamiltonian is given by $\hat{H}_{SB}=\sum_{j}\hat{a}_{j}^{\dagger}\hat{B}_{j}+\hat{B}_{j}^{\dagger}\hat{a}_{j},\>\hat{B}_{j}=\int d\omega\ g_{j}(\omega)\hat{c}_{\omega}$ ($\hat{a}_{j}^{\dagger}$ is the creation operator of the impurity with index $j$, 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: $C^{>B}_{jj^{\prime}}$ and $C^{<B}_{jj^{\prime}}$. For any fermionic Gaussian state $\hat{\rho}_{B}(0)$, $C^{>B}_{jj^{\prime}}$ and $C^{<B}_{jj^{\prime}}$ are defined as,

$$\begin{gathered}C^{>B}_{jj^{\prime}}(t,t^{\prime})=\operatorname{Tr}_{B}\left[\hat{B}_{j}(t)\hat{B}^{{\dagger}}_{j^{\prime}}(t^{\prime})\hat{\rho}_{B}(0)\right],\\ C^{<B}_{jj^{\prime}}(t,t^{\prime})=\operatorname{Tr}_{B}\left[\hat{B}^{\dagger}_{j}(t)\hat{B}_{j^{\prime}}(t^{\prime})\hat{\rho}_{B}(0)\right].\end{gathered}$$

(14)

Corresponding to these two BCFs, we use two different auxiliary baths with an initial reduced density operator $\hat{\rho}_{A}(0)=\hat{\rho}_{A_{1}}(0)\otimes\hat{\rho}_{A_{2}}(0)=\mathinner{|{\bm{0}}\rangle}\!\!\mathinner{\langle{\bm{0}}|}\otimes\mathinner{|{\bm{1}}\rangle}\!\!\mathinner{\langle{\bm{1}}|}$. The bath-only Hamiltonian and dissipator are set to be diagonal without loss of generality, as discussed in the bosonic case. The Hamiltonian is $\hat{H}_{A}=\sum_{k_{1}}E_{k_{1}}\hat{d}^{\dagger}_{k_{1}}\hat{d}_{k_{1}}+\sum_{k_{2}}E_{k_{2}}\hat{d}^{\dagger}_{k_{2}}\hat{d}_{k_{2}}$ where $k_{1}$ and $k_{2}$ refer to the pseudomode index in each bath, respectively, and the dissipator is set to be

$$\begin{gathered}\mathcal{D}_{A_{1}}\bullet=\sum_{k_{1}}2\gamma_{k_{1}}\left(\hat{d}_{k_{1}}\bullet\hat{d}^{{\dagger}}_{k_{1}}-\frac{1}{2}\{\hat{d}^{{\dagger}}_{k_{1}}\hat{d}_{k_{1}},\bullet\}\right),\\ \mathcal{D}_{A_{2}}\bullet=\sum_{k_{2}}2\gamma_{k_{2}}\left(\hat{d}^{\dagger}_{k_{2}}\bullet\hat{d}_{k_{2}}-\frac{1}{2}\{\hat{d}_{k_{2}}\hat{d}^{\dagger}_{k_{2}},\bullet\}\right),\end{gathered}$$

(15)

which satisfies $\mathcal{D}_{A_{i}}\hat{\rho}_{A_{i}}(0)=0$ for both $i=1,2$. The system-bath Hamiltonian coupling is $\hat{H}_{SA}=\hat{a}_{j}^{\dagger}\hat{A}_{j}+\hat{A}_{j}^{\dagger}\hat{a}_{j}$, $\hat{A}_{j}=\sum_{k_{1}}(\bm{V}_{1})_{jk_{1}}\hat{d}_{k_{1}}+\sum_{k_{2}}(\bm{V}_{2})_{jk_{2}}\hat{d}_{k_{2}}$ and the system-bath Lindblad coupling is set to be

$$\begin{gathered}\mathcal{D}_{SA_{1}}\bullet=\sum_{j}\hat{a}_{j}\bullet\hat{L}_{1j}^{\dagger}+\hat{L}_{1j}\bullet\hat{a}_{j}^{\dagger}-\frac{1}{2}\{\hat{L}_{1j}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{L}_{1j},\bullet\},\\ \mathcal{D}_{SA_{2}}\bullet=\sum_{j}\hat{a}^{\dagger}_{j}\bullet\hat{L}_{2j}+\hat{L}^{\dagger}_{2j}\bullet\hat{a}_{j}-\frac{1}{2}\{\hat{L}_{2j}\hat{a}^{\dagger}_{j}+\hat{a}_{j}\hat{L}^{\dagger}_{2j},\bullet\},\end{gathered}$$

(16)

where $\hat{L}_{1j}=\sum_{k_{1}}2(\bm{M}_{1})_{jk_{1}}\hat{d}_{k_{1}}$ and $\hat{L}_{2j}=\sum_{k_{2}}2(\bm{M}_{2})_{jk_{2}}\hat{d}_{k_{2}}$. The BCF with this ansatz becomes

$$\begin{gathered}\bm{C}^{>A}(\Delta t)=(\bm{V}_{1}-i\bm{M}_{1})e^{-\bm{Z}_{1}\Delta t}(\bm{V}_{1}+i\bm{M}_{1})^{\dagger},\\ \bm{C}^{<A}(\Delta t)=(\bm{V}_{2}-i\bm{M}_{2})^{*}e^{-\bm{Z}_{2}\Delta t}(\bm{V}_{2}+i\bm{M}_{2})^{T},\end{gathered}$$

(17)

where $\bm{Z}_{1}$ and $\bm{Z}_{2}$ are diagonal matrices with elements $z_{k_{1}}=\gamma_{k_{1}}+iE_{k_{1}}$ and $z_{k_{2}}=\gamma_{k_{2}}-iE_{k_{2}}$. We note that the baths $A_{1}$ and $A_{2}$ only contribute to the BCF $\bm{C}^{>A}$ and $\bm{C}^{<A}$, 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 $j=1$). In this case, the BCF simplifies to a scalar function without indices:

$$C^{A}(\Delta t)=\sum_{k}w_{k}e^{-z_{k}\Delta t},$$

(18)

where $w_{k}=(V_{k}-iM_{k})(V_{k}^{*}-iM_{k}^{*})$. This expression is a sum of complex exponential functions with complex weights. In contrast, the conventional pseudomode description corresponds to $M_{k}=0$, resulting in positive weights $w_{k}=|V_{k}|^{2}$. 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 pseudomode. 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 [37, 38], which has previously been reported as one of the most competitive methods for BCF fitting [44], to determine the complex exponents, $z_{k}$. We then perform a least-squares fitting of the complex weights $w_{k}$. We note that this simple unconstrained least-squares approach is possible because we do not restrict $w_{k}\geq 0$ as would be required in the conventional pseudomode approach. Once $w_{k}$ is determined, we choose the gauge for $V_{k}$ and $M_{k}$. In this limit, the gauge choice is parameterized as $V_{k}-iM_{k}=\kappa_{k}\sqrt{w_{k}},V_{k}+iM_{k}=\kappa_{k}^{*-1}\sqrt{w^{*}_{k}}$ where $\kappa_{k}\in\mathbb{C}$. Here, we follow a simple heuristic of minimizing the magnitude of $M_{k}$, as this term determines the extent of CP condition violation. We describe the result for the bosonic bath expression in Eq. 13 for $\bm{V}$ and $\bm{M}$, and its adaptation to the fermionic BCF expression is straightforward. The solution of $V_{k}$ and $M_{k}$ with minimum $|M_{k}|$ is then given by $V_{k}-iM_{k}=\sqrt{w_{k}}$ with $V_{k},M_{k}\in\mathbb{R}$, ($\kappa_{k}=1$) which we provide its proof in the SM [52]. We will examine other choices of $V_{k}$ and $M_{k}$ and their effect on the stability of the system-bath dynamics in Sec. IV.

There are several approaches to relax the CPTP condition in pseudomode representations. One approach involves a non-Hermitian pseudomode formulation [32, 33, 34, 35], where the system-bath coupling is given by a non-Hermitian Hamiltonian (without Lindblad coupling). This formulation results in real weights, $w_{k}\in\mathbb{R}$ [32, 35], or complex weights from a pair of pseudomodes [35]. Another approach, presented in [41], 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 [41] represents one gauge choice for $\kappa_{k}$ ³³3The corresponding $\kappa_{k}$ for the Lindblad-like equation in [41] is $\kappa_{k}=\sqrt{\operatorname{Re}[w_{k}]/w_{k}}$. within our quasi-Lindblad pseudomode framework.

For general couplings with multiple $j$, a similar strategy can be adopted. In this case, we first find a sum of exponential functions with a matrix-valued weight $\bm{C}^{A}(\Delta t)\approx\sum_{k}\bm{W}_{k}e^{-z_{k}\Delta t}$. We determine $z_{k}$ by fitting a scalar quantity, such as the sum of the BCF entries, $\sum_{jj^{\prime}}C^{A}_{jj^{\prime}}(\Delta t)\approx\sum_{k}\left(\sum_{jj^{\prime}}(\bm{W}_{k})_{jj^{\prime}}\right)e^{-z_{k}\Delta t}$, or the trace $\sum_{j}C^{A}_{jj}(\Delta t)\approx\sum_{k}\left(\sum_{j}(\bm{W}_{k})_{jj}\right)e^{-z_{k}\Delta t}$. The matrix-valued complex weights $\bm{W}_{k}$ are then obtained through the least-squares fitting.

Finding $\bm{V}$ and $\bm{M}$ from $\bm{W}_{k}$ involves additional matrix decompositions. However, the representation from Eq. 13 has a rank-1 matrix factorization form, $(\bm{W}_{k})_{jj^{\prime}}=(\bm{V}-i\bm{M})_{jk}(\bm{V}+i\bm{M})^{*}_{j^{\prime}k}$, and $\bm{W}_{k}$ from the least-squares fitting is generally not a rank-1 matrix, but of rank $N_{S}$, after denoting the range of $j$ as $N_{S}$. To address this, we represent the rank-$N_{S}$ $\bm{W}_{k}$ by indexing the pseudomode with $n=(k,s)$ where the index $s$ is an additional index with a range of $N_{S}$. Then, we assign the diagonal elements of $\bm{Z}$ as $z_{n}=z_{(k,s)}=z_{k}$. This gives the following expression:

$$C^{A}_{jj^{\prime}}(\Delta t)=\sum_{n=(k,s)}(\bm{V}-i\bm{M})_{jn}(\bm{V}+i\bm{M})^{*}_{j^{\prime}n}e^{-z_{k}\Delta t},$$

(19)

where $(\bm{W}_{k})_{jj^{\prime}}=\sum_{s}(\bm{V}-i\bm{M})_{jn}(\bm{V}+i\bm{M})^{*}_{j^{\prime}n}$. Thus, given $N_{\exp}$ exponents $z_{k}$ and $N_{S}$ coupling operators, $N_{S}N_{\exp}$ pseudomodes are used in the BCF fitting. With this setting, we decompose $\bm{W}_{k}$ with a singular value decomposition (SVD), $\bm{W}_{k}=\bm{U}_{k}\bm{\Sigma}_{k}\widetilde{\bm{U}}_{k}^{\dagger}$, where $\bm{\Sigma}_{k}$ is a diagonal matrix. One possible choice for $\bm{V}$ and $\bm{M}$ is $(\bm{V}-i\bm{M})_{jn}=(\bm{U}_{k}\bm{\Sigma}^{1/2}_{k})_{js},(\bm{V}+i\bm{M})_{jn}=(\widetilde{\bm{U}}_{k}{\bm{\Sigma}_{k}}^{1/2})_{js}$, which reduces to the gauge choice above for the single coupling limit ⁴⁴4The SVD representation has a phase redundancy in $\bm{U}_{k}$, $\bm{U}_{k}\rightarrow\bm{U}_{k}e^{i\bm{\Theta}_{k}}$, and also in $\bm{M}_{k}$, $\bm{M}_{k}\rightarrow\bm{M}_{k}e^{i\bm{\Theta}_{k}}$, where $\bm{\Theta}_{k}$ is a real diagonal matrix. However, it does not affect the overall norm of $\bm{M}_{k}$..

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 [40, 42, 41, 43, 44] and tensor network influence functionals [45, 46]. 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 $J_{jj^{\prime}}(\omega)=g_{j}(\omega)g_{j^{\prime}}^{*}(\omega)$, we will fit the associated zero-temperature BCF [2],

$$C_{jj^{\prime}}(t)=\int J_{jj^{\prime}}(\omega)e^{-i\omega t}d\omega,$$

(20)

using the three different ansatz. We examine the performance of these ansatz for three spectral densities, as illustrated in Fig. 2. The first case is a subohmic spectral density,

$$J(\omega)=\frac{\alpha}{2}\omega_{c}^{1-s}\omega^{s}e^{-\omega/\omega_{c}},\quad\omega\in[0,\infty),$$

(21)

the second, a semicircular spectral density,

$$J(\omega)=\frac{\Gamma}{\pi}\sqrt{1-\frac{\omega^{2}}{W^{2}}},\quad\omega\in[0,W],$$

(22)

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

$$J_{jj^{\prime}}(\omega)=\begin{pmatrix}1&r(\omega)\\ r^{*}(\omega)&1\end{pmatrix}J(\omega),$$

(23)

where $r(\omega)=\exp(-i\omega/2W)$ and $J(\omega)$ is the semicircular spectral density specified above. For the subohmic case, we use $s=0.5$ and $\alpha=\omega_{c}=1$, and for the semicircular case, we set $\Gamma=1$ and $W=10\Gamma$.

Figs. 2 (a–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 multi-site case) weights (blue dash-dotted), and (3) a discrete Hamiltonian (green dotted). Fitting scheme (2) represents the Lorentzian pseudomode. Here, the parameters are optimized through numerical minimization of the fitting error. For fitting (3), we used a discretization scheme based on Gaussian quadrature using Legendre polynomials [14]. Details of these fitting schemes are in the SM [52]. The fits to the quasi-Lindblad ansatz with $N_{\exp}=4$ ($N_{\exp}=6$) for Fig. 2a (Fig. 2b and Fig. 2c), respectively, show high accuracy already with a small number of exponentials. In Fig. 2d, Fig. 2e, and 2f, the maximum fitting error as a function of $N_{\text{exp}}$ is illustrated. This comparison shows the clear advantage of the quasi-Lindblad pseudomode compared to the Lorentzian pseudomode and direct discretization schemes and illustrates an exponential decrease of error with respect to $N_{\text{exp}}$ in both cases.

We then test the accuracy of the quasi-Lindblad pseudomode 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 $\hat{H}_{S}=\omega_{0}\hat{\sigma}_{z}/2$ with $\omega_{0}=4$ and $\hat{S}_{j=1}=\hat{\sigma}_{z}$ and the initial system-reduced density operator $\hat{\rho}_{S}(0)=\mathinner{|{+}\rangle}\!\!\mathinner{\langle{+}|}$, with $\mathinner{|{+}\rangle}=1/\sqrt{2}(\mathinner{|{0}\rangle}+\mathinner{|{1}\rangle})$. While the analytical solution for this setup is known [2, 4], obtaining accurate numerical results remains nontrivial [56, 26, 57], particularly due to the requirement for a precise description of the BCF.

We choose the $s=0.5$ subohmic spectral density as above ($\alpha=\omega_{c}=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 pseudomode, and a discrete Hamiltonian from Gaussian quadrature. Fig. 3 shows the time evolution of the coherence $\langle\hat{\sigma}_{x}(t)\rangle$ using these three different bath discretizations. We see that the time evolution from the discrete Hamiltonian shows unphysical oscillations in the long-time limit, while the Lorentzian pseudomode qualitatively captures the dissipative dynamics but at a much lower accuracy than the quasi-Lindblad pseudomode representation.

We next study the dynamics of the spinless noninteracting single-site and two-site models with $\hat{H}_{S}=\varepsilon\sum_{j}\hat{n}_{j}$ ($\hat{n}_{j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$). Here, we can carry out efficient Gaussian calculations, which we do for both a continuous unitary bath and for the Lindblad pseudomode dynamics [22, 58]. We choose the semicircular spectral density in Eq. 22 and Eq. 23 with $W=10\Gamma$, but with $\omega\in[-W,W]$, $\varepsilon=-4\Gamma$, and consider the initial bath as a fully occupied (empty) state for $\omega<0$ ($\omega>0$). We fit the two BCFs, $C^{>}_{jj^{\prime}}(t)$ and $C^{<}_{jj^{\prime}}(t)$, with $N_{\text{exp}}$ exponential functions each and represent the bath with in total $2N_{\text{exp}}$ ($4N_{\text{exp}}$) pseudomodes for the single-site (two-site) impurity model, respectively. The initial impurity state is assumed to be an empty state, $\hat{\rho}_{S}(0)=\mathinner{|{0}\rangle}\!\!\mathinner{\langle{0}|}$, for the single-site model, and $\hat{\rho}_{S}(0)=\mathinner{|{10}\rangle}\!\!\mathinner{\langle{10}|}$ for the two-site model. Fig. 4a and 4b show the time-dependence of the impurity occupation, $\langle\hat{n}_{j}(t)\rangle$, with $N_{\text{exp}}=4$ for the exact dynamics, the quasi-Lindblad pseudomode, Lorentzian pseudomode, and discrete Hamiltonian dynamics. These results clearly show the improved accuracy of the quasi-Lindblad pseudomode relative to the Lorentzian pseudomode, which reaches the wrong steady state. The inset in Fig. 4a and 4b shows the maximum error of the impurity occupation as a function of $N_{\text{exp}}$, which also shows an exponential convergence in $N_{\text{exp}}$.

Next, we carry out a test with interacting models - a spinful single-site impurity model, known as the single-impurity Anderson model (SIAM), with $\hat{H}_{S}=U\hat{n}_{\uparrow}\hat{n}_{\downarrow}+\varepsilon(\hat{n}_{\uparrow}+\hat{n}_{\downarrow})$, and a spinless two-site impurity model with $\hat{H}_{S}=U\hat{n}_{0}\hat{n}_{1}+\varepsilon(\hat{n}_{0}+\hat{n}_{1})$ - and choose $U=-2\varepsilon=8\Gamma$. For the SIAM model, we assume that there is the same bath for each spin with a semicircular spectral density without any coupling between the baths. The initial bath state is assumed to be the same as in the noninteracting case. As a reference for the exact unitary dynamics, we use a large bath discretization obtained from Gaussian quadrature with Legendre polynomials [14] - 80 (200) bath modes for the single-site (two-site) impurity model, respectively. For both the unitary reference and pseudomode Lindblad dynamics, we use the time-dependent density matrix renormalization group (tdDMRG) to propagate the dynamics [59, 60] with bond dimensions of up to 250, using the Block2 [61, 62] package. Fig. 4c and 4d show the time-dependence of $\langle\hat{n}_{j}(t)\rangle$ ($j=\uparrow$ for the SIAM, $j=0,1$ for the two-site impurity). The quasi-Lindblad ansatz was constructed from $N_{\text{exp}}=4$ ($2N_{\text{exp}}=8$ pseudomodes per each spin) for SIAM and $N_{\text{exp}}=8$ ($4N_{\text{exp}}=32$ pseudomodes) for the two-site impurity model. We see that the quasi-Lindblad pseudomode accurately reproduces the original unitary dynamics, including the long-lived oscillation in the two-site impurity model.