---
title: Hierarchy of Pure States (HOPS)
url: https://www.emergentmind.com/topics/hierarchy-of-pure-states-hops
type: topic
---

# Hierarchy of Pure States (HOPS)

The Hierarchy of Pure States (HOPS) is a wavefunction-based method for open quantum system dynamics in non-Markovian structured environments. Introduced by D. Suess, A. Eisfeld, and W. T. Strunz, it reformulates the non-Markovian stochastic Schrödinger equation as a hierarchy of stochastic evolution equations for pure states, or quantum trajectories, such that the exact reduced density operator is recovered by ensemble averaging the lowest-tier state [1402.4647]. In its original formulation, HOPS was demonstrated for the Spin–Boson model, the calculation of absorption spectra of molecular aggregates, and energy transfer in a photosynthetic pigment-protein complex [1402.4647].

## 1. Origin in non-Markovian quantum state diffusion

HOPS starts from a standard system–bath Hamiltonian,
$$
H_{\rm tot}=H_{\rm sys}\otimes \mathbb{1}_B+\mathbb{1}_{\rm sys}\otimes H_B+H_{\rm int},
$$
with
$$
H_B=\sum_\lambda \omega_\lambda a^\dagger_\lambda a_\lambda,\qquad
H_{\rm int}=\sum_\lambda (g_\lambda L\otimes a^\dagger_\lambda+g_\lambda^* L^\dagger\otimes a_\lambda).
$$
The bath is characterized by a spectral density \(J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)\) and, at temperature \(T\), by the two-point correlation function
$$
\alpha(t-s)=\int_0^\infty d\omega\,J(\omega)\Big[\coth\!\big(\omega/2k_BT\big)\cos[\omega(t-s)]-i\sin[\omega(t-s)]\Big].
$$
Assuming an initially factorized state \(|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B\), the reduced density operator can be obtained exactly from stochastic pure states \(|\psi_t(z^*)\rangle\) as
$$
\rho(t)=\mathbb{E}\big[\,|\psi_t(z^*)\rangle\langle\psi_t(z^*)|\,\big],
$$
where \(z_t\) is a complex Gaussian process satisfying
$$
\mathbb{E}[z_t]=0,\qquad \mathbb{E}[z_t z_s]=0,\qquad \mathbb{E}[z_t z_s^*]=\alpha(t-s).
$$
The corresponding non-Markovian stochastic Schrödinger equation is [1402.4647]
$$
\partial_t|\psi_t\rangle=
\Bigl(-iH_{\rm sys}+Lz_t^*-L^\dagger\int_0^t ds\,\alpha(t-s)\frac{\delta}{\delta z_s^*}\Bigr)|\psi_t\rangle.
$$

The central obstruction to direct simulation is the functional derivative \(\delta/\delta z_s^*\). HOPS is the procedure that removes this obstacle by replacing the functional-derivative structure with a closed hierarchy of auxiliary pure states. This places HOPS within the NMQSD framework while turning the formal exactness of the unraveling into a numerically tractable hierarchy [1402.4647].

## 2. Construction of the hierarchy and the representation of memory

For a single-exponential bath kernel,
$$
\alpha(\tau)=g\,e^{-w\tau}\quad (\tau\ge 0),\qquad w=\gamma+i\Omega,
$$
one defines
$$
D_t=\int_0^t ds\,\alpha(t-s)\frac{\delta}{\delta z_s^*},
\qquad
|\psi_t^{(k)}\rangle:=(D_t)^k|\psi_t\rangle,\quad k=0,1,2,\dots .
$$
The hierarchy then becomes
$$
\partial_t|\psi_t^{(k)}\rangle
=
\Bigl(-iH_{\rm sys}-k\,w+Lz_t^*\Bigr)|\psi_t^{(k)}\rangle
+k\,\alpha(0)\,L\,|\psi_t^{(k-1)}\rangle
-L^\dagger |\psi_t^{(k+1)}\rangle,
$$
with initial conditions \(|\psi_0^{(0)}\rangle=|\psi_0\rangle\) and \(|\psi_0^{(k>0)}\rangle=0\). The physical state is \(|\psi_t\rangle=|\psi_t^{(0)}\rangle\) [1402.4647].

