The 3s Rydberg state as a doorway state in the ultrafast dynamics of 1,1-difluoroethylene

The deactivation dynamics of 1,1-difluoroethylene after light excitation is studied within the surface hopping formalism in the presence of 3s and 3p Rydberg states using multi-state second order perturbation theory (MS-CASPT2). Due to the proximity of the Rydberg p -3s state with the pp * state, the states are mixed favoring ultrafast exchange of population via a conical intersection that closely resembles the equilibrium structure. After excitation, it is found that the p -3s state acts as a doorway state, trapping the population and delaying internal conversion to the pp * state, from which deactivation to the closed-shell ground state takes place. Besides the conical intersection between the p -3s and pp * states, five additional conical intersections between the pp * state and the ground state are found, indicating that after the system is excited, it stretches the C Q C bond before it twists and pyramidalizes at any of the carbon atoms, in the spirit of a hula-twist mechanism.


Introduction
Ethylene is a prototype for conjugated p systems in a large number of interesting systems, including polymers, 1 organic dye sensitized solar cells, 2,3 non linear optical chromophores, 4 or organic light-emitting diodes using low molecular weight compounds with high fluorescence quantum yields. 5 Additionally, its light-induced cis-trans isomerization plays an important role in biological processes, such as vision, [6][7][8] and it can be exploited to operate molecular rotors such as those pioneered by Feringa based on overcrowded alkenes. 9 Given its many applications, it has been a recurrent topic of interest for theory [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28] as well as for experiment. [29][30][31][32][33][34][35][36][37][38] As a result, the electronic states of ethylene are well characterized, with an intense pp* or valence band dominating the absorption spectrum and surrounded by Rydberg series where the p-3s lies below the absorption maximum, well separated from the other p-3p states, and a mostly dark double excited or zwitterionic state. 32,33 From the electronic ground state, an excitation to the bright pp* state triggers an ultrafast deactivation in less than 50 fs, 33,34,37,38 which has been interpreted as internal conversion to the electronic ground state through several conical intersections (CI), including twisted, twisted-pyramidalized, hydrogen-migration and ethylidenelike structures. 14,20,22,23,[39][40][41] An interesting question is whether the Rydberg states play an active role in the deactivation of ethylene. Theoretically, the 3s Rydberg orbital has been included in a number of stationary 42,43 as well as dynamical calculations. 40,41 The work by Mori et al. 40 using multiple spawning dynamics showed that although the p-3s state lies below the bright pp* state, little population arrives in this state during the first 10 fs and the excited state lifetime is only slightly longer (90 fs) than when the p-3s state is excluded (60 fs), thus concluding that the Rydberg state is merely a spectator. Another dynamical study by Sellner et al., 41 including the 3s orbital, also found the contribution of the p-3s state to be minor, although the details of the deactivation dynamics depend slightly on the method employed for the calculation of the excited states. Regardless of the methods employed, both studies agree that the inclusion of Rydberg states has a modest effect in ethylene and its photodynamics can be theoretically described including valence and zwitterionic states exclusively.
Adding halogen atoms can substantially alter the photochemical properties and the deactivation mechanism of ethylene. Experiments indicate that 1,1-difluoroethylene and fluoroethylene show a very similar absorption spectra, 29,44,45 only slightly blueshifted with respect to ethylene. However, while in ethylene, internal conversion from the pp* to the ground state requires a CQC stretching followed by a twisting around the double bond and a pyramidalization at one of the methylene units, 46 in fluoroethylene, the torsion over the double bond directly drives to a CI between the pp* state and the ground state. 39 It is also curious that depending on which of the two H atoms are substituted by F in ethylene, CIs similar to those predicted in either ethylene or fluoroethylene can be located in 1,1-and 1,2difluoroethylene, 47 respectively. Also the CI corresponding to H-migration present in ethylene and 1,2-difluoroethylene was never found in 1,1-difluoroethylene. 47 The aim of the present paper is to investigate the deactivation dynamics of 1,1-difluoroethylene (1,1-DFE hereafter) after light irradiation using non-adiabatic surface hopping simulations. Up to our knowledge, no excited state dynamics study is available in ethylene-fluoro-derivatives, beyond our previous one-dimension wave packet simulations in 1,1-DFE. 48,49 Our intention is twofold. First, we aim at comparing the internal conversion mechanism of deactivation of 1,1-DFE against bare ethylene. Second, we want to investigate whether Rydberg orbitals play a role in the nonadiabatic dynamics of 1,1-DFE. Generally, we are interested to learn which main coordinates play a main role in the deactivation of 1,1-DFE. To this purpose, multiconfigurational calculations including the 3s and the three 3p orbitals have been carried out and used to perform non-adiabatic simulations at the multistate second order perturbation (MS-CASPT2) level of theory.

