As we approach the large N a limit with N a ∼ 10 6 , which is considered the experimentally relevant regime, we note that one should take advantage of the wellknown mean-field solution 101,102 or the recently proposed approach 103,104 (for the zero-temperature case and with a single high-frequency vibration mode per molecule) that effectively only treats one or a few molecules. We can also exploit the sparsity of the HTC Hamiltonian to enable direct simulations of N ∼ 10 6 molecules coupled to a single cavity mode in the single excitation subspace with trajectory-based methods. One could also explore using the relatively short memory kernel time compared to the density matrix dynamics time to facilitate the simulations. These directions remain to be explored in the future. the benchmark simulations, and polariton relaxation dynamics for models 5 and 6.