where ρS is the system reduced density operator, β = 1/k B T is the inverse temperature, and Z B = Tr[e -β ĥB ] is the bare-bath partition function. In Eq. 39, we assume that the bath is in thermal equilibrium so that the bath reduced density matrix operator takes the form of a Boltzmann distribution. Also, the system is initially placed in the upper polariton state |Ψ(0)⟩ = |+⟩ and the associated system reduced density operator is given by ρS = |+⟩ ⟨+|. This initial state is chosen so that we can observe the relaxation from the upper polariton state to the dark states and lower polariton state because of exciton-phonon coupling. Model Parameters. We construct six HTC benchmark models to investigate the performance of the trajectory-based non-adiabatic methods to simulate polariton relaxation dynamics. These models scan a range of the following physical parameters: (1) number of molecules (N a ), (2) light-matter detuning (ω c -ω a -λ), (3) single-molecule coupling strength (g c ), (4) reorganization energy (λ), (5) bath cutoff frequency (ω ν ), and ( ) temperature (β = 1/k B T ).