and when there is no dipole orientation disorders, simply set cos θ n = 1 in the above expression. The operation in Eq. 28 has a linear scaling in the number of states O(K), and thus ∝ O(N ) for a single excitation space and ∝ O(N 2 ) double excitation space when having N -molecules for the polaritonic systems. The action of ĤQ on |Ψ⟩ as outlined in Eq. 28 can be used to accurately calculate electronic evolution via a simple RK4 solution of Eq. 17 or through polynomial expansions of exponential propagators such as the Chebyshev series, which will be discussed next.