Low-Lying ππ* States of Heteroaromatic Molecules: A Challenge for Excited State Methods.

The description of low-lying ππ* states of linear acenes by standard electronic structure methods is known to be challenging. Here, we broaden the framework of this problem by considering a set of fused heteroaromatic rings and demonstrate that standard electronic structure methods do not provide a balanced description of the two (typically) lowest singlet state (La and Lb) excitations. While the Lb state is highly sensitive to correlation effects, La suffers from the same drawbacks as charge transfer excitations. We show that the comparison between CIS/CIS(D) can serve as a diagnostic for detecting the two problematic excited states. Standard TD-DFT and even its spin-flip variant lead to inaccurate excitation energies and interstate gaps, with only a double hybrid functional performing somewhat better. The complication inherent to a balanced description of these states is so important that even CC2 and ADC(2) do not necessarily match the ADC(3) reference.

Introduction

The extensive computational investigations of (hetero)aromatic systems were prompted by the importance and broad applications of these compounds in organic electronics. In particular, the investigations of electronically excited states with theoretical tools should allow not only an in-depth understanding of the properties of known molecules but also the efficient design of new compounds. In this context, time-dependent density functional theory (TD-DFT)[1,2] within its standard approximations (i.e., the linear-response, adiabatic approximation) has become the primary framework,[3] mainly due to its good compromise between accuracy and computational efficiency. Out of the many distinct types of molecular excitations present in conjugated molecules, local ππ* states (with prevailing single excitation character) are typically considered as the least problematic for TD-DFT. In π-conjugated systems, these states are of great relevance for both absorption and emission properties,[4−7] and they play a major role in determining available decay channels.[8−11]

Despite the general reliability of TD-DFT, several investigations uncovered sizable errors in the description of low-lying ππ* states of fairly simple organic compounds. For instance, in 2001 Grimme et al.[12] reported an imbalanced description of the two lowest singlet states of oligoacenes, La and Lb. The La and Lb notation[13] was originally introduced for polycyclic alternant hydrocarbons with La corresponding to the bright state of dominant HOMO → LUMO character and Lb corresponding to the dark state encompassing nearly equal contributions from the HOMO–1 → LUMO and HOMO → LUMO+1 transitions. The low oscillator strength of the Lb state was explained by the cancellation of the transition dipole moments associated with these two contributions.[14] As emphasized in ref (12), the La excitations are significantly underestimated by standard local and semilocal functionals (such as the popular PBE[15] and B3LYP[16]), with a state order inversion in the case of naphthalene and a large excitation energy downshift for larger acenes. Due to the fundamental importance of oligoacenes, the conundrum has gained significant interest in the literature.[5,17−31] Large improvements of the La excitation energies were later reported with the use of range-separated hybrid functionals,[19−21,24] however at the expense of deteriorating the Lb excitation energy values.[21] The difficulty to provide a balanced description was attributed to the significant impact of contributions from double-excitation (mainly for Lb), that cannot be properly described with the standard adiabatic TD-DFT implementations.[25] Therefore, the La-Lb problem originates from the description of both states. Indole (a building block of tryptophan amino acid) and several structurally related compounds were found to behave similarly.[32−34] Indeed, it was shown[33] that both hybrid and meta-GGA functionals predict a wrong ordering of the La and Lb ππ* states, whereas range-separated hybrid functionals, despite providing a qualitatively correct ordering of the states, predict much too small energy gaps compared to the experimental values. The inversion of the La and Lb states was also spotted for 9H-adenine[35] by comparing TD-DFT estimates to high level reference values obtained with wave function-based approaches, such as EOM-CCSD(T)[36] (equation of motion-coupled cluster singles and doubles with perturbative triples) or CASPT2[37] (complete active space second-order perturbation theory), even though the ambiguity still remains in this case.[35] More recently, two of us have unraveled similar discrepancies for the low-lying ππ* states of thiophene and thienoacenes,[38] which constitute popular building blocks in organic electronics. Regardless of the exchange-correlation functional used, we found not only a spurious state inversion but also a wrong distribution of oscillator strengths and erroneous potential energy surfaces. In contrast to TD-DFT, the performances of the several correlated single reference methods including contributions from double excitations such as CC2[39] (approximate coupled cluster singles and doubles) and ADC(2)[40] (algebraic diagrammatic construction up to second-order) were found rather satisfying.[10,38]

