Pseudomode Method in Open Quantum Systems

• Pseudomode Method is a non-perturbative framework that replaces a complex bath with discrete damped modes encoding memory and information backflow.
• It converts non-Markovian, integro-differential dynamics into a time-local Lindblad master equation via analytic pole decomposition of the bath correlation function.
• Extensions like coupled Lindblad, quasi-Lindblad, and purified models link the method to HEOM, enhancing simulation efficiency for strong-coupling and memory effects.

Searching arXiv for recent and foundational pseudomode-method papers to ground the article. arXiv search query: pseudomode open quantum systems non-Markovian site:arxiv.org The pseudomode method is a non-perturbative framework for open quantum dynamics in which a structured environment is replaced by a finite set of auxiliary discrete modes, each damped into a Markovian residual reservoir, so that the enlarged dynamics becomes time-local while the reduced dynamics of the original system retains its non-Markovian character. In its standard form, the method is built from the analytic structure of the bath spectral density or, equivalently, from an exponential representation of the bath correlation function; the resulting pseudomodes encode memory, information backflow, and strong-coupling effects in a compact auxiliary space (Pleasance et al., 2021). Later developments extended the construction to coupled Lindblad embeddings, quasi-Lindblad and non-Hermitian formulations, purified input-output models, and a constructive one-to-one correspondence with HEOM for exponential bath correlation functions (Huang et al., 12 Jun 2025, Müller et al., 7 Apr 2026).

1. Open-system formulation

In the bosonic setting emphasized by Pleasance and Petruccione, one considers a system $S$ with Hamiltonian $H_S$ coupled linearly to a bosonic bath $E$,

$$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$

with

$B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$

The reduced state

$$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$

is fully characterized, for a Gaussian, stationary bath with $\langle B\rangle=0$, by the two-point correlation function

$$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$

For a thermal bath one often writes

$$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$

with one-sided spectral density $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$ (Pleasance et al., 2021).

This formulation makes explicit that the continuum bath enters the reduced dynamics only through $H_S$0. The pseudomode strategy therefore targets the bath correlation function rather than the microscopic bath itself. In the rotating-wave two-level setting, the same logic appears in the non-local amplitude equation

$H_S$1

with a structured reservoir approximated in practice by a sum of Lorentzian peaks,

$H_S$2

where each pole sits at $H_S$3; here $H_S$4 is the reservoir’s correlation (memory) time and $H_S$5 sets its coupling strength (Ohyama et al., 2017).

2. Pole decomposition and the basic pseudomode mapping

The standard exact construction assumes that the two-sided spectral density $H_S$6 is meromorphic in the lower-half complex $H_S$7-plane and decays faster than $H_S$8 as $H_S$9. If its poles are

$E$0

with residues $E$1, then contour integration yields

$E$2

Each pole is then mapped to a pseudomode $E$3 with frequency $E$4, decay rate $E$5, and coupling $E$6. The enlarged Hamiltonian on $E$7 is

$E$8

and, after attaching each pseudomode to its own zero-temperature Markovian reservoir, the joint density matrix obeys

$E$9

Matching the two-point bath correlation functions $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$0 implies

$$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$1

so that all non-Markovian memory effects are encoded exactly in the discrete pseudomodes (Pleasance et al., 2021).

In the Lorentzian case this mapping reduces to the familiar one-pseudomode picture. For a single-peak structure function

$$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$2

there is one simple pole $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$3 with residue $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$4. The corresponding extended dynamics has one pseudomode $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$5, effective Hamiltonian

$$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$6

and jump operator

$$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$7

thereby converting the original integro-differential dynamics into a Markovian Lindblad equation on the enlarged Hilbert space (Behzadi et al., 2017).

3. Formulations and generalizations

Several distinct pseudomode formulations are now in use.

Formulation	Distinctive feature	Representative paper
Generalized theory	Complex residues can produce non-Hermitian couplings and non-Lindblad generators	(Pleasance et al., 2020)
Quasi-Lindblad pseudomodes	Trace preserving, but not completely positive on $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$8	(Park et al., 2024)
Coupled Lindblad pseudomodes	Pseudomodes may be mutually coupled through $$H = H_S + H_E + H_I,\qquad H_E = \sum_k \omega_k a_k^\dagger a_k,\qquad H_I = A\otimes B,$$9 and $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$0	(Huang et al., 12 Jun 2025)
Purified input-output model	Pure auxiliary modes for system and bath observables, including non-Gaussian initial states	(Liang et al., 2024)
Purified nonlinear model	Extension from linear $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$1 couplings to general $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$2	(Zhang et al., 25 Feb 2025)

