Pseudomode Elimination in Open Quantum Systems

• Pseudomode elimination is a method that reformulates non-Markovian dynamics into a finite Markovian framework by introducing auxiliary damped modes corresponding to spectral-density poles.
• It compresses complex bath correlations into a finite sum of exponential terms, enabling efficient simulation and analytical treatment of open-system dynamics.
• The approach has practical applications in quantum thermodynamics, dissipative cavity QED, and photosynthetic systems, while addressing challenges in rational approximation and gauge stability.

Pseudomode elimination is the operation by which non-Markovian open-system dynamics is recovered from, or encoded through, an enlarged Markovian model containing a finite set of auxiliary damped modes called pseudomodes. In the canonical Garraway-type construction, a continuum bosonic environment is replaced by discrete harmonic oscillators associated with poles of the analytically continued spectral density, each coupled to a Markovian reservoir; tracing out those pseudomodes reproduces the reduced dynamics of the original system (Pleasance et al., 2021). In later work, the same expression is used more broadly for closely related reductions: replacing any eliminated sector whose influence admits a rational self-energy by damped auxiliary modes, compressing exponential bath representations to near-minimal pseudomode sets, or replacing unphysical pseudomodes by physical ensembles together with extrapolation (G. et al., 5 May 2026, Huang et al., 12 Jun 2025, Cirio et al., 2023).

1. Canonical construction from spectral-density poles

In the standard bosonic setting, the system $S$ is coupled linearly to an environment $E$ through

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

with $A=A^\dagger$. The environment is characterized by the two-time correlation function

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

where $\gamma(\omega)$ is the full spectral density associated with $B$ (Pleasance et al., 2021).

If $\gamma(\omega)$ is meromorphic in the complex plane and decays faster than $O(1/|\omega|)$ as $|\omega|\to\infty$, and if its poles in the lower half-plane are

$E$0

with residues $E$1, then contour integration yields

$E$2

Each exponential term defines a pseudomode: a discrete bosonic oscillator $E$3 with frequency $E$4, damping rate $E$5, and system coupling

$E$6

Thus, in the canonical formulation, pseudomodes are in one-to-one correspondence with poles of the analytically continued spectral density (Pleasance et al., 2021).

This pole-based interpretation has a direct physical reading. The factor $E$7 represents an oscillatory memory contribution with decay time $E$8. Slow pseudomodes, with small $E$9, encode long-lived memory; fast pseudomodes encode short-time correlations and approach Markovian behavior. A closely related formulation uses the frequency-domain self-energy $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$0: if $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$1 is rational,

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

then the corresponding memory kernel is a finite sum of exponentials and can likewise be represented by a finite pseudomode set (G. et al., 5 May 2026).

2. Markovian embedding and recovery of non-Markovian dynamics

Once the pseudomodes have been introduced, the enlarged system $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$3 is governed by a time-independent GKSL generator. In the formulation of Ref. (Pleasance et al., 2021), the auxiliary Hamiltonian is

$$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$4

and the density matrix obeys

$$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$5

Because $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$6 may be complex, $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$7 can be complex and $$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$8 is generally non-Hermitian; the dissipator cancels this non-Hermiticity and yields a completely positive Markovian semigroup on the extended space (Pleasance et al., 2021).

Pseudomode elimination, in the strict sense, is then

$$H = H_S + H_E + H_I,\qquad H_I=A\otimes B,\qquad B=\sum_k g_k(a_k+a_k^\dagger),$$9

For Gaussian environments, the reduced dynamics is completely determined by the two-time correlation function. The core equivalence argument is therefore correlation matching: in an extended auxiliary model $A=A^\dagger$0, each pseudomode is coupled to an independent Markovian reservoir $A=A^\dagger$1, and the auxiliary bath operator $A=A^\dagger$2 is constructed so that its correlation function satisfies

$A=A^\dagger$3

Under appropriate initial states, this implies

$A=A^\dagger$4

so tracing out the pseudomodes reproduces exactly the physical non-Markovian dynamics of the original model (Pleasance et al., 2021).

The same logic appears in the older Friedrichs and Jaynes–Cummings context. There, a non-local integro-differential equation for an excited-state amplitude is converted into a finite system of ordinary differential equations by introducing pseudomode amplitudes, and pseudomode elimination means solving for those amplitudes and substituting them back into the system equation, or equivalently tracing them out from the enlarged GKSL dynamics (Teretenkov, 2019). In the single-Lorentzian two-level example, the pseudomode amplitude provides a direct measure of temporarily stored excitation, and the expected time the pseudomode remains excited is

$A=A^\dagger$5

which was proposed as a simple physical meaning of the reservoir memory time (Ohyama et al., 2017).

3. Rational self-energies, generalized embeddings, and compressed mode sets