Computational details
The non-adiabatic dynamics simulations on 1,1-DFE were carried out using the Surface Hopping including ARbitrary Couplings (SHARC) approach. [50][51][52] The trajectories were propagated during 200 fs with a nuclear time step of 0.5 fs, where the integration of the nuclear motion was done with the velocity-Verlet algorithm. 53,54 Decoherence corrections were taken into account following the energy-based method of Granucci and Persico 55 with a parameter of a = 0.1 Hartree. As a hopping scheme, the ''Standard SHARC surface hopping probabilities'' scheme was used. 51 Direct dynamical simulations including Rydberg and valence states are challenging as both type of states need to be described on the same footing with similar accuracy. The description of Rydberg states requires extensive diffuse functions and valence states profit from the inclusion of large amounts of dynamical correlation. 56 A balanced description can be achieved using the multi-state complete active space second-order perturbation theory 57,58 (MS-CASPT2) method and at least the aug-cc-pVDZ basis set. 59 No IPEA shift 60,61 and an imaginary level shift of 0.3 a.u. were employed. Accordingly, the electronic energies and gradients were calculated on-the-fly for each nuclear time step using MS-CASPT2, as implemented in MOLCAS8.0. 62 To avoid the expensive calculation of non-adiabatic couplings at every time step, the propagation was done using the local diabatization scheme, that allows the use of (the full) wave function overlaps instead, as described in ref. 63.
Two different active spaces (see corresponding orbitals in Fig. 1) were employed in the underlying Complete Active Space Self-Consistent Field 64 (CASSCF) calculation to investigate the effect of the Rydberg states on the deactivation dynamics. The smallest active space includes only the two electrons and the two active p and p* orbitals; accordingly, three roots were stateaveraged (SA) with equal weights, resulting in the neutral (N) closed-shell ground state p 2 , the valence (V) pp* state and the dark zwitterionic (Z) p* 2 state -this calculation will be refereed as SA(3)-CAS(2,2)PT2. The larger active space includes the Rydberg 3s, 3p x , 3p y , 3p z orbitals besides the p and p* pair and it requires 6 or 11 roots to include all the possible single or single and double excitations from the p to the p* and Rydberg orbitals, respectively; these calculations are labelled SA(6)-or SA(11)-CAS(2,6)PT2.
In order to generate initial conditions for the SHARC dynamical simulations, the ground state equilibrium geometry was optimized at both the MP2/aug-cc-pVDZ and SA(11)-CAS(2,6)PT2/aug-cc-pVDZ levels of theory. As it can be seen from Table 1, the geometrical parameters are very similar with both methods. Thus, the MP2 geometry was further employed to calculate harmonic frequencies that were used to generate a Wigner distribution of 1000 uncorrelated velocities and geometries. From each geometry, vertical excitations were calculated at MS-CASPT2 level of theory with both SA(3)-CAS(2,2) and SA(11)-CAS(2,6) active spaces to simulate an absorption spectrum from the oscillator strengths, as explained in ref. 65. The oscillator strengths served as a condition to select initial trajectories for the dynamics 66 which were excited assuming an infinitely short d-pulse. Since the dynamical results will be given in the adiabatic picture (i.e. governed by states strictly ordered by their energy that never cross), a transformation to the diabatic or spectroscopic picture 51 was carried out a posteriori, directly analyzing the configurational interaction vectors resulting from subsequent MS-CASPT2 single point calculations carried out for the geometries taken every 2 fs of every trajectory. The state character was assigned to the largest reference weight. Note here that the diabatization performed does not consider the full configuration interaction vector, only the largest contribution to it. This however, allows us to follow in an approximate way the electronic character of the states, even though the unitary transformation diabatic-adiabatic is lost.    To cluster the hopping geometries to a minimum energy CI (MECI), the root-mean-square deviation (RMSD) was minimized using quaternions. 67 The structures were assigned according to their minimum value of the RMSD with respect to every MECI. The optimization of MECIs was done using MOLCAS in combination with the the external optimizer of ORCA, 68 as described in ref. 69. The evolution of the relevant coordinates during the dynamics was performed using a Gaussian convolution of all trajectories. Trajectories were stopped after hopping to the lowest energy state (S 0 ) and remaining there during ten femtoseconds.

