Quantum-classical hierarchical equations of motion

We develop a quantum-classical hierarchical equations of motion (QC-HEOM) approach for simulating non-Markovian open quantum systems. The method combines the ensemble-averaged classical path reference of the quantum-classical path integral formalism with a hierarchy of auxiliary quantum influence functionals. By incorporating thermal fluctuations through an ensemble average over reference trajectories, the hierarchy is required to represent only the residual quantum memory associated with the imaginary part of the bath response function. Consequently, unlike conventional hierarchical equations of motion, QC-HEOM does not require Matsubara or Padé expansions of the thermal kernel and exhibits only weak temperature dependence of the hierarchy size. Furthermore, because thermal fluctuations are supplied through reference classical trajectories, the framework provides a natural route toward extensions beyond harmonic baths and enables the incorporation of anharmonic and molecular environments through externally generated trajectories. We derive the QC-HEOM formalism and demonstrate its exactness for a harmonic bath. Applications to an asymmetric spin-boson model and an eight-site model of the Fenna-Matthews-Olson complex illustrate the accuracy of QC-HEOM. It reproduces benchmark exact quantum dynamics results while requiring substantially fewer auxiliary objects, particularly at low temperatures. These results establish QC-HEOM as an efficient framework for treating residual quantum memory in quantum-classical descriptions of open-system dynamics. The separation of thermal fluctuations from residual quantum memory through the use of Wigner trajectories suggests an approximate route toward hierarchical treatments of complex anharmonic environments beyond harmonic bath simulations, which shall be developed in the future.