For a sum of exponentials,
$$
\alpha(\tau)=\sum_{j=1}^J g_j e^{-w_j \tau},
$$
HOPS introduces a multi-index \(\mathbf{k}=(k_1,\dots,k_J)\) and auxiliary states
$$
|\psi_t^{\mathbf{k}}\rangle:=\prod_{j=1}^J (D_{j,t})^{k_j}|\psi_t\rangle,
$$
leading to
$$
\partial_t|\psi_t^{\mathbf{k}}\rangle
=
\Bigl(-iH_{\rm sys}-\sum_{j=1}^J k_j w_j+\sum_{j=1}^J L z_{j,t}^*(t)\Bigr)|\psi_t^{\mathbf{k}}\rangle
+\sum_{j=1}^J k_j \alpha_j(0)L|\psi_t^{\mathbf{k}-\mathbf{e}_j}\rangle
-\sum_{j=1}^J L^\dagger |\psi_t^{\mathbf{k}+\mathbf{e}_j}\rangle .
$$
Here \(\gamma_j={\rm Re}\,w_j\), \(\Omega_j={\rm Im}\,w_j\), and the noises satisfy \(\mathbb{E}[z_j(t)z_k(s)]=0\) and \(\mathbb{E}[z_j^*(t)z_k(s)]=\delta_{jk}\alpha_j(t-s)\) [1402.4647].

The multi-index carries the environmental memory. In the original formulation, each \(k_j\) counts how many times the corresponding “history operator” \(D_{j,t}\) has acted; retaining larger weights in \(\mathbf{k}\) means that the stochastic state probes deeper into the bath memory kernel [1402.4647]. This makes HOPS an explicit memory hierarchy rather than a Markovian embedding in density-operator space.

A non-linear, normalized version is obtained via a Girsanov transform. In that form, the hierarchy retains the same coupling structure but acquires drift terms containing the instantaneous expectation \(\langle L\rangle_t\), and it is designed to improve Monte Carlo convergence [1402.4647].

## 3. Exact reconstruction, normalization, and observables

In the linear formulation, the reduced state is reconstructed as
$$
\rho(t)=\mathbb{E}\big[\,|\psi_t^{\mathbf{0}}\rangle\langle\psi_t^{\mathbf{0}}|\,\big].
$$
In the non-linear formulation, one evolves normalized states \(|\tilde\psi_t^{\mathbf{0}}\rangle\) and still has
$$
\rho(t)=\mathbb{E}\big[\,|\tilde\psi_t^{\mathbf{0}}\rangle\langle\tilde\psi_t^{\mathbf{0}}|\,\big].
$$
Thus, in either case, the exact non-Markovian reduced dynamics is obtained from the ensemble average over sufficiently many noise realizations [1402.4647].

This exactness should be interpreted carefully. HOPS is formally exact for the infinite hierarchy; numerical implementations require a finite hierarchy depth, a finite exponential representation of the bath correlation function, and Monte Carlo sampling. A common misconception is therefore to identify a specific truncated calculation with exact dynamics without convergence checks. The literature instead treats HOPS as systematically improvable, with truncation order and trajectory number increased until observables stabilize [1402.4647; 1802.04530].

The normalized, non-linear hierarchy is particularly important in strong-coupling regimes. For the spin-boson problem with sub-Ohmic spectral density, the non-linear variant functions as an importance-sampling mechanism and is described as necessary in the strong-coupling regime; for non-zero temperature, it is reported to be highly favorable to use the zero temperature bath correlation function and include temperature via a stochastic Hermitian contribution to the system Hamiltonian [1710.08268].