Vertical excitations and absorption spectra
The vertical excitations of 1,1-DFE computed at SA(3)-CAS(2,2)PT2, SA(6)-and SA(11)-CAS(2,6)PT2 levels of theory are collected in Table 2. The Rydberg states originate from excitations from the p orbital to every one of the Rydberg orbitals. As it can be observed, the inclusion of the Rydberg orbitals leads to closer agreement of the energy of the valence state V with the experimental value, 44,45 but only when all single and double excitations are included in the calculation. If only single excitations are allowed (6 roots, SA (6)), the V state is mixed with the p-3p z state, deteriorating considerably the energy of the pp* state. The most extensive SA(11)-CAS(2,6)PT2 calculation predicts that the spectrum of 1,1-DFE is composed of a main band due to the excitation to the bright V state surrounded by Rydberg series. The p-3s state lies slightly below the V state (E0.5 eV) at the equilibrium geometry whilst the three p-3p states are more than 1 eV above.
The absorption spectra calculated from a Wigner distribution including 1000 geometries at the SA(3)-CAS(2,2)PT2 and SA(11)-CAS(2,6)PT2 levels of theory are superimposed to the experimental spectrum in Fig. 2. Again, the spectrum obtained with the smallest active space is slightly red-shifted with respect to that calculated with the largest active space and the experimental one from Bélanger and coworkers. 44 The overall width of the spectrum is well described with both methods, but it is noticeable that the fine structure of the main experimental band is very well reproduced with the (2,6) active space. The band at energies larger than 8 eV corresponds to excitations to 3d and higher order Rydberg orbitals that are not included in any of our calculations.

