Using these two electronic states as well as two photon number states, |0⟩ and |1⟩, as a basis, the four polaritonic Born-Oppenheimer surfaces were constructed at a light-matter coupling strength A 0 = 0.01 (Fig. ) with the curves colored by their average photon number, ⟨ â † â ⟩ (see colorbar). The cavity frequency is ω c = 0.1 a.u. After vertical excitation from the R = -4.0 a.u. to the |e, 0⟩ basis state (i.e., a linear combination of polaritonic states), the quantum dynamics of the system was performed using the (Fig. ) MQC spin-LSC, (Fig. ) MQC Ehrenfest, and (Fig. ) MQC FSSH approaches. The MQC results are shown in circles while the exact results are shown as solid curves. For these MQC simulations, the exact nuclear gradient was implemented (as discussed in the previous section, Sec. 4.3) in order to compute the exact polaritonic nuclear force. At short times (< 20 fs), the nuclear wavepacket starts in and quickly leaves the region of strong light-matter coupling, showcasing a fast decay in coherent oscillations from the light-matter coupling. At ∼20 fs, the nuclear wavepacket crosses the electronic non-adiabatic region near R = 2.0 a.u. where much of the population decays to the total ground state, |g, 0⟩, before the simulation ends. In this example, spin-LSC MQC dynamics performs better than either Ehrenfest or FSSH. It should be noted here that FSSH can generally be improved by using various decoherence corrections. Specifically, spin-LSC is able to capture the weak population of the |e, 1⟩ state (green) while Ehrenfest (and to a lesser extent FSSH) dramatically over-populates this state. Furthermore, Ehrenfest and FSSH do not accurately predict the electronic non-adiabtic crossing time, instead showing a delay of 5-10 fs.