is the classical nuclear velocity. The matrix elements ϵ α,n , μαβ , η nm , and D αβ were defined previously in Eq. 18. To achieve this expression, one assumes that the second-derivative couplings are vanishingly small and that the nuclear distribution is sharply peaked around its average position. Additionally, we have made use of the orthogonality of the vacuum Fock states, ⟨n |m ⟩ = δ nm as well as noting that the nuclear gradient does not act on the photonic basis states (i.e., ⟨n | R |m ⟩ = 0) since they do not depend on the nuclear positions R. This non-dependency of the photonic basis states leads to a large amount of sparsity in the matrices, and as noted earlier in Section 3.1.1, sparse matrix methods could be employed to achieve more efficient calculations. Note that this is not true, for example, when employing the PFS photonic basis (Section 3.1.2), where the definition of the photonic basis states relies on the nuclear positions through the use of the molecular dipole operator. In order to propagate the classical nuclear DOFs, the classical-like force is simply the sum of all contributions from the elements of the reduced density matrix, ρ αn,β m (t ) = c * α,n (t )c β ,m (t ). In this way, the force can be written as,