Of course, when only single-and double-excitation subspaces are considered, these results are (almost) analytically available. When only considering the single excitation subspace |j⟩ = {|g⟩⊗|1⟩, |e n ⟩⊗|0⟩}, the expression of S kj (∆t) ≡ α ⟨k|α⟩e -i ℏ ϵα∆t ⟨α|j⟩ is analytically available and can be precomputed and stored. One just needs to multiply the Skj ≡ e -i ℏ λ k ∆t/2 S kj (∆t)e -i ℏ λj ∆t/2 matrix onto the c j (t) column matrix to get the result of c k (t + ∆t), which in principle still requires a cost of O(N 2 ). As such, for the quantum dynamics problem of the HTC system, the main computational bottleneck is the ĤQ |Ψ⟩ evaluation. One would also want to avoid directly diagonalizing ĤQ , as well as increase the time step ∆t for the dynamics propagation as much as possible (because each step requires at least one ĤQ |Ψ⟩ evaluation) while keeping the propagation stable. In the following sections, we introduce several techniques that can significantly reduce the computational costs of propagating Eq. 17. These techniques can be broadly applied to any MQC trajectory-based methods for polariton quantum dynamics simulations.