The time evolution in open quantum systems, such as a molecular aggregate in contact with a thermal bath, still poses a complex and challenging problem. The influence of the thermal noise can be treated using a plethora of schemes, several of which decompose the corresponding correlation functions in terms of weighted sums of exponential functions. One such scheme is based on the hierarchical equations of motion (HEOM), which is built using only certain forms of bath correlation functions. In the case where the environment is described by a complex spectral density or is at a very low temperature, approaches utilizing the exponential decomposition become very inefficient. Here we utilize an alternative decomposition scheme for the bath correlation function based on Chebyshev polynomials and Bessel functions to derive a hierarchical equations of motion approach up to an arbitrary order in the environmental coupling. These hierarchical equations are similar in structure to the popular exponential HEOM scheme, but are formulated using the derivatives of the Bessel functions. The proposed scheme is tested up to the fourth order in perturbation theory for a two-level system and compared to benchmark calculations for the case of zero-temperature quantum Ohmic and super-Ohmic noise. Furthermore, the benefits and shortcomings of the present Chebyshev-based hierarchical equations are discussed.