The present contribution explores La- and Lb-like excitations in a large and diverse set of fused aromatic and heteroaromatic compounds (Scheme ). These are typically the lowest ππ* excited states in the spectrum and are consequently of huge chemical and physical relevance. In contrast to earlier case studies, dealing mostly with oligoacenes and occasionally with specific compounds relevant to applications, here we generalize the problem to a broader class of heteroaromatic molecules and propose a simple diagnostic for identifying these challenging excited states. Our objective is to pinpoint the excited state methods providing a properly balanced description of the two states. We critically examine the performance of standard TD-DFT, using several functionals and wave function based approaches (ADC(2) and CC2) as well as nonstandard TD-DFT based (spin flip (SF),[41] double hybrid[42]) methods, that are all likely to be used for “real-life” applications due to their reasonable computational costs. The paper is organized as follows. In Section , we provide computational details. In Section , the naphthalene example is used as an illustrative case, followed by the examination of three exemplary heteroaromatic systems and the overall analysis of the excitation energy trends for ten different compounds. Concluding remarks are given in Section .

Scheme 1

Investigated Heteroaromatic Compounds

Computational Details

The optimized structures and corresponding transition energies are listed in the Supporting Information (SI). If not stated otherwise, the aug-cc-pVTZ[43] atomic basis set was used throughout.

Ground state geometries were optimized at the MP2 level (employing the resolution of identity, RI)[44] using the Turbomole 6.5 package.[45] Excited state computations with TD-DFT (PBE and PBE0[46] functionals), CC2, and ADC(2) (which can be seen as an approximation to CC2)[47] were performed with Turbomole 6.5. The latter two methods were employed using the frozen core approximation and the RI approach (with an auxiliary aug-cc-pVTZ basis set taken from the Turbomole library).[48]

TD-DFT computations with the M06-2X,[49] M06-HF,[50] BHHLYP,[51,52] ωB97X-D,[53] and LC-PBE*[54] functionals, as well as TD-HF calculations,[55] were performed with Gaussian09 (version D.01).[56] For both M06-2X and M06-HF, the ultrafine integration grid was employed to ensure numerical stability. Note that the long-range corrected LC-PBE functional was optimally tuned (here notation LC-PBE*) according to the nonempirical procedure described in ref (57). As such, the range separation parameter γ was optimized to minimize the function |εHγ(N)+IPγ(N)| + |εHγ(N+1)+IPγ(N+1)|, where εHγ is the energy of the HOMO orbital and IPγ is the vertical ionization potential of the neutral (N) and anionic (N+1) system, N being the number of electrons. For those systems (I, III, VI, VIII, and IX in Scheme ) where the HOMO level of the anion was close to zero or positive (indicating an unbound electron), the tuning was solely based on the HOMO of the neutral system, i.e., the function |εHγ(N)+IPγ(N)| was minimized by varying γ. For some systems, the default SCF convergence parameters led to higher energy solutions for the cation, typically resulting in large IPs and large optimal γ values. The Stable = Opt approach implemented in Gaussian was then used to ensure the convergence to the lower energy solution.

CIS/CIS(D)[58] and B2LYP/B2PLYP[42] (within the Tamm-Dancoff approximation[59]) computations were performed with the Orca 3.0.2 software.[60] Here B2LYP denotes a global hybrid functional (53% of exact exchange) that is underlying the B2PLYP double hybrid.

