Below, we briefly summarize the details of performing the polariton dynamics. First, we use Eq. 53 to generate the initial nuclear condition, including the initial positions and momenta. Then, we use Eq. 27a and Eq. 27b to propagate the quantum dynamics of the electronic-photonic DOFs, while the nuclear DOFs are propagated with force from the active polaritonic state calculated using Eq. 33. In the construction of the system's Hamiltonian and calculation of the nuclear forces, we use the machine learning model to get the dipoles and their derivative using Eq. 49 and Eq. 52, respectively, whereas all of the other quantities, including all energy, gradients, and derivative couplings are computed on-the-fly using CASSCF. During the TSH dynamics, we use Eq. 39 to calculate the hopping probability from a current active state to any other states, and use the Monte Carlo algorithm presented in Eq. 40 to stochastically determine which state the system will hop. Finally, we use Eq. 44 and Eq. 46 to obtain the population of the adiabatic Fock state and polaritonic state, respectively. All the dynamics calculations were performed with a development version of the JADE-NAMD package. presents the electronic adiabatic PESs for |g⟩ (red curve) and |e⟩ (blue curve) along the reaction path, with geometries obtained using procedure described in Sec. 3.1. The inset depicts the molecular structures of azomethane. Two isomers, cis and trans, exist for this molecule, defined by the CNNC dihedral angle at ∼ 0 • and ∼ 180 • , respectively. These two isomers can convert to each other upon photoexcitation. To explore how the photoisomerization process of azomethane can be influenced by the light-matter coupling, we first perform the non-adiabatic dynamics of azomethane outside the cavity. Although the trans isomer is energetically more stable (with a lower energy of ∼ 0.5 eV) than the cis isomer in the ground state (|g⟩), we start the dynamics from the cis isomer upon a Franck-Condon excitation, due to a larger transition dipole between the ground state (|g⟩) and the first excited state (|e⟩) at the cis configuration compared to that of the trans isomer. Fig. presents the time-dependent adiabatic population in the TSH dynamics process. After being excited to the |e⟩ state, the system decays fast and nearly half of the trajectories hop to the |g⟩ state within ∼ 50 fs. The fast excited-state population decay process can be explained by the shape of the PES of the |e⟩ state (see Fig. ), where no energy barrier exits along the reaction pathway from the cis isomer to the CI. After 50 fs, the decay of the |e⟩ state becomes slower and almost all the trajectories are at the |g⟩ state after ∼ 330 fs. Fig. presents the time-dependent population of the trans configuration during the photo-excitation dynamics. In the beginning, all trajectories have the configurations of the cis isomer, and thus the trans population is zero. This is correlated with the CNNC dihedral angles presented in Fig. , where all trajectories are ∼ 0 • at t = 0. Then the system starts to rotate through the central NN bond and isomerize, with the CNNC dihedral angle passing by 90 • or -90 • , and the population of trans configuration almost increases to one at ∼ 70 fs. Many trajectories hop to the ground state in the vicinity of CI, where the CNNC dihedral angles are ∼ 90 • or ∼ -90 • , as shown in Fig. . The trajectories hop to the ground state and then get trapped in the potential minimum of the trans isomer configuration (see Fig. ). These trajectories will oscillate around 180 • or -180 • for the long-time dynamics. Other trajectories will keep rotating, and go back to the cis configuration with the CNNC dihedral angles passing by 270 • or -270 • . Similarly, they are also trapped in the potential minimum of the cis isomer and oscillate around 360 • or -360 • . The rotation from the trans isomer back to the cis isomer of some trajectories explains why the population in Fig. decreases after ∼ 70 fs. The population of trans configuration does not change too much after ∼ 300 fs, where most of the trajectories are already at the |g⟩ state, as shown in Fig. .