Garraway's Pseudomode Construction

• Garraway’s pseudomode construction is a method that represents non-Markovian dynamics by replacing a continuous bosonic environment with a finite set of damped harmonic oscillators.
• It maps the pole structure of the bath's spectral density onto auxiliary modes so that the enlarged system obeys an exact Lindblad master equation.
• Modern generalizations extend the framework to interacting, fermionic, and nonlinear cases while addressing representational limits and design freedom.

Searching arXiv for recent pseudomode papers and Garraway-related developments. arxiv_search(query="pseudomode Garraway open quantum systems", max_results=10, sort_by="relevance") arxiv_search(query="pseudomode Garraway open quantum systems", max_results=10, sort_by="relevance") is not available in this interface, so proceeding with the provided arXiv records as the source set. Garraway’s pseudomode construction is a method for representing non-Markovian open-system dynamics by embedding the original system in a larger, but Markovian, auxiliary model. In its canonical form, a continuous structured bosonic environment is replaced by a finite set of damped harmonic oscillators—the pseudomodes—chosen so that the enlarged system obeys an exact Lindblad master equation while reproducing the same reduced dynamics as the original continuum model. In the modern formulation, the construction is exact whenever the bath correlation function is a finite sum of complex exponentials, and more general spectral densities are treated by pole extraction or exponential fitting (Pleasance et al., 2021).

1. Historical setting and core idea

Barry M. Garraway introduced the pseudomode method in the context of strongly coupled atom-reservoir decay, where a structured reservoir cannot be reduced to a memoryless decay channel. The central observation is that the non-Markovian memory kernel of the reservoir can be encoded in a small number of auxiliary damped modes, each associated with a pole of the analytically continued spectral density. This converts an integro-differential reduced equation into a time-local dynamics on an enlarged Hilbert space.

In the zero-temperature rotating-wave setting, the construction is closely connected to both the Friedrichs model and the dissipative Jaynes-Cummings model. Teretenkov emphasized that the one-excitation continuum problem can be re-expressed in a GKSL form for the system plus pseudomodes, with the continuum poles becoming discrete damped levels (Teretenkov, 2019). In a complementary single-excitation analysis, Ohyama and Tokura interpreted the pseudomode amplitudes as carriers of reservoir memory: excitation is temporarily stored in a pseudomode and can later be fed back to the system, while the inverse linewidth of the corresponding spectral feature sets the characteristic memory time (Ohyama et al., 2017).

A common misconception is that pseudomodes are literal physical modes of the environment. The later literature instead treats them as effective auxiliary degrees of freedom whose role is to reproduce the environmental correlation structure seen by the system. This suggests that the construction is best understood as an embedding principle rather than a microscopic identification of specific bath oscillators.

2. Pole structure of the spectral density

In the general bosonic setting, the environment is characterized by a one-sided spectral density

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

which is analytically continued to a meromorphic function $J(z)$ in the complex plane, with no poles in the upper half-plane and sufficiently fast decay at infinity. If its simple poles in the lower half-plane are

$z_k=\xi_k-i\lambda_k,$

with residues

$R_k=\operatorname{Res}[J(z)]_{z=z_k},$

then one may write a pole decomposition

$$J(z)=\sum_k \frac{R_k}{z-z_k} + (\text{regular terms}).$$

At zero temperature, contour integration yields the bath two-time correlation function

$$C(\tau)=\frac{1}{\pi}\int_0^\infty d\omega\,J(\omega)e^{-i\omega\tau} =-i\sum_k R_k e^{-i z_k \tau},\qquad \tau\ge 0.$$

Each pole therefore contributes one exponential memory component. Garraway’s construction assigns one bosonic pseudomode $a_k$ to each such pole, with frequency $\Omega_k=\xi_k$ and local decay rate $\gamma_k=2\lambda_k$. The enlarged Hamiltonian may be written as

$$H_0=H_S+\sum_k \Omega_k a_k^\dagger a_k + S\otimes B',$$

with

$J(z)$0

In an RWA-like Hermitian limit, one may choose the residues purely imaginary so that the pseudomode couplings are real and the non-Hermitian contributions are shifted into the dissipator (Pleasance et al., 2021).