Later work extended HOPS beyond reduced-system observables. Bath energy change and interaction energy can be computed directly from HOPS by expressing them through stochastic averages of matrix elements involving first-tier auxiliary states \(|\psi^{\mathbf{e}_\mu}\rangle\), which enables a fully quantum dynamical treatment of energetic contributions in a strongly coupled quantum heat engine [2402.06039]. Multitime correlation functions can likewise be formulated through pure-state decompositions propagated under the same noise realization, and this has been used for absorption and resonance fluorescence spectra in quantum-dot models with phonon coupling [2412.20598].

## 4. Truncation schemes, convergence control, and computational scaling

Practical HOPS calculations truncate the hierarchy at total order \(|\mathbf{k}|=K\). For a single exponential, a common terminator is
$$
|\psi_t^{(K+1)}\rangle \simeq \frac{\alpha(0)}{w}L^\dagger |\psi_t^{(K)}\rangle,
$$
and analogous expressions are used for multi-exponential triangular truncation [1402.4647]. Convergence is tested by increasing \(K\) until observables such as populations, coherences, or spectra change negligibly. A posteriori, comparison between truncation orders \(K\) and \(K+1\) provides an upper bound on the truncation error because the method is systematically improvable [1402.4647].

The combinatorics of the auxiliary set is central to the numerical cost. For \(|\mathbf{k}|\le K\), the number of auxiliary states is
$$
N_{\rm aux}=C(K+J,J),
$$
and the overall dense cost per time step scales as \(\sim O(M\,N_{\rm aux}\,d^2)\), where \(M\) is the number of independent noise realizations and \(d\) is the system dimension. The non-linear hierarchy typically requires far fewer trajectories \(M\) because of importance sampling [1402.4647].

The limitations of standard triangular truncation motivated alternative truncation strategies. The \(n\)-particle approximation (\(n\)PA) retains only those multi-indices for which the number of distinct sites \(\ell\) with \(\sum_j k_{\ell j}>0\) is at most \(n\), while the \(n\)-mode approximation (\(n\)MA) retains only those multi-indices for which the total number of modes \((\ell,j)\) with \(k_{\ell j}>0\) is at most \(n\) [1802.04530]. These schemes were introduced to make convergence checks numerically feasible, since the jump in equation count from one triangular depth to the next can be very large [1802.04530].

Benchmark calculations in absorption and excitation transfer showed that \(n\)MA and \(n\)PA can be combined with moderate depth \(D\) to generate a ladder of hierarchy sizes that grows slowly enough for practical convergence checks [1802.04530]. This suggests that truncation in HOPS is not a single prescription but a family of controllable approximations adapted to bath structure and localization.

## 5. Adaptive, dyadic, and tensor-network variants

Adaptive HOPS (adHOPS) exploits the locality of each trajectory by constructing time-dependent reduced bases of molecular site states \(S_t\) and auxiliary-index wavefunctions \(A_t\), chosen such that the local truncation error in the time derivative satisfies
$$
\|\partial \Phi_{\rm full}-\tilde\partial \Phi_{\rm reduced}\|<\delta,
\qquad
\delta^2=\delta_S^2+\delta_A^2.
$$
In dyadic adaptive HOPS, relative bounds are used,
$$
\Delta_S(t)=\delta_S\|\psi^{(0)}(t)\|,\qquad
\Delta_A(t)=\delta_A\|\psi^{(0)}(t)\|.
$$
For sufficiently large aggregates, this adaptive strategy yields size-invariant, \(O(1)\), scaling because the average basis size saturates once the system size exceeds the delocalization length [2301.03718].

