  The hierarchical equations of motion (HEOM) approach is an accurate method to
simulate open system quantum dynamics, which allows for systematic convergence
to numerically exact results. To represent the effects of the bath, the
reservoir correlation functions are usually decomposed into the summation of
multiple exponential terms in the HEOM method. Since the reservoir correlation
functions become highly non-Markovian at low temperatures or when the bath has
complex band structures, a present challenge is to obtain accurate exponential
decompositions that allow efficient simulation with the HEOM. In this work, we
employ the barycentric representation to approximate the Fermi function and
hybridization functions in the frequency domain. The new method, by
approximating these functions with optimized rational decomposition, greatly
reduces the number of basis functions in decomposing the reservoir correlation
functions, which further allows the HEOM method to be applied to ultra-low
temperature and general bath structures. We demonstrate the efficiency,
accuracy, and long-time stability of the new decomposition scheme by applying
it to the Anderson impurity model (AIM) in the low-temperature regime with the
Lorentzian and tight-binding hybridization functions.
