  Understanding the dynamics of open quantum many-body systems is a major
problem in quantum matter. Specifically, efficiently solving the spectrum of
the Liouvillian superoperator governing such dynamics remains a critical open
challenge. Here, we put forward a method for solving the many-body Liouvillian
spectrum and dynamics based on the non-Hermitian kernel polynomial method and
tensor-network techniques. We demonstrate the faithfulness of our method by
computing the dynamics of the dephasing quantum compass model with a gradient
magnetic field and comparing it with exact results. In particular, we show that
our method allows us to characterize the quantum Zeno crossover and the
reduction of relaxation rate due to Stark localization in this model. We
further demonstrate the ability of our method to go beyond exact results by
exploring nearest-neighbor interaction effects on the Liouvillian dynamics,
elucidating the interplay between Stark localization and many-body
interactions. Our method provides an efficient solution to many-body
Liouvillian spectrum and dynamics, establishing a methodology to explore large
open quantum many-body systems.