Dyadic HOPS was developed for spectroscopic response functions. For linear absorption, the dipole autocorrelation can be written in terms of pure states in the one-exciton manifold, and a normalized dyadic equation propagates ket and bra states in different electronic Hilbert spaces [2111.01089]. DadHOPS combines this dyadic construction with adHOPS and introduces an initial-state decomposition that reconstructs the linear absorption spectrum from a sum over locally excited initial conditions. The method is reported to allow trivial inclusion of static disorder in the Hamiltonian and to achieve size-invariant scaling for sufficiently large aggregates [2301.03718].

A different line of development replaces the explicit hierarchy of vectors by a tensor-network representation. The hierarchy of matrix product states (HOMPS) rewrites HOPS using formal creation and annihilation operators and then formulates the resulting stochastic first-order differential equation in terms of matrix product states and matrix product operators. In this way, the exponential complexity of HOPS can be reduced to scale polynomial with the number of particles [2109.06393].

These variants do not alter the foundational HOPS idea—ensemble reconstruction from stochastic pure states and auxiliary memory tiers—but they target the principal computational bottlenecks: combinatorial hierarchy growth, large Hilbert-space dimension, and slow Monte Carlo convergence.

## 6. Applications, performance regimes, and later generalizations

The original HOPS paper demonstrated the method on the Spin–Boson model, on absorption spectra of molecular aggregates, and on excitation-energy transfer in a photosynthetic pigment-protein complex [1402.4647]. For the Spin–Boson model with a single-exponential correlation function,
$$
H_{\rm sys}=-(\Delta/2)\sigma_x+(\epsilon/2)\sigma_z,\qquad L=\sigma_z,
$$
the hierarchy reduces to
$$
\partial_t|\psi_t^{(k)}\rangle=
\Bigl[-iH_{\rm sys}-k(\gamma+i\Omega)+\sigma_z z_t^*\Bigr]|\psi_t^{(k)}\rangle
+k\,g\,\sigma_z |\psi_t^{(k-1)}\rangle
-\sigma_z |\psi_t^{(k+1)}\rangle.
$$
In practice, the non-linear HOPS version was found to converge with \(M\sim 10^2\)–\(10^3\) trajectories where the linear version may need \(M\sim 10^4\), and truncation order \(K\simeq 4\)–\(8\) was sufficient even for strong coupling [1402.4647].

Subsequent work broadened the performance envelope. For sub-Ohmic environments with algebraically decaying bath correlation functions, HOPS was reported to show perfect agreement with other methods from weak to strong coupling and for zero and non-zero temperature [1710.08268]. For linear absorption of large molecular aggregates, DadHOPS was applied to the photosystem I core complex and to perylene bis-imide aggregates; the former showed Dyadic HOPS in quantitative agreement with DM-HEOM, while the latter displayed basis-size saturation and the onset of size-invariant cost [2301.03718].

HOPS has also been extended into domains where bath observables or multitime observables are essential. In a strongly coupled quantum heat engine, HOPS was used to compute bath energy and interaction energy during both transient and periodic steady-state operation [2402.06039]. In quantum-dot spectroscopy with a super-Ohmic phonon bath, HOPS was used for multitime correlation functions underlying absorption and resonance fluorescence spectra, including temperature- and coupling-dependent Mollow triplets [2412.20598]. In singlet fission, adaptive HOPS was generalized to include both Holstein and Peierls vibrations and applied to EP–PDI, where Peierls vibrations were found to accelerate singlet fission and to support singlet-mediated triplet transport on the 100-nm scale [2505.02292].

A recurring theme across these applications is that HOPS remains formally exact at the level of the infinite hierarchy, while practical success depends on how effectively one manages bath-correlation fitting, hierarchy truncation, stochastic sampling, and system-space compression. This suggests that the contemporary significance of HOPS lies not only in its original derivation but also in the family of adaptive, dyadic, tensor-network, and observable-specific extensions that preserve the stochastic pure-state architecture while enlarging the accessible class of non-Markovian problems [1402.4647; 2301.03718; 2109.06393].

Source: https://www.emergentmind.com/topics/hierarchy-of-pure-states-hops