Hierarchical Equations of Motion

Hierarchical equations of motion, usually abbreviated HEOM, are a family of auxiliary-density-operator methods for simulating open quantum systems coupled to structured Gaussian environments. The method rewrites bath memory as a hierarchy of coupled equations. The physical reduced density matrix is the lowest member of the hierarchy, and higher auxiliary density operators encode progressively deeper system-bath memory.

HEOM is widely used in chemical physics, molecular exciton dynamics, condensed-phase spectroscopy, spin-boson models, quantum thermodynamics, and other regimes where weak-coupling Markovian master equations can be unreliable.

The phrase “numerically exact” is common in this literature, but it has a precise meaning: exact within a specified system-bath Hamiltonian, bath correlation expansion, hierarchy truncation, Hilbert-space cutoff, and time-integration tolerance. It does not mean exact for arbitrary environments or automatically converged without checks.

System-Bath Setting

A typical HEOM starting point is a system linearly coupled to one or more harmonic baths:

$$H = H_S + \sum_\alpha V_\alpha\otimes B_\alpha + H_B.$$

For a Gaussian bath, the influence of $B_\alpha$ is determined by two-point correlation functions,

$$C_{\alpha\beta}(t) = \langle B_\alpha(t)B_\beta(0) \rangle_B.$$

The essential HEOM step is to represent each relevant bath correlation as a sum of exponentials:

$$C(t) \approx \sum_{k=0}^{K} c_k e^{-\nu_k t}, \qquad t\ge0.$$

The coefficients $c_k$ may be complex. Their real parts encode fluctuations, and their imaginary parts encode dissipative response and bath-induced phase shifts. The decay constants $\nu_k$ set memory times.

Auxiliary Density Operators

For one coupling operator $V$ and one correlation expansion, introduce a multi-index

$$\mathbf n = (n_0,n_1,\ldots,n_K), \qquad n_k\in\{0,1,2,\ldots\}.$$

The physical reduced density matrix is

$$\rho_{\mathbf 0}(t), \qquad \mathbf 0=(0,0,\ldots,0).$$

The other objects

$$\rho_{\mathbf n}(t)$$

are auxiliary density operators, often called ADOs. They are not ordinary density matrices. They need not be positive or normalized. Their role is to store the bath-memory information required to evolve $\rho_{\mathbf 0}$.

The tier of an ADO is

$$|\mathbf n| = \sum_{k=0}^{K} n_k.$$

Higher tiers represent higher-order memory contributions.

Canonical One-Bath Form

Conventions vary, but a common compact form is easiest to see with $\hbar=1$ and

$$\mathcal L_S X = -i[H_S,X].$$

Let $\mathbf e_k$ be the multi-index with $1$ in the $k$th position and $0$ elsewhere. Then a standard HEOM form is

$$\begin{aligned} \dot\rho_{\mathbf n} =& \left( \mathcal L_S - \sum_k n_k\nu_k \right) \rho_{\mathbf n} \\ & -i \left[ V, \sum_k \rho_{\mathbf n+\mathbf e_k} \right] \\ & -i \sum_k n_k \left( c_k V\rho_{\mathbf n-\mathbf e_k} - c_k^*\rho_{\mathbf n-\mathbf e_k}V \right). \end{aligned}$$

Terms with negative multi-index entries are omitted. The equation for $\rho_{\mathbf 0}$ couples upward to first-tier ADOs. Those ADOs couple upward and downward, creating a hierarchy.

This equation is best read structurally:

• $\mathcal L_S$ gives coherent system evolution;
• $\sum_k n_k\nu_k$ damps higher-tier memory variables;
• upward coupling creates memory variables driven by the system operator $V$;
• downward coupling feeds stored bath correlations back into lower tiers.

Different HEOM variants include multiple baths, multiple system coupling operators, fermionic reservoirs, scaled ADOs, Matsubara or Padé expansions, terminators, and correlated initial states. The core idea remains the same.

Drude–Lorentz Bath

A common bosonic spectral density is the Drude–Lorentz form,

$$J(\omega) = \frac{ 2\lambda\gamma\omega } { \omega^2+\gamma^2 }, \qquad \omega>0.$$