In addition, low-lying transitions were computed with the spin-flip (SF) version of TD-DFT (SF-DFT)[41] in combination with BHHLYP,[51,52] i.e., SF-BHHLYP.[41] Excitation energies were also computed with the algebraic diagrammatic construction up to third-order (ADC(3)).[61] Due to the steep computational scaling of this method (M6) and large memory requirements (M4), where M is the number of basis functions, the ADC(3) computations were converged with a smaller aug-cc-pVDZ[43] atomic basis set. To obtain our ADC(3) best estimates, basis set corrections based on ADC(2) computations (i.e., E(aug-cc-pVTZ)-E(aug-cc-pVDZ)) were added to the ADC(3)/aug-cc-pVDZ values. Due to the generally weak basis set dependence of the ππ* states, these ADC(3) best estimates are expected to be close to the actual ADC(3)/aug-cc-pVTZ values. SF-BHHLYP and ADC(3) computations were performed with the Q-Chem 4.3 package.[62]

The spectral simulations of acridine (compound VII in Scheme ) were performed with the Newton-X package.[63] The nuclear configurations used for the spectral simulations were sampled by an uncorrelated Wigner distribution[64,65] in the ground state (the Hessian was obtained by reoptimizing the structure at the PBE0/aug-cc-pVDZ[43] level with the Turbomole package). 200 structures were taken, and the vertical excitation energies and oscillator strengths were computed at both the TD-PBE0 and ADC(2) levels using the aug-cc-pVDZ atomic basis set. The transitions were broadened by a Lorentzian using a phenomenological width of 0.05 eV.

Results and Discussion

The Case of Naphthalene La and Lb States

We selected naphthalene, an intensively studied example, to serve as a prototype example for the excitation energy trends found in fused (hetero)aromatic compounds. Naphthalene allows for the illustration of the major issues regarding the imbalanced description of the two lowest ππ* states. First, let us provide an overview of the main conclusions raised in the literature. Although most of the qualitative results shown in Figure have been described previously,[12,19−21,24,25] they were recomputed here to minimize the impact of using different ground state geometries and diverse atomic basis sets. Additional insights are also obtained from the CIS/CIS(D), SF-BHHLYP, and ADC(3) results.

Figure 2

Computed excitation energies of the La (red) and Lb (blue) excited states of a) acridine, b) 3,1,3-benzooxadiazole, and c) benzo[2,1-b:3,4-b′]dithiophene with the aug-cc-pVTZ basis set; see Computational Details for ADC(3).

As pointed out in Grimme’s seminal study,[12] local and semilocal functionals such as PBE (generalized gradient approximation, GGA functional, 0% of exact exchange) and PBE0 (global hybrid functional, 25% of exact exchange) severely underestimate the La excitation energies and provide incorrect state ordering (left, Figure ). As a side note, we remind that some improvements were reported when Tamm-Dancoff approximation (TDA) was used (it fixes the state order when combined with PBE0 although the energy gap remains rather inaccurate).[66] While TDA was found beneficial in some studies,[26,66,67] in others the improvements were attributed to a fortuitous cancellation of errors.[21] We will return to TDA later in the text. The inclusion of a larger portion of exact exchange as in M06-2X (meta-GGA hybrid, 54% of exact exchange) or in range-separated hybrid functionals such as ωB97X-D and LC-PBE* not only upshifts the HOMO → LUMO (La) state toward the reference value but also overshoots the energy of the Lb state. Consistently with its accurate description by range-separated functionals, the La state shows some similarities with charge transfer states and was called “charge transfer in disguise”[21] or “charge transfer-like excitation”.[24] Nevertheless, according to standard analysis tools, such as the Tozer Λ diagnostic based on the overlap between the MOs,[19] there is no net charge transfer,[21] and both states can be characterized as local ππ* excitations. Alternatively, the valence bond picture describes La and Lb as ionic and covalent states, respectively.[17,68] The CIS/CIS(D) excitation energies bring up a relevant trend. Unlike La, the Lb state is highly sensitive to the differential correlation effects introduced by the perturbative correction for contribution from double excitations. We note that, in contrast, the CASSCF analysis of the Lb state wave function shows the dominant contributions from single excitations,[68] but it still misses important contributions from the dynamic correlation. A more detailed analysis with a high level post-Hartree–Fock method (CC3, coupled cluster singles, doubles, and triples[69]) reveals 15% of nonsingles, compared to the 10% in La.[5] Therefore, it is not surprising that the description of the Lb state is rather problematic at the TD-DFT level. This issue was already recognized by Grimme et al.[25] who applied double hybrid functionals to the series of linear and nonlinear acenes, obtaining significant improvements over standard TD-DFT computations. As shown in Figure , B2PLYP indeed produces excitation energies comparable to correlated single reference methods with explicit contributions from the doubles, such as CC2 and ADC(2).

