where T is the maximum time of the recorded time-domain spectra and cos(πt/2T ) is a smoothing function. As can be seen, to compute R (1) (t 1 ) and R (ω), one needs to act the dipole operator, μ, on the density matrix, which in principle also scales as O(K 4 ). However, μ also takes the form of a very sparse matrix in the single and double excitation subspace, as illustrated in Fig. . The sparsity of the dipole operator scales as 1 -1/N 2 with N molecules coupled to a cavity mode. The action of μ on a ket vector (following the same arrangement of the basis as in Eq. 25) can be written as The above becomes a very simplified modification of |Ψ⟩,