The main challenge of solving Eq. 18 for the HTC Hamiltonian under the collective coupling regime is the large number of N that one has to incorporate. The right-hand side of the differential equation, ĤQ |Ψ⟩, requires the operation of ĤQ matrix on a state vector c j (t). For the HTC Hamiltonian with N molecules and one cavity mode, there are K = N + 1 ∝ O(N ) states in the single excitation subspace and K = N (N -1)/2 + N + 1 ∝ O(N 2 ) states in the double excitation subspace. The force evaluation for a mean-field-like approach of the HTC model (Eq. 19), on the other hand, only scales with O(K) due to the purely diagonal structure of Ĥsb , in the diabatic site basis. The typical experimental condition requires N = 10 6 -10 12 , resulting in a very large K, and the computational work is often limited by N = 10 1 -10 2 due to the computational cost.