  The hierarchical equations of motion technique has found widespread success
as a tool to generate the numerically exact dynamics of non-Markovian open
quantum systems. However, its application to low temperature environments
remains a serious challenge due to the need for a deep hierarchy that arises
from the Matsubara expansion of the bath correlation function. Here we present
a hybrid stochastic hierarchical equation of motion (sHEOM) approach that
alleviates this bottleneck and leads to a numerical cost that is nearly
independent of temperature. Additionally, the sHEOM method generally converges
with fewer hierarchy tiers allowing for the treatment of larger systems.
Benchmark calculations are presented on the dynamics of two level systems at
both high and low temperatures to demonstrate the efficacy of the approach.
Then the hybrid method is used to generate the exact dynamics of systems that
are nearly impossible to treat by the standard hierarchy. First, exact energy
transfer rates are calculated across a broad range of temperatures revealing
the deviations from the Forster rates. This is followed by computations of the
entanglement dynamics in a system of two qubits at low temperature spanning the
weak to strong system-bath coupling regimes.
