Next, we present the ab initio polariton dynamics simulation results obtained using the TSH method described in Sec. 3.3. In this section, the field is polarized along the Y-axis, which means (ϵ x , ϵ y , ϵ z ) in Eq. 16 is set to be (0, 1, 0). In this case, only the component of dipoles along the Y-axis will contribute to the polariton dynamics. Two different light-matter coupling strengths g c = 0.005 a.u. and g c = 0.05 a.u., are used here. These light-matter couplings cause the Rabi splitting between the polariton states |E 2 ⟩ (commonly referred to as the Upper Polariton state) and |E 1 ⟩ (commonly referred to as the Lower Polariton state), with hΩ R = 68 meV (for g c = 0.005 a.u.) and hΩ R = 680 meV (for g c = 0.05 a.u.). In our simulations, we will consider a range of cavity loss, with the largest one being Γ = 64 meV. Note that strong coupling in cavity QED refers to hΩ R ≪ Γ (as well as larger than the molecular excitation decay rate). As such, g c = 0.05 a.u. always satisfies the strong coupling condition and for g c = 0.005 a.u., one needs to consider the cavity loss rate Γ. Fig. presents the PESs of the polaritonic states (see Eq. 17) along the reaction pathway from the cis isomer to the trans isomer, with the cavity frequency hω c = 2.72 eV. We use two adiabatic electronic states (|g⟩ and |e⟩) and three Fock states (|0⟩, |1⟩ and |2⟩) in the construction of the adiabatic-Fock state basis, which produces six polaritonic states, all of which are included in the dynamics process. The PESs are obtained with the light-matter coupling strength g c = 0.005 a.u. (depicted in solid lines) and g c = 0.05 a.u. (depicted in open circles). Fig. provides a zoom-up of the square region highlighted in Fig. . An avoided crossing exists between the polaritonic states |E 2 ⟩ and |E 1 ⟩ due to the light-matter coupling, and the energy gap (so-called the Rabi splitting) increases along with the increase of the coupling strength. Fig. presents the nuclear gradients at the |E 1 ⟩ and |E 2 ⟩ states, and Fig. presents NAC between them along the reaction pathway. The one-dimension gradients and NAC for each geometry along the reaction pathway are obtained by projecting the gradients and NACs in the whole space (including all nuclear DOFs) to the vector defined by the difference between the Cartesian coordinates of geometries before and after the current geometry. The gradients consist well with the change of the PESs (Fig. ), and the NAC between the |E 1 ⟩ and |E 2 ⟩ states decreases when the light-matter coupling strength increases as expected.