A popular approach to simulate the non-adiabatic dynamics of systems with coupled electronic-nuclear degree of freedom is the mixed quantum-classical (MQC) methods , which has been widely used to model polariton relaxation dynamics. Two commonly used MQC methods are the mean-field Ehrenfest method (MFE) and the surface hopping method. These methods treat the electronic DOFs quantum mechanically while propagating the nuclear DOFs classically (while often sampling the initial condition through quantum Wigner distributions). To simulate the dynamics of molecular polaritons, one extends the MQC approaches to treat both the electronic and photonic DOFs quantum mechanically while propagating the nuclear DOFs classically . However, the mixed quantumclassical approximations built into these methods have been known to produce unphysical results, such as the breakdown of detailed balance (that is, the long-term populations) for the MFE method and the introduction of artificial electronic coherence or incorrect chemical kinetics for the surface hopping method .