Equilibrium Geometries, Adiabatic Excitation Energies and Intrinsic C=C/C–H Bond Strengths of Ethylene in Lowest Singlet Excited States Described by TDDFT

Seventeen singlet excited states of ethylene have been calculated via time-dependent density functional theory (TDDFT) with the CAM-B3LYP functional and the geometries of 11 excited states were optimized successfully. The local vibrational mode theory was employed to examine the intrinsic C=C/C–H bond strengths and their change upon excitation. The natural transition orbital (NTO) analysis was used to further analyze the C=C/C–H bond strength change in excited states versus the ground state. For the first time, three excited states including p0 y ! 3s, p0 y ! 3py and p0 y ! 3pz were identified with stronger C=C ethylene double bonds than in the ground state.


Introduction
Computational modeling of excited states is nowadays gaining more and more attention due to its increasing application in photochemistry and the rapid development of accurate and affordable computational methods [1][2][3][4][5][6]. While chemistry is majorly about the breaking/formation of chemical bonds, an important topic in modeling excited states is to analyze the chemical bonding of molecules in excited states.
The earliest theoretical tool in this regard dates back to 1983 when Dick and Freund developed a Cohen-type bond order [7] and also applied Wiberg's index [8] to excited-state organic molecules. In the past two decades, various methods have been proposed to analyze excited-state chemical bonding including (1) the time-dependent electron localization function (ELF) [9] providing the time-resolved monitoring of chemical bonding dynamics, (2) the parity function [10] based on the spin-traced one-particle density matrix, (3) the quantum theory of atoms in molecules (QTAIM) [11] with its derived interacting quantum atoms (IQA) energy partition scheme [12], (4) the density overlap region indicator (DORI) [13] and (5) the adaptive natural density partitioning (AdNDP) method [14]. However, none of these aforementioned methods can directly quantify whether a chemical bond in a molecule is strengthened or weakened when this molecule is excited to a different electronic state.
One may think it is trivial to investigate chemical bonding of excited states because molecules with higher electronic energies should have generally lower bond energies implying that an excited-state molecule will have weaker chemical bonds (with the exception of excimers [15]). However, one needs to be aware that such a simplified picture overlooks the different types of electronic excitations and the heterogeneity in the molecular charge distribution Therefore, to understand the chemical bond strength and its change caused by electronic excitation is non-trivial. A bond weakening effect in excited states has been discovered in cyclopropyl ketones [16], double-bond species [17], metal-ligand bonds [18] and ligands within a manganese-nitride complex [19]. Han and co-workers pioneered the investigation of the so-called excited state hydrogen bond between a chromophore and protic solvents. They found that these hydrogen bonds can be strengthened or weakened upon excitation [20][21][22] and characterized the strengthening/weakening effect by calculating the hydrogen bond binding energy between chromophore and solvent molecules or by the red-/blue-shift of the O-H bond vibrational frequency.
Inspired by these studies on the excited state hydrogen bonding, we fostered the idea to apply the local vibrational mode theory [23] to excited-state molecules in order to obtain the intrinsic bond strength. The local vibrational mode theory was first developed by Konkoli and Cremer in 1998 [24] and its major idea is to extract from normal vibrational modes the local counterparts in terms of individual internal coordinates. For example, the local vibrational mode of a chemical bond corresponds to the motion initiated by an infinitesimal change in that bond length followed by the relaxation of all other atoms in this molecule. This local stretching mode in terms of bond length can then be used to characterize the intrinsic strength of the bond via the corresponding local stretching force constant k a or frequency w a [23]. In the past decade, the local vibrational mode theory has been applied to investigate the chemical bonding including covalent bonds, metal-ligand bonds and non-covalent interactions of molecules and solids in the ground state [23]. However, it has not been exclusively applied to study the bonding properties in the excited states, although this is theoretically feasible, because to conduct the local mode analysis for either ground state or excited states of a molecule requires only the equilibrium geometry of a molecule and its second-order energy derivative matrix (Hessian).
In this pilot study, we applied the local vibrational mode theory to quantify the intrinsic bond strength in the ethylene molecule in its excited states. The reason why we chose ethylene as a model system is because it is a rather simple molecule with one C=C double bond and four C-H single bonds. In addition, there are plenty of theoretical studies on the excited states of this molecule [17,[25][26][27] to be used as a reference. We aimed to answer the following questions through this work: The paper is structured in the following way. First, the computational details are given. In Section 3, we present the vertical excitation energies, adiabatic excitation energies, local C=C/C-H stretching force constants of ethylene in its excited states. The natural transition orbitals (NTOs) of the electronic excitations are shown to rationalize the intrinsic bond strength of ethylene. Important findings are summarized in Section 4.

Computational Details
The geometry of ethylene was optimized using the coupled-cluster singles and doubles (CCSD) [28] method with Dunning's aug-cc-pVDZ basis set [29,30]. The optimized structure was then used in time-dependent density functional theory (TDDFT) and equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) [31] calculations with Dunning's aug-cc-pVTZ basis set [29,30] to obtain the vertical excitation energies (VEEs) of the lowest 17 singlet excited states (by specifying 20 excited state roots to find). The natural transition orbitals (NTOs) [32] analysis implemented in the Multiwfn 3.7 package [33] was employed to determine the character of the excited states where each electronic excitation is transformed into one or two pairs of compact orbital representations.
Eight out of the 17 calculated VEEs were compared with the experimentally measured excitation energies, then the root mean square (rms) and the mean signed average (MSA) deviations were calculated as following.
Based on the performance of different methods in reproducing the experimental excitation energies, the density functional CAM-B3LYP [34] was chosen for the (TD)DFT calculations to obtain optimized geometries of the singlet excited states as well as of the ground state. Noteworthy is that no symmetry constraint was imposed during the geometry optimizations. The vibrational analysis based on analytic second-order energy derivatives of TDDFT [35,36] was then carried out to confirm that the optimized geometries were indeed local minima on the potential energy surfaces (PES). In addition, the NTO analysis was applied to each optimized excited-state geometry to reassure that it is staying on the desired excited-state PES.
All calculations of VEEs, optimized geometries, adiabatic excitation energies (AEEs), Hessian matrices and normal vibrational analyses were carried out with the Gaussian 16 package [37].
The intrinsic C=C and C-H bond strength quantified by their local stretching force constants k a was calculated with the standalone LModeA package [23,38] based on the Hessian matrices.

Results and Discussion
Before characterizing the bonding properties of ethylene in its excited states, we benchmarked several levels of theory in terms of VEEs against 8 experimental excitation energies [25,26] in order to find the optimal method for theoretical calculations. As most low-lying excited states of ethylene are Rydberg states, we tested several frequently used density functionals with long-range corrections including LC-wHPBE [39], whPBE0 [40], wB97X-D [41], CAMh-B3LYP [40], CAM-B3LYP [34] and M06-2X [42] for TDDFT calculations. In addition, the EOM-CCSD theory which is believed to outperform TDDFT was considered in our benchmark.
As shown in Figure 1, most density functionals have comparable rms deviations below 0.4 eV while LC-wHPBE has both large rms and MSA deviations. Surprisingly, the EOM-CCSD has larger errors (rms = 0.57, MSA = 0.47 eV) than TDDFT in calculating VEEs for ethylene. Among the tested density functionals, three of them have rms deviations equal or less than 0.3 eV including whPBE0 (0.31 eV), CAM-B3LYP (0.31 eV) and wB97X-D (0.27 eV). The latter two functionals have MSA deviations less than 0.1 eV. Based on this result, we chose the CAM-B3LYP functional for the remaining calculations although wB97X-D would have been the other ideal option.
We calculated 17 lowest singlet excited states of ethylene at the CAM-B3LYP/aug-cc-pVTZ level with their irreducible representations, assignment and vertical excitation energies collected in Table 1. In order to quantify the intrinsic bond strength of ethylene in its excited states with the local vibrational mode theory, it is required to first optimize the ethylene geometry accordingly. The equilibrium geometries of 11 among the 17 excited states were successfully determined while the remaining 6 states could not be optimized due to two possible reasons: (1) the target excited-state PES has no local minima (with vanishing gradients) as it may form a crossing seam or canonical intersection with another electronic state; (2) the geometry optimizer fails to stay on the desired excited-state PES and crosses to the PES of another undesired electronic state [43]. The C=C/C-H bond lengths, point group and adiabatic excitation energies of the optimized geometries of ethylene in different excited states are complied in Tables 1 and 2, while the Cartesian coordinates are available in the Supplementary Materials. Due to the pseudo-Jahn-Teller effect [44], there exists some connection between the symmetry of the excited state and the point group of its equilibrium geometry. Among the 11 states investigated in this work, the B 3u states always leads to D 2 -symmetry non-planar structures. The A g states leads to either D 2 or D 2h structures. The C 2h -symmetry structures can arise from the B 1u and B 1g states. The B 3g state can lead to very unique C 2v structure where two carbon atoms are no longer equivalent.
By employing the local vibrational mode theory, the local stretching force constants k a were determined for C=C and C-H bonds of these optimized geometries of ethylene and the results are collected in Table 2. Figures 2-4 show the distribution of the local stretching force constants as intrinsic bond strength measure for C=C and C-H bonds of ethylene in its excited-state equilibrium geometries. For both types of bonds, the bond strength is overall weakened as expected for the 11 excited states investigated in this work. However, three and two excited states exhibit stronger C=C and C-H bonds than ground state (S 0 ) respectively.
By correlating the local stretching force constant k a and bond length r for C=C bonds as shown in Figure 2, we found a strong linear correlation for 9 out of 11 excited states with the coefficient of determination (R 2 ) of 0.966. The generalized Badger's rule [45] implying that shorter bonds are stronger is valid. However, the ethane C-C single bond is out of trend in this k a -r relationship, which is in line with the fact that single and double CC bonds have different bonding nature.
In order to translate the local stretching force constant k a for C=C bonds of ethylene into a more chemically intuitive quantity, we propose to use the bond strength order (BSO) [46] in reminiscence of Pauling's bond order (BO) [47]. As shown in Figure 3, we first used the local stretching force constants of ethane C-C, ethylene C=C and ethyne C⌘C bonds to define the reference BSO values of 1.0, 2.0 and 3.0 respectively. Then a cubic function between BSO and k a was used to fit these three points and the origin. We found the C=C bonds of most excited states have a BSO of 1. In addition, Figure 4 shows the relationship between k a and r for C-H bonds of ethylene and we observed a weaker correlation between these two quantities (R 2 = 0.902). The C=C bond axis is defined as the z direction while the x direction is defined to be perpendicular to the geometrical plane of ethylene (see Figure 2). b Vertical excitation energies in eV. Experimental excitation energy is recorded in brackets. c Adiabatic excitation energies in eV. d Oscillator strength.  CC bond CAM-B3LYP/aug-cc-pVTZ S1 S16 S9 S14 S15 S4 S11 S8 S6 S7  In order to better understand the variation in the intrinsic C=C/C-H bond strength of ethylene in excited states, we calculated the natural transition orbitals (NTOs) [32] of all 17 singlet excited states investigated in this work. The underlying idea of the NTO analysis is a unitary transformation of the occupied and virtual orbitals with no change towards the transition density. By correlating transformed occupied/virtual orbitals (i.e., NTOs) into pairs, one can easily obtain the straightforward physical picture of the excitation of "particle" into the empty "hole". Figures 5-7 illustrate the NTO diagrams of the 11 excited states whose equilibrium geometries and bonding properties are investigated in this work. The NTO diagrams for the remaining 6 singlet excited states can be found in the Supplementary Materials.
Most of the 11 singlet excited states of ethylene investigated in this work have the electronic excitation from the p bonding orbital, except S 5 , S 9 , S 14 and S 15 . The "particle" NTO of these 4 excited states is p 0 y orbital which consists of the 2p y atomic orbitals of C and the s orbitals of H. One interesting observation is that S 9 , S 14 and S 15 have C=C bonds stronger than the ground state (S 0 ) ethylene, while the S 5 state has weaker C=C bond. The p 0 y orbital has a similar shape as the p ⇤ anti-bonding orbital except its orientation and therefore this orbital has anti-bonding character. When the electron is excited from p 0 y orbital, the depletion in this anti-bonding orbital can lead to stronger C=C bonds for S 9 , S 14 and S 15 . However, in the case of the S 5 state the electron is excited to the real p ⇤ anti-bonding orbital from the p 0 y anti-bonding orbital, thus the C=C bond is weakened compared to the ground state. Noteworthy is that the C=C bond of S 15 state is slightly stronger than for the S 9 and S 14 states. Unlike S 9 and S 14 , the "hole" NTO of S 15 has some s ⇤ anti-bonding character due to the presence of a node between two carbon atoms, it is not quite clear why such excitation could lead to the strongest C=C bond at this moment.  Table 1. The left and right NTO diagrams shown as mesh isosurfaces depict the excited particle (occupied) and empty hole (unoccupied) respectively. The red and blue colors are chosen arbitrarily to indicate two different phases. Below each isosurface is the corresponding isovalue in atomic units (a.u.). Each NTO pair accounts for more than 90% contribution to the target electronic excitation unless otherwise specified.
For those excited states whose particle NTO is the p bonding orbital, the C=C bonds are weakened to different extent. The S 6 and S 7 are relatively stronger than others. In the case of S 6 , the electron is excited from the 2p x p bonding orbital into 3p x p bonding orbital, however, the C=C bond is weakened.
For the C-H bonds of ethylene in its excited states, the S 16 state has almost the same C-H strength as the ground state and only S 7 has a marginally stronger C-H bond. The NTO diagrams for the S 7 state show that the electron is excited from the p bonding orbital into the 4d yz orbital which contains 4 s(C-H) bonding orbitals and this explains why the C-H bonds in S 7 state are strengthened.
When the electron is excited from the p 0 y orbital in the cases of S 5 , S 9 , S 14 and S 15 , the C-H bonds are unanimously weakened because their particle NTO consists of 4 s(C-H) bonding orbitals and electron depletion from this orbital can lead to weaker C-H bonds. The C-H bonds in S 15 state are greatly weakened and this can be related to the existence of the anti-bonding nodes of the 3p z orbital as shown in Figure 6. Interesting is that the S 4 state has also quite weak C-H bonds (see Figure 4). Its hole NTO is the same as the hole NTO of S 15 where 4 anti-bonding nodes exist for the C-H bonds. Due to symmetry, not all C-H bonds are equivalent in ethylene for some excited states. S 5 , S 9 and S 14 have two sets of C-H bonds and the difference in the local stretching force constant between two sets is 0.86, 2.34 and 1.55 mdyn/Å, respectively. In the case of S 9 where the electron is excited from p 0 y orbital to 3s orbital leading to unbalanced charge distribution of the ethylene molecule, the C-H bonds at one side have their local stretching force constant being 1.6 mdyn/Å as the lowest value among all C-H bonds.
The S 1 and S 11 states have very similar NTO character in their excitations as the former has the 3s orbital as the hole NTO and the latter has 4s orbital as the hole NTO while their particle NTO is the same. Both states have weakened C-H bond strength compared to the ground states. However, the S 11 state has its C-H bonds weakened to a much larger extent and this can be related to its more pronounced C-H anti-bonding character in the hole orbital as shown in Figure 5.

Conclusions
In this work, we calculated the 17 lowest singlet excited states of ethylene in gas phase with TDDFT and analyzed the C=C/C-H bonding properties by combining the local vibrational mode theory and natural bond orbitals (NTOs). This led to the following important findings:

1.
Among the several commonly used density functionals for describing the Rydberg states, wB97X-D, CAM-B3LYP and whPBE0 give relatively low rms deviations (⇠0.3 eV) against the experimental excitation energies for ethylene.

2.
The equilibrium geometries for 11 out of 17 excited states of ethylene were successfully modeled at the CAM-B3LYP/aug-cc-pVTZ level. Six of them are twisted in contrast to the planar ground-state structure.

3.
For the first time, the local vibrational mode theory has been applied to molecules in excited states to characterize the intrinsic bond strength which can be directly compared with the ground-state counterpart. In this study on ethylene, the majority of the 11 excited states with optimized geometries exhibit weakened C=C/C-H bonds as characterized by the local stretching force constant k a , which is consistent with the assumption that a molecule in its excited states with higher electronic energy has in general smaller bond (dissociation) energies. However, the results from the local mode analysis show that 3 excited states have stronger C=C bonds with BSO values of 2.4, which almost reaches the midpoint between ethylene C=C (BSO = 2.0) and ethyne C⌘C (BSO = 3.0). Meanwhile, 2 excited states show only marginally stronger C-H bond strength the ground state. Overall, the local stretching force constants for these 11 excited states of ethylene show a larger range of variation for the C=C bond compared to the ground state (variation of 48 % to +27 %) than for C-H bonds (variation of 72 % to +1 %). This leads to the important conclusion that the C=C bond of ethylene is more tunable by electronic excitation, while this does not hold for the C-H bond.

4.
The NTO analysis plays an important role in helping interpret the result of local stretching force constants for C=C/C-H bonds of ethylenes in excited states.
• When the particle NTO is the p bonding orbital, the C=C bond is always weakened in terms of local stretching force constant. • When the particle NTO is the p 0 y orbital arising from the s-backbone of ethylene, the C=C bond is in general strengthened except when the hole NTO is the p ⇤ anti-bonding orbital. Surprisingly, when the hole NTO orbital is 3p z the C=C bond is the strongest and this hole NTO has s ⇤ (C=C) anti-bonding character. • The excited state leading to slightly increased C-H bond strength has its hole NTO as 4d yz which consists of 4 C-H bonding orbitals. • When the hole NTO has significant C-H anti-bonding character, the C-H bonds are greatly weakened.

5.
Computational modeling of molecular excited states is generally more challenging than for the ground state. Therefore, we used ethylene as a simple theoretical model system to provide the proof of concept that fine-tuning chemical bond strength via electronic excitation is feasible, although the calculations performed in this work can be further improved by more advanced wavefunction-based correlation methods other than TDDFT [5,27]. It is also possible to find the equilibrium geometries of more excited states with dedicated algorithms [43], which we will test in the future.
We hope our work will spur a similar interest in exploring chemical bonding properties of excited states among both experimental and theoretical chemists.
Funding: This work was financially supported by the National Science Foundation (Grant CHE 1464906).

Acknowledgments:
We thank SMU for providing generous computational resources. Y.T. thanks Tian Lu for helpful discussions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: