Note that the delta-function in Eq. 14 is usually broadened with by a finite-width Gaussian or Lorentzian function to account for excitonic, photonic, and/or environmental relaxation/broadening processes present in realistic experimental conditions. Here, we chose to use a Lorentzian function where σ = 0.1 eV is the broadening parameter. In Eq. 15, µ M 0j is the molecular part of the polaritonic dipole calculated from Eq. 6 with µ M 0j = Â B = μ ⊗ 1ph . The matrix elements of μ are computed according to Eq. 3.