The generalized theory of Pleasance, Garraway, and Petruccione makes explicit that exact pseudomode embeddings need not always be Lindbladian. When the analytically continued spectral density has multiple poles or complex residues, some couplings become intrinsically complex, the interaction can be non-Hermitian, and the resulting exact generator may violate complete positivity. For two discrete modes, however, the theory demonstrates how to convert between pathological and Lindblad forms of the master equation (Pleasance et al., 2020). The quasi-Lindblad framework sharpens this point: it represents the bath correlation function as

$B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$3

and emphasizes that the pseudomode representation is not unique but differs by a gauge choice; the physical reduced dynamics of $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$4 remains exact if the enlarged dynamics is propagated exactly, while the gauge can affect numerical stability (Park et al., 2024).

The coupled Lindblad formulation replaces independent pseudomodes by a finite set of damped modes with

$B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$5

with $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$6 and $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$7, thereby enlarging the admissible bath-correlation structures beyond diagonal sums of Lorentzians (Huang et al., 12 Jun 2025). The purified input-output construction, by contrast, analytically continues each exponential into one-sided components realized by pure pseudomodes $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$8 and $B=\sum_k g_k(a_k+a_k^\dagger),\qquad A=A^\dagger.$9; this permits pure-state propagation, access to bath observables, and treatment of bosonic baths prepared in non-Gaussian initial states (Liang et al., 2024). The subsequent purified nonlinear extension preserves the same philosophy for general nonlinear system-bath interactions $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$0 (Zhang et al., 25 Feb 2025).

A recurring subtlety is that the auxiliary model may itself be unphysical. In particular, optimal pseudomode fits can require complex or negative $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$1, imaginary $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$2, or classical fields with nonstandard statistics. A separate line of work therefore proposes reproducing the effects of an unphysical pseudomode model through measurement results over an ensemble of physical systems involving ancillary harmonic modes and an optional stochastic driving field, combined with an extrapolation technique (Cirio et al., 2023).

4. Exactness, memory, and relation to HEOM

The core exactness statement of the standard method is correlation matching: once the pseudomode operator $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$3 is chosen so that its two-point function $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$4 equals the physical bath correlation $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$5, the reduced system state produced by tracing out the continuum bath is identical to that produced by tracing out the pseudomodes after Markovian evolution on the enlarged space (Pleasance et al., 2021). The method therefore trades temporal nonlocality for an enlarged but time-local description.

This enlarged description also gives a direct interpretation of reservoir memory. In the two-level atom setting, the no-jump state can be expanded as

$$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$6

where $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$7 satisfies exactly the original integrodifferential equation, while the amplitudes $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$8 encode the temporary storage of excitation in pseudomode $$\rho_S(t)=\mathrm{Tr}_E\!\left\{e^{-iHt}\bigl[\rho_S(0)\otimes \rho_E(0)\bigr]e^{+iHt}\right\}$$9. Defining

$\langle B\rangle=0$0

one finds

$\langle B\rangle=0$1

and $\langle B\rangle=0$2 is exactly the average time the pseudomodes remain excited, i.e. the reservoir’s memory time. In the Markovian limit $\langle B\rangle=0$3, $\langle B\rangle=0$4, $\langle B\rangle=0$5, and no backflow occurs (Ohyama et al., 2017).

A major structural advance is the proof of a one-to-one correspondence between HEOM and pseudomodes for exponential bath correlation functions. For every physical bath correlation function that can be written as a sum of $\langle B\rangle=0$6 exponential terms, there exists a physical model with $\langle B\rangle=0$7 interacting pseudomodes damped in Lindblad form; conversely, a non-unitary linear transformation mirrors the evolution of the system-pseudomode state onto the HEOM hierarchy, and vice versa (Müller et al., 7 Apr 2026). This correspondence also gives elegant derivations of HOPS and nuHOPS. Related work on non-Hermitian pseudomodes further shows that generalized pseudo-Lindblad equations admit a quantum-jump unraveling into trajectories of non-Hermitian states, enabling easily parallelizable Monte Carlo simulations (Menczel et al., 2024).

5. Numerical construction, scaling, and design subtleties

Recent work has shifted part of the emphasis from analytic pole identification to constructive numerical design. In coupled Lindblad pseudomode theory, one starts from a sampled spectral density $\langle B\rangle=0$8, computes the Fourier-domain bath correlation $\langle B\rangle=0$9, forms the Loewner matrix

$$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$0

takes a rank-$$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$1 truncated SVD, extracts a state-space realization $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$2, and then enforces the coupled Lindblad constraints by solving the convex SDP

$$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$3

