For the NE-FGR approach, the population dynamics are computed by solving the EOM for a 3-element array P = [P + , P -, P D ], where P + (P -) is the upper (lower) polariton population, and P D is the sum of the N -1 identical dark mode populations. The time local (TL) EOM in Eq. 16a (which only depends on P j (t)) is solved numerically, using the memory Kernal in Eq. 12 computed from the correlation function C i→f (t, s), with expressions in Eq. 23 and Eq. 26a. The EOM is numerically solved using a 4th-order Runge-Kutta algorithm on a quadrature, following the same procedure in Ref. 56. The step size dt for integration is taken to be dt = 0.1 fs. Note that solving TL-QME in Eq. 16a with memory kernel (Eq. 12) is equivalent to solving the rate equation in Eq. 16b with the NE-FGR rate (Eq. 17). For N = 10 6 under the collective coupling limit, numerically, C 0 ±→D = C 0 D→± → 1 and F ±→∓ (t, s) ∝ O(N -2 ) → 0, it is safe to ignore them in the calculation. A similar numerical procedure is used to solve the TNL-QME in Eq. 11 (note that it depends on P j (t -s)) with the same memory kernel in Eq. 12.