This pole-based construction makes explicit what the pseudomodes encode: not the full bath, but the analytic structure of the bath correlation function relevant to the reduced dynamics. Later formulations retain this logic even when the pseudomodes interact, become non-Hermitian, or are generalized to fermionic or nonlinear settings.

3. Markovian embedding and conditions for exactness

Once the pseudomodes are introduced, each is coupled to a local Markovian decay channel. The density operator $J(z)$1 of the original system plus pseudomodes then obeys

$J(z)$2

The construction is exact when the pseudomode bath operator reproduces the original bath correlation function,

$J(z)$3

because, for Gaussian environments, matching the two-time correlation function guarantees the same reduced system dynamics (Pleasance et al., 2021).

The simplest exact case is therefore a bath correlation function that is already a finite sum of exponentials. More general spectral densities are handled by fitting the correlation function or a thermalized spectral density to such a sum. This is the principal route by which the original Lorentzian, zero-temperature construction has been extended to arbitrary spectral-density shapes, non-RWA couplings, multi-operator couplings, and thermal baths in principle (Pleasance et al., 2021).

Finite temperature introduces additional terms. For the underdamped Brownian oscillator,

$J(z)$4

the exact correlation function splits into an analytic part $J(z)$5 associated with poles

$J(z)$6

and an infinite Matsubara sum. In the Pleasance-Petruccione construction, $J(z)$7 is mapped exactly with two pseudomodes, while the Matsubara contribution is approximated by one or two fitted exponentials or by a local $J(z)$8 dephasing term. Benchmarks against HEOM show excellent agreement for $J(z)$9 even at strong coupling, with only a handful of pseudomodes and local truncation such as 5–10 Fock states each (Pleasance et al., 2021).

Several later works relax the auxiliary dynamics further. Menczel et al. showed that the auxiliary master equation need not preserve Hermiticity or positivity on the pseudomode sector, provided the relevant advanced and retarded auxiliary correlations match the physical bath correlations. In that generalized setting, fewer pseudomodes can suffice for underdamped finite-temperature environments, and the dynamics admits a double-Hilbert quantum-jump unraveling (Menczel et al., 2024). Cirio et al. likewise formulated thermal pseudomode models with Lindblad terms containing both $z_k=\xi_k-i\lambda_k,$0 and $z_k=\xi_k-i\lambda_k,$1, and showed that temperatures may also be handled by extra zero-frequency pseudomodes or an equivalent classical colored noise (Cirio et al., 2023).

4. Representative models and applications

The pseudomode method has been used as a concrete simulation tool across several non-Markovian settings. In the spin-boson model, the construction turns a structured bosonic bath into a finite damped-mode problem whose reduced dynamics can be benchmarked directly against numerically exact methods. In transport models, it yields a Lindblad description of excitation transfer through spin channels coupled to Lorentzian reservoirs. In strong-coupling thermodynamics, it supplies finite Markovian embeddings from which heat, work, and system-bath interaction energy can be computed. In spectroscopy and cavity QED, it provides explicit finite-state realizations of memory effects and Fano interference (Pleasance et al., 2021).

Domain	Pseudomode structure	Reported role
Spin-chain transport	One pseudomode for a Lorentzian reservoir	Exact Lindblad embedding and transport enhancement via auxiliary chains
Strong-coupling thermodynamics	Small networks of damped, possibly interacting, modes	Efficient evaluation of heat, work, and interaction energy
Multidimensional spectroscopy	One pseudomode per Lorentzian or Drude–Lorentz component	Markovian embedding for linear and nonlinear response simulation
Dissipative cavity QED	Single pseudomode plus constant background	Re-derivation of a Lindblad master equation with Fano interference

