A new method to improve the numerical stability of the hierarchical equations of motion for discrete harmonic oscillator modes.

The hierarchical equations of motion (HEOMs) have developed into an important tool in simulating quantum dynamics in condensed phases. Yet, it has recently been found that the HEOM may become numerically unstable in simulations using discrete harmonic oscillator modes [I. S. Dunn, et al., J. Chem. Phys. 150, 184109 (2019)]. In this paper, a new set of equations of motion are obtained based on the equivalence between the HEOM for discrete harmonic oscillator modes and the mixed quantum-classical Liouville equation. The new set of equations can thus be regarded as the expansion of the same phase space partial differential equation using different basis sets. It is shown that they have similar structures as the original HEOM but are free from the problem of numerical instability. The new set of equations are also incorporated into the matrix product state method, where it is found that the trace of the reduced density operator is not well conserved during the propagation. A modified time-dependent variational principle is then proposed to achieve better trace conservation.