The hierarchical equations of motion formalism (HEOM) is a numerically exact, nonperturbative method to study the dynamics of open quantum systems interacting with harmonic baths. Pioneered by works of Kubo and Tanimura,1 further developed by Ishizaki, Tanimura, and Fleming,2 HEOM has found applications in a variety of contexts, such as in spectroscopy and excitation energy and charge transfer studies.2–8 In particular, the method has found numerous uses in studying exciton transfer and coherence dynamics in biological systems, such as in the Fenna-Mathews-Olson (FMO) complexes.9–11 The HEOM calculations allowed accurate calculations at low temperatures, whereas perturbative approaches such as Redfield12,13 or Forster12,14fail.15–22 The HEOM has often been used to obtain numerically exact results for spin-Boson problems across several regimes, and it can be used as a reference methodology for other approximate methods8,23–25 applied to modeling open quantum systems.
Despite the great utility and potential of the method, it has been a rather niche approach, likely due to its computational complexity. As a result, only a handful of open-source software implementing HEOM dynamics and spectra calculations has been made available. Notably, these are the Fortran codes from the Tanimura group,26including a GPU implementation,27 the PHI code from the Schulten group,5,6 which emphasizes an efficient CPU parallelization, as well the nanoHub implementation of Kreisbeck and Kramer.28 Although these codes have been utilized in the past, they are designed as standalone programs with fixed user interaction protocol. A modular all-Python implementation of HEOM is also available within the QuTiP2 program for quantum optics.29 Although it is flexible in its design, it is primarily meant to be used within Python programs and may not be easily used within C++ codes.
In this work, we present our own implementation of the HEOM method within the Libra library.30 Libra is a comprehensive, open-source, quantum dynamics package. Adding HEOM to the Libra repertoire not only provides researchers with easy access to the method, but also allows them to integrate the method easily into already established workflows. The HEOM implementation reported in this work emphasizes modularity, transferability, and user friendliness via intuitive and flexible Python interface. Our code includes multiple distinctive features, such as the ability to compute absorption spectra, general system bath coupling, and more. The essential components of the HEOM code are written in C++ in a very modular fashion, so that they can be used within other C++ codes. The openmp parallelization and scaled HEOM are enabled for better efficiency. The low-level components are exposed to Python for their use in Python programs, including in Jupyter notebooks. They can be used in various independent contexts or organized into pre-defined workflows. Our higher-level Python modules provide implementation for the corresponding workflows for computing population dynamics and absorption spectra as well as provide auxiliary functions for storing and analyzing the data produced.
This paper aims to establish the required vocabulary and provide a brief overview of the underlying concepts and methodologies for users and developers to be able make a connection between the implementation and the underlying theory. We explain the package organization and some key algorithmic details. We illustrate the use cases and explain how users can setup the calculations with the HEOM module of the Libra package. Some qualitative features of the dynamics in various types of environments are demonstrated.