Figure 1

Computed excitation energies of the La (red) and Lb (blue) excited states of naphthalene compared to the experimental estimates taken from ref (12). The aug-cc-pVTZ basis set was used; see Computational Details for ADC(3).

While SF-BHHLYP provides a rather balanced description of the La/Lb states in naphthalene, these results deserve a closer analysis. The transition energy to La is slightly overestimated. On the other hand, SF-BHHLYP underestimates the Lb energy, i.e., it shows the behavior opposite compared to TD-DFT when similar exchange-correlation functionals are used. To understand this difference, one should stress that the SF excitation scheme is best suited for the study of small HOMO to LUMO gaps and for the computation of states with the main contributions coming from HOMO → LUMO and (HOMO)2→ (LUMO)2 transitions. Although other electronic transitions can be computed with SF-DFT, the final expression for these states is not spin complete.[70] This is precisely the situation for the Lb state, for which SF-BHHLYP generates a broken-symmetry solution, that is, the mixing between the singlet and triplet Lb states.

Regarding the experimental values, it is important to note that as the computed vertical excitation energies are not experimental observables, they should not be compared to the experimental band maxima directly but preferably to the results obtained with higher levels of theory.[71] The experimental estimates of vertical excitation energies shown in Figure were back corrected[12] from the accurate measures of adiabatic excitation energies. The resulting energies values (4.13 eV for Lb and 4.66 eV for La) compare very well with our ADC(3) estimates (4.11 eV for Lb, 4.70 eV for La) as well as with earlier CASPT2 computations (4.03[72]/4.24[5]/4.06[73]eV for Lb, 4.56[72]/4.77[5]/4.49[73] eV for La). Predictions from other high level methods include CR-EOM-CCSD(T) (4.13 eV for Lb and 4.79 eV for La)[23] and CC3 (4.27 eV for Lb and 5.03 eV for La; a triple-ζ atomic basis set with no diffuse functions was used).[5] Despite some spread of excitation energies, each of these methods predict relatively large interstate gaps, which is not the case for the lower level methods. Overall, TD-DFT energies are clearly dependent upon the extent of exact exchange, but none of the functionals provides a simultaneous good description of both states. Moderate improvements are achieved with wave function based (ADC(2), CC2) and more sophisticated TD-DFT based (double hybrid, spin flip) methods, but even these results suffer from significant errors. In the upcoming sections, we rely on ADC(3) reference values to evaluate the systematic shortcomings of different excited state methods on a larger number of heteroaromatic compounds. However, given the variations observed among high-level methods (as evident from the naptahlene example), we will restrict the forthcoming discussions to large quantitative deviations (i.e., > 0.2 eV) as well as clear trends.

Criteria for La- and Lb-like States