Here $\lambda$ is a reorganization-energy scale and $\gamma^{-1}$ is a bath correlation time. In one common convention, the bath correlation can be written

$$C(t) = \sum_{k=0}^{\infty} c_k e^{-\nu_k t}.$$

The first decay rate is

$$\nu_0=\gamma,$$

and the Matsubara rates are

$$\nu_k = \frac{2\pi k}{\beta\hbar}, \qquad k\ge1.$$

The corresponding coefficients depend on the normalization of $J(\omega)$. A typical convention gives

$$c_0 = \lambda\gamma \left[ \cot \left( \frac{\beta\hbar\gamma}{2} \right) -i \right],$$

with thermal correction terms for $k\ge1$. At low temperature, the Matsubara terms decay slowly and many terms may be needed unless a Padé or other optimized decomposition is used.

The important point is not this one convention. HEOM needs a controlled expansion of the bath correlation function.

Why the Hierarchy Works

A non-Markovian reduced equation often involves a memory integral. HEOM avoids evaluating a growing history integral directly by promoting each exponential memory component to a dynamical variable.

For a scalar analogy, if

$$y_k(t) = \int_0^t ds\, e^{-\nu_k(t-s)}x(s),$$

then

$$\dot y_k(t) = x(t)-\nu_k y_k(t).$$

The memory integral has become a time-local auxiliary equation. HEOM does the quantum operator version of this idea, with additional tiers required because system operators do not commute and bath fluctuations can act repeatedly.

Truncation

The exact hierarchy is infinite. In computation one chooses a maximum tier $N_{\max}$ and discards, approximates, or terminates higher tiers.

If there are $M$ exponential components in the hierarchy and the maximum tier is $N_{\max}$, the number of multi-indices through that tier is

$$N_{\mathrm{ADO}} = \binom{M+N_{\max}}{N_{\max}}.$$

Each ADO is a system operator. If the system Hilbert space dimension is $d$, storing all ADOs scales like

$$N_{\mathrm{ADO}}\,d^2$$

complex numbers, before accounting for multiple baths, sparse structure, symmetries, or parallelization.

This combinatorial growth is the main cost of HEOM. Low temperatures, slow baths, strong coupling, and many system-bath coupling channels all increase the required hierarchy size.

Convergence Checklist

A credible HEOM calculation should report or test:

• the system-bath Hamiltonian and coupling operators;
• the spectral density and temperature;
• the exponential decomposition of $C(t)$;
• the number of Matsubara, Padé, or fitted exponential terms;
• the hierarchy depth $N_{\max}$;
• the truncation or terminator scheme;
• the system Hilbert-space truncation, if any;
• time-step and integrator tolerances;
• convergence of observables under increasing $K$ and $N_{\max}$;
• whether the initial state is factorized, correlated, or prepared by imaginary-time HEOM.

Without these checks, “HEOM result” is not yet a reproducible statement.

Initial States

Many introductory derivations assume a factorized initial state,

$$\rho_{\mathrm{tot}}(0) = \rho_S(0)\otimes\rho_B.$$

That assumption is often useful for nonequilibrium preparation, but it is not always physical at strong coupling. HEOM can also be formulated in imaginary time or extended schemes to prepare correlated thermal states. This matters for thermodynamics and equilibrium comparisons, where the correct reduced equilibrium may involve the Hamiltonian of mean force rather than a bare Gibbs state.

For the conceptual issue, see Initial Correlations and Reaction-Coordinate Mapping.

Relation to Other Methods

HEOM, pseudomodes, and reaction-coordinate mappings all make memory explicit, but they do it differently.

Method	Memory carrier	Typical strength
HEOM	auxiliary density operators	Gaussian baths with exponential correlation expansions
Pseudomodes	damped effective modes	Lorentzian spectra and pole-dominated reservoirs
Reaction coordinate	explicit collective bath coordinate	strong coupling, structured spectra, thermodynamics
Memory kernels	time-history integral	formal reduced descriptions and projection methods

HEOM is often more systematic for Gaussian harmonic baths, while reaction-coordinate and pseudomode methods may give more physical intuition in terms of explicit modes.

Common Mistakes

Treating ADOs as density matrices

