, where N el is the number of electronic adiabatic states (i.e., N el = 2 in this example), as defined in Eq. 18. For additional details regarding the nuclear gradients, we refer to interested reader to Ref. 116 and Ref. 154. These nuclear gradients have been the hardest hurdle to overcome in performing on-the-fly simulations with realistic, ab initio systems, since the nuclear gradients on the molecular dipole are not widely available analytically. However, some recent work has been done in achieving MQC dynamics at the Jaynes-Cummings level (i.e., without the dipole self-energy term and ignoring the highly oscillatory terms in the light-matter interaction Hamiltonian) through iterative schemes or through the approximation of the gradient via a Taylor expansion. However, for further theoretical advancement, the community requires work toward the efficient and analytic evaluation of these nuclear gradients to enable another theoretical leap in the simulation of polariton dynamics for realistic systems. Fig. presents an example of the effects of the nuclear gradients on the quantum dynamics of the system. As before, the one-electron Shin-Metiu model is used and provides access to a realistic ab initio system, albeit simplified. Fig. show the bare electronic properties, lowest two potential energy surfaces E α (R ), the non-adiabatic coupling d ge (R ) (Eq. 4), and the molecular dipole matrix elements µ αβ (R ). cavity photon in the four-state basis of |g, 0⟩, |e, 0⟩, |g, 1⟩, and |e, 1⟩. The colorbar indicates the average photon number ⟨ â † â ⟩. The gradients in this basis take the form of Eq. 62. Fig. shows only a few of these matrix elements as a function of the nuclear coordinate. For example, the [V ] e0,g1 ≡ [H pl ] e0,g1 is the off-diagonal gradient of the Hamiltonian shown in Eq. 64, and the other elements are of the quantity