🔬 mesoscale physics

Improving the Stability of the Hierarchical Equations of Motion for Open Quantum Systems with Strong Coupling to Structured Bosonic Baths

This paper addresses the numerical instability of the Hierarchical Equations of Motion (HEOM) method for open quantum systems with strong coupling to bosonic baths by applying a non-unitary similarity transformation that balances hierarchy-raising and lowering terms, thereby enabling stable simulations of much stronger system-bath interactions.

Original authors: Salvatore Gatto, Samuel L. Rudge, Bokang Hou, Eran Rabani, Michael Thoss

Published 2026-09-29

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Salvatore Gatto, Samuel L. Rudge, Bokang Hou, Eran Rabani, Michael Thoss

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Quantum systems are rarely isolated islands. Whether it is an electron moving through a crystal, a molecule reacting in a solution, or a qubit in a quantum computer, these tiny entities are constantly interacting with their surroundings. This environment, often made of vibrating atoms or electromagnetic waves, acts as a heat bath that can drain energy from the system or scramble its delicate quantum properties. Understanding how a quantum system behaves while constantly exchanging energy and information with this noisy background is one of the central challenges in modern physics. It is essential for designing better quantum computers, creating more efficient solar cells, and understanding how energy moves through biological molecules.

To study these interactions, scientists use mathematical models that describe the system and its environment as a single, evolving whole. One of the most powerful tools for this is a method called the hierarchical equations of motion. This approach breaks down the complex influence of the environment into a series of layers, or a hierarchy, where each layer accounts for a deeper level of interaction between the system and the bath. By solving these layers step by step, researchers can simulate how the system changes over time with extreme precision. However, because the environment is theoretically infinite, scientists must stop the calculation at a certain depth. This necessary cut-off has long been a source of trouble, causing the calculations to become unstable and produce nonsensical results, especially when the system interacts strongly with the environment or when the environment has a complex, structured nature.

In a recent study, researchers set out to fix this instability. They focused on a specific type of environment made of bosons, which are particles that can occupy the same state in large numbers, such as the vibrations in a solid material. The team realized that the mathematical structure of the standard method contained an imbalance. In the hierarchy of layers, there were terms that constantly pushed the calculation to higher, more complex levels, but the terms that pulled them back down were not strong enough to counteract this push. When the calculation is forced to stop at a finite depth, this imbalance causes the numbers to grow uncontrollably, leading to a numerical explosion that ruins the simulation. The researchers identified this as a non-normal amplification, a phenomenon where the mathematical operator driving the system behaves in a way that allows small errors to swell into massive distortions, even if the underlying physics is stable.

To solve this, the team applied a mathematical transformation directly to the auxiliary layers of the hierarchy. Instead of changing the physical model of the system or the environment, they changed the way the layers were represented. They introduced a weighting scheme that effectively rebalanced the equation. This process converted the unbalanced, one-way push toward higher layers into a balanced combination of pushing up and pulling down. By doing so, they suppressed the artificial amplification that caused the instability. The result was a new, stabilized version of the equations that could handle much stronger interactions and more complex environments without breaking down.

The researchers tested their new method using a model system known as a spin-boson model, which describes a simple two-state system interacting with a bath of vibrations. First, they used a standard, smooth type of environment called a Brownian oscillator. In these simulations, the traditional method failed as the interaction strength increased, with the calculated population of the system states diverging wildly after a short time. The new, stabilized method, however, produced smooth, physically realistic results that remained stable even when the interaction was very strong. They also tested the method on a more difficult case: a structured environment with a complex pattern of frequencies, similar to what is found in real quantum dots. Here, the traditional method became unstable even at moderate interaction strengths, while the stabilized approach continued to work perfectly, accurately capturing the complex, memory-driven dynamics of the system.

A key finding of the study was that simply making the calculation deeper did not solve the problem. In the traditional method, increasing the number of layers often made the instability worse, causing the simulation to fail even sooner. This counterintuitive result showed that the issue was not just a lack of precision, but a fundamental flaw in how the finite calculation was structured. The new method, by contrast, showed clear convergence, meaning that as the calculation depth increased, the results became more accurate and stable. This allowed the researchers to simulate systems with reorganization energies as high as 100 millielectronvolts, a regime where the standard method would have failed completely.

The work demonstrates that the stability of a numerical simulation can depend heavily on how the equations are written, even if the underlying physics remains the same. By reorganizing the mathematical representation of the hierarchy, the team created a tool that is robust enough to handle the strong couplings and complex environments found in real-world quantum materials. This advancement opens the door to more reliable simulations of quantum dynamics in condensed-phase chemistry and nanoscale materials, where the interaction between a system and its environment is often strong and highly structured. The stabilized approach provides a way to explore these regimes with confidence, ensuring that the results reflect the true behavior of the quantum world rather than artifacts of the calculation method.

Technical Summary: Improving the Stability of the Hierarchical Equations of Motion for Open Quantum Systems with Strong Coupling to Structured Bosonic Baths