Only $\rho_{\mathbf 0}$ is the physical reduced density matrix. Higher-tier ADOs are bookkeeping operators for bath memory and need not be positive or normalized.

Calling an unconverged run numerically exact

The method may be exact in the limit of complete correlation expansion and infinite hierarchy, but a finite calculation must be converged.

Mixing spectral-density conventions

Factors of $\pi$, $\hbar$, and $2$ differ across communities. The coefficients $c_k$ must match the chosen definition of $J(\omega)$ and $C(t)$.

Using too few low-temperature terms

Low temperature makes Matsubara or equivalent correction terms more important. A hierarchy that works at high temperature may fail at low temperature.

Ignoring initial correlations

Strong system-bath coupling can make factorized initial states inappropriate for equilibrium questions. Correlated preparation must be specified.

Applying HEOM to the wrong bath class

Standard bosonic HEOM assumes Gaussian harmonic environments. Non-Gaussian spin baths, nonlinear environments, and strongly driven reservoirs may require different methods or generalized hierarchies.

Exercises

Counting ADOs

Suppose a hierarchy has $M$ exponential components and is truncated at tier $N_{\max}$. Show that the number of ADOs through that tier is

$$\binom{M+N_{\max}}{N_{\max}}.$$

Solution

An ADO is labeled by a multi-index

$$\mathbf n=(n_1,\ldots,n_M)$$

with nonnegative integer entries. Keeping tiers through $N_{\max}$ means

$$n_1+\cdots+n_M\le N_{\max}.$$

Introduce one slack variable

$$n_{M+1} = N_{\max} - \sum_{j=1}^M n_j.$$

Then the problem is to count nonnegative integer solutions of

$$n_1+\cdots+n_M+n_{M+1} = N_{\max}.$$

By the stars-and-bars count, this is

$$\binom{M+N_{\max}}{N_{\max}}.$$

Two-Component Hierarchy

For two exponential components, list the multi-indices through tier $2$.

Solution

The multi-index is $(n_0,n_1)$. Tier $0$ gives

$$(0,0).$$

Tier $1$ gives

$$(1,0), \qquad (0,1).$$

Tier $2$ gives

$$(2,0), \qquad (1,1), \qquad (0,2).$$

There are $6$ ADOs total, matching

$$\binom{2+2}{2}=6.$$

Why Exponentials Help

Let

$$y(t) = \int_0^t ds\, e^{-\nu(t-s)}x(s).$$

Derive the time-local equation for $y(t)$.

Solution

Differentiate using the upper limit and the derivative of the kernel:

$$\dot y(t) = x(t) + \int_0^t ds\, \left[ -\nu e^{-\nu(t-s)} \right]x(s).$$

The integral is $y(t)$, so

$$\dot y(t) = x(t)-\nu y(t).$$

This is the scalar version of the HEOM idea: exponential memory kernels can be represented by auxiliary time-local variables.

ADO Interpretation

Why should one not demand that a first-tier ADO be positive semidefinite?

Solution

A first-tier ADO is not a physical state of the system. It is an auxiliary operator encoding system-bath correlation information associated with one exponential component of the bath memory. Positivity is required for the physical reduced density matrix $\rho_{\mathbf 0}$, not for every bookkeeping operator in the hierarchy.

Cross-Links

• Non-Markovian Dynamics
• Memory Kernels
• Noise Spectra
• Caldeira–Leggett Model
• Pseudomode Methods
• Reaction-Coordinate Mapping
• System-Bath Hamiltonians
• Initial Correlations
• Approximation Checklist

References

• Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath,” Journal of the Physical Society of Japan 58, 101-114 (1989).
• Y. Tanimura, “Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath,” Physical Review A 41, 6676-6687 (1990).
• Y. Tanimura, “Stochastic Liouville, Langevin, Fokker–Planck, and master equation approaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006).
• A. Ishizaki and G. R. Fleming, “Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach,” Journal of Chemical Physics 130, 234111 (2009).
• J. Strümpfer and K. Schulten, “Open quantum dynamics calculations with the hierarchy equations of motion on parallel computers,” Journal of Chemical Theory and Computation 8, 2808-2816 (2012).
• Y. Tanimura, “Numerically exact approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” Journal of Chemical Physics 153, 020901 (2020).