To overcome this unfeasible scaling, many works have shown that excited polaritonic states at weak-to-moderate light-matter coupling strengths can be simulated in the so-called first-excited subspace of the total Hamiltonian (Eq. 17). In singly excited subspace, the overall scaling can be written as dim[ Ĥ ] ∝ N el N mol N F N mode , which is now linear scaling in all quantities. The zeroth S 0 and singly excited S 1 subspaces can be explicitly written as,