Photochemistry in the strong coupling regime: A trajectory surface hopping scheme

Abstract The strong coupling regime between confined light and organic molecules turned out to be promising in modifying both the ground state and the excited states properties. Under this peculiar condition, the electronic states of the molecule are mixed with the quantum states of light. The dynamical processes occurring on such hybrid states undergo several modifications accordingly. Hence, the dynamical description of chemical reactivity in polaritonic systems needs to explicitly take into account the photon degrees of freedom and nonadiabatic events. With the aim of describing photochemical polaritonic processes, in the present work, we extend the direct trajectory surface hopping scheme to investigate photochemistry under strong coupling between light and matter.


| INTRODUCTION
The coherent interaction between light and matter in confined systems offers an alternative pathway to tailor optical and chemical properties of molecules. While the spectroscopy of atoms and molecules in resonant cavities is well established, the possibility to manipulate the molecular reactivity through quantum coupling with light has only recently been addressed. By devising microcavities [1,2] and nanocavities, [3,4] the experimental efforts [5,6] to bring molecules in the strong coupling regime down to the single molecule level have driven an increasing theoretical interest. [7][8][9][10] Yet, the modeling of such complex systems experiences limitations both theoretical and computational.
Understanding which approximations can hold for a correlated nuclear-photonic-electronic system is indeed challenging. [11,12] Even more, an important option is whether to couple the photonic degrees of freedom to the nuclear ones or to the electronic ones. [8,13,14] Within the first approach, the photonic degrees of freedom are treated so as the nuclear ones, allowing to study the effect of the electron-nuclei-photon coupling on adiabatic potential energy surfaces. Such approach provides insightful tools of analysis for phenomena like Raman Scattering, [15,16] modified molecular properties, [14,17] and ground state reactivity. [18][19][20] Instead, the second approach, in which electronic and photonic states are mixed, is suitable to describe the modified photochemical properties [8,21] and reactivity, [22][23][24] provided that nonadiabatic couplings are taken into account.
A full quantum approach has been developed by Rubio's group in the DFT framework. The method is based on rewriting the DFT formulation in terms of a current density functional which allows to include the photonic degrees of freedom [10,25] (QEDFT). Later on, the same group reformulated the Born-Oppenheimer approximation to partially decouple the nuclear-photonic-electronic problem with the so-called Cavity Born-Oppenheimer approximation. [13,14] Such works opened a way to a full ab-initio investigation of strongly coupled light-matter systems, [10,17,26] with successful applications in strongcoupling modified properties of single and many molecules.
Aiming to investigate polaritonic photochemical reactions, the complexity of the system can quickly become cumbersome. The correct computation of excited states is mandatory, together with the treatment of the photonic degree of freedom. [27][28][29] Further complexity to the problem is added by interfacing a propagation scheme for the nuclei [30] and by accounting for environmental effects. In addition, the common problems encountered in photochemical simulations [31] are directly transposed to the study of polaritonic photochemical reactions.
A pioneering conceptual display of novel photochemical events in the strong coupling regime is offered by the works of Mukamel and collaborators [7,32,33] and Feist and collaborators [8,22,23] on model molecules. Such works collect a plethora of insights for a novel chemical reactivity ranging from single-molecule strong coupling to collective strong coupling effects. We have also recently shown how moving beyond model treatments to investigate polaritonic chemistry can also reveal noteworthy effects like enhanced photoisomerization quantum yields. [34] To simulate mixed light-molecule systems, a toolbox of strongcoupling techniques for photochemistry has been developed in the last three years. [24,35,36] Among them, we mention the multiscale MD approach devised by Groenhof and collaborators, which allows to investigate the collective polariton behavior in biological environments through a QM/MM approach. [35,37,38] For events occurring in small ensembles of realistic molecules in cavities on a shorter timescale, the extension of the MCTDH technique to polaritonic systems [36,39] is also remarkable. In recent works, [24,34] we showed how the Surface Hopping scheme in a semiempirical framework can be used to describe the photochemistry of molecules in a strong coupling environment with a high level of realism. A similar Surface Hopping scheme has been used by Tretiak et al., [40] which studied the stilbene photoisomerization under strong coupling, employing a single reference quantum chemical approach for the electronic states. Comparatively, in our scheme, we also include cavity losses, and our semiempirical multireference FOMO-CI scheme allows for a qualitatively correct description of potential energy surfaces and couplings, which is also quantitatively accurate since we reparametrize the semiempirical Hamiltonian.
In the present contribution, we show in detail the theoretical approximations and the implementation techniques of our approach to polaritonic chemistry. To this aim, we show first how polaritonic states are built on top of the semiempirical FOMO-CI [41][42][43][44] technique for the computation of electronic states. Then, we derive the analytical gradients for the strong coupling contribution to the CI energy via the Z-vector [45,46] algorithm. We also discuss the interface with the on-the-fly Direct Trajectory Surface Hopping (DTSH), with emphasis on the method we have adopted to include the effect of cavity losses on the dynamics.
We want to stress that the hereby presented gradients and Surface Hopping interface have a general applicability to multiconfigurational wavefunction methods. The choice of a semiempirical approach to solve the electronic problem resides in the good compromise between efficiency and accuracy. [31] We also mention that such approach has been successfully applied to deal with the molecular complexity of polaritons when all the degrees of freedom are taken into account, [24] and also in the presence of an environment [34] inspired to a realistic setup. [5] We also stress that, while our method carries the potential to treat a few chromophores, the study of a large ensemble of molecular systems is beyond the aim of the present work.

