Benchmarking semi-classical simulation schemes for molecular dynamics in the strong coupling regime

Benchmarking semi-classical simulation schemes for molecular dynamics in the strong coupling regime Arun Kumar Kanakati1and Gerrit Groenhof1∗ 1Department of Chemistry and Nanoscience Center, University of Jyv¨askyl¨a, 40014 Jyv¨askyl¨a, Finland Abstract Experiments suggest that collectively coupling molecules to optical cavities can alter their photochemistry. To understand such cavity effects in atomic detail, semi-classical molecular dynamics have been introduced that can model the collective interaction of many molecules with chemical accuracy within the Tavis-Cummings framework of quantum optics. Here, we verify the validity of three such semi-classical approaches, namely Ehrenfest, Fewest-Switches Surface Hopping (FSSH) and the multi-state mapping approach to surface hopping (MS-MASH), for modeling the dynamics of strongly-coupled carbon-monoxide molecules against full quantum dynamics simulations with the multi-configuration time-dependent Hartree (MCTDH) method. The trajectories obtained with the semi-classical schemes agree qualitatively with the results of the quantum dynamics simulations, with FSSH providing the best quantitative agreement when a decoherence correction is applied. Our results thus suggest that semi-classical approaches can provide a valid alternative to quantum dynamics methods for modeling photochemistry in the collective strong coupling regime, for systems that are too large or complex for a fully quantum mechanical treatment. ∗Electronic mail: arun.k.kanak[email protected] gerrit.x.gro[email protected] 1 I. INTRODUCTION Over the past decades, significant efforts have focused on engineering materials capable of precisely controlling the properties of light. The reverse process of controlling the fundamental properties of materials using light represents a far greater challenge but also a profound opportunity. Recent experiments suggest that optical resonators, such as cavities and plasmonic lattices, may provide such control. Indeed, embedding materials inside optical cavities has been demonstrated to reshape their physico-chemical properties, influencing energy transport [1?–9], charge mobility [10–13], lasing thresholds [14, 15], and even photo-chemistry [16–20]. The possibility of harnessing cavity effects to steer photochemical reactions could open the door to transformative applications in artificial light harvesting, energy storage, and quantum technologies. Yet, despite its promise, progress in this new field of polaritonic chemistry is significantly hampered by a lack of theoretical understanding. Optical resonators, such as a Fabry-P´erot microcavity, enhance light-matter interactions by confining the electromagnetic field into very small volumes [21]. If the strength of the light-matter interaction becomes sufficiently high, which can be achieved by increasing confinement or by increasing the collective oscillator strength of the material with more molecules, the system can enter the strong coupling regime, where quantum transitions in the material (e.g., electronic or vibrational transitions) hybridize with the confined photon modes of the optical resonator [22, 23]. These light-matter hybrid states are called polaritons and are characterized by a Rabi splitting of the coupled material’s absorption spectrum into an upper and lower polariton. Because these polaritonic states are coherent superpositions of molecular transitions and cavity modes, an excitation is delocalized over many molecules. In addition to these bright and delocalized polaritonic states, also “dark” states form, which are the remaining superpositions that lack cavity mode contribution. Because changes in material properties have so far only been observed in combination with Rabi splitting, these changes have been attributed to polaritons [23]. However, there is currently no consensus on why polariton formation would affect photochemistry, in particular, because the large majority of states are dark and hence similar to the uncoupled molecular states [24, 25]. The paradox of polaritonic chemistry is that under the collective strong coupling conditions in optical cavities, an excitation is coherently shared by many molecules, whereas 2 chemistry involves only a few. Although the underlying physics of light-matter interactions in the strong coupling regime is known [?], the challenge facing theoretical efforts at resolving this paradox is to model such interactions for a sufficiently large number of molecules with chemical accuracy. To circumvent modeling large numbers of molecules, most theoretical approaches focus on single molecules instead, using much stronger cavity vacuum fields [26–36]. Because the Rabi splitting scales with the number of molecules, N, as √N[22], these fields are enhanced by the same factor to achieve strong coupling with only a single molecule. However, such scaled fields can easily become nonphysical [37], and hence induce changes to the photo-chemistry that are not real. On the other hand, models from quantum optics that focus on describing collective strong coupling in large Nlimit [38–40], tend to lack chemical details that are needed to unravel how the coherent coupling impacts the chemistry locally. To overcome the limitations of these two modeling extremes for modeling electronic strong coupling, a divide-and-conquer strategy was proposed. Combining the Born-Oppenheimer and long-wavelength approximations, polaritonic states are obtained within the TavisCummmings (TC) framework of quantum optics using the adiabatic electronic ground and excited states of the molecules evaluated at a suitable level of quantum chemistry, as the basis [41, 42]. By propagating the nuclear degrees of freedom classically under the influence of the polaritonic wavefunction, while simultaneously propagating that wavefunction as a time-dependent superposition of the TC eigenstates along the classical trajectory, the molecular dynamics (MD) in the collective strong coupling regime can be modeled with chemical accuracy [43, 44]. Through extensive parallelization, such semi-classical MD simulations of thousands of molecules in Fabry-P´erot micro-cavities helped resolve important questions, such as: (i) why, despite the short lifetimes of cavity modes, the polariton appears longlived [43, 45]; (ii) why polariton transport is not ballistic but diffusive [46, 47] and can even reverse on longer timescales [48]; and (iii) what is the role of molecular disorder on polaritonic effects [25, 45]. However, with Newton’s equations of motion governing the propagation of the nuclear degrees of freedom, semi-classical ans¨atze neglect nuclear quantum effects. To understand the impact of this approximation on the computed dynamics under electronic strong coupling, we here compare semi-classical trajectories of polaritonic systems to the results of fully quantum dynamics simulations. To keep the latter simulations tractable, we focus on up to 3 four carbon-monoxide (CO) molecules strongly coupled to a single-mode optical cavity that is resonant with the optical transition into the electronically excited 1Π state. The paper is organized as follows: In section II, we present our model of a CO in a single-mode optical cavity and share the details of our simulations on this system. Then, in section III we compare the observables obtained from semi-classical and quantum dynamics simulations for a varying number of CO molecules coupled to the cavity. Finally, we conclude in section IV with a summary and outlook. II. COMPUTATIONAL DETAILS A. Potential energy surfaces The equilibrium geometry of the electronic ground state (S0) of CO was optimized with the Gaussian09 program [49] at the coupled-cluster singles and doubles (CCSD) level of ab initio theory employing the augmented correlation-consistent polarized valence double zeta (aug-cc-pVDZ) basis set. The S0minimum energy geometry has a C∞vpoint group symmetry with an equilibrium bond length of 1.1405 ˚ A and a fundamental anharmonic frequency of 2172 cm−1. The potential energy curves and dipole moments as a function of the CO bond length were calculated with the MOLPRO program package [50] at the CASSCF(14,10)/aug-ccpVDZ level of multi-configuration self-consistent field theory, with state-averaging over the six lowest-energy singlet electronic states. Energies, dipole moments, and transition dipole moments (TDM) were calculated for eighty inter nuclear distances between 0.7 ˚ A (∼1.3 a.u.) and 2.3 ˚ A (∼4.3 a.u.). The energy profiles for electric ground state (1Σ) and for one of the doubly degenerate electronic excited states (1Π) are shown in Fig. 1a, along with the transition dipole moment between these states [cf. Fig. 1b]. A Morse potential was fitted to each potential energy profile [cf. Fig. 1(a)]: Vi(R) = Ei+Di[1 −exp (−αi(R−Req i))]2(1) with Req iis the equilibrium C-O distance, and Dithe dissociation energy in the ground (i= 0) and excited state (i= 1), respectively. The values of these fitting parameters are 4 FIG. 1. Adiabatic potential energy curves of the bare CO molecule (a), and transition dipole moment curve between the ground 1Σ and the first excited 1Π electronic states (b). PESs of one CO molecule in a cavity: The coupling strength is equal to zero (c) and 0.050 (d). The PES without coupling shows the energy of the ground state shifted to the energy mode of the cavity (i.e. cavity mode excitation, V0+ωc), as well as the ground and excited states of the CO molecule. Panel (b) shows the hybrid states as result of the coupling between the cavity and the molecule. These cuts are obtained by the diagonalization of the cavity and electronic Hamiltonian [cf. Eq. 9]. provided in Table I. TABLE I. Morse potential fitting parameters (in a.u.) for both 1Σ+and 1Π electronic states Parameters 1Σ1Π Ei0.0000 0.3311 Di0.4013 0.1086 αi1.2710 1.4144 Req i2.15 2.42 5 B. Quantum dynamics simulations of strongly coupled CO molecules 1. Molecule-cavity Hamiltonian We consider a one-dimensional array of Nidentical and non-interacting CO molecules inside an optical cavity. The molecular ensemble-cavity Hamiltonian contains a molecular part, cavity part and their interaction: ˆ H= N X k=1 ˆ H(k) mol +ˆ Hcav +ˆ Hcav-mol (2) with ˆ H(k) mol =ˆ T(k) n+ˆ H(k) ethe Hamiltonian of the kth CO molecule, including the kinetic and potential energy operators of the bare molecule: ˆ Hmol =−ℏ2 2µ ∂2 ∂R21+  V0(R) 0 0V1(R) (3) where µis the reduced nuclear mass, Ris the C-O internuclear distance, and V0(R) and V1(R) the potential energy curves in the electronic ground and excited states, respectively (Fig. 1a), fitted to Morse potentials (Equation 1). Within the long-wavelength approximation, the Hamiltonian describing the cavity mode and its interaction with the molecules is [26, 27, 41, 51–54] ˆ Hcav +ˆ Hcav-mol =ℏωc1 2+ ˆa†ˆa+g−→ ˆ D·−→ ϵcˆa†+ ˆa(4) Here, ωcthe cavity mode frequency, ˆa†and ˆathe photon creation and annihilation operators, respectively, gthe coupling strength defined as g=qℏωc 2V ϵ0,−→ ˆ Dthe molecular dipole operator, and −→ ϵcthe cavity mode polarization. For simplicity, we only consider a single polarization direction. 2. Quantum dynamics propagation The time evolution of the molecular ensemble-cavity wave function is computed using the multi-configuration time dependent hartree (MCTDH) approach [55, 56] implemented in the 6 Heidelberg MCTDH program version 86.2 [57]. In this approach, the wave function is approximated as a product of time-dependent coefficients with time-dependent basis functions for each nuclear, electronic, and cavity degree of freedom. Within this formalism, the cavity mode is treated as an harmonic oscillator in terms of the position (ˆqc=pℏ/2ωc[ˆa†+ ˆa]) and momentum operators (ˆpc=ipℏωc/2[ˆa†−ˆa]) rather than the annihilation and creation operators in Equation 4 [53, 54, 58, 59]. To facilitate comparison with semi-classical molecular dynamics simulations, in which the photo-excitation is modeled as an instantaneous population transfer into one of the polaritonic states, we generate the initial state for the MCTDH simulations through application of the operator ˆ T±, which directly excites the cavity-molecule system to a 1:1 light-matter superposition state (with the minus sign in the index for LP and plus for UP) [58] ˆ T±=1 √2ˆa†+ ˆa∓ N X κ 1 √2N(|0κ⟩⟨1κ|+ h.c.) .(5) 3. Wavefunction analysis The populations of the UP and LP states are the key observables characterizing the time-evolution of the hybrid system. A convenient way to extract such information from the MCTDH wavefunction is to compute the expectation value of the light-matter interaction term in the Hamiltonian divided by the interaction strength:[58, 59] Vint,r =⟨Ψ(t)|ˆ O|Ψ(t)⟩(6) with ˆ O=−ˆa†+ ˆaX κ (|0κ⟩⟨1κ|+|1κ⟩⟨0κ|) (7) where |0κ⟩and |1κ⟩are the two electronic states of the κth CO molecule coupled by the optical cavity mode. Because the LP and UP states correspond to positive and negative linear combinations of the electronic and photonic excitations, Vint,r approaches -1 for the LP and 1 for the UP. 7 C. Semi-classical molecular dynamics simulations of strongly coupled CO molecules Starting from the Born-Oppenheimer approximation, we separate the electronic plus cavity mode degrees of freedom from the nuclear degrees of freedom [41]. Neglecting the dipole-self energy, which is very small for realistic cavity setups [37], and adopting the rotating wave approximation (RWA), valid when the cavity mode and molecular excitation are resonant, we can recast the light-matter interaction in Equation 4 into the Tavis-Cummings Hamiltonian [60, 61]: ˆ HTC =ωc1 2+ ˆa†ˆa+g−→ µ·−→ ϵc(ˆa†+ ˆa) (8) Within the single excitation manifold, valid for the weak driving typically employed in experiments, this Hamiltonian can be represented in the basis of product states formed from the adiabatic electronic excitations of the molecules and the cavity mode [42]: H= N X i=1 V0(Ri)1N+         ℏωcγ(R1)γ(R2)··· γ(R1) ∆(R1) 0 ··· γ(R2) 0 ∆(R2)··· . . .. . .. . ....         (9) where 1is a unit matrix of dimension N,V0(Ri) is the ground state potential energy of the ith molecule, ∆(Ri) = V1(Ri)−V0(Ri) is the vertical change in potential energy upon excitation, and the dipole coupling of the ith molecule to the cavity mode is γ(Ri) = g0d01(Ri), d01 being the dipole transition moment of the molecule along the cavity polarization axis. Diagonalizing the matrix representation in the basis of the single-excitation product states (Equation 9) of the multi-scale Tavis-Cummings Hamiltonian (Equation 2) yields the N+ nmode adiabatic hybrid light-matter eigenstates, |ψm⟩: |ψm⟩= N X j βm jˆσ+ j+ nmode X p αm pˆa† p!|ϕ0⟩(10) with eigenenergies Em. The βm jand αm pexpansion coefficients reflect the contribution of the molecular excitations, |Sj 1(Rj)⟩, and the cavity mode excitations, |1p⟩, to polariton |ψm⟩. Due to their parametric dependence on the nuclear degrees of freedom, these eigenenergies 8 form adiabatic potential energy surfaces Em(R), with Rthe 3 ×Nmol ×Natoms coordinates of all atoms in the system [41, 62]. Polaritonic surfaces for a single CO molecule coupled to the cavity are shown in Fig. 1d. In semi-classical dynamics simulations, nuclear trajectories are evolved by numerically integrating Newton’s equations of motion under the influence of the forces due to the quantum degrees of freedom, for which the wave function, |Ψ(t)⟩, is propagated along the classical trajectory. To model the non-adiabatic dynamics in the manifold of eigenstates (Equation 10), we used three popular approaches: (i) Ehrenfest; (ii) fewest-switches surface hopping [63, 64] with and (iii) without decoherence correction [65] and (iv) the Mapping approach to Surface Hopping (MASH) [66, 67]. The details of their implementation for semi-classical molecular dynamics in the collective strong coupling regime are described in Sokolovskii and Groenhof [44]. All semi-classical simulations were performed with GROMACS version 4.5.4, using a timestep of 0.1 fs. For each semi-classical method and single-molecule coupling strength, g, 200 non-adiabatic trajectories were computed with an integration time step of 0.1 fs. The starting coordinates and velocities were randomly sampled from a Wigner distribution at 0 K and observables were averaged over the trajectories. III. RESULTS AND DISCUSSION A. Bare CO dynamics Before focusing on the dynamics of CO molecules in a cavity, we first simulate and compare the dynamics of a single CO molecule in vacuum. In Figure 2a we show the linear absorption spectrum of a bare CO molecule obtained with MCTDH and classical molecular dynamics simulations. In the MCTDH simulations, the ground-state vibrational wavefunction, |ψ(0)⟩was promoted from the electronic ground to excited state and propagated for 200 fs. The overlap, S(t) = ⟨ψ(0)|ψ(t)⟩was computed, multiplied with an empirical damping function e−t/τ (τ= 3 fs) and Fourier transformed into the linear absorption spectrum shown in Figure 2a. The classical absorption spectrum was obtained a sum of the vertical excitation energies along MD trajectories, convolved with either a gaussian (width, σ=0.2 eV) or a Lorentzian (half-width at half maximum, γ=0.2 eV) functions. Whereas both con9 FIG. 6. Same as Fig. 5, when the four CO molecules are interacting with the cavity mode. Finally, we increase the number of CO molecules in the cavity to four. In Fig. 6 we plot the interaction potential, Vint,r for a vertical excitation into the LP r UP states. The overall dynamical behavior follows similar trends to the two-molecules case, as discussed above. However, the collective coupling leads to enhanced stability of the LP state and a faster decay of the UP population, consistent with the emergence of additional dark states in the multi-molecule configuration. These results underline the cooperative nature of molecularcavity coupling and its impact on the redistribution of excitation energy among polaritonic and dark manifolds. While Ehrenfest and MASH FSSH without decoherence corrections now capture the decay from the UP into the dark states, including decoherence corrections to FSSH suppresses this dephasing process. We attribute this to an over-correction of the decoherence, in particular when the non-adiabatic coupling is small. More problematic is the rapid relaxation of population from the LP into the dark states, which is happening at a much faster rate for the semi-classical approaches compared to MCTDH, even for the largest coupling strenghts. 16 This is a major discrepancy, and can lead to qualitatively incorrect results. IV. SUMMARY AND OUTLOOK In this work, we systematically investigated the nonadiabatic dynamics of CO molecules strongly coupled to an optical cavity, combining full quantum and semi-classical molecular dynamics approaches. The light–matter interaction was modeled through a molecule–cavity Hamiltonian, and the MCTDH method was employed to obtain reference quantum dynamics results. The time-dependent populations of the upper and lower polaritonic states, as well as the cavity–molecule interaction potentials, were analyzed to characterize the coherent energy exchange between photonic and molecular degrees of freedom. To assess the performance of computationally efficient semi-classical methods, three widely used approaches Ehrenfest dynamics, FSSH, and MS-MASH were benchmarked against the MCTDH results. The comparison revealed that, while the Ehrenfest approach captures the overall oscillatory features of polaritonic energy exchange, it tends to overestimate coherence loss due to its mean-field nature. The FSSH method, particularly when augmented with decoherence corrections, provided the closest quantitative agreement with the fully quantum simulations, accurately reproducing both the Rabi oscillation period and the population transfer between polaritonic states. The MS-MASH approach also performed well, maintaining a balance between computational efficiency and physical accuracy, though with slightly reduced fidelity in describing population transfer from the UP state. The results demonstrate that semi-classical molecular dynamics methods can reliably reproduce the key dynamical signatures of strong light–matter coupling, including Rabi oscillations, and population exchange, when appropriately parameterized. Importantly, the inclusion of decoherence effects is crucial for achieving quantitative consistency with quantum benchmarks. These findings suggest that semi-classical approaches, especially FSSH with decoherence, can serve as viable alternatives to fully quantum treatments like MCTDH, particularly for larger or more complex systems where full quantum propagation becomes computationally intractable. Looking forward, the methodological insights gained here provide a foundation for exploring collective polaritonic effects in multi-molecule systems, as demonstrated by our simulations of two and four CO molecules interacting with a single cavity mode. The observed 17 scaling of Rabi splitting with the number of molecules confirms the cooperative nature of the coupling and highlights the potential of these approaches for studying collective strongcoupling phenomena. Future extensions of this work will focus on including dissipative environments, and multi-mode cavities, enabling a more comprehensive understanding of polaritonic chemistry and light-induced reactivity in realistic condensed-phase environments. ACKNOWLEDGMENTS AUTHOR DECLARATIONS Conflict of Interest The authors have no conflict of interest to disclose. DATA AVAILABILITY 18 [1] D. M. Coles et al., Nat. Mater. 13, 712 (2014). [2] G. Lerario et al., Light: Science & Applications 6, e16212 (2017). [3] G. G. Rozenman, K. Akulov, A. Golombek, and T. Schwartz, ACS Photonics 5, 105 (2018). [4] K. Georgiou, R. Jayaprakash, A. Othonos, and D. G. Lidzey, Angewandte Chemie 133, 16797 (2021). [5] A. M. Berghuis et al., ACS photonics 9, 2263 (2022). [6] S. Hou et al., Advanced Materials 32, 2002127 (2020). [7] R. Pandya et al., Advanced Science 9, 2105569 (2022). [8] M. Balasubrahmaniyam et al., Nature Materials 22, 338 (2023). [9] N. Krupp, G. Groenhof, and O. Vendrell, Nature Communications 16, 5431 (2025). [10] E. Orgiu et al., Nature Materials 14, 1123 (2015). [11] N. Krainova, A. J. Grede, D. Tsokkou, N. Banerji, and N. C. Giebink, Physical review letters 124, 177401 (2020). [12] K. Nagarajan et al., ACS nano 14, 10219 (2020). [13] P. Bhatt, K. Kaur, and J. George, ACS nano 15, 13616 (2021). [14] S. K´ena-Cohen and S. Forrest, Nature Photonics 4, 371 (2010). [15] T. K. Hakala et al., Nature Physics 14, 739 (2018). [16] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Angewandte Chemie International Edition 51, 1592 (2012). [17] B. Munkhbat, M. Wers¨all, D. G. Baranov, T. J. Antosiewicz, and T. Shegai, Science Advances 4, eaas9552 (2018). [18] K. Stranius, M. Hertzog, and K. B¨orjesson, Nature Communications 9, 2273 (2018). [19] J. Mony et al., Advanced Functional Materials 31, 2010737 (2021). [20] Y. Yu, S. Mallick, M. Wang, and K. B¨orjesson, Nature communications 12, 3255 (2021). [21] K. J. Vahala, Nature 424, 839 (2003). [22] P. T¨orm¨a and W. L. Barnes, Rep. Prog. Phys. 78, 013901 (2015). [23] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, eabd0336 (2021). [24] G. D. Scholes, C. A. DelPo, and B. Kudisch, The journal of physical chemistry letters 11, 6389 (2020). 19 [25] A. Dutta et al., Nat. Commun. 15, 6600 (2024). [26] M. Kowalewski, K. Bennett, and S. Mukamel, J. Phys. Chem. Lett. 7, 2050 (2016). [27] J. Flick, H. Appel, M. Ruggenthaler, and A. Rubio, Journal of chemical theory and computation 13, 1616 (2017). [28] J. Fregoni, G. Granucci, E. Coccia, M. Persico, and S. Corni, Nat. Commun. 9, 4688 (2018). [29] T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio, and H. Koch, Phys. Rev. X 10, 041043 (2020). [30] C. F´abri, G. J. Hal´asz, L. S. Cederbaum, and ´ A. Vib´ok, Chemical Science 12, 1251 (2021). [31] C. Sch¨afer, S. Hultmark, Y. Yang, C. M¨uller, and K. B¨orjesson, Chemistry of Materials 34, 9294 (2022). [32] T. E. Li, A. Nitzan, and J. E. Subotnik, Angewandte Chemie International Edition 60, 15533 (2021). [33] T. E. Li, A. Nitzan, and J. E. Subotnik, The Journal of Chemical Physics 154 (2021). [34] J. Sun and O. Vendrell, The Journal of Physical Chemistry Letters 13, 4441 (2022). [35] T. E. Li, Z. Tao, and S. Hammes-Schiffer, J. Chem. Theory Comput. 18, 2774 (2022). [36] I. S. Lee, M. Filatov, and S. K. Min, Nat. Commun. 16, 4554 (2025). [37] D. F. de la Pradilla, E. Moreno, and J. Feist, arXiv preprint arXiv:2508.00702 (2025). [38] F. Herrera and F. C. Spano, Physical Review Letters 116, 238301 (2016). [39] W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkins, Science 380, 1165 (2023). [40] J. B. P´erez-S´anchez, F. Mellini, N. C. Giebink, and J. Yuen-Zhou, Phys. Rev. Res. 6, 013222 (2024). [41] J. Galego, F. J. Garcia-Vidal, and J. Feist, Physical Review X 5, 041022 (2015). [42] H. L. Luk, J. Feist, J. J. Toppari, and G. Groenhof, Journal of chemical theory and computation 13, 4324 (2017). [43] R. H. Tichauer, J. Feist, and G. Groenhof, The Journal of Chemical Physics 154 (2021). [44] I. Sokolovskii and G. Goenhof, J. Chem. Phys. 161, 134106 (2024). [45] G. Groenhof, C. Climent, J. Feist, D. Morozov, and J. J. Toppari, J. Chem. Phys. Lett. 10, 5476 (2019). [46] I. Sokolovskii, R. H. Tichauer, D. Morozov, J. Feist, and G. Groenhof, Nature Communications 14, 6613 (2023). 20 [47] I. Sokolovskii, Y. Luo, and G. Groenhof, The Journal of Physical Chemistry Letters 16, 6719 (2025). [48] R. H. Tichauer, I. Sokolovskii, and G. Groenhof, Advanced Science 10, 2302650 (2023). [49] M. J. Frisch et al., Gaussian 09 Revision E.01, Gaussian Inc. Wallingford CT 2009. [50] H.-J. Werner et al., Molpro, version 2010.1, a package of ab initio programs, http://www.molpro.net. [51] F. H. Faisal, Theory of Multiphoton Processes, Springer Science & Business Media, 1987. [52] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 114, 3026 (2017). [53] O. Vendrell, Chemical Physics 509, 55 (2018). [54] O. Vendrell, Physical review letters 121, 253001 (2018). [55] H.-D. Meyer, U. Manthe, and L. Cederbaum, Chem. Phys. Lett. 165, 73 (1990). [56] M. Beck, A. J¨ackle, G. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000). [57] G. A. Worth, M. H. Beck, A. Jackle, and H.-D. Meyer, The MCTDH package, version 8.2, (2000), University of Heidelberg, Heidelberg, Germany. H.-D. Meyer, version 8.3 (2002), version 8.4 (2007). O. Vendrell and H.-D. Meyer, version 8.5 (2011)., See http://mctdh.unihd.de. [58] I. S. Ulusoy, J. A. Gomez, and O. Vendrell, J. Phys. Chem. A 123, 8832 (2019). [59] I. S. Ulusoy and O. Vendrell, J. Chem. Phys. 153 (2020). [60] E. Jaynes and F. Cummings, Proceedings of the IEEE 51, 89 (1963). [61] M. Tavis and F. W. Cummings, Phys. Rev. 188, 692 (1969). [62] J. Galego, F. J. Garcia-Vidal, and J. Feist, Nat. Comm. 7, 13841 (2016). [63] J. C. Tully, J. Chem. Phys. 93, 1061 (1990). [64] J. C. Tully, Int.J. Quant. Chem. 25, 299 (1991). [65] G. Granucci and M. Persico, J. Chem. Phys. 126 (2007). [66] J. R. Mannouch and J. O. Richardson, The Journal of Chemical Physics 158 (2023). [67] J. E. Runeson and D. E. Manolopoulos, The Journal of Chemical Physics 159 (2023). 21 Supporting Material: Benchmarking semi-classical simulation schemes for molecular dynamics in the strong coupling regime FIG. S1. Linear absorption spectrum for a resonant excitation into the 1Π electronic state, when the molecule is without the cavity and a single molecule inside an optical cavity at different coupling strengths. TABLE S1. Rabi splitting energy ℏΩR(at the FC geometry), Rabi oscillation period τR, and polaritonic states energies (ELP ,EUP ), when the two CO molecules are coupled with a single cavity mode. g(a.u.) ΩR(eV) τR(fs) ELP (eV) EUP (eV) 0.001 0.0327 126.47 9.6095 9.6422 0.002 0.0653 63.33 9.5932 9.6585 0.003 0.0980 42.20 9.5768 9.6748 0.004 0.1307 31.64 9.5604 9.6911 0.005 0.1634 25.31 9.5440 9.7074 0.006 0.1960 21.10 9.5276 9.7237 0.008 0.2614 15.82 9.4948 9.7562 0.010 0.3267 12.66 9.4619 9.7886 0.020 0.6536 6.33 9.2967 9.9503 0.030 0.9809 4.22 9.1302 10.1110 0.040 1.3086 3.16 8.9623 10.2709 0.050 1.6369 2.53 8.7930 10.4299 22 FIG. S2. Expectation value of the cavity photon number ⟨Nph⟩=⟨Ψ(t)|ˆa†ˆa|Ψ(t)⟩as a function of time. FIG. S3. Linear absorption spectrum for an excitation into the molecular 1Π electronic excited state (i.e., bare CO molecule, orange color) and the LP/UP state, when a single molecule inside an optical cavity at different coupling strengths, g. 23 FIG. S4. PESs of two CO molecules in a single cavity mode. The coupling strength considered here is 0.050. The PESs without coupling shows the energy of the ground state shifted to the energy mode of the cavity (i.e. cavity mode excitation, V0+ωc), as well as the ground and excited states of the CO molecule. The hybrid states as result of the coupling between the cavity and the molecule. These cuts are obtained by the diagonalization of the cavity and electronic Hamiltonian [cf. Eq. 9]. 24