A preliminary step necessary for the assessment of the performances of the different excited state methods is to distinguish the two states (i.e., La- and Lb-like states) on the set of small and middle sized fused heteroaromatic systems. However, the definition of such states in terms of quasiparticle levels appears somewhat arbitrary, as in practice orbitals obtained from Hartree–Fock or generalized Kohn–Sham methods may be largely distorted when large and diffuse basis sets are used, giving rise to multiple contributions of orbital excitations with sizable coefficients.[74] This is why we relied on natural transition orbitals (NTO) to distinguish between the states (see the SI; also note that the NTO analysis is rather qualitative due to the neglect of correlation effects, the proper treatment of which is important for excitation energies and excited state properties). La is typically well described by a single pair of NTO, while Lb consists of two major configurations, which generally do not have equal weights. Also, in contrast to the oligoacenes, the Lb-like state can have oscillator strength as large as, or even larger than, the corresponding La state. Compounds with permanent dipole moment are characterized by a Lb state presenting a dipole of a magnitude similar to its ground state counterpart, while the values for the “ionic” La are typically larger (the exceptions are compounds VII and IX due to the more symmetric charge distribution).

The patterns of the La and Lb states for the set of Scheme are less systematic than in the oligoacene series. The presence of heteroatom(s) induces some variations on the nature of the ππ* transitions, such as character mixing with nearby πσ* and nπ* states, or between the La and Lb states.[34] For compound VIII, a moiety frequently found in organic electronics, there is even an ambiguity in determining the character of La and Lb since both states belong to the same irreducible representation and, therefore, mix (see the SI for the assignment used herein).

Three Illustrative Compounds

Figure provides the detailed analysis of three individual cases, e.g., acridine, 2,1,3-benzooxadiazole, and benzo[2,1-b:3,4-b′]dithiophene (respectively VII, V, and IX in Scheme ), that manifest the problem illustrated for naphthalene. The reference ADC(3) relative and absolute excitation energies of La and Lb differ significantly in these three compounds, which makes them interesting study cases. In particular, IX and V possess their lowest lying states at a similar energy and a similar energy gap, but the state ordering is reversed. In VII, both states are computed to be nearly degenerate at the ADC(3) level. As can be seen in Figure , the three distinct functionals, PBE, PBE0, and ωB97X-D, fail to reproduce both the absolute excitation energies and the excitation energy gaps predicted by ADC(3). There is a characteristic dependence of excitation energies on the amount of exact exchange, with a similar upshift for both La and Lb. For the three systems, ωB97X-D predicts the La state very close to ADC(3), illustrating its remarkable performance for states of dominant HOMO → LUMO character. On the other hand, the same functional severely overestimates Lb (0.3 to 0.5 eV). The comparison between CIS and CIS(D) uncovers an essential trend: the perturbative double correction has a large impact on on Lb but a much smaller impact on La. Considering the significant double excitation character of Lb, the apparent “good” performance of PBE and PBE0 for this state is most probably fortuitous. The double hybrid, B2PLYP, which incorporates a CIS(D)-like correction, improves upon standard TD-DFT and provides excitation energies comparable to ADC(2). Alternatively, the SF-BHHLYP results are rather unsatisfactory due to the dramatic underestimation of the Lb excitation energy and the inconsistent description of La. Finally, ADC(2) (and similarly CC2), although being the closest to the reference, shows some lack of systematic behavior. Both La and Lb are slightly overestimated, though the trends are not perfectly equivalent for all compounds. While Lb is overestimated for VII, La is too high in V, whereas in IX the energy gap of ADC(3) is well reproduced due to the similar upshift for both states. In fact, this clearly shows that ADC(2)/CC2 might not be a sufficiently accurate benchmark to assess the quality of TD-DFT as already discerned in previous benchmark studies.[5,75] Along this line, the performance of range-separated hybrid functionals might be even superior for excited states of HOMO → LUMO character.

