Coupling molecules inside an optical cavity generates polaritons, which have been the target of interest in recent years. These hybrid polariton states arise from strong coupling between molecular states (electronic or vibrational excitations) and quantized Fock states of a cavity-confined radiation mode inside. The formation and subsequent dynamics of the polariton states are hypothesized to be the key to various important phenomena, such as cavity-induced alternation of chemical reactivity , as well as the enhanced transport of excitonpolaritons. In this work, we investigate the excitonic polariton relaxation process, where the excited states of N molecules (hereon referred to as sites) are coupled to the photonic excitation of the cavity. The system consists of two polaritonic states which are highly delocalized mixtures of the photon state and excited states of a large number of sites (for N ∼ 10 6 -10 12 ), as well as a large number of dark states (a total of N -1 of them) with purely excitonic character and carry zero transition dipole from the ground state (hence optically dark). The coupling, and corresponding population transfer as a result, in between the polaritonic as well as the dark state is facilitated by the nuclear/phonon (hereon referred to as the bath) degrees of freedom (DOF). This is crucial for an accurate description of the polariton relaxation process, as demonstrated by the collective scaling of the transfer rates between the polariton and dark states as well as the polaron decoupling effect observed in the collective limit, the interesting energy gap law, and phonon and polariton bottleneck effect. While the relaxation process of cavity polaria) Electronic mail: ylai19@ur.rochester.edu b) Electronic mail: pengfei.huo@rochester.edu ton has been the subject of various experimental studies, accurate simulation of the population dynamics is essential to the theoretical studies of an excitonic polaritonic system. Numerically exact methods such as multi-configurational time-dependent Hartree (MCTDH) and its multilayer extension, and hierarchical equations of motion (HEOM) adopt full quantum mechanical treatment of every DOF, often facing the challenge of a rapidly increasing computational cost when increasing the number of molecules N . Reported simulation of collective dynamics using truncated equations (CUT-E) takes advantage of permutational symmetry of the material states and can significantly reduce the computational cost through an effective 1/N expansions. Mixed-quantum-classical (MQC) methods treat bath DOF classically and are feasible for stimulating a large number of molecules N collectively coupled to the cavity. These methods have been applied previously to study the collective behavior of the polaritonic dynamics. Reduced population dynamics based on transfer rates represent an alternative way to investigate the polariton population dynamics. These approaches often use the well-known Fermi's golden rule (FGR) to the polariton dynamics to describe transition rate and simulate the population dynamics. The population transfer rates between polariton states and dark states can be obtained without any explicit propagation of the bath equation of motion (EOM), and in principle can be performed with an arbitrarily large number of sites N . Another strength of the rate-based approach is the decomposition of the EOM into explicit transfer rates, which provides valuable physical intuition to the underlying dynamics. However, the widely-known conventional frequency-domain formulation of FGR commonly referred to as the equilibrium Fermi's golden rule (E-FGR) assumes that the bath DOF initial condition is at equilibrium with the initial system state, which is not necessarily valid upon photoexcitation. Additionally, E-FGR leads to a fully Markovian EOM, which means a different time-dependent rate equation is required to capture any non-Markovian characteristics of the dynamics.