A few important properties of this matrix are as follows: (I) The N F block diagonals M n are composed of the diagonal energies and diagonal DSE elements ϵ α,n + D αα in addition to the off-diagonal DSE elements D αβ . (II) The super-and sub-diagonal blocks are composed only of the matter dipole operator matrix elements μαβ weighted by the photon number of the larger diagonal Fock state label n. Note here that for clarity, we have neglected the zero-point energy of the photonic mode and R dependence. One can easily see from the right-most side of Eq. 19 that this matrix is extremely sparse, especially for larger number of included Fock basis states (needed for convergence), so many approximate diagonalization schemes can be used which rely on the properties of sparse matrices, such as the Lanczos and Krylov subspace techniques, and are able to return the lowest eigenvalues and eigenvectors without loss of physics but with a large computational speed-up.