The hierarchical equations of motion (HEOM) provide a numerically exact approach for simulating the dynamics of open quantum systems coupled to a harmonic bath. However, its applicability has traditionally been limited to specific spectral forms and relatively high temperatures. Recently, an extended version called Free-Pole HEOM (FP-HEOM) has been developed to overcome these limitations. In this study, we demonstrate that the FP-HEOM method can be systematically employed to investigate higher order master equations by truncating the FP-HEOM hierarchy at a desired tier. We focus on the challenging scenario of the spin-boson problem with a sub-Ohmic spectral distribution at zero temperature and analyze the performance of the corresponding master equations. Furthermore, we compare the memory kernel for population dynamics obtained from the exact FP-HEOM dynamics with that of the approximate Non-Interacting-Blip Approximation (NIBA).