We combine matrix product operator techniques with Chebyshev polynomial expansions and present a method that is able to explore spectral properties of quantum many-body Hamiltonians. In particular, we show how this method can be used to probe thermalization of large spin chains without explicitly simulating their time evolution, as well as to compute full and local densities of states. The performance is illustrated with the examples of the Ising and PXP spin chains. For the non-integrable Ising chain, our findings corroborate the presence of thermalization for several initial states, well beyond what direct time-dependent simulations have been able to achieve so far.The study of one-dimensional quantum many-body systems has motivated the emergence of a number of techniques, based on tensor network states (TNS). More concretely, they use matrix product states (MPS) and matrix product density operators (MPDO) [1-5] to approximate the ground states, low-lying excitations, thermal states, as well as time evolution. These methods have enabled the in-depth study of a multitude of models and the analysis of relevant physical phenomena.The success of such techniques is rooted in the ability of MPS and related ansatzes to accurately describe states that fulfill an area law of entanglement [6,7], satisfied (or only slightly violated) by many of the problems mentioned above [8]. There are, however, important open questions that such techniques cannot easily solve. In particular, excited states at finite energy density are difficult to approximate, except in very particular cases [9, 10], as they generically display volume law entanglement and, additionally, are embedded in highly dense spectral regions, which severely hinders the convergence of the algorithms. Out-of-equilibrium dynamics is also problematic: under time evolution a volume law often emerges, that makes an MPS approximation inadequate, except for short times. As a consequence, it is virtually impossible for standard MPS techniques to address the fundamental questions of equilibration and thermalization of relatively large closed quantum systems.A few alternative tensor network algorithms have tried to overcome these problems by avoiding the explicit representation of the states [11][12][13][14]. Although they extend the applicability of the toolbox and allow access to additional dynamical quantities in some scenarios, the fundamental goal of accessing the long time behavior in a general case, and thus deciding the appearance of equilibration or thermalization, has not been achieved.Here we introduce a new powerful tool to fill in these gaps. Our method is based on the use of MPO to approximate a family of generalized densities of states, and provides a means to directly address thermalization. More concretely, we combine TNS and the kernel polynomial method (KPM) [15] in a general scheme that provides access not only to the full density of states (DOS) of a given many-body Hamiltonian [16], but also to energy functions that are intimately related to the out-o...