  The hierarchical equations of motion (HEOM) provide a numerically exact
approach for computing the reduced dynamics of a quantum system linearly
coupled to a bath. We have found that HEOM contains temperature-dependent
instabilities that grow exponentially in time. In the case of continuous-bath
models, these instabilities may be delayed to later times by increasing the
hierarchy dimension; however, for systems coupled to discrete, non-dispersive
modes, increasing the hierarchy dimension does little to alleviate the problem.
We show that these instabilities can also be removed completely at a
potentially much lower cost via projection onto the space of stable eigenmodes;
furthermore, we find that for discrete-bath models at zero temperature, the
remaining projected dynamics computed with few hierarchy levels are essentially
identical to the exact dynamics that otherwise might require an intractably
large number of hierarchy levels for convergence. Recognizing that computation
of the eigenmodes might be prohibitive, e.g. for large or strongly-coupled
models, we present a Prony filtration algorithm that may be useful as an
alternative for accomplishing this projection when diagonalization is too
costly. We present results demonstrating the efficacy of HEOM projected via
diagonalization and Prony filtration. We also discuss issues associated with
the nonnormality of HEOM.
