  In this work, we propose the Prony fitting decomposition (PFD) as an accurate
and efficient exponential series method, applicable to arbitrary interacting
bath correlation functions. The resulting hierarchical equations of motion
(HEOM) formalism is greatly optimized, especially in extremely low temperature
regimes that would be inaccessible with other methods. For demonstration, we
calibrate the present PFD against the celebrated Pad\'e spectrum decomposition
method, followed by converged HEOM evaluations on the single-impurity Anderson
model system.
