In Fig. , we present the population dynamics of the polariton states and the dark states for model 1, which varies the number of molecules N a from 5 molecules to 15 molecules while keeping the other parameters fixed. The trajectory-based methods are depicted using open circles and are compared to numerically exact results (HEOM) depicted using solid lines. We can see that in Figs. 2a-2l, all trajectory-based methods are able to semi-quantitatively account for the relaxation of the upper polariton (UP) state into the dark states (DS), and eventually to the lower polariton (LP) state. However, both mixed quantum-classical (MQC) methods, Ehrenfest, and global flux surface hopping (GFSH) predict a slower relaxation rate and a larger steady-state population (see Fig. -Fig. ) for the upper polariton state compared to the HEOM results. Focusing on the transitions into and out of the dark states, we observe that both MQC methods are only able to qualitatively capture the increase in the total dark state populations. After the dark state populations have reached a maximum value, both MQC methods predict little changes in the dark state population, which is in contrast to numerically exact results from the HEOM simulation. As a result, the increase in lower polariton population calculated from both MQC methods, which comes through the transitions from the dark states, is smaller than that predicted by the HEOM method.