Non-adiabatic SHARC dynamics
SHARC simulations were initiated from the initial geometries and velocities contained in an excitation window of 0.5 eV centered at the absorption maxima (color boxes in Fig. 2). For the SA(3)-CAS(2,2)PT2 calculation, this provides 152 geometries excited to the S 1 state, which here is the pp* state. Likewise, for the SA(11)-CAS(2,6)PT2 simulations, the same energetic window corresponds to 121 trajectories, from which 108 (89%) are excited to the S 2 , 11 (9%) to the S 1 and 1 to the S 3 as well as the S 4 states. Note that in this case, the states are labelled by increasing energy from S 0 to S 10 . Initially, besides the brightest pp* state, also the p-3s state is excited. Note however, that during the dynamics the character of the states can change so that a one-to-one correlation with the energy ordering and character of the states at the Franck-Condon (FC) point is not straightforward.
The evolution of the adiabatic populations for both sets of simulations is shown in Fig. 3a and b, together with exponential fits (thinner lines). For the small active space (Fig. 3a), after 20 fs the population of the S 1 state decays very fast to the ground state, which then recovers to 93% at the end of the simulation. An exponential fit from the time interval between 20 and 200 fs returns a time scale of 30 fs for the decay from the S 1 to the S 0 state. Table 2 Vertical excitations (DE, in eV) of 1,1-DFE at SA(3)-CAS(2,2)PT2, SA(6)-and SA(11)-CAS(2,6)PT2 levels of theory compared with experimental 44,45 data, with associated oscillator strengths f and reference weights. The electronic states are labelled following the Mulliken notation for ethylene, 71   Interestingly, the dynamics including Rydberg states (Fig. 3b) shows during the first 20 fs that population is mainly transferred from the S 2 to the S 1 state, with some negligible amount of population in the S 3 and S 4 states. The ground state recovery starts after 26 fs, and at first glance one can see that the slope of the associated curve is smaller than when Rydberg states are excluded (compare with Fig. 3a). After 200 fs, 72% of the trajectories have decayed to the electronic ground state S 0 . In this case, the time constants for the exponential decay of the population of the respective states were fitted separately in two time intervals: from 0 to 26 fs and from 26.5 to 200 fs. Within the first 26 fs, population decays from the S 2 to the S 1 with a time constant of 36 fs. After 26 fs the S 1 state decays with a time constant of 50 fs. The time constant for the decay from the S 1 to the S 0 was only fitted in the interval between 26 and 200 fs and was found to be 32 fs. This time constant appears to be the same as when Rydberg states are excluded (30 fs); however, one should notice that when Rydberg states are present, the dynamics starts with population mainly in the S 2 state (which for some geometries is the V state and for others the Rydberg p-3s) and this population needs to be transferred first to the S 1 before decaying to the S 0 ; therefore, the overall effective decay is slower.
In order to compare the kinetic models obtained with both active spaces in the same way, another fit was then performed considering the sum of the population of all the electronic excited states versus that of S 0 . The evolution of these populations as well as the corresponding exponential fit is shown in Fig. 3c. In this case, a lifetime of 60 fs was obtained; this is the time that all the excited state population takes to decay to the electronic ground state when the Rydberg states are included. Thus, the population transfer to the ground state when Rydberg states are excluded occurs approximately twice as fast than when Rydberg orbitals were included.
For both sets of trajectories, the amount of hops between the two lowest energy states in time intervals of 5 fs was analyzed. Fig. 4 displays an histogram with the frequency of hops between S 1 and S 0 per time interval over all the trajectories. In the (2,2) active space the hopping frequency is more localized, with the majority of the hops occurring around 26 fs. The (2,6) active space produces a broader and generally slower behavior, mainly due to the transfer of population to the Rydberg states that at later times need to find their way back to the valence pp* state before deactivating to the S 0 . The correspondence of the geometries where the above mentioned hops occurred and its relationship to the optimized MECIs will be discussed in the next section.
As it is clear that Rydberg states have an influence in the deactivation mechanism of 1,1-DFE, henceforth, only the ensemble of trajectories from the SA(11)-CAS(2,6)PT2 will be further analyzed and discussed. Fig. 5 shows the variation of the coordinates expected to play an important role before deactivation to the S 0 takes place. Taking into account the initial Wigner distribution as a reference, during the first 10 fs of the dynamics the CQC distance starts elongating while the other coordinates remain localized around the equilibrium geometry. This CQC stretching oscillates

