Simulation of vibrationally resolved absorption spectra of neutral and cationic polyaromatic hydrocarbons

The identification of the carriers of the absorption features associated with the diffuse interstellar bands (DIBs) is a long-standing problem in astronomical spectroscopy. Computational simulations can contribute to the assignment of the carriers of DIBs since variations in molecular structure and charge state can be studied more readily than through experimental measurements. Polyaromatic hydrocarbons have been proposed as potential carriers of these bands, and it is shown that simulations based upon density functional theory and time-dependent density functional theory calculations can describe the vibrational structure observed in experiment for neutral and cationic naphthalene and pyrene. The vibrational structure arises from a small number of vibrational modes involving in-plane atomic motions, and the Franck–Condon–Herzberg–Teller approximation improves the predicted spectra in comparison with the Franck–Condon approximation. The study also highlights the challenges for the calculations to enable the assignment in the absence of experimental data, namely prediction of the energy separation between the different electronic states to a sufficient level of accuracy and performing vibrational analysis for higher-lying electronic states.


Introduction
The diffuse interstellar bands (DIBs) are absorption features associated with the interstellar medium [1][2][3]. DIBs were first reported in 1922 [1], and the identification of the carriers of the DIBs has proved to be elusive, with their assignment being referred to as the longest standing challenge in astronomical spectroscopy [4]. It has long been believed that these bands have a molecular origin, and many potential candidates have been proposed [5][6][7]. It has been postulated that DIBs are associated with neutral and charged large carbon containing molecules that are present within the interstellar dust [7,8]. One example is C + 60 which has been confirmed as the source of two lines in the DIBs [9,10]. Polyaromatic hydrocarbons (PAHs) have also been proposed as potential carriers of DIBs [11,12], and there has been considerable effort toward confirming the association of PAHs with DIBs [13][14][15].
The unambiguous assignment of lines in the DIBs is extremely challenging and a detailed characterisation of the relevant spectra of potential carriers is a vital component to achieve a definitive assignment. Computational calculations can make a valuable contribution to this process through enabling the spectra of structural variations of potential carriers with different charge states to be determined much more readily than would be possible via experimental measurements. There has been a number of studies that have characterised the infrared (IR) spectra of PAH-based molecules through predominantly density functional theory (DFT) calculations within the context of assigning features in the DIBs [16][17][18][19][20][21][22]. The focus of this study is the UV/ Vis region of the spectrum, and the excited states of neutral and cationic PAH molecules have been studied using a variety of computational approaches. Several authors have used semi-empirical approaches, such as ZINDO/S [23], to determine the excited states of PAH molecules [24][25][26][27][28][29], with the low computational cost of these methods allowing large PAH molecules to be studied. At the other end of the computational cost spectrum, the excited states of naphthalene [30] and the ground-and excited-state potential energy surfaces of the naphthalene radical cation have been characterised at the complete active space self-consistent field (CASSCF) level with multiconfigurational perturbation theory (CASPT2) [31]. Second-order coupled cluster theory calculations of the excited states of neutral naphthalene have also been reported [32]. These types of approaches have also been applied to larger PAH molecules; for example, pyrene has been studied with multireference perturbation theory [33] and multireference configuration interaction methods [34]. Time-dependent density functional theory (TDDFT) lies at an intermediate computational cost, and the excited states of neutral and cationic forms of a wide range of PAH molecules have been studied using TDDFT [14,[35][36][37][38][39].
With an appropriate choice of exchange-correlation functional, accurate excitation energies can be determined [36,40], while for some functionals an incorrect energy ordering of the states for naphthalene is predicted [32]. A recent study of the radical cations of naphthalene and pyrene assessed the potential energy surfaces predicted by TDDFT methods through comparison with CASPT2 calculations; it was shown that TDDFT approaches were sufficiently accurate to provide a foundation for the description of the photophysical behaviour that could be applied to study much larger PAH molecules [41].
In order to enable a direct comparison between observational data and computational simulations, it is necessary to go beyond a calculation of the just the excitation energies and the associated oscillator strengths and incorporate vibronic effects to simulate vibrationally resolved spectra. This represents a significantly greater challenge since, in addition to the excitation energies and transition dipole moments, it requires the vibrational frequencies and normal modes of the ground and excited states to be described accurately. Simulated vibrationally resolved spectra for neutral acenes, including naphthalene, and pyrene have been reported [42,43]. Dierksen and Grimme reported calculations of the vibronic structure in electronic spectra for the low-energy bands of range of organic molecules, including pyrene, with the Franck-Condon (FC) and Franck-Condon-Herzberg-Teller (FCHT) approximations based upon TDDFT calculations with the B3LYP exchange-correlation functional [42]. The second derivatives of the energy required to determine harmonic frequencies were evaluated by finite difference of analytical gradients. The calculations provided a good description of the vibrational structure although an energy shift needed to be applied to the vertical excitation energy. More recently, the calculation of vibrationally resolved absorption spectra of acenes and pyrene has been reported using a real-time generating function method in combination with TDDFT and the PBE0 functional [43]. For naphthalene, the second-order approximate coupled-cluster level was also used. The vibronic effects were simulated using the Franck-Condon approximation, and subject to an energy shift to the excitation energy, the vibrational structure was predicted accurately.
In this study, we focus on the UV/Vis region of the spectrum and explore how DFT-based methods can be used to accurately determine the vibrationally resolved UV/Vis spectra of representative PAH molecules, namely naphthalene and pyrene. These methods are then applied to predict the corresponding spectra for the cationic species.

Computational details
The ground-state structures were optimised with DFT using the B3LYP [44] exchange-correlation functional and the 6-31G* basis set. Vertical excitation energies and transition dipole moments were determined with TDDFT using the CAM-B3LYP functional [45] and 6-311+G* basis set, which represents a higher level of theory that is more appropriate for describing the excited states of these systems. The major limiting factor for these simulations is determining the vibrational frequencies for the excited states. The ability to compute vibrational frequencies for excited states has been enhanced greatly through the implementation of analytical second derivatives within TDDFT [46]. We note that these derivatives have been implemented for TDDFT with the Tamm-Dancoff approximation (TDA) [47] and are not available for all functional types. For this reason, harmonic vibrational frequencies and the normal modes were determined with TDDFT/TDA using B3LYP/6-31G*, and the vibrational frequencies were scaled by the factor of 0.96 to account for anharmonicity [48]. Even with the availability of analytical derivatives, determining vibrational frequencies can remain problematic, particularly for the higher-energy states. Despite accounting for changes in the energy ordering of the states during the geometry optimisation process and relaxing symmetry constraints, we were unable to determine minimum energy structures with no imaginary frequencies for some higher-energy states. A definitive reason for this has not been established, but we observe that this appeared to be prevalent when other electronic states were close in energy and the eigenvector describing the transition had significant contributions from several single excitations.
The spectra were simulated using the FCclasses program [49] at a temperature of 10 K. At this temperature the origin of all transitions is the vibrational ground state of the lower-energy (ground) state. Spectra were simulated based upon FC and FCHT approximations. The difference between these two approximations is that the FCHT calculations incorporate the effects of the transition dipole moment changing with the nuclear coordinates. The spectra were generated through convolution with a Lorentzian function with full width half maximum of 0.015 eV.

Results and discussion
The calculated CAM-B3LYP/6-311+G* vertical excitation energies for the low-energy states of naphthalene are given in Table 1 along with a description of the electronic transitions in terms of the molecular orbitals which are shown in Fig. 1. The transitions have → * character with the exception of the 1A u state which arises from a →Rydberg 3s excitation. In the context of the UV/Vis spectrum, only the excitations to form the 1B 2u and 2B 3u states have significant intensity. These transitions are calculated to lie at 4.65 eV and 6.05 eV, which are close to the values from experiment of 4.60 eV and 5.90 eV, and compare well with the previous work that reported values of 4.89 eV and 6.22 eV from second-order approximate coupled cluster (CC2) calculations [43].
For the calculations of the vibrationally resolved spectrum, we focus on the lower-energy 1B 2u state. We were unable to successfully optimise the structure of the 2B 3u state so that it did not give an imaginary frequency using TDDFT. The 1B 2u state arises from the 1a u → 2b 2g transition and the structural parameters associated with the ground and 1B 2u states are shown in Fig. 2. There is little variation in the bond angles between the two states, and the most significant change in the structure is an increase in Table 1 Calculated CAM-B3LYP/6-311+G* vertical transition energies and oscillator strengths for naphthalene and naphthalene +• The orbitals are shown in Fig. 1, and in naphthalene +• the unpaired electron has spin

Molecule
State the bond length between the carbon atoms labelled B and C (and D and E) combined with a reduction in the bond length between carbon atoms C and D. These changes are consistent with the form of the molecular orbitals associated with the transition. Examination of the orbitals in Fig. 1 shows the 1a u orbital to be bonding between atoms B and C, while the 2b 2g orbital has a node between these two atoms. Similarly, the 1a u orbital has a node between atoms C and D, while the 2b 2g orbital has bonding character between these atoms. The vibrationally resolved absorption band for the 1a u → 2b 2g transition leading to the 1B 2u state is shown in Fig. 3. The calculated spectra capture all of the features observed in experiment and allow the vibrational structure to be assigned with some confidence. There is no significant difference between the FC and FCHT calculations for the first two vibrational lines; however, for the remainder of the band the intensity from the FCHT calculation is greater, and the relative heights of the peaks are closer to experiment. The transition involves an excitation from a orbital to a orbital with more antibonding character. This is consistent with the observed notable changes in the structure and the less dominant 0-0 transition. While the calculations with the FCHT approximation lead to an increase in the intensity of the higher-energy lines relative to the 0-0 transition, the 0-0 line does remain the most intense, whereas in experiment the most intense line corresponds to a higher-energy transition. The assignment of the vibrational structure to the vibrational modes of the excited state is given in Table 2. The calculations find the most intense line occurs for the 0-0 transition, which is consistent with there being little change of the structure in the excited state relative to the ground state. The remaining vibrational structure arises predominantly from just four vibrational modes which have A g symmetry. These modes are shown in Fig. 4 and comprise in-plane stretching of carbon-carbon bonds and wagging of hydrogen atoms. The assignment of the vibrational structure is consistent with a previous study that used a different approach to simulate the spectrum [43]. It is also possible to examine the effects of the FCHT approximation in more detail. A further comparison of the FC and FCHT approximations is presented in Supporting Information ( Figure S1), which shows that the FCHT approximation leads to a greater number of transitions with nonzero intensity since transitions to modes without A g symmetry gain intensity. However, the difference in the resulting spectra arises predominantly from an increase in the intensity for the transitions associated with the 33 Fig. 3 Calculated and experimental spectra for the 1a u → 2b 2g transition leading to the 1B 2u state in naphthalene. The computed spectrum has been shifted by −0.08 eV to align with the experiment. The lower panel shows individual transitions for the FCHT spectrum, and the assignment of the dominant transitions is given in Table 2. The experimental data are adapted from reference [50]  The accurate description of the spectrum for neutral naphthalene provides a good platform to study the cation. The calculated vertical excitation energies for naphthalene +• are also shown in Table 1. In the cation, the 1a u orbital is singly occupied, and the resulting low-energy part of the spectrum is richer with four transitions with nonzero oscillator strength. The first two of these involve excitation of a electron from the orbitals of lower energy to the singly occupied orbital. The remaining two transitions include excitation of an electron to higher-energy unoccupied orbitals with * character. The available experimental data show that the transition energy for the 1B 3u state lies at about 1.85 eV which is lower than the calculated value of 2.09 eV, but within the typical error that is associated with TDDFT calculations.
The structure of naphthalene +• in the different states is shown in Figure S2. For one of the states (1B 2u ), there is a distinct deviation from planarity with the hydrogens bonded to the end carbon atoms distorted from the plane of the molecule. Figure 5 shows the calculated vibrationally resolved spectrum for naphthalene +• and includes the four transitions with nonzero oscillator strength, and the decomposition of the spectrum into the individual contributions from the four states is also shown. This spectrum has been computed using the FCHT approximation. We note that we were unable to determine the vibrational frequencies for the 2B 2u state with TDDFT and have used frequencies determined using a selfconsistent field method [51,52] to enable a full prediction of the lower-energy region of the spectrum. The experimental data include the lowest energy state, but do not fully capture the higher-energy states.
Once a constant energy shift of -0.18 eV is applied to the calculated spectrum, it provides an excellent description of the low-energy band for which experimental data are available, with the distinct vibrational structure reproduced accurately. The assignment of the observed vibrational structure is given in Table 3. For the 1B 3u state, the most intense line corresponds to the 0-0 transition and the remaining vibrational structure can be assigned to contributions from the two vibrational modes 9 and 35 . These two modes are illustrated in supporting information ( Figure S3) and are analogous to the 9 and 33 modes of the 1B 2u state of the neutral molecule.
The band corresponding to the 1B 2u state is much weaker and its intensity is gained through the FCHT approximation, and within the FC approximation, this band is not observed. The different nature of this band is reflected in the vibrational structure which arises from the 5 mode which has B 3g symmetry and involves atomic displacements out of the plane of the molecule. There is some evidence of this band being present in the experimental spectrum, which indicates the importance of going beyond a purely FC approximation-based description of the spectra to describe all of the  Table 3, and the experimental data are adapted from reference [15] features that are observed in experiment. The two bands of higher energy that arise from the 2B 2u and 2B 3u states are not covered by the experimental data. One of the common features of these two bands is that the 0-0 transition is not the most intense line. As indicated in Table 1, these two transitions involve excitation of an -spin electron in a transition with → * character in contrast to the two bands at lower energy that arise from -spin electron in transitions with → character. As a consequence, a greater change in the geometry is to be expected which is consistent with the 0-0 transition not being dominant. The relevant vibrational modes ( Figure S3) involve in-plane motions which have common features with the those leading to the vibrational structure in the neutral molecule and the 1B 3u state of the cation. The analysis is now extended to the larger molecule pyrene. The excited states of pyrene have been studied in earlier work using multiconfigurational wave function-based methods and TDDFT [33,34,36,43] and the vibrational resolved spectra for the lowest allowed transition (1B 3u ) have been reported [42,43]. The computed vertical excitation energies along with a description of the transition in terms of the molecular orbitals, shown in Fig. 6, are given in Table 4. For neutral pyrene, the calculations predict three transitions

Fig. 6
Relevant molecular orbitals for the low-energy transitions in pyrene and pyrene +• . In pyrene +• the 2b 3g orbital is singly occupied Table 4 Calculated CAM-B3LYP/6-311+G* vertical transition energies and oscillator strengths for pyrene and pyrene +• The orbitals are shown in Fig. 1 3.73 0.061 2b 2g → 2a u with significant oscillator strength, all of which can be characterised as → * transitions. This is consistent with the spectrum measured in experiment [15], where the bands are found to lie at 3.8 eV, 4.7 eV and 5.3 eV. The corresponding calculated values of 3.95 eV, 5.01 eV and 5.61 eV lie within about 0.3 eV of the observed values which is consistent with the typical error associated with TDDFT. While the energy separation between the 2B 2u and 2B 3u states is predicted accurately, the separation between the 1B 2u and 2B 2u states is overestimated. Although it does not affect this study directly, we note that the calculated energy for the 1B 3u state is much too high, compared with an experimental value of 3.4 eV [53]. The TDDFT excitation energies can be compared with those from the CC2 method which gives values of 4.07 eV, 4.87 eV and 5.67 eV with the def2-TZVP basis set for the three observed transitions [43]. While these calculations predict the energy separation between the two lower states well, the separation between the higher-energy states is too large. This highlights the challenge of predicting the vertical excitation energies with sufficient accuracy that allows direct comparison with experimental spectra that involve contributions from a range of electronic transitions. Multiconfigurational wave function-based methods have considered the lowest two energy states (1B 2u and 1B 3u ) [33,34], so while they predict accurate values for these states, their performance for the higher-energy states has not been established. Figure 7 shows a comparison between the experimental spectrum and the calculated spectrum using the FCHT approximation that incorporates the lowest two energy allowed transitions and the bond lengths for the different excited states are shown in Figure S4. We were unable to eliminate all imaginary vibrational frequencies for the TDDFT calculations for the 2B 3u state. The consequence of the overestimation of the vertical transition energy for the 2B 2u state is apparent with the calculated band clearly displaced relative to experiment. Within each band, the vibrational structure is described well, and the distinct structure observed in experiment can be identified in the simulated spectrum. The assignment of the vibrational lines is given in Table 5 with the vibrational modes depicted in Supporting Information ( Figure S5). Similar to naphthalene, the vibrational structure arises from symmetric vibrational modes that involve atomic displacements within the plane of the molecule, for example expansion and compression of the molecule in the x and y directions. A comparison between the FC and FCHT approximations is shown in Figure S6, which shows that the FCHT approximation leads to an increase in intensity for some of the higher-lying vibrational lines, particularly those associated with 62 (1B 3u ) state and 61 (2B 2u ) state.
The computed and experimental spectra for pyrene +• are shown in Fig. 8 Table 5. Experimental data adapted from reference [15]  excited states are shown in Figure S7. The spectrum is dominated by the 2B 3u state, and for this state, the observed vibrational structure is reproduced well by the calculations. The assignments of the vibrational structure (Table 6) and the visualisation of the vibrational modes ( Figure S8) show that the modes responsible for the vibrational structure are similar to those for neutral pyrene. The lower-energy transitions are much weaker in comparison, and although they have a rich vibrational structure, it is not possible to assess the accuracy of the calculations. There is a notable discrepancy for the higher-energy band 2B 2u . In the experimental spectra, structure is observed, while the calculated band is much weaker and is difficult to distinguish.

Conclusion
The simulation of vibrationally resolved UV spectra is an important component for enabling the identification of the carriers of DIBs. The focus of this study was neutral and cationic PAH molecules that are widely believed to the carriers of DIBs. It has been shown that DFT-based methods are capable of providing an accurate description of the available experimental data for the vibrational structure of neutral and cationic naphthalene and pyrene. The calculations show that the vibrational structure is associated with a small number of symmetric vibrational modes that involve in-plane atomic displacements. Furthermore, the FCHT approximation can lead to significant changes in the calculated spectra that improve the agreement with experiment compared with the FC approximation. More generally, the study has highlighted a number of challenges that need to be addressed if simulations of vibrationally resolved spectra are to achieve a level of accuracy to be used to assign molecular carriers of complex spectroscopic data in the absence of experimental measurements. The typical error associated with TDDFT calculations is usually considered to be approximately 0.3 eV, but this is not sufficient to predict the energy separation between different states to the required level of accuracy. High-quality multiconfigurational wave function-based methods, such as CASPT2 or multiconfigurational configuration interaction (MRCI), have the potential to improve upon this, but these methods are difficult to apply widely. Determining the vibrational frequencies for the excited states is also problematic.   Table 6, and the experimental data are adapted from reference [15]  In this study, it was not possible to eliminate all imaginary frequencies in the vibrational analysis for some of the higher-lying vibrational states. This was particularly evident when other states where close in energy. Other methods for obtaining vibrational modes for the excited states such as those based upon SCF methods provide another approach; however, these approaches are not well suited for describing excited states that are a mixture of single excitations which is the case for many of the transitions in PAHs. Alternatively, this problem can be avoided through different methods for determining the vibrationally resolved spectra such as the real-time generating function approach [32].