  The proliferation of quantum fluctuations and long-range entanglement
presents an outstanding challenge for the numerical simulation of interacting
spin systems with exotic ground states. Here, we present a toolset of Chebyshev
polynomial-based iterative methods that provides a unified framework to study
the thermodynamical properties, critical behavior and dynamics of frustrated
quantum spin models with controlled accuracy. Similar to previous applications
of the Chebyshev spectral methods to condensed matter systems, the algorithmic
complexity scales linearly with the Hilbert space dimension and the Chebyshev
truncation order. Using this approach, we study two paradigmatic quantum spin
models on the honeycomb lattice: the Kitaev-Heisenberg (K-H) and the
Kitaev-Ising (K-I) models. We start by applying the Chebyshev toolset to
compute nearest-neighbor spin correlations, specific heat and entropy of the
K-H model on a 24-spin cluster. Our results are benchmarked against exact
diagonalization and a popular iterative method based on thermal pure quantum
states. The transitions between a variety of magnetic phases, namely
ferromagnetic, N\'eel, zigzag and stripy antiferromagnetic and quantum spin
liquid phases are obtained accurately and efficiently. We also accurately
obtain the temperature dependence of the spin correlations, over more than
three decades in temperature, by means of a finite temperature Chebyshev
polynomial method introduced here. Finally, we report novel dynamical
signatures of the quantum phase transitions in the K-I model. Our findings
suggest that the efficiency, versatility and low-temperature stability of the
Chebyshev framework developed here could pave the way for previously
unattainable studies of quantum spin models in two dimensions.
