First, we present the population dynamics of the SM1 model in a cavity compared between 2-level's and 4level's results of the exact quantum simulation in Fig. ) along with the polaritonic surface in Fig. . For all simulations, we start with |e0 state and a light-matter coupling strength of g c = 0.005. Fig. presents the first four polaritonic surfaces of the SM1 model coupled to a resonance optical cavity, where the potential is color coded based on the expectation value of â † â indicated on top of this panel. Note that this should not be viewed as a "photon number" operator under the dipole gauge used in the PF Hamiltonian, becasue the rigorous photon number operator should be obtained by applying the Power-Zienau-Woolley (PZW) Gauge transformation on the photon number operator â † â. Nevertheless, it can be viewed an approximate estimation of the photon number when the lightmatter couplings are not enter into the ultra-strong coupling regime. Fig. presents the exact quantum dynamics within the JC subspace {|e0 , |g1 } (dashed lines) or in the larger {|g0 , |e0 , |g1 , |e1 } (solid lines) sub-space , under which the polaritonic eigenstates and dynamics are converged. Clearly, additional states beyond the JC subspace will be explored by the quantum dynamics of the hybrid system (b) due to the NACs among these states. The population of the |e1 state is mainly contributed from the |e0 state due to the light-matter coupling carried from the permanent dipole µ ee . In addition, the population transition between |e0 and |g1 is mediated by the cavity-induced coupling near R = 0, where the PES exhibits a strong mixing between these two states as shown in Fig. at R=0. Fig. presents several X ij components of the general gradient expression (Eq. 39) in the larger subspace {|g0 , |e0 , |g1 , |e1 }. Comparing the gradient expression in the 2-state (Eq. 35) and 4-state subspaces (Eq. 38 and Eq. 39), the [X] terms (Eq. 39) do not show up in the 2-level JC case (Eq. 35). The JC gradient component ∇ [V] e0,g1 = g c • ∇µ ge and the regular gradient from the electronic derivative couplings X g0,e0 = d ge (V e -V g + D 2 ee -D 2 gg ) are compared with the rest of the non-JC type of gradients, such as X g0,g0 = 2d ge D 2 eg (where D 2 eg is the matrix of the DSE operator), X g0,e1 = d ge g c • (µ eeµ gg ), and X g0,g1 = 2d ge g c • µ ge . As we show in Fig. , all of the gradients we discussed above have a similar magnitude in this Shin-Metiu-cavity model. Some elements can be even larger than ∇[V ] e0,g1 , such as X g0,e0 and X g0,g1 . All of these non-JC type of gradient terms, unfortunately, are not included in the previous MQC studies of polariton quantum dynamics based upon the JC (or the TC ) models. When the system starts to explore all of the states and generate sizable population and coherences (ρ ij = c * i c j , see Eq. 16), in the nuclear forces (for example, in Eq. 23) from these new gradient terms will make a non-negligible contribution which is required to be explicitly and correctly included.