This numerical challenge can be easily resolved by decomposing the Hamiltonian as Ĥ = ε1 + ∆ Ĥ, where ε is a constant reference energy level, and ∆ Ĥ is the fluctuation around ε. Because ε1 commute with ∆ Ĥ, the quantum evolution becomes e -i ℏ Ĥ∆t = e -i ℏ ε∆t • e -i ℏ ∆ Ĥ•∆t , where the highly oscillatory part of the phase e -i ℏ ε∆t will be evaluated analytically, and the rest of the propagation e -i ℏ ∆ Ĥ•∆t will be numerically updated using the Chebyshev expansion. Further, in the HTC Hamiltonian, these subspaces are not directly connected by Hamiltonian terms. Following this logic, we rewrite the Hamiltonian in Eq. 37a as