where X 0α vo is the transition density matrix similar to that found in Eq. 142 between the ground and α th excited electronic state in the CISapproximation (see additional discussion in Sec. 3.2) and µ ov is the transition dipole between atomic orbitals o and v. From a computational perspective, obtaining both terms in Eq. 161 is not always trivial and may require additional methods such as iterative optimization algorithms (e.g., bi-conjugate gradient optimization) to acquire the individual terms themselves, which adds an additional layer of complexity and consideration when performing on-the-fly NAMD simulations inside the cavity. 243 When generalizing beyond the Jaynes-Cummings model, one needs to additionally account for the excited state permanent and transition dipoles matrix elements. For more complicated excited state methods (e.g., EOM-CC, CISD, etc., see Sec. 3.2), acquiring analytic gradients is not trivial and extremely expensive. However, the analytic expression for the nuclear gradients on the atomic orbital dipoles ∇ R µ and transition density ∇ R X 0α (as well as the bare electronic non-adiabatic couplings d αβ and excited state PES gradients ∇E α ) can, in principle, be achieved analytically in any TD-SCF method and has been shown possible in similar works. However, the implementation of such quantities in commercial or open-source electronic structure packages are few and far between.