with the polariton wavefunction |Ψ(t)⟩ expressed in Eq. 25, and reduced density ρ ij (t) in Eq. 26. Since most of the photonic character is contributed from the |g1⟩ state, the evolution of the photon number is very similar to that of the population of the |g1⟩ state as shown in Fig. . For a larger light-matter coupling strength (g c = 0.05 a.u.), because of a larger population for |g1⟩ state (Fig. ), the resulting photon number inside the cavity is also larger. Fig. presents the time-dependent population of the trans configurations when coupling the molecule to the cavity. In particular, we compare the case of the trans population dynamics outside the cavity (green solid line) with the case of coupling the molecule inside the cavity, with the coupling strength of g c = 0.005 a.u. (red solid lines) and g c = 0.05 a.u. (red open circles). We have observed a significant suppression of the trans population when the molecule is strongly coupled to the cavity. As such, this simulation demonstrates the suppression of the isomerization reaction from cisto-trans configuration, when considering realistic cavity loss. The reason for the suppression of the cis-to-trans reaction is that more trajectories are staying on the |E 2 ⟩ state after leaving the Franck-Condon region due to the decrease of the polaritonic NAC between the |E 2 ⟩ state and |E 1 ⟩ state, as clearly shown in Fig. . The above results clearly demonstrate that the chemical reaction of the molecular system can be controlled by the lightmatter coupling, and even effectively controlled by tuning the light-matter coupling strength, agreeing with what we have observed in our previous work using a simple model system. Fig. presents the CNNC dihedral angles during the dynamics process to illustrate how light-matter interactions affect the reaction ratio of the azomethane molecule, with the light-matter coupling strength g c = 0.005 a.u. In this figure, we print out the event of the hopping between |E 1 ⟩ polariton state to |E 0 ⟩ polariton state. Here, we randomly selected 100 trajectories out of the 500 trajectories we generated in the simulation. As opposed to the case of outside the cavity (Fig. ) where essentially no trajectory appears around the 0 • CNNC dihedral angle, when coupling to the cavity, more trajectories are trapped around the cis isomer (Fig. ) located at the 0 • CNNC dihedral angle, and the CNNC dihedral angles of these trajectories oscillate centered at 0 • . This indicates that fewer trajectories will go to the trans configuration, resulting in the suppression of the reaction ratio inside the cavity. Fig. presents the same results with a stronger light-matter coupling strength g c = 0.05 a.u., where the fraction of trajectories trapped in the cis configuration region is further enhanced. The reason for the restricted motion of these trajectories is clear. After the trajectory leaves the Franck-Condon region, a large fraction of trajectories will evolve on the |E 2 ⟩ state, as opposed to all moving to the electronic conical intersection between |e⟩ and |g⟩. Note that although the NAC between the |E 2 ⟩ state and |E 1 ⟩ state is large (Fig. ), not all trajectories will hop to the |E 1 ⟩ state immediately when they approach the avoided crossing region between the |E 2 ⟩ and |E 1 ⟩ state. As shown in Fig. , there is a large energy barrier on the |E 2 ⟩ polariton surface (violet curve), with the peak of the potential barrier located at the CNNC dihedral angle ∼ 90 • (because it is largely the |g1⟩ state in that region). Hence, the rotation of these trajectories on the |E 2 ⟩ surfaces from cis to the trans configuration will be suppressed, and these trajectories are eventually forced back to the cis configuration region on the |E 2 ⟩ surface. These trajectories will pass by the |E 2 ⟩ and |E 1 ⟩ avoided crossing region again and either hop to the |E 1 ⟩ state or stay on the |E 2 ⟩ surface. Eventually, they will decay to the |E 0 ⟩ state around the cis isomer due to the cavity loss, giving rise to the hopping event indicated as the open circles in Fig. . Note that there are many number of hops between |E 2 ⟩ and |E 1 ⟩ surfaces as well. For clarity, we have not indicated these hopping events in the current plot. (left column) and g c = 0.05 a.u. (right column). For realistic molecule-cavity hybrid systems, due to the interaction between the cavity mode and the far-field photon modes outside the cavity, the lifetime of the cavity is always finite. In this work, we use the Lindblad dynamics to incorporate cavity loss in the TSH dynamics. To investigate how cavity loss can impact the chemical reaction, we perform the dynamics inside the cavity with different cavity loss rates Γ. For two coupling cases where g c = 0.005 a.u. and g c = 0.05 a.u. the average photon number decrease when the cavity loss rate Γ increases, as clearly shown in Fig. and.