The initial nuclear distribution of MQC simulations are generated by sampling the Wigner density [ R|χ ] w = 1 π e -M (P 2 +ω 2 (R-R0) 2 )/ω , which is the Wigner transformation of the nuclear wavefunction χ(R) = R|χ in the initial state (see Eq. 53). Here, R and P are the nuclear coordinate and momentum, respectively. The initial state for the electronic-photonic subsystem is set to be |e0 . The population dynamics were averaged over 2000 trajectories, although 500 trajectories was good enough to produce the basic trend of the polariton dynamics. The light-matter coupling strength g c was chosen to be 0.001, 0.005 and 0.01.