In Figure , we present the population dynamics of the polariton states and the dark states for model 4, which varies the reorganization energy λ from 10 meV to 50 meV while keeping the other parameters fixed. All trajectory-based methods are able to qualitatively account for the relaxation of the upper polariton state into the other states. However, both MQC methods (Ehrenfest and GFSH) predict a slower relaxation rate and a larger steady-state population for the upper polariton state compared to the HEOM results. For λ = 10 meV and λ = 30 meV (Fig. , Fig. , Fig. and Fig. ), both MQC methods predict little changes in the dark state population after the dark state populations have reached a maximum value, which is in contrast to numerically exact results from the HEOM method. On the other hand, for λ = 50 meV (Fig. and Fig. ), the MQC methods underestimate the populations of the dark states as compared to the HEOM method. As a result, the increase in lower polariton population calculated from both MQC methods, which comes through the transitions from the dark states, is always smaller than that predicted by the HEOM method. The mapping-based methods (γ-SQC and spin-LSC) show more reliable relaxation rates for the upper polariton state compared to the HEOM method and the steady-state population for the upper polariton is also similar to the HEOM method for λ = 10 meV and λ = 30 meV (Fig. , Fig. , Fig. , and Fig. ). We also observe that the transitions into the dark states are better captured with the mapping-based methods, and the γ-SQC method outperforms the spin-LSC method in predicting the longer time transitions from the dark states to the lower polariton states. Consequently, the population of the lower polariton state computed from the mapping-based method is also comparable to the populations predicted from the HEOM method.