Note that the polariton state defined in Eq. 14 and the dark states in Eq. 17 are diabatic in nature, because there is no nuclear configuration dependence. In this case, the exciton-phonon couplings ĤSB will couple these polariton states and make transitions among them. On the other hand, one can also choose to define ĤS + ĤSB as the polariton Hamiltonian, whose eigenvector will be adiabatic polariton and dark states, and both of their characters will parametrically depend on the nuclear configuration. In this case, the phonon fluctuation caused by ĤSB on polariton and dark states will be counted as the adiabatic polariton energy fluctuations, and the transitions among these adiabatic polariton and dark states are caused by the nuclear kinetic energy operators (as the derivative couplings). This is the picture used in the theoretical simulations in Ref. 26 (see Fig. in that reference). Nevertheless, rigorous quantum mechanical descriptions of the dynamics (for all DOFs) will generate identical results in both pictures. Ehrenfest dynamics, as well as the mean-field-like approach (such as the mapping-based methods used in this work), are representation-independent, and thus generate identical results for both representations. The trajectory surface hopping approach, on the other hand, is not representation independent, and often performs the best in the adiabatic representation, so the adiabatic polariton Hamiltonian Ĥpl = ĤS + ĤSB should be used in these surface hopping simulations.