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: convergence to the dynamics of a specified system–bath Hamiltonian as the correlation expansion, hierarchy depth, Hilbert-space cutoff, and integration accuracy are improved. A finite truncation is an approximation, not a new sense of exactness. 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.$$

Assume centered bath forces, $\langle B_\alpha\rangle=0$, with any mean absorbed into $H_S$. For a Gaussian bath with linear coupling and a specified Gaussian initial preparation, the bath influence 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].$$

The following unscaled form assumes real decay constants $\nu_k>0$, so the coefficients of $C(t)^*$ are $c_k^*$ with the same exponents. Complex-conjugate exponent pairs require the correspondingly paired coefficients in the downward term, rather than blindly conjugating $c_k$ at a fixed index.

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

$$\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).$$

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;
• the term containing $\rho_{\mathbf n+\mathbf e_k}$ feeds higher-tier correlations into the equation at tier $\mathbf n$;
• the term containing $\rho_{\mathbf n-\mathbf e_k}$ drives a higher-tier variable from a lower-tier one.

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 with a Fixed Convention

Continue to use $\hbar=1$ in this example, so energies and angular frequencies have the same units. Define

$$J(\omega)=\pi\sum_k g_k^2\delta(\omega-\omega_k),\qquad C(t)=\frac1\pi\int_0^\infty d\omega\,J(\omega) \left[\coth\!\left(\frac{\beta\omega}{2}\right)\cos\omega t-i\sin\omega t\right].$$

For the Drude–Lorentz spectrum,

$$J(\omega)=\frac{2\lambda\gamma\omega}{\omega^2+\gamma^2},\qquad \lambda=\frac1\pi\int_0^\infty\frac{J(\omega)}\omega\,d\omega,$$

$\lambda$ is the reorganization energy and $\gamma^{-1}$ is the cutoff time. The correlation expansion for $t>0$ has

$$\begin{aligned} \nu_0&=\gamma,& c_0&=\lambda\gamma[\cot(\beta\gamma/2)-i],\\ \nu_k&=\frac{2\pi k}{\beta},& c_k&=\frac{4\lambda\gamma}{\beta}\frac{\nu_k}{\nu_k^2-\gamma^2},\quad k\ge1. \end{aligned}$$

Each $c_k$ has units of energy squared. Coincident poles are evaluated by the combined limiting expression, not by separate divergent coefficients. The ideal Drude spectrum also has a short-time ultraviolet singularity; the finite-temperature expansion and any short-time correction must be handled consistently. At low temperature, the Matsubara rates $\nu_k$ are small and many terms may be needed. Padé decompositions can represent the same thermal correlation more efficiently.

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.

Converging a Calculation

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.

Vary the correlation expansion and hierarchy depth separately: convergence in one does not establish convergence in the other.

Initial States

Many introductory derivations assume a factorized initial state,

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

For the unscaled hierarchy above, a factorized preparation gives $\rho_{\mathbf n\ne0}(0)=0$; correlated preparations generally require nonzero auxiliary operators. 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.

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 auxiliary operator in the hierarchy.

References

• 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 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, “Numerically exact approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” Journal of Chemical Physics 153, 020901 (2020), doi:10.1063/5.0011599.
• 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).