Minimum energy conical intersections
In order to understand better the internal conversion mechanism of 1,1-DFE to the S 0 , all the hopping geometries encountered during the SA(11)-CAS(2,6)PT2 dynamics were optimized as statecrossings. Starting from the S 1 /S 0 hopping events, 78 geometries were optimized, which collapsed into the five MECIs, CI 1 -CI 5 , shown in Fig. 6 (see ESI † for Cartesian coordinates). All structures are substantially different from the equilibrium (FC) geometry, see Table 3, but have a similar torsional angle around 60-751. According to their relative energies, after the pp* state is excited, all MECIs are energetically accessible. The MECIs with the highest energy are CI 1 and CI 3 , characterized by a very elongated CQC bond and a very elongated C-F bond, respectively. While the CI 1 and CI 2 both show the large CF 2 pyramidalization angle of about 801, the latter is found at lower energy. In contrast, CI 4 and CI 5 have the lowest energy. The former is pyramidalized at the CH 2 carbon whilst the latter is the most twisted one, showing also a small pyramidalization at the CH 2 carbon of about 201. In terms of the one-bond-flip (OBF) or hula-twist (HT) mechanisms, 39 the CI 1-4 allow for the HT deactivation pathway while the CI 5 operates for the OBF mechanism. In the CI 1-2 , it is the CF 2 which hula-hoops and the CH 2 that twists; in CI 3-4 the CF 2 twists and the CH 2 hula-hoops.
Starting from the S 2 /S 1 hopping events, 26 geometries were optimized, all of which collapsed to the same geometry, the CI a . This structure enables population transfer between the p-3s Rydberg and the pp* state. The geometry is very similar to that of the equilibrium, only with a slightly larger CQC bond distance; accordingly, both structures are energetically very close, allowing for a very fast radiationless transfer between the pp* state and the p-3s Rydberg state.
As all the MECIs are energetically accessible from the V state, their relative importance can only be assessed from the dynamical simulations. By minimizing the RMSD between the structures, it is possible to classify the hopping geometries to find that the majority of the trajectories went to the S 0 through CI 2 (44.9%) followed by CI 5 and CI 1 (21.8% and 20.5%). The less active MECIs are CI 4 and CI 3 with 9.0% and 3.8% of the trajectories, respectively. Considering their geometries, one can conclude that deactivation operates after the system elongates its CQC bond, twists and pyramidalizes the CF 2 carbon atom, the latter motion only being present in two out of the three predominant MECIs. This means that 65% of the CIs correspond to HT of the CF 2 , 13% to HT of the CH 2 and 22% can be assigned to the OBF mechanism.
In comparison to the ethylene, 40,41 no ethylidene-like CI was found here. Interestingly, in the dynamics of Sellner et al. 41 on ethylene, the hops to the S 0 required pyramidalized geometries, while in 1,1-DFE, 22% of the deactivation is possible without pyramidalization, i.e. following a OBF mechanism. Our optimized CI a between the pp* and p-3s is energetically in agreement with the one found by Mori et al. 40 for ethylene, close to the FC region. However, in 1,1-DFE, the CI a is closer to the equilibrium geometry, triggering a faster population transfer to the p-3s state than in ethylene. In contrast to the study of fluoroethylene done by Barbatti et al., 39 in 1,1-DFE neither H-migration CI or ethylidene-like were found within the dynamics.