Problem Statement
The Hierarchical Equations of Motion (HEOM) approach is a leading method for simulating the non-Markovian dynamics of open quantum systems, offering a numerically exact representation of the Feynman-Vernon influence functional. However, the method relies on a finite truncation of an inherently infinite hierarchy of auxiliary density operators (ADOs). For bosonic environments, this finite truncation can induce fundamental numerical instabilities, particularly in regimes of strong system-bath coupling or when dealing with structured spectral densities. These instabilities manifest as unphysical divergences in the dynamics, often preventing simulations even when the truncation tier ( $n_{max}$ ) is increased. Previous attempts to resolve this, such as increasing the truncation depth or using projection methods, have proven insufficient for certain spectral densities, and sophisticated pole decomposition schemes for structured environments can exacerbate the issue.

Methodology
The authors address these instabilities by analyzing the algebraic structure of the HEOM generator within an auxiliary Fock-space picture.

Fock-Space Representation: The hierarchy is reformulated using bosonic creation and annihilation operators acting on an auxiliary Fock space. This representation reveals that the HEOM generator is a non-normal operator.
Identification of Instability Mechanism: The analysis identifies the source of instability as an imbalance between damping terms (which suppress auxiliary occupations) and "hierarchy-raising" terms (which transport weight to higher tiers). In a finite truncation, the removal of couplings to tiers beyond $n_{max}$ leaves unpaired hierarchy-raising terms. This creates a non-normal algebraic structure capable of generating strong transient amplification (pseudospectral growth), even when the eigenvalues of the generator have non-positive real parts.
Stabilization via Similarity Transformation: To suppress this amplification, the authors apply a non-unitary similarity transformation directly in the auxiliary Fock space. This transformation utilizes a Gaussian weighting operator, $e^{-S}$ , which contracts the relevant part of the Fock space.
Bogoliubov Rotation: Following the Gaussian transformation, a Bogoliubov-like rotation is applied to the transformed hierarchy. This reorganizes the operators such that the previously unbalanced hierarchy-raising terms are converted into balanced combinations of hierarchy-raising and hierarchy-lowering terms.
Algorithmic Implementation: The resulting stabilized equations of motion are derived for the transformed ADOs. The physical reduced density matrix is recovered by transforming the results back to the original representation using analytically evaluated matrix elements.

Key Contributions

Mechanistic Insight: The paper provides a rigorous explanation for finite-truncation instabilities in bosonic HEOM, attributing them to the non-normal nature of the generator and the specific imbalance of raising/lowering terms at the truncation boundary.
Stabilized Formalism: The authors introduce a stabilized HEOM formulation that employs a Gaussian similarity transformation and a subsequent Bogoliubov rotation. This approach fundamentally alters the finite-dimensional representation of the generator without changing the exact dynamics of the untruncated system.
Algorithmic Framework: A complete algorithmic procedure is presented, including the definition of initial conditions in the transformed space, the propagation of the stabilized hierarchy, and the reconstruction of physical observables.

Results
The efficacy of the stabilized approach is demonstrated through two case studies using a spin-boson model:

Brownian Oscillator Spectral Density: Simulations show that while the standard HEOM becomes unstable at moderate to strong reorganization energies ( $\Lambda$ ), the transformed HEOM remains stable and converges. Notably, increasing the truncation tier ( $n_{max}$ ) in the standard approach does not guarantee stability and can sometimes lead to earlier divergence due to enhanced pseudospectral sensitivity. In contrast, the transformed approach exhibits clear convergence with increasing $n_{max}$ .
Structured Spectral Density: For a more complex spectral density featuring both broad low-frequency and sharp high-frequency components (representing nonradiative loss pathways in quantum dots), the standard HEOM fails at weak-to-moderate coupling strengths. The structured nature of the bath introduces long-lived memory effects that severely destabilize the standard truncation. The transformed HEOM, however, successfully simulates the dynamics across all tested coupling strengths, maintaining stability where the standard method diverges.

Significance and Claims
The paper claims that the proposed stabilization technique substantially broadens the applicability of the HEOM approach. By suppressing transient non-normal amplification, the method enables numerically exact simulations of open quantum systems in regimes characterized by:

Strong system-bath coupling.
Highly structured spectral densities.
Pronounced non-Markovian effects.

The authors emphasize that this improvement is particularly critical for modeling experimentally relevant systems in quantum technologies, condensed-phase chemistry, and nanoscale materials, where standard truncation schemes often fail. The work supports the broader theoretical view that non-Markovian dynamics can admit multiple equivalent extended-space embeddings, where the choice of representation significantly impacts numerical stability under finite truncation. The stabilized hierarchy offers a robust representation specifically designed to mitigate the transient amplification issues inherent in standard bosonic HEOM implementations.

Technical Summary: Improving the Stability of the Hierarchical Equations of Motion for Open Quantum Systems with Strong Coupling to Structured Bosonic Baths

More like this