To further demonstrate that the discrepancies shown in Figure have a major impact on the theoretical prediction of absorption properties, we computed the absorption spectra of acridine at two illustrative levels (Figure ). Unlike anthracene, which has an optically dark Lb state, the oscillator strength of Lb in the widely used acridine dye is similar to La. Because of the large gap between the two states, PBE0 predicts two distinct peaks, while ADC(2) predicts two overlapping peaks, which is consistent with the experiment.[76]

Figure 3

Photoabsorption spectra of acridine computed from the Wigner distribution, employing PBE0 and ADC(2) combined with aug-cc-pVDZ.

Statistical Analysis of Excitation Energy Trends

The mean signed deviation for each excited state method (Figure ) best illustrates the overall performances and general trends associated with the set of compounds represented in Scheme . The corresponding mean absolute deviations can be found in the SI. In line with the individual molecular cases, the most striking feature is the difference between CIS and CIS(D). CIS severely overestimates the excitation energies of both La and Lb, and the Lb state energy goes down by a large amount (ca. 0.8 eV) when including the correction for the contribution of the doubles. In contrast, La is rather constant. Given that this effect is characteristic for all the investigated compounds, the CIS/CIS(D) computations ideally serve as a simple diagnostic for identifying the La- and Lb-like states in real life applications: the excitation energies of Lb states are much more sensitive to the dynamical correlation effects, that are absent in CIS. In addition, the overestimation of the ionic La state is rationalized by the well-documented CIS large positive bias for charge transfer states.[77] The introduction of the second order perturbative corrections clearly improves the description but remains insufficient to provide well-balanced excitation energies. A more balanced picture is achieved with CC2 and ADC(2). The two methods give practically the same energies in line with earlier studies.[78,79] The averaged overestimations of both the La and Lb states are around 0.15 and 0.20 eV, respectively, indicating the reliability of both methods for practical applications on medium-sized organic molecules. Nevertheless, the shifts of each state are not systematic (see for instance the naphthalene example in Figure , where La is well positioned but Lb is overestimated, and the irregular deviations for the individual compounds in Figure ), and the similar mean signed deviation for the two states are somewhat misleading. These scattered results illustrate that the treatment of correlation is still incomplete in ADC(2) and CC2. The correlation effects are albeit crucial (as it is apparent from CIS results) to obtain accurate energies for both states.

Figure 4

Mean signed deviations of the La (red) and Lb (blue) states of compounds in Scheme obtained by comparison of a) wave function and b) TD-DFT based methods to ADC(3) reference. Bars with dashed lines correspond to the results with TDA, while dotted lines on top of the CIS correspond to TD-HF. The aug-cc-pVTZ atomic basis set was systematically used.

From the TD-DFT perspective, Figure b confirms that none of the tested functionals provides a balanced picture of the two relevant states. PBE severely undershoots the excitation energy of La, and the more accurate description of Lb is fortuitous. PBE0, which contains a moderate fraction of exact exchange, underestimates the excitation energies of the La states but to a smaller extent. Global hybrids with large amount of exact exchange (M06-2X and BHHLYP) as well as range-separated hybrids (ωB97X-D and LC-PBE*) improve the description of La, owning to its charge transfer-like character, but overshoot significantly the energy of Lb. This is valid for the optimally tuned variant, LC-PBE*, which does not improve the results in comparison with the range-separated hybrid functional with fixed γ, such as ωB97X-D (0.20 bohr–1). The improvement of the HOMO → LUMO excitation within the range-separation framework goes along with the more accurate quasiparticle energies. However, the amount of exact exchange that is optimal for La is not necessarily optimal for multiconfigurational Lb. It is of course possible to tune parameters to specifically minimize the Lb errors, but this is neither a practical nor satisfying solution. As noted earlier, frequency independent (i.e., adiabatic) TD-DFT does not perform well for excited states with strong configuration mixing.[80] Alternatively a functional with 100% of exact exchange, M06-HF, severely overestimates the transition energies of both states. These results generally indicate that the TD-DFT errors for the La and Lb excitation energies are rather systematic and depend mainly upon the treatment of exact exchange.