Analysis of the dynamics in the diabatic representation
As during the dynamics the character of the states mix, it is convenient to transform the adiabatic populations into diabatic ones. Diabatic states can be described as those where the   Table 3. Table 3 Geometrical parameters of the critical points optimized for 1,1-DFE (see also Fig. 6). Distances (r) in Å and pyramidalization (pyrCF 2 and pyrCH 2 ) and torsion angles in degrees electronic character does not change with respect to the nuclear coordinates. A transformation from adiabatic to diabatic states is often a tedious task, as no general black-box recipe exists.
Here, the diabatization is based on the configuration interaction vectors of the MS-CASPT2 method taken every 2 fs for each trajectory; the diabatic states are defined according to the state character of the largest squared coefficient of this vector. The resulting diabatic populations for the N, V, Z and p-3s states are depicted in Fig. 3d.
Comparing to the adiabatic picture (Fig. 3a), only 60% of the trajectories deactivated to the ground state after 200 fs. This decrease evidences a mixing of the pp* state and the closedshell ground state during the dynamics, as expected for torsional angles close to 901, where both orbitals p and p* are degenerated and thus indistinguishable. At the beginning of the dynamics, the population is well distributed between the valence state and the energetically close 3s-Rydberg in a ratio of 3 : 1. Note that the initial distribution between both states is different from that obtained for the adiabatic S 2 and S 1 states (9 : 1), as state inversions take place for some geometries of the initial Wigner distribution. None of p-3p Rydberg states is populated during the propagation time, meaning that their contributions to the configuration interaction state vector are always minor. The population of the double excited p* 2 state (Z) is always low, but not negligible, despite this state lies ca 7 eV above the V state at the FC region. This is the result of the ultrafast stabilization of the Z state, which after about 10 fs is already the adiabatic S 2 state. In passing we note that this emphasizes the importance to include double excitations in the electronic wavefunction during the dynamics. As in Fig. 3a, the diabatic populations show a ground state recovery after a short time delay, in this case 24 fs. Thus, in this case the populations are fitted within two time intervals, one from 0 to 24 fs, where the V state transfers population to the p-3s and another from 24.5 fs onwards, in which the deactivation from V to N takes place. A scheme showing the kinetic constants obtained in Fig. 7. Note that the fact that the final populations on the ground state are different in the adiabatic and diabatic representations clearly evidences the importance of the diabatization, as it emphasizes the strong mixing between pp* and pp closed shell contributions.
From the FC structure, 75% of the population is vertically excited to the V state and 25% to the Rydberg p-3s state. The population transferred to the V state leaks non-adiabatically to the p-3s state with a time constant of 9 fs. After this initial ultrafast exchange to the p-3s state, the population is transferred back to the V state with a time constant of 51 fs, from which it deactivates to the close-shell N state, in another 51 fs. The population directly promoted to the p-3s state follows the same route, undergoing via the V state.
The initial ultrafast population inversion between the valence and the Rydberg state is in good agreement with the fact that the CI a structure found for this process is geometrically very close to the FC geometry. After the V state is populated, 1,1-DFE twists around the CQC bond, stretches the CQC bond and depending, on which of the CI is used for deactivation, pyramidalizes at any of the carbon atoms to undergo a radiationless transition to the N ground state.
Notably, the dynamics of 1,1-DFE is rather different from that reported by Sellner et al. 41 for ethylene, where the Rydberg states, slightly populated at the beginning of the dynamics, decay very fast to the V state, becoming depopulated at already 10 fs. One should nevertheless note that in that study 41 only three adiabatic states (S 0 , S 1 and S 2 ) were included in the dynamics and this might not be sufficient, in the light of the present work. In contrast, the population dynamics obtained here is similar to that obtained by Mori et al. 40 for bare ethylene, where ca. 30% of the population goes initially from the valence state to the Rydberg 3s state, before deactivating to the ground state. However, in 1,1-DFE, almost all the population makes the detour via the Rydberg 3s; this makes the ground state recovery slower than in ethylene, as the system needs to find its way back to the V state to continue its deactivation pathway to the ground state.

Conclusions
In this paper, the non-adiabatic dynamics of 1,1-DFE including the 3s and 3p Rydberg manifold was investigated using surfacehopping simulations at the MS-CASPT2 level of theory. It is found that the inclusion of Rydberg states, (i) considerably improves the description of its absorption spectrum, and (ii) has an influence on the deactivation time scale, with respect to calculations that exclude Rydberg states. In particular, the p-3s state acts as a doorway state for the ultrafast deactivation of 1,1-DFE as it traps all the population before it gets transferred to the V state and deactivates through internal conversion to the closed-shell N electronic ground state.
The p-3s and the V states are excited in a ratio 1 : 3, as they lie rather close in energy and there is intensity borrowing from the valence state. A conical intersection very similar to the equilibrium geometry allows for ultrafast exchange between these states within less than 10 fs. Once the population is finally in the V state, five structures have been located which enable radiationless decay from the V to the N state. From them, the most important three channels for internal conversion involve CQC stretching, torsion, and pyramidalization on the CF 2 carbon atom, following a HT mechanism. One of the other two less important pathways involves slight pyramidalization on the CH 2 carbon atom and CQC stretching, while the least important of these three conical intersections is strongly twisted. The black numbers are the kinetic constants in fs. Deactivation to the closed-shell N ground state occurs mainly via three conical intersections (CI) between the V and the N state. The Rydberg p-3s state slows down the deactivation as all the population excited to the V state makes a detour via the CI a to the p-3s state.
A sequential kinetic fit indicates that the ultrafast deactivation of 1,1-DFE in the presence of the 3s Rydberg state is delayed due to the Rydberg detour experienced by the population excited to the V state. In either case, it takes about 50 fs to empty the p-3s state to the V state and another 50 fs to deactivate from the V to the N state. For comparison, when the Rydberg states are neglected, direct excitation to the V state deactivates to the N state within 30 fs.

Conflicts of interest
The authors declare no conflicts of interest.