Abstract: A hierarchical equations of motion (HEOM) approach is developed for general open quantum systems coupled to fermionic environment. The HEOM method is in principle formally exact, as it resolves nonperturbatively the combined effects of many-body interaction, system-bath dissipation, and non-Markovian memory 1-4. In practice, the HEOM approach is highly accurate and efficient for the characterization of strongly correlated quantum impurity systems. The HEOM approach has been successfully applied to characterize a wide range of equilibrium and nonequilibrium properties of quantum dot systems. These include the studies on dynamical Coulomb blockade, transient current response, and thermoelectric properties of Kondo-correlated quantum dots 5,6. The numerical results of the HEOM approach are considered to be quantitatively accurate, as long as they converge with respect to the truncation of the hierarchy. It has been demonstrated that the HEOM approach is capable of addressing Kondo resonance and Fermi-liquid characteristics, achieving same level of accuracy as the latest high-level numerical renormalization group method. In this talk, the accuracy and practicality of the HEOM approach will be exemplified by a number of examples. These include the use of HEOM as an impurity solver for dynamical mean-field theory, the characterization of thermopower of multi-level quantum dots, and the understanding of Kondo resonance in the d-CoPc/Au(111) adsorption system.