We now consider the case where the upper polariton is initially populated, followed by population relaxation. Fig. presents the population dynamics for N = 5 case with zero detuning ∆ = 0, where E-FGR (short dashed), NE-FGR (long dashed), and TNL (thin solid) are used to generate the population dynamics, and are compared against to the numerically exact solution from HEOM (dotted line). Fig. and Fig. present dynamics for N = 5, with the positive and negative detuning cases, respectively. Fig. -f presents the corresponding dynamics, with N = 10 6 , without the HEOM results. Among all of the approximate population dynamics, we consider the TNL (Eq. 11) as the most accurate one, because it only makes the second-order perturbation approximation of the system bath coupling. The NE-FGR theory (Eq. 16a) carries an additional approximation by using time local populations P k (t -s) → P k (t) compared to TNL QME. The E-FGR population dynamics (Eq. 33) has the most approximations, which only uses the E-FGR rate in the population dynamics, compared to the NE-FGR population dynamics (Eq. 16b). For all of the dynamics presented here, the three perturbative methods resulted in virtually identical dynamics. This means that the transient (∼ 50 fs) difference between the three meth- ods is not significant enough to influence the dynamics. The population dynamics therefore can be analyzed in the context of E-FGR, where in the case of N = 5 the lower polariton and dark states reach an equilibrium corresponding to the relative strength of k D→-and k -→D . In the case of N = 10 6 , the equilibrium populations consist of mostly dark states, consistent with predictions from the density of states argument and entropy-based analysis. Nevertheless, we must emphasize that for other model systems with different spectral densities J(ω), it is possible for NE-FGR and E-FGR to differ significantly, in which case NE-FGR is much preferable to use.