A major generalization replaces the pole expansion of a spectral density by a representability criterion in terms of response functions. In the Heisenberg-picture framework of Ref. (G. et al., 5 May 2026), one partitions the full problem into retained operators $A=A^\dagger$6 and an eliminated sector, derives a generalized Langevin equation with memory kernel $A=A^\dagger$7, and defines the self-energy

$A=A^\dagger$8

If $A=A^\dagger$9 is rational, the eliminated sector can be replaced by a finite set of damped auxiliary modes. The paper’s central claim is that pseudomode elimination is governed by rational self-energy, not by linearity: Kerr nonlinearities, cross-Kerr terms, bilinear exchange, and three-wave mixing can all be treated provided the relevant sector reduces to a rational propagator (G. et al., 5 May 2026).

This perspective supports a wider class of auxiliary constructions. In quasi-Lindblad pseudomode theory, the bath correlation function is represented directly as a complex weighted sum of complex exponentials,

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

so the global generator is trace preserving but not necessarily completely positive (Park et al., 2024). The representation is not unique: different gauge choices give different pseudomode couplings while leaving the reduced system dynamics unchanged when the global dynamics is simulated exactly. This non-uniqueness makes pseudomode elimination partly a problem of representation theory: pseudomodes can be rearranged, mixed, or effectively absorbed into one another without changing the reduced dynamics, although gauge choice can strongly affect numerical stability (Park et al., 2024).

Coupled Lindblad pseudomode theory restores CPTP dynamics while retaining compact exponential representations. Its bath correlation function has the form

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

and the main theorem states that if a quasi-Lindblad and a coupled-Lindblad representation generate the same bath correlation function, then they generate the same reduced system dynamics. When a suitable feasibility condition is satisfied, the coupled-Lindblad model uses the same number of pseudomodes as the quasi-Lindblad one, and the number of coupled pseudomodes only needs to scale as $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$2 in simulation time $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$3 and precision $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$4 (Huang et al., 12 Jun 2025). In the fermionic setting, a complementary construction based on analytic pole placement or AAA rational approximation followed by interpolative matrix decomposition yields a compressed pseudomode count

$$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$5

so pseudomode elimination becomes a controlled compression of redundant exponentials while preserving impurity observables up to the target error (Thoenniss et al., 2024).

4. Representative models and benchmark systems

The spin-boson model is a standard benchmark because its bath correlation function combines a structured analytic part with thermal Matsubara tails. For the underdamped Brownian oscillator spectral density

$$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$6

the correlation function splits as $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$7, where $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$8 is already a finite sum of exponentials and $$C(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty} d\omega\,\gamma(\omega)e^{-i\omega\tau},\qquad \tau\ge 0,$$9 is an infinite Matsubara series. In Ref. (Pleasance et al., 2021), $\gamma(\omega)$0 is captured exactly with only two pseudomodes because complex couplings are allowed, while $\gamma(\omega)$1 is treated either by two additional real exponentials or by keeping only the first Matsubara term plus a local dephasing term. Benchmarking against numerically exact HEOM for $\gamma(\omega)$2, $\gamma(\omega)$3, $\gamma(\omega)$4, $\gamma(\omega)$5, and $\gamma(\omega)$6 shows excellent agreement over both short and long times (Pleasance et al., 2021).

In photosynthetic light harvesting, the pseudomode approach has been used to model vibronic non-Markovianity in simplified Fenna–Matthews–Olson dimers. With a single dominant vibronic mode at $\gamma(\omega)$7, Huang–Rhys factor $\gamma(\omega)$8, and effective coupling $\gamma(\omega)$9, the model predicts coherence lifetimes on the $B$0 scale, compared with $B$1 in naive Markovian models, and oscillation frequencies of order $B$2 (Teretenkov, 2019). In a closely related two-level setting, the pseudomode amplitude was interpreted as the temporary storage of information and excitation in the reservoir, making pseudomode elimination a concrete picture of backflow rather than only an operator identity (Ohyama et al., 2017).

In dissipative cavity QED, a single pseudomode can encode Fano interference between direct emission and cavity-mediated emission. Ref. (Kobayashi et al., 15 Jan 2026) shows that a previously known Lindblad master equation for a two-level system plus cavity with Fano interference can be rederived from a structured-reservoir model by introducing one pseudomode. The corresponding spectral function contains a constant part and a non-Lorentzian contribution forming the Fano profile, and the constant term is essential for obtaining a Lindblad master equation (Kobayashi et al., 15 Jan 2026). This example illustrates a recurrent theme: pseudomode elimination can convert what appears as non-Markovian reservoir structure into a compact explicit mode plus simple Markovian damping.

5. Purified, nonlinear, and thermodynamic extensions

