whereas the δ (ξ) E (as expressed in Eq. 39a and Eq. 40a) remains unchanged. Because different subspaces are not directly coupled by ĤQ , when propagating the quantum dynamics governed by H Q , the above Chebyshev scheme will be used independently for each subspace. On the other hand, different subspaces are connected through the cavity loss process (such as |G 1 ⟩ → |G 0 ⟩ in the single excitation subspace, and |G 2 ⟩ → |G 1 ⟩ and |E 1 n ⟩ → |E 0 nm ⟩ in the double excitation space), and when computing the response function, the instantaneous action of the dipole operator will connect them. For those two types of operations, one will apply the stochastic Lindblad algorithm (Appendix A) and the rule of dipole on state vector directly, without shifting the reference energies. This simple technique outlined in Eq. 42 significantly reduces the noise of the numerical propagation and increases the accuracy of the Chebyshev expansion method, as shown in Appendix B.