With $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$4, one sets $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$5, $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$6, and decomposes $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$7, with $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$8. Under suitable analyticity conditions on $$C(\tau)=\langle B(t+\tau)B(t)\rangle =\frac{1}{2\pi}\int_{-\infty}^{\infty}\gamma(\omega)e^{-i\omega\tau}\,d\omega,\qquad \tau\ge 0.$$9, the minimal number $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$0 of pseudomodes required to fit $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$1 on $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$2 up to error $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$3 scales as

$$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$4

and the construction avoids the non-convex optimization required by existing approaches (Huang et al., 12 Jun 2025).

For fermionic quantum impurity models, a complementary route approximates the hybridization kernels $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$5 by sums of complex exponentials, each defining a fermionic pseudomode. Two numerical constructions are described: an analytic contour-rotation scheme and a rational approximation of the bath spectral density using the AAA algorithm, followed in both cases by interpolative decomposition. Their combination yields a pseudomode count scaling as

$$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$6

and the agreement between the two approaches suggests that the result is close to optimal (Thoenniss et al., 2024).

Pseudomode design is not merely a fitting problem. In particular, once pseudomodes couple to each other, the effective spectral density is no longer a sum of Lorentzians; non-diagonalizability of the effective single-particle non-Hermitian Hamiltonian can generate terms that cannot be obtained by diagonalizable non-Hermitian Hamiltonians. Moreover, for many uncoupled pseudomodes, the effective spectral density does not necessarily converge in the limit of an infinite number of pseudomodes, a phenomenon attributed to the non-completeness of Lorentzians as basis functions (Alford et al., 19 Sep 2025). A common misconception is therefore that arbitrarily refining an uncoupled Lorentzian grid must recover an arbitrary target bath; current analysis indicates that the issue is subtler.

6. Applications and physical interpretation

The spin-boson model remains the canonical benchmark. For

$$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$7

and the underdamped Brownian oscillator spectral density

$$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$8

the thermal correlation function splits into an analytic part $$C(\tau)=\frac{1}{\pi}\int_0^\infty J(\omega)\bigl[\coth(\beta\omega/2)\cos(\omega\tau)-i\sin(\omega\tau)\bigr]\,d\omega,$$9, which gives two pseudomodes exactly, plus an infinite Matsubara sum $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$0. One may either fit $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$1 by two additional exponentials, or map the leading Matsubara term to one pseudomode and absorb the remainder as a local $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$2-dephasing. Numerical comparisons of $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$3 against numerically exact HEOM data show that, with as few as two pseudomodes plus optional dephasing, the pseudomode master equation reproduces population dynamics with high accuracy, with errors $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$4 over many oscillation periods across regimes $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$5 and strong coupling $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$6, while reducing computational cost sharply relative to HEOM at low temperature and strong coupling (Pleasance et al., 2021).

The framework has also become a vehicle for strong-coupling thermodynamics. For bosonic baths linearly coupled to a system, exact expressions for heat, work, and average system-bath interaction energy can be rewritten in the pseudomode picture as one-time expectation values of the extended state $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$7, for example

$J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$8

This reformulation was applied to the entropy production of a driven two-level system in an Ohmic bath and to a two-bath thermal machine in which an appropriate sinusoidal modulation of the coupling with the cold bath only is enough to obtain work extraction (Albarelli et al., 2024).

Beyond reduced system dynamics, purified input-output pseudomodes provide access to environmental observables and multi-time bath correlations. They have been used to simulate non-Markovian multi-photon transfer processes on a coupled cavity waveguide system in the large time delay regime, with four purified pseudomodes exactly fitting the delayed cavity-waveguide correlation structure (Liang et al., 2024). The purified nonlinear extension carries the same logic to general nonlinear couplings and has been demonstrated for spontaneous decay inside a single-mode lossy cavity and for the resonance fluorescence spectrum of a quantum dot in the presence of a phonon environment (Zhang et al., 25 Feb 2025).

The method has also entered quantum simulation and many-body contexts. A collision-model implementation has been integrated into a quantum algorithm for simulating linear and nonlinear response functions in multidimensional electronic spectroscopy, where the key step is the pseudomode embedding of the system-environment problem into a finite Markovian master equation (Gallina et al., 2024). In a different direction, a Heisenberg-recursion construction produces a pseudomode expansion of infinite-temperature many-body autocorrelation functions, and the first few pseudomodes already give a good approximation in the quantum Ising and $J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k)$9 spin-$H_S$00 models on the square lattice (Teretenkov et al., 2024).

In this broader perspective, the pseudomode method is best understood not as a single master equation ansatz but as a family of correlation-matching embeddings. Its central claim is modest and precise: when the bath correlation structure can be represented by poles or exponential terms, a small auxiliary set of damped modes can reproduce the original non-Markovian reduced dynamics exactly or to controlled accuracy, while often exposing memory storage, backflow, thermodynamic bookkeeping, and numerical complexity in a form that is unavailable in the original continuum description.