Purified pseudomode theory reinterprets pseudomodes so that positive- and negative-time bath correlations are carried by separate auxiliary modes acting only on one side of the density matrix. In the purified input-output model, the original non-Markovian bath is viewed as already having been obtained by mathematically eliminating pseudomodes, and the construction effectively undoes that elimination in a way that preserves access to environmental observables and non-Gaussian initial bath states (Liang et al., 2024). The enlarged dynamics is time-local, while tracing out the purified pseudomodes recovers the same reduced system dynamics; the same formalism also gives access to field observables such as $B$3, cavity occupations, and emission spectra (Liang et al., 2024).

The nonlinear extension replaces linear couplings $B$4 by

$B$5

For such interactions, reproducing the reduced dynamics requires not only matching the two-time function $B$6 but also the equal-time moment

$B$7

Ref. (Zhang et al., 25 Feb 2025) therefore augments the pseudomode construction with an extra zero-frequency pseudomode and then performs a Liouvillian Schrieffer–Wolff elimination of fast modes to derive a purified pseudomode master equation. The method is validated for spontaneous decay of a two-level atom in a lossy cavity and for the resonance fluorescence spectrum of a quantum dot with a phonon environment (Zhang et al., 25 Feb 2025).

In strong-coupling quantum thermodynamics, pseudomodes are retained rather than immediately eliminated, because they provide a finite Markovian representation from which bath-dependent observables can be recovered. For bosonic baths linearly coupled to the system, exact expressions for heat, work, and average system–bath interaction energy can be written in terms of bath autocorrelation functions and two-time system correlators; after the pseudomode replacement, these quantities can be evaluated efficiently in terms of one-time expectation values of the system and the pseudomodes (Albarelli et al., 2024). For example, the heat associated with bath $B$8 becomes

$B$9

so the physical continuum bath has been eliminated in favor of a finite pseudomode network and Lindblad generator (Albarelli et al., 2024).

6. Alternative meanings, numerical reduction, and limitations

Outside the canonical open-system setting, the phrase “pseudomode elimination” acquires additional technical meanings. In the many-body recursion method, a dissipatively deformed Krylov-space Liouvillian yields an autocorrelation expansion

$\gamma(\omega)$0

and elimination means truncating this expansion to the first few dominant pseudomodes. In the 2D XX and 2D Ising examples, only a few low-lying pseudomodes are needed to reproduce the correlation function accurately over the accessible time range (Teretenkov et al., 2024).

A different reduction appears in the stochastic pseudomode model, which decomposes a Gaussian bosonic environment into a finite set of zero-temperature ancillary bosonic modes plus classical stochastic fields. For rational spectral densities, all parameters can be specified analytically, and the number of ancillary quantum degrees of freedom is reduced because Matsubara and other symmetric contributions are moved into the classical field sector (Luo et al., 2023). The classical fields can even be imaginary-valued, which the authors show can decrease the entropy of the system, in contrast to real-valued fields (Luo et al., 2023). Ref. (Cirio et al., 2023) pushes this further by replacing explicitly unphysical pseudomodes with measurement results over ensembles of physical systems and an extrapolation procedure, thereby eliminating unphysical auxiliary modes as direct dynamical degrees of freedom while preserving the target reduced dynamics (Cirio et al., 2023).

The main limitations are structural. Exactness typically requires that the bath correlation function be representable as a finite sum of exponentials, equivalently that the relevant spectral density or self-energy be rational or meromorphic with finitely many relevant poles (Pleasance et al., 2021, G. et al., 5 May 2026). Thermal Matsubara tails, low-frequency singular structure, or long-time algebraic behavior generally require truncation, fitting, or hybrid stochastic treatments (Pleasance et al., 2021, Luo et al., 2023). In quasi-Lindblad formulations, the pseudomode representation is more compact but the global dynamics may lose positivity, and gauge choice becomes a stability problem rather than a merely formal one (Park et al., 2024). In coupled-Lindblad and fermionic compression schemes, the asymptotic mode count can be polylogarithmic in time and accuracy, but the construction still depends on the quality of rational approximation and on feasibility conditions or interpolation tolerances (Huang et al., 12 Jun 2025, Thoenniss et al., 2024).

Taken together, these developments establish pseudomode elimination as both a structural theorem and a practical methodology. In its narrow sense it is the tracing out of auxiliary damped modes from an enlarged GKSL evolution to recover non-Markovian reduced dynamics. In its broader modern sense it also denotes the systematic reduction, compression, purification, or physical re-realization of auxiliary-mode descriptions so that environmental memory is represented with as few effective degrees of freedom as possible, while retaining control over reduced dynamics, observables, and numerical stability (Pleasance et al., 2021, Huang et al., 12 Jun 2025).