The action of Ô on |Ψ⟩ is the same as ĤQ on |Ψ⟩ as described in Eq. 28, with the coefficients now scaled according to Eq. 36. Fig. presents the magnitude of the Chebyshev coefficients as a function of the expansion term index for various z-values. With an increase in z, more terms are needed for the expansion to converge. For z ≪ 1, only the first three or four terms in the expansion are required. For comparison, a typical propagation scheme such as RK-4 or a symplectic integrator (velocity verlet) typically needs a much smaller substep, δt ≈ ∆t/10-∆t/100 for a stable integration of the dynamics. The Chebyshev expansion approach thus provides a significant speedup of the dynamics propagation compared to solving Eq. 31 using RK4 because one can use a much larger ∆t (up to 100 times), with only a slight increase of the prefactor of the computational cost associated with how many terms are needed to form a converged Chebyshev expansion (each associated with the Ô|Ψ⟩ operations). In all simulations, we are in the regime of z ≪ 1 and we found that only n = 4 terms of the Chebyshev polynomial are required to converge the results.