| Polaritonic wavefunction in a semiclassical framework
To build polaritonic states, we consider a generalized correlated photon-electron-nuclear system: where the electronic degrees of freedom are described by the e subscript (r coordinates), the nuclear ones by the n subscript (R coordinates), and the photon one by the ph subscript (q coordinates). The total wavefunction of the correlated electron-photon-nuclei system is Ψ(r, R, q). Two approaches to approximate the eigenstates and the time-evolution of strongly coupled systems have been applied so far: the first is to embed the photon degrees of freedom into the nuclear wavefunction [14] while the second is to embed the photon into the electronic wavefunction. [8] Such two different approaches provide different insight on two classes of processes. In fact, the molecular properties and the dynamics in the Cavity Born-Oppenheimer approximation [13,17] are optimally described by incorporating the photon in the nuclear wavefunction (Ψ n + ph, e ). Instead, the processes involving nuclear dynamics on polaritonic states, that is, photochemical processes, are accurately described by considering hybrid electron-photon states (polaritons, Ψ n, e + ph ). [7,24,32,34,47] The Born-Huang factorizations of the wavefunction in these cases respectively correspond to: Equation (2) represents the case where the photon degrees of freedom are considered slow. Hence, they are treated alike to the nuclear degrees of freedom in the Cavity Born-Oppenheimer framework. [13,17] Based on this assumption, the purely electronic wavefunction and the related electronic potential energy surfaces show a parametric dependence on both the nuclear and photonic coordinates. This framework explicitly requires to compute the quantum nuclear wavefunction to include the photon effects, hence it is not properly interfaced with semiclassical methods developed treating the whole nuclear dynamics as classical.
In the factorization presented in Equation ( In that case, the nuclei are moving according to a classical trajectory R (t), and the polaritonic nonadiabatic couplings can be included as for the purely electronic case. In the semiclassical case, we define by analogy with Equation (3) the polaritonic wavefunction: where jAi are the semiclassical analogous of the polaritonic states ϕ ph + e k of Equation (3). We choose jAi to label such states to directly refer to their adiabatic behavior, as they are the eigenstates of the Here, b H el is the standard electronic Hamiltonian and b H ph is the Hamiltonian of the quantized electromagnetic field (we consider here a single mode for the field), where b b † , b b are the creation and annihilation operators for the electromagnetic field. In principle, the approach considered in this work could be extended in a straightforward way to consider several cavity modes. However, it is uncommon that many modes can reach a coupling strength large enough to require a strong coupling treatment, not to mention that they may also be well separated in energy. As interaction Hamiltonian H int , we take the dipolar light-matter interaction in the Coulomb gauge and long wavelength approximation: In the light-matter interaction, we refer to ℰ 1ph as the 1-photon field strength, with the electromagnetic field polarization λ. b μ tr is the electronic transition dipole moment between the electronic states.
Notice that b H pol is parametrically dependent on the nuclear coordinates through b H el and b H int . As we have numerically verified in previous works, [24,34] for the case of strong coupling with optical frequencies, it is enough to restrain to the transition dipole operator. In the next section we focus term-by-term on the two individual subcomponents of the polaritonic Hamiltonian, namely b H el and b H ph .

| FOMO-CI wavefunction and uncoupled states
The method for the computation of electronic states, that is, the eigenstates of b H el , is based on floating occupation of molecular orbitals (FOMO). [42,43] This variant relies on the optimization of a single determinant wavefunction with fractional variational occupation of the molecular orbitals through a self-consistent field calculation (SCF). The single-determinant SCF calculation is formally closedshell. Here, the energy of the i-th orbital (φ i ) is the Fock eigenvalue ε i corresponding to that orbital, while the occupation number O i of φ i is obtained integrating a function f i (ε) normally distributed along the energy axis around ε i : Here, σ is an arbitrary parameter and the Fermi energy ε F is determined by imposing that the sum of the orbital occupation numbers O i equals the total number of electrons. The Fock operator b F is obtained from the density (orbitals are considered as real functions in the pre- Through this procedure, the lower virtual orbitals can be popu- Hamiltonian. [41] Notice that the standard semiempirical parameters are normally determined to reproduce ground state properties, with a SCF wavefunction. Therefore, to deal with excited states, a reparametrization is often mandatory, as what has been done in Reference [41]. As we adopt a CI-type wavefunction, the (approximated) eigenstates jni of b H el and the corresponding eigenenergies U n are obtained by diagonalizing the electronic Hamiltonian matrix: b on the basis of a chosen set of N CI Slater's determinants {Φ}, so that Similarly to the electronic states, the photon states are the The meaning of p is a photon occupation state number, for the single electromagnetic mode of frequency ω ph considered here.
The product states between the electronic and photonic eigenstates jn, pi are then the eigenstates of the light-matter non- We shall address to them as uncoupled states through all the present work. Such set of uncoupled states jn, pi are the polaritonic equivalent of, for example, the spindiabatic states for the purely electronic case with spin-orbit coupling. [48]

| Polaritonic states evolution and energies
The time evolution of the wavefunction is performed in terms of the polaritonic adiabatic states jAi, which are obtained by diagonalization of the matrix of b H pol (Equation (6)) on a selected subspace of N × (p max + 1) uncoupled states {n, p}, where N ≤ N CI is the number (usually small) of electronic states considered, and p max is the maximum value of the photon occupation number. The set of adiabatic states jAi is used to perform the time evolution as the surface hopping approach is representation-dependent, and usually performs better in the adiabatic basis. However, the set of uncoupled states jn, pi is still useful, mainly in order to ease the interpretation of the results.
Within the described framework, the polaritonic wavefunction evolves according to the "polaritonic TDSE" iℏ _ where and G AB is the derivative coupling vector between the polaritonic states jAi and jBi, namely According to Equation (13), a polaritonic state can be written as and its energy is where the contribution of the uncoupled part can be extracted by exploiting Equations (14) and (18), resulting in: The interaction term E A int is given by where we used the shorthand Notice that D(A j m, n) = D(A j n, m).
When m < n the process described is the molecule exchanging the photon of frequency ω ph with the cavity. The rate of such exchange is the Rabi splitting (Jaynes-Cummings Hamiltonian). [49,50] In this regime, the emission rate and efficiency is greatly enhanced through the Purcell effect [47,51] and the energy is coherently exchanged between matter and cavity. Such energy contribution is the Rabi splitting. Instead, when m > n, the so-called counter rotating terms account for the simultaneous creation/annihilation of two off-resonant excitations within the cavity.
Such terms become non-negligible, together with the dipolar self-energy of the molecule, in the ultrastrong coupling regime. [11,12] From now on, we will use i, j, … to label CI-active molecular orbitals (MO) and a, b for any kind of MO. A more appealing expression, from the computational point of view, of the interaction energy E A int is obtained by using the spinless electronic density matrix, suitably modified, of the polaritonic state jAi considered. In particular, we have where μ ij = i h jλÁb μ j j i and Δ el ij m, n ð Þ is the spinless transition density matrix between the electronic states m and n, expanded on the molecular orbital basis.
The action of the bosonic creators and annihilators of Equation (8) is embedded into the D(A j m, n) coefficients. Therefore, Δ el ij m, n ð Þ is purely electronic: Within our method, we are able to compute the Polaritonic Potential Energy Surfaces (PoPESs) up to an arbitrary occupation number of the photonic mode involved in strong coupling.
We shall now briefly discuss the dependence of the PoPESs on the molecular transition dipole moments. Upon diagonalization, a crossing of the uncoupled states PES is converted to a polaritonic avoided crossing. The magnitude of the splitting (Rabi splitting) is proportional to the transition dipole moment between the crossing states, potentially reaching zero for vanishing transition dipole moments (polaritonic conical intersection). The strong dependence of the Rabi splitting on the transition dipole moment also embodies a strong dependence on the nuclear geometry at which the crossing between uncoupled states occurs, as the transition dipole moments variation with nuclear geometry may be large. Note that the orientation of the molecule here plays the same role as the internal coordinates, because it affects the projection of the transition dipole on the field polarization vector.
While the polaritonic conical intersection and avoided crossings have been reported in previous works, [8,24,39,52,53] here we stress that they are an easy-to-predict feature only when limited to two level strong coupling models, that is, Jaynes-Cummings like. Two-level models imply a linear dependence of the Rabi splitting on the coupling constant  Figure 1a,b for the azobenzene molecule.
We stick to azobenzene as a test case, since the phenomenology of polaritonic photochemistry has been investigated in recent works. [24,34] In the present work, instead, we focus on discussing the change of shape of the polaritonic avoided crossing regions, com- nanocavities, [5] the extreme limit of $1 nm 3 has been accessed experimentally via single-atom hotspots. [54,55] In all the conditions examined in this work, the mode volume is enough to fully embed the molecule (the molecular volume being $0.25 nm 3 ). A few works pioneer the interaction beyond the dipolar approximation for small mode volumes for TERS experiments, [56,57]

| Analytical gradients for CI-expanded polaritonic states
After showing the strong coupling contribution to the energy in the previous section (Equation (24)), here we derive the analytical energy gradient with respect to the nuclear coordinates R α for a generic FOMO-CI expanded polaritonic state. The present approach is based on previous works, [43,58] where the Z-vector method has been applied.
In particular, here we adapt to the polaritonic case the "contracted" strategy that was developed in a spin-orbit framework. [44] As in References [44,58], only the active MOs are allowed to have floating occupation numbers. The gradient of the energy can be partitioned in a response term, containing the derivatives of CI and MO coefficients, and a static term. The static contribution specific to the present case is given by the derivative of the molecular dipole operator matrix elements in terms of atomic orbitals (AO). As for the response terms, notice that the derivatives of the expansion coefficients D A n,p of the polaritonic adiabatic state jAi give a null contribution to n,p = 0 by construction. As a consequence, since E A ph does not involve geometry dependent quantities other than the D A coefficients (in the long wavelength approximation), it does not contribute to the gradient and will not be considered further here. At variance, the derivatives of the electronic CI coefficients C I, n have to be considered.
We have then The gradients for the electronic energies U n entering ∂E A el ∂Rα are known. [43,44,46,58] Hence, here we only show explicitly the evaluation of ∂E A int ∂Rα . By making use of Equation (24) one gets MO set (c is real orthogonal in the semiempirical framework considered here). We have then μ = c † μ AO c, where μ AO is the matrix of b μ λ in the AO basis. Therefore, the derivatives of μ can be expressed as which can be decomposed into a static part and a response part, [43,44,58] ∂μ The static term (Equation (30)) is easily evaluated as follows. Let The dipole matrix elements μ στ are then, in the semiempirical framework where −e is the electronic charge and the Kronecker delta δ αβ is due to the semiempirical NDDO approximation. Moreover, the term is independent on the nuclear coordinates. Therefore, the derivative of μ ! στ with respect to R ! α vanishes unless the two atomic orbitals σ and τ are both centered on the nucleus α, and in that case it simply In an ab initio context, one would have to compute also the derivative of the dipole matrix elements between atomic orbitals centered on different atoms, which has a more involved expression with respect to the term considered here. However, that would not be expected to have a large impact on the computational cost, which is mainly influenced by the response part of the gradient.
The contribution of the response term of Equation (31) to the derivative of E A int (Equation (28)) can be recast in this way, following Patchkovskii and Thiel [45,46] where ε i is the energy of MO i and As q int ii = 0, the term containing the derivative of the orbital energy ε i gives a null contribution to the sum of Equation (34). Such term has been included to recover the same formalism of previous works. [43,44,58] We now turn to the derivative of ρ int ij A ð Þ, which is a response term (i.e., the CI response contribution to the polaritonic energy), evaluated by taking the derivative of Δ el ij m, n ð Þ. Such derivative is obtained by following the same procedure outlined for the MOs response terms (Equation (31)) Notice that the sum in Equation (36) is extended to N CI rather than to the number N of states selected: in principle, the evaluation of the CI response contribution requires the full diagonalization of the CI space considered. While this may be too demanding in an ab initio context, normally it does not represent a problem in a semiempirical framework, where N CI is usually small. The antisymmetric matrix d α nm , expressing the response of the CI coefficients, represents the CI contribution to the derivative couplings. We have then where According to Reference [44], we evaluate the derivative coupling terms d α mn by exploiting the Hellmann-Feynman theorem for m 6 ¼ n, and d α are the twoelectron density matrices and the terms ε + ij are defined in equation (36) of Reference [44].
Inserting Equation (41) into (38), we obtain the following expression for the CI response term induced by the strong coupling interaction: where Here, Δ symm ij k, n ð Þ is the symmetric part of Δ el ij k, n ð Þ. Notice that it is symmetric with respect to both i, j and k, n indices, since Δ el ij k, n ð Þ= Δ el ji n, k ð Þ.
To obtain the final expression for the gradient of E A pol , we have also to consider the contribution given by the derivative of the naked electronic state energy U n (see Reference [58]) By putting all the terms together we arrive at where we made use of the modified electronic density matrices The evaluation of the gradient of E A pol can proceed in the way outlined in Reference [58], using the modified density matrices Δ pol ij A ð Þ and Γ pol ijkl A ð Þ. In particular, the response term is where q el , defined as done in [58] (see also [44] ), is explicitly reported below for reader's convenience In the above equations, we used the shorthand hijkkli = 2hij| kli − hik| jli, and f i is the Gaussian function defined in Equation (10).
Finally, for the static part, one has just to add the last term of Equation (48), representing the static dipole derivative (see above).

| Surface hopping
In the framework of Direct Trajectory Surface Hopping, the formulation of strong coupling given in this work allows to include the decoherence corrections [59] and environmental effects through the QM/MM interface previously devised. [60][61][62] For the time evolution of the polaritonic wavefunction, we adopt the local diabatization technique, [48,63] with a recently improved evaluation of transition probabilities. Such probabilities are compliant with Tully's Fewest Switches prescription and particularly effective when many states are involved in nonadiabatic events, [64] as commonly happens in single-molecule polaritonic systems (see Figure 1b).
As a test case, we examine the azobenzene strong coupling dynamics with ω ph = 2.7 eV and ℰ 1ph = 0.004 au ($12 nm 3 ) in the absence of the cavity losses. The initial conditions are sampled on the ground state via a 20 ps dynamics, thermostated at room temperature (with a Bussi-Parrinello thermostat [65] ). In particular, 230 starting structures and velocities are extracted from the sampling dynamics, and the system is initially vertically excited to the jR 8 i state, that is mostly jS 2 , 2i state. Rather than the simulation of a realistic excitation (the transition jS 0 , 0i ! j S 2 , 2i would require a multiphoton pumping), this is a test case to investigate the effect of photon occupation numbers greater than 1 (up to p = 3). In Figure 3a, we show the behavior of the population of the photon states during the dynamics, in the absence of cavity losses. The blue line with circle markers (right y scale) shows the total photon number within the cavity, The full lines (left y scale) show the populations of each photon state, that is, P n |D n, p | 2 , with p = 0, …, 3. While no cavity loss is explicitly included in the dynamics, still the total photon number in the system decreases. Through strong coupling, a photon is continuously exchanged between the molecule and the cavity. However, the electronic component keeps decaying via internal conversion, meaning that when the photon is absorbed, its energy can be redistributed to the nuclear degrees of freedom. While the total number of photons decreases in the ongoing dynamics, the energy of the system is still conserved (Figure 3b). The initial conditions are chosen such that the resonant region between the p = 3 states (namely S 0 , 3) and the p = 2 states jS 1 , 2i is interested, so both the subspaces jn, 3i and jn, 2i are populated. This results into an average photon occupation number greater than j2i, namely 2.23.

| Cavity losses
Aiming to provide a realistic model, we deal with the issue of lossy cavities. The strong coupling regime for single molecules is usually reached by exploiting a nanocavity setup of the system. [5,6,66,67] The typical lifetime of nanocavities is few tens of femtoseconds. However, we have recently shown that the overall photon lifetime of the system is way longer than the individual cavity lifetime, [34] reaching a time scale comparable to several ultrafast photochemical processes. This effect is due to the transient passage of the wavepacket through strongly coupled regions, so that the composition of the polaritonic state keeps interchanging between electronic and photonic. As a consequence of the mixing, the lifetime of states with the photon partially absorbed is extended up to hundreds of fs, depending on the strong coupling conditions. Within our model, we adapt a quantum jump algorithm [68][69][70][71] already exploited in the Stochastic Schrödinger Equation (SSE) framework [72] to account for relaxation and dephasing channels. Stochastic methods in the framework of SSE are also commonly exploited as an equivalent alternative to master equations in treating cavity losses. [73][74][75] We then follow a standard implementation of this approach, similar to others already present in Quantum Optics simulation packages like QuTip. [76] The quantum jump is a natural choice as it fully exploits the trajectory-based machinery of the surface hopping. Having to deal with semiclassical trajectories, both the polaritonic wavefunction and the "current state" (i.e., the adiabatic state on which PES the nuclei are evolving) must be taken into account whenever a photon loss occurs.
We start with the expression of the polaritonic wavefunction in terms of the uncoupled states basis: where d n,p = P A D A n,p C A . Only states with free photons can decay via cavity losses, namely states with p ≥ 1. We evaluate the photon loss probability P dec by taking the squared modulus of the uncoupled states coefficients with p ≥ 1 in the total wavefunction Ψ pol , that is, d n, p ≥ 1 : Here, Δt is the integration time step and τ is the cavity lifetime, namely the inverse of the cavity decay rate κ. We generate a uniform random number within the interval [0, 1]. If the random number falls in [P dec , 1], the photon is retained and the cavity loss does not occur.
If not, the photon is lost. Upon photon loss, the photon occupation number is lowered by 1 via application of the projector b P which includes the photon annihilation operator b b: F I G U R E 3 Photon statistics and energy conservation. (a) Dynamics of the photons in the cavity during the internal relaxation of the stronglycoupled azobenzene molecule, in absence of cavity losses. The dynamics is started from the jR 8 i state and runs for 1 ps, with ℰ 1ph = 0.004 au ($12 nm 3 mode volume), ω ph = 2.7 eV and longitudinal field polarization. The molecule is in gas phase. [24] The curves with full lines show the dynamics of each p subspace, while the light blue line with circle markers (with the scale on the right) represents the total photon number within the cavity. Error bars, represented as lighter bands, are also shown. Even in absence of cavity losses, the average photon number decreases during the dynamics. While the photon is in its absorbed state, the energy stored within the molecule is redistributed via internal conversion to nuclear kinetic energy. The overall process is still conserving the energy, as shown in panel (b) [Color figure can be viewed at wileyonlinelibrary.com] Here, the projector b P preserves the electronic coherence within each p subspace, apart of course for p = 0 that is annihilated. The photonic annihilation operator b b is applied to mimic the loss of the photon from the cavity, resulting in a manifold of states with photon number lowered by one. The wavefunction is normalized after application of P.
To reinitialize the dynamics after the photon loss has occurred, we need both the wavefunction to propagate and a polaritonic surface to resume the nuclear trajectory integration, jFi. The wavefunction is simply a linear combination of polaritonic states jAi: where the ' symbol denotes the quantities after the jump. As a polaritonic energy surface to resume the nuclear trajectory propagation, we choose the polaritonic surface jFi that has the maximum overlap with the polaritonic wavefunction after the jump: If the quantum jump does not occur, the wavefunction is propagated with the non-Hermitian Hamiltonian [68,72,77] : For each timestep, this is accomplished by first propagating according to H pol , and then modifying Ψ pol in the following way The polaritonic wavefunction Ψ pol is then renormalized. The propagation between each attempted jump should be performed with the non-Hermitian b H eff of Equation (60), leading to un-normalized wavefunction. Anyway, in our algorithm, the jump is attempted at each time step and so the wavefunction is always normalized, one way or the other. A consequence of the photon loss is that the total energy of the system is not conserved.
In Figure 4, we replicate the dynamics performed for the lossless case ( Figure 3) in presence of a cavity lifetime τ of 65 fs. The same color scheme and notation is applied. While the relaxation dynamics is of course quicker (Figure 4a) due to the presence of an extra relaxation channel (cavity loss), the decay dynamics is better described by a biexponential function, rather than a simple exponential (Figure 4b).
The main reason is that photons can be exchanged back and forth between the cavity and the molecule, via transitions jn, pi ! j n + 1, p − 1i and vice versa, slowing down the cavity loss rate. This is especially important when p = 1, as there is no way to lose the photon from a state jn 0 , 0i with zero free photons in the cavity. In particular, this happens for the system considered here, which shows transitions back and forth from jS 1 , 0i to jS 0 , 1i. Here, the single photon remaining appears to decay with a lifetime which is 20 fs longer than the nominal decay time of the cavity. Notice that, if the single photon remaining is adsorbed by the molecule due to strong coupling, the lifetime of the system is ascribable to that of the pure electronic states.
Conversely, when the photon is free within the cavity, the lifetime of the system becomes that of the nominal cavity lifetime.
The consequence of the cavity losses becomes also evident in the energy conservation plot (Figure 4c), where the initial part of the dynamics is characterized by a quick drop of the total energy due to the photon losses with no kinetic energy compensation. As a last remark, we stress that the current implementation takes advantage of dressing the chemical quantities for the strong coupling effect. Consequently, it directly supports the interface with the TINKER package to F I G U R E 4 Cavity losses in strong coupling. Same conditions of Figure 3, with the same notation and color scheme. A cavity lifetime τ = 65 fs is considered. (a) The overall population dynamics is definitely shorter in this case, with a transient population of the jn, 1i subspace. (b) Photon number in the cavity at each time step. Remarkably, the kinetics is not simply dissipative. While p ≥ 1, the photon loss occurs at a faster rate than the cavity lifetime (circle markers fit). After only one photon remains, the loss dynamics slows down, as the only photon remaining is partially absorbed by the molecule and cannot be lost. (c) Breakdown of the energy conservation, due to cavity losses [Color figure can be viewed at wileyonlinelibrary.com] perform QM/MM simulations with electrostatic embedding, as described and applied in References [34,61,62].

| CONCLUSIONS
In the present work, we describe a scheme we have implemented to perform direct nonadiabatic molecular dynamics simulations for semiclassical molecules in strong coupling, based on classical nuclear trajectories and on multiconfigurational wavefunctions. We build polaritonic states and present the evaluation of analytical gradients for polaritonic CI energies, extending the DTSH machinery to the polaritonic systems. Among the DTSH [42,43,58] exploitable features, we count the decoherence corrections, [59] the QM/MM interface with electrostatic embedding [60,62] and the local diabatization scheme [48,63] for wavefunction propagation. Cavity losses are included in the simulations through quantum jumps, relying on the stochastic nature of Surface Hopping. We choose the test case to highlight the complex features of the potential energy surfaces arising when moving beyond the one-dimensional 2-level molecular models. The results presented for the test dynamics highlight the delicate interplay between radiative and nonradiative emissions, both impacting the relaxation dynamics of strongly coupled systems. Especially, we show that losses are competitive with usual nonadiabatic events and that the outcoming dynamics cannot be described as simply dissipative, the photon actually living longer than the nominal lifetime of the cavity. The content of this work provides both formal and conceptual tools to approach the polaritonic photochemical simulations within a semiclassical ansatz, allowing to simulate complete photochemical reactions with a trivially parallelizable technique.