In the spin-channel model of Haddadi, Salimi, and Khorashad, diagonalization of the chain Hamiltonian reduces the system to an $z_k=\xi_k-i\lambda_k,$2-fold V-type subsystem coupled to a Lorentzian reservoir. Because the Lorentzian spectral density has a single simple pole, one introduces a single pseudomode $z_k=\xi_k-i\lambda_k,$3 with frequency $z_k=\xi_k-i\lambda_k,$4 and decay width $z_k=\xi_k-i\lambda_k,$5. The resulting master equation is of Lindblad type, and the inclusion of $z_k=\xi_k-i\lambda_k,$6 auxiliary chains enlarges a dark subspace of dimension $z_k=\xi_k-i\lambda_k,$7, with numerics showing that the sink population $z_k=\xi_k-i\lambda_k,$8 rises substantially as $z_k=\xi_k-i\lambda_k,$9 increases, especially for longer chains (Behzadi et al., 2017).

In strong-coupling thermodynamics, Trushechkin and collaborators considered linearly coupled thermal bosonic baths and replaced each bath by a small network

$R_k=\operatorname{Res}[J(z)]_{z=z_k},$0

with local Lindblad damping. They reported that a Lorentzian peak requires one or two pseudomodes, while a moderate cut-off Ohmic spectral density can be captured by 3–5 pseudomodes to $R_k=\operatorname{Res}[J(z)]_{z=z_k},$1 on times $R_k=\operatorname{Res}[J(z)]_{z=z_k},$2, and that for an Ohmic bath with $R_k=\operatorname{Res}[J(z)]_{z=z_k},$3 the entropy production of the spin-boson model converges already at $R_k=\operatorname{Res}[J(z)]_{z=z_k},$4 interacting pseudomodes (Albarelli et al., 2024).

In multidimensional electronic spectroscopy, pseudomodes are used as a preprocessing step before collision-model simulation. The bath spectral function is approximated by a finite sum of Lorentzian or Drude–Lorentz components, one pseudomode is introduced for each component, and the enlarged system follows a strictly Markovian master equation. This embedding is then integrated into a quantum algorithm for linear and nonlinear optical response functions, including spectral diffusion and relaxation along delay times (Gallina et al., 2024).

In dissipative cavity QED, Takeuchi and collaborators showed that a master equation previously derived within the Born-Markov approximation can be rederived by introducing a single pseudomode for a spectral function of the form

$R_k=\operatorname{Res}[J(z)]_{z=z_k},$5

where $R_k=\operatorname{Res}[J(z)]_{z=z_k},$6 is a constant background and $R_k=\operatorname{Res}[J(z)]_{z=z_k},$7 has one simple pole in the lower half-plane. In that construction, the constant term yields a direct Markovian decay channel, while the pole part yields the pseudomode. The final Lindblad generator includes off-diagonal dissipative terms proportional to $R_k=\operatorname{Res}[J(z)]_{z=z_k},$8, which encode Fano interference between direct radiation and cavity-mediated emission (Kobayashi et al., 15 Jan 2026).

5. Modern generalizations and formal correspondences

Subsequent work has broadened Garraway’s original framework in several directions. One major development is the use of interacting pseudomodes. When the exponential expansion of the bath correlation function contains complex weights, a set of uncoupled damped oscillators is not always sufficient, and a tridiagonal interacting network can reproduce the required interference among poles (Albarelli et al., 2024). This already moves the formalism beyond the original “one pole, one independent damped mode” picture.

A second development is the rigorous relation to hierarchical methods. Hartmann, Mascherpa, and collaborators proved that every physical bath-correlation function that can be written as a sum of $R_k=\operatorname{Res}[J(z)]_{z=z_k},$9 exponential terms can be obtained from a physical model with $$J(z)=\sum_k \frac{R_k}{z-z_k} + (\text{regular terms}).$$0 interacting pseudomodes damped in Lindblad form, and that there exists a non-unitary linear transformation mirroring the system-pseudomode state onto the HEOM hierarchy, and vice versa. In the single-Lorentzian case, the construction reduces to Garraway’s original one-pseudomode model (Müller et al., 7 Apr 2026). This places pseudomodes and HEOM on the same structural footing for exponential bath decompositions.

