Hierarchy of Pure States (HOPS)

• HOPS is a wavefunction-based method that reformulates non-Markovian stochastic Schrödinger equations into a hierarchy of pure-state trajectories to recover the exact reduced density operator.
• It systematically manages environmental memory and truncation challenges, allowing accurate simulations of systems like the spin-boson model and photosynthetic complexes.
• Adaptive, dyadic, and tensor-network variants enhance HOPS by improving convergence, scalability, and efficiency in modeling complex open quantum dynamics.

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 (Süß et al., 2014). 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 (Süß et al., 2014).

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 (Süß et al., 2014)

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

The central obstruction to direct simulation is the functional derivative $$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).$$1. 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 (Süß et al., 2014).

2. Construction of the hierarchy and the representation of memory

For a single-exponential bath kernel,

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

one defines

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

The hierarchy then becomes

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

with initial conditions $$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).$$5 and $$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).$$6. The physical state is $$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).$$7 (Süß et al., 2014).

For a sum of exponentials,

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

HOPS introduces a multi-index $$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).$$9 and auxiliary states

$$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$0

leading to

$$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$1

Here $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$2, $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$3, and the noises satisfy $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$4 and $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$5 (Süß et al., 2014).

The multi-index carries the environmental memory. In the original formulation, each $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$6 counts how many times the corresponding “history operator” $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$7 has acted; retaining larger weights in $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$8 means that the stochastic state probes deeper into the bath memory kernel (Süß et al., 2014). 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 $$J(\omega)=\sum_\lambda |g_\lambda|^2\delta(\omega-\omega_\lambda)$$9, and it is designed to improve Monte Carlo convergence (Süß et al., 2014).

3. Exact reconstruction, normalization, and observables

In the linear formulation, the reduced state is reconstructed as

$T$0

In the non-linear formulation, one evolves normalized states $T$1 and still has

$T$2

Thus, in either case, the exact non-Markovian reduced dynamics is obtained from the ensemble average over sufficiently many noise realizations (Süß et al., 2014).

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 [(Süß et al., 2014); (Zhang et al., 2018)].

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 (Hartmann et al., 2017).

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 $T$3, which enables a fully quantum dynamical treatment of energetic contributions in a strongly coupled quantum heat engine (Boettcher et al., 2024). 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 (Toivonen et al., 2024).

4. Truncation schemes, convergence control, and computational scaling

Practical HOPS calculations truncate the hierarchy at total order $T$4. For a single exponential, a common terminator is

$T$5

and analogous expressions are used for multi-exponential triangular truncation (Süß et al., 2014). Convergence is tested by increasing $T$6 until observables such as populations, coherences, or spectra change negligibly. A posteriori, comparison between truncation orders $T$7 and $T$8 provides an upper bound on the truncation error because the method is systematically improvable (Süß et al., 2014).

The combinatorics of the auxiliary set is central to the numerical cost. For $T$9, the number of auxiliary states is

$$\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].$$0

and the overall dense cost per time step scales as $$\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].$$1, where $$\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].$$2 is the number of independent noise realizations and $$\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].$$3 is the system dimension. The non-linear hierarchy typically requires far fewer trajectories $$\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].$$4 because of importance sampling (Süß et al., 2014).

The limitations of standard triangular truncation motivated alternative truncation strategies. The $$\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].$$5-particle approximation ($$\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].$$6PA) retains only those multi-indices for which the number of distinct sites $$\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].$$7 with $$\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].$$8 is at most $$\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].$$9, while the $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$0-mode approximation ($|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$1MA) retains only those multi-indices for which the total number of modes $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$2 with $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$3 is at most $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$4 (Zhang et al., 2018). 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 (Zhang et al., 2018).

Benchmark calculations in absorption and excitation transfer showed that $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$5MA and $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$6PA can be combined with moderate depth $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$7 to generate a ladder of hierarchy sizes that grows slowly enough for practical convergence checks (Zhang et al., 2018). 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 $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$8 and auxiliary-index wavefunctions $|\Psi(0)\rangle=|\psi_0\rangle\otimes|0\rangle_B$9, chosen such that the local truncation error in the time derivative satisfies

$|\psi_t(z^*)\rangle$0

In dyadic adaptive HOPS, relative bounds are used,

$|\psi_t(z^*)\rangle$1

For sufficiently large aggregates, this adaptive strategy yields size-invariant, $|\psi_t(z^*)\rangle$2, scaling because the average basis size saturates once the system size exceeds the delocalization length (Gera et al., 2023).

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 (Chen et al., 2021). 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 (Gera et al., 2023).

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 (Gao et al., 2021).

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 (Süß et al., 2014). For the Spin–Boson model with a single-exponential correlation function,

$|\psi_t(z^*)\rangle$3

the hierarchy reduces to

$|\psi_t(z^*)\rangle$4

In practice, the non-linear HOPS version was found to converge with $|\psi_t(z^*)\rangle$5–$|\psi_t(z^*)\rangle$6 trajectories where the linear version may need $|\psi_t(z^*)\rangle$7, and truncation order $|\psi_t(z^*)\rangle$8–$|\psi_t(z^*)\rangle$9 was sufficient even for strong coupling (Süß et al., 2014).

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 (Hartmann et al., 2017). 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 (Gera et al., 2023).

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 (Boettcher et al., 2024). 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 (Toivonen et al., 2024). 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 (Lynd et al., 4 May 2025).

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 [(Süß et al., 2014); (Gera et al., 2023); (Gao et al., 2021)].