Excitation energies obtained with TDA are systematically blue-shifted with respect to full TD-DFT (Figure b). The shift is similar for different functionals but somewhat larger when higher fractions of exact exchange are included. La is typically shifted more (∼0.2 eV) than Lb (∼0.1 eV). As seen earlier, CIS, being a Tamm-Dancoff approximation of the TD-HF scheme (dotted line in Figure a), also leads to an upshift of the computed excitation energies. In short, no general improvement originating from TDA can be identified as the accuracy is improved for the functionals underestimating the excitation energies but is deteriorated for the others (such as range-separated hybrids for La and Lb).

In comparison with standard TD-DFT, the double hybrid approach certainly delivers a more balanced treatment of La and Lb in heteroaromatic molecules. B2PLYP gives results similar to both CC2 and ADC(2), highlighting once again the importance of accounting for double excitations, for the Lb state. The SF-BHHLYP energies follow the behavior already observed for the individual compounds, that is a slight overestimation of La for the same magnitude as non-SF TD-DFT energies with a similar amount of exact exchange and an underestimation of the Lb excitation energy.

The overall performance of the different approximations for the calculation of the energy gaps between the two states (see the mean absolute deviations in Figure ) is also relevant given that the relative position of the excited states is sometimes more important than the absolute transition energies. Despite the uncertainty intrinsic to any excited state method, and also to the ADC(3) reference, the gap obtained from TD-DFT is systematically away from the reference values (∼0.4 eV). Unbalanced gaps result, for various cases, in qualitatively incorrect state ordering. TDA reduces the errors for the gap between the two states, essentially because the shifts for La and Lb are not equivalent. CIS provides incorrect gaps, and significant improvements are achieved by CIS(D). The smallest (although still relatively large) deviations from the reference are obtained with (TDA-)B2PLYP, ADC(2), and CC2 methods.

Figure 5

Mean absolute deviations of energy gaps (La-Lb) compared to the ADC(3) reference. Dashed lines correspond to the results with TDA, while dotted lines on top of the CIS correspond to TDHF. All results are obtained with the aug-cc-pVTZ atomic basis set.

Conclusions

We demonstrated and generalized the problem of the imbalanced description of La and Lb, which is well-known for oligoacenes, to a set of heteroaromatic-fused systems. A CIS/CIS(D) computational check was proposed as a simple diagnostic for identifying these two problematic (typically lowest) excited states, with Lb being highly sensitive to the correlation effects introduced by the doubles (D) correction. A pronounced difference between the CIS/CIS(D) excitation energies is expected to foretell significant difficulties when employing the most widely used excited state methods. As a matter of fact, TD-DFT (within its standard approximations) does not provide balanced excitation energies nor accurate interstate gaps, which occasionally results in the spurious inversion of the states. Yet, TD-DFT outperforms CIS, thanks to the approximate treatment of correlation, absent in CIS. Pure DFT functionals and those with a small amount of exact exchange tend to underestimate the excitation energies of La. Functionals with a larger amount of Hartree–Fock-like exchange as well as range-separated hybrid functionals describe La very accurately but overestimate the energy of Lb. The benchmarking of such functionals exclusively on excited states with dominant HOMO → LUMO character is therefore somewhat biased. Since changing the functional does not solve the overall issue, it is likely that the approximations used in standard TD-DFT (i.e., adiabatic approximation) are at the origin of the problem. Some improvements over standard TD-DFT are achieved by using a double hybrid functional in which part of the correlation is described by a posteriori (D)-like correction. Spin-flip DFT, which is generally a very good approach for the description of low-lying energy states in molecules with diradical or triradical character, is not able to accurately reproduce the relative energies between La and Lb in heteroaromatic molecules. Better performances are obtained with ADC(2) and CC2 albeit higher levels of theory are necessary to reach high and robust accuracy.