A third extension concerns fermionic baths. Thoenniss et al. reformulated the pseudomode problem for out-of-equilibrium fermionic quantum impurity models by approximating the bath influence kernel by sums of complex exponentials, each defining a fermionic pseudomode. Under analyticity assumptions on the bath spectral density, they derived a construction with pseudomode number scaling polylogarithmically with the maximum time $$J(z)=\sum_k \frac{R_k}{z-z_k} + (\text{regular terms}).$$1 and the approximation error $$J(z)=\sum_k \frac{R_k}{z-z_k} + (\text{regular terms}).$$2, and then reduced the mode count further using interpolative decomposition and an AAA-based rational-approximation route (Thoenniss et al., 2024). This is a direct fermionic analogue of the bosonic correlation-function expansion.

A fourth extension is the purified input-output formulation. Tamascelli and collaborators introduced “right” and “left” purified pseudomodes for the causal and anticausal parts of the bath correlation, allowing access not only to reduced system dynamics but also to environmental observables and non-Gaussian initial bath states. In their formulation, negative-time correlations are handled by a second set of purified pseudomodes with frequencies $$J(z)=\sum_k \frac{R_k}{z-z_k} + (\text{regular terms}).$$3, and the enlarged evolution is governed by a Lindblad-like equation acting one-sidedly on pure pseudomode states (Liang et al., 2024).

A fifth extension abandons linear auxiliary structure as the defining principle. In a circuit-QED setting, Chia and coauthors argued that pseudomode elimination is not fundamentally tied to linearity but to representability: any eliminated sector whose influence on the retained subsystem admits a rational self-energy can be replaced by a finite set of damped auxiliary modes. They formulated this through a Heisenberg-picture Dyson equation and demonstrated closed-form elimination for Kerr-coupled and three-wave-mixing systems (G. et al., 5 May 2026). This suggests a broader “rational self-energy” interpretation of the construction.

6. Subtleties, limitations, and recurrent misconceptions

The main limitation of the pseudomode method is representational rather than conceptual. Exactness requires that the bath correlation function be representable as a finite sum of exponentials, or at least be well approximated by such a sum. For environments that cannot be captured by a handful of poles or exponentials, one needs many pseudomodes or alternative fitting strategies, and the computational cost rises correspondingly (Pleasance et al., 2021).

Recent work has shown that pseudomode design itself contains substantial freedom and several nontrivial pitfalls. de Vega, Zanardi, and collaborators studied coupled pseudomodes and showed that once pseudomodes couple to each other, the effective spectral density is generally no longer a simple sum of Lorentzians. If the effective single-particle non-Hermitian pseudomode Hamiltonian is defective, Jordan blocks generate polynomially weighted exponentials in time and corresponding spectral terms such as squared Lorentzians. They also gave an algebraic inversion procedure for matching a fitted spectral density exactly, and emphasized that the design freedom is enormous (Alford et al., 19 Sep 2025).

The same work identifies a widely held but incorrect expectation: increasing the number of uncoupled Lorentzian pseudomodes on a finer and finer grid does not necessarily produce convergence to the target spectral density. In their analysis, the effective spectral density oscillates about the target because Lorentzians are not a complete basis in the relevant sense; even the infinite-uncoupled-modes limit can fail to converge pointwise (Alford et al., 19 Sep 2025). This is a technical warning against equating “more Lorentzians” with systematic convergence.

Another recurrent misconception is that the auxiliary model must itself be directly physical. Menczel et al. explicitly described pseudomodes as unphysical in general, and Cirio et al. showed that a direct experimental implementation of the unphysical pseudomode model is, in general, impossible. They nevertheless demonstrated that its effects can be reproduced from measurement results over ensembles of physical systems involving ancillary harmonic modes and, optionally, a stochastic driving field (Menczel et al., 2024, Cirio et al., 2023). The same point appears in purified and non-Hermitian constructions: the auxiliary modes need not correspond to literal bath excitations for the reduced dynamics to be exact.

Taken together, these developments recast Garraway’s construction as a general embedding framework for non-Markovian dynamics. In its original form, each pole of a meromorphic spectral density becomes one damped harmonic oscillator. In its modern form, the same principle supports interacting, non-Hermitian, purified, fermionic, and nonlinear auxiliary models, provided the relevant bath influence admits a finite rational or exponential representation.