To simulate the 2DES spectra for the polaritonic system, we follow the previous procedure outlined in Ref. 19. Furthermore, we explicitly compute the contribution from individual Liouville pathways (Feynman diagrams) through Eq. 4a to Eq. 4d. The quantum DOFs, which are described as the exciton-Fock states {|a⟩, |b⟩} ≡ {|G⟩ ⊗ |k⟩, |m⟩ ⊗ |k ′ ⟩} where |k⟩ and |k ′ ⟩ are the Fock states (eigenstates of ĤPh in Eq. 36). The nuclear DOFs are sampled from the initial Wigner distribution based on Eq. 10. The total initial condition for the molecule-cavity hybrid system (the system reduced density operator, see Eq. 15) is