# Classification of doubly excited molecular electronic states

Mariana T Do Casal, Josene M Toldo, Mario Barbatti, Felix Plasser

► To cite this version:

Mariana T Do Casal, Josene M Toldo, Mario Barbatti, Felix Plasser. Classification of doubly excited molecular electronic states. Chemical Science, 2023, 14, pp.4012 - 4026. 10.1039/d2sc06990c. hal-04102333

HAL Id: hal-04102333

https://hal.science/hal-04102333

Submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
EDGE ARTICLE

View Article Online

View Journal | View Issue

<figure><img src="image_2.png" /></figure>

Check for updates

Cite this: *Chem. Sci.*, 2023, 14, 4012

All publication charges for this article have been paid for by the Royal Society of Chemistry

# Classification of doubly excited molecular electronic states†

Mariana T. do Casal, ID a Josene M. Toldo, ID a Mario Barbatti ID ab and Felix Plasser ID *c

Electronic states with partial or complete doubly excited character play a crucial role in many areas, such as singlet fission and non-linear optical spectroscopy. Although doubly excited states have been studied in polyenes and related systems for many years, the assignment as singly vs. doubly excited, even in the simplest case of butadiene, has sparked controversies. So far, no well-defined framework for classifying doubly excited states has been developed, and even more, there is not even a well-accepted definition of doubly excited character as such. Here, we present a solution: a physically motivated definition of doubly excited character based on operator expectation values and density matrices, which works independently of the underlying orbital representation, avoiding ambiguities that have plagued earlier studies. Furthermore, we propose a classification scheme to differentiate three cases: (i) two single excitations occurring within two independent pairs of orbitals leaving four open shells ($D_{OS}$), (ii) the promotion of both electrons to the same orbital, producing a closed-shell determinant ($D_{CS}$), and (iii) a mixture of singly and doubly excited configurations not aligning with either one of the previous cases ($D_{mix}$). We highlight their differences in underlying energy terms and explain their signatures in practical computations. The three cases are illustrated through various high-level computational methods using dimers for $D_{OS}$, polyenes for $D_{mix}$, and cyclobutane and tetrazine for $D_{CS}$. The conversion between $D_{OS}$ and $D_{CS}$ is investigated using a well-known photochemical reaction, the photodimerization of ethylene. This work provides a deeper understanding of doubly excited states and may guide more rigorous discussions toward improving their computational description while also giving insight into their fundamental photophysics.

Received 21st December 2022

Accepted 14th March 2023

DOI: 10.1039/d2sc06990c

rsc.li/chemical-science

## 1. Introduction

Electronic states with doubly excited characters have aroused interest and generated lively debate in recent years. Their involvement in singlet fission¹⁻³ can provide a promising route towards highly efficient photovoltaic devices, but they are also of particular interest in other technological applications, such as non-linear optical spectroscopy,⁴⁻⁶ and thermally activated delayed fluorescence.⁷ In the photochemistry of polyenes⁸⁻¹⁰ and derived systems, such as carotenoids,¹¹,¹² states with partially doubly excited character play a crucial role and have

been investigated for over 50 years, yet, still inciting fiery discussions. Despite this substantial interest, there is no well-defined framework for classifying doubly excited states or even a well-accepted definition of doubly excited character. Different authors use different definitions; consequently, even the simple example of the $A_g$ state of butadiene, and its assignment as singly or doubly excited, has recently sparked considerable controversy.¹³,¹⁴

Indeed, the description of doubly excited states is still a challenge for computational chemistry.¹⁵ Many commonly used methods, such as linear-response time-dependent density functional theory (TDDFT),¹⁶ second-order approximate singles and doubles coupled cluster (CC2), or the second-order algebraic diagrammatic construction (ADC(2))¹⁷ fail in their description. Describing double excitations within a single-reference framework requires going up hierarchies and using more intricate and expensive methods. For example, the CC3 (ref. 18) and ADC(4) methods¹⁹ are considered appropriate, whereas even ADC(3)²⁰ can be problematic.¹⁹,²¹ Moreover, in a multireference framework,²² one can produce accurate descriptions of doubly excited states. However, these methods are accompanied by the ever-present problems of choosing an

<footnote>ªAix-Marseille University, CNRS, Marseille, France</footnote>

<footnote>ᵇInstitut Universitaire de France, 75231, Paris, France</footnote>

<footnote>Department of Chemistry, Loughborough University, Loughborough, LE11 3TU, UK.
E-mail: f.plasser@lboro.ac.uk</footnote>

<footnote>† Electronic supplementary information (ESI) available: An alternative discussion of $\Omega$ (S1); $CO_2$ example with 1TDM norm larger than 1 (S2); generalization of the excitation number (S3); DFT/MRCI parametrization comparison (S4, S5); formaldehyde dimer orbitals and discussion of doubly excited states within ADC(3) (S6); density descriptors of the polyenes' series (S7); density descriptors of ethylene's cycloaddition (S8); density descriptors of tetracene, diketopyrrolopyrrole derivatives and s-tetrazine (S9). See DOI: https://doi.org/10.1039/d2sc06990c</footnote>
Edge Article

Chemical Science

appropriate active space and related parameters. Spin-flip methods<sup>23,24</sup> and state-specific orbital-optimized DFT<sup>25,26</sup> present themselves as interesting alternatives but require special care in their applications as well. The problem is exacerbated by the fact that there is no clear rule to indicate when such more involved methods are required, and related questions are heavily contested in the literature. A typical example is an ongoing discussion of which methods are suitable to describe the lowest A<sub>g</sub> state of butadiene.<sup>8,13,14</sup> Therefore, having a well-defined and method-independent quantifier for double excitation character could greatly help for issues of this type.

Aside from methodological questions, it is also desirable to gain a deeper understanding of the underlying physics of doubly excited states with the eventual goal of designing optimised molecules for specific tasks. At this point, it is particularly interesting to contrast doubly excited states on individual molecules with intermolecular doubly excited states in terms of their energies and wave function properties. However, no rigorous and transferable classification scheme exists that would allow comparing these cases meaningfully. Note that discussions of classification schemes for doubly excited states in the literature are restricted to two-electron atoms where notably different physics is at play.<sup>27,28</sup>

To obtain a basic definition of singly or doubly excited character, one might sum over the weights of all singly excited configurations (denoted %T₁ henceforth). However, such an assignment is only meaningful within a given wave function model and set of reference orbitals. As a consequence, the assignment may vary if a different computational method is chosen, and it is not even immediately clear whether the concept of a doubly excited state possesses intrinsic physical meaning at all (cf. ref. 29). More specifically, it is unclear whether any given doubly substituted Slater determinant should be interpreted as contributing to correlation or as an actual double excitation.<sup>13,14</sup> Furthermore, one should realise from a fundamental physical viewpoint that a double excitation is represented by four correlated particles—two excitation holes and two excited electrons. Formally, such two-body processes should not be described by orbitals but by geminals.<sup>30</sup> Therefore, significantly enhanced complexity can be expected compared to singly excited states. As a consequence, doubly excited states have remained quite elusive and ambiguous in the discussions so far.

We propose solving the first problem, the assignment of doubly excited states, by using density matrices. They are well-defined independently of the wave function model, thus, allowing us to extract molecular orbital pictures and numerical descriptors from correlated wave functions.<sup>31-33</sup> More specifically, we base our analysis on transition and difference density matrices, cancelling out correlation contributions also present in the ground state. The use of density matrices provides meaning to these descriptors *via* their connection to physical observables. Specifically, we elaborate on the viewpoint that a doubly (or higher) excited state is a state that cannot be coupled to the ground state with any conceivable one-electron operator (cf. ref. 31 and 34). Secondly, to deal with the enhanced complexity of doubly excited states, we combine and

contrast several analysis methods to obtain a well-defined and simple yet comprehensive picture. These methods amount to the squared norm of the 1-electron transition density matrix (1TDM),<sup>34,35</sup> the promotion number<sup>36</sup> based on the attachment and detachment densities, and the occupation of the natural orbitals (NO), collectively, through the number of unpaired electrons<sup>37</sup> or, individually, through the occupation of the lowest unoccupied NOs (LUNO and LUNO+1).<sup>38</sup> In addition, we apply an extension of the excitation number as defined by Barca *et al.*<sup>14</sup>

This work aims to comprehensively describe molecular doubly excited states in quantum chemistry computations. We start with a Theory section presenting various definitions of doubly excited character, contrasting different limiting cases for doubly excited states, and discussing the underlying energy contributions. Three illustrative examples follow (Fig. 1). We use the formaldehyde dimer to illustrate the limiting case of a doubly excited state involving two independent pairs of open-shell orbitals (denoted the D<sub>OS</sub> case). Subsequently, we study polyenes highlighting the complexity of their A<sub>g</sub> excited states involving a mix of partial doubly and singly excited character (denoted D<sub>mix</sub>). To examine the interconversion between the open- (D<sub>OS</sub>) and closed-shell (D<sub>CS</sub>) limiting cases, we investigate the dimerization of ethylene. Finally, the three different archetypes of doubly excited states are reviewed in a more extended set of molecules.

## 2. Theory

### 2.1 Definition of doubly excited character *via* transition density matrices

Double excitations are traditionally defined *via* the %T₁ values, which reflect the total weight of single excitations. However, the challenge in using %T₁ values is that they are only defined within a given computational method, and it is unclear, for example, how to compare results from single- and multi-reference computations. Therefore, we choose a different route here. We start with a method-independent definition of doubly excited states based only on physical observables without any explicit reference to orbitals or wave functions. Subsequently,

<figure><img src="image_3.png" /><figcaption>Fig. 1 Structures of the molecules studied within this work: the series of all-*trans*-polyenes up to 6 alternating double bounds, ethylene dimer and formaldehyde dimer.</figcaption></figure>

<figure><img src="image_4.png" /></figure>
we use this starting point to derive the squared 1TDM norm $\Omega$ as a rigorous and method-independent substitute for %T$_1$.

Within the 1TDM picture, we define a state as being doubly (or higher) excited via the condition that it is impossible to couple it to the ground state with any conceivable one-electron operator. We can turn this definition into a practical rule if we first realize that an arbitrary transition property of a one-electron operator between wave functions $\Psi_i$ and $\Psi_f$ is given as

$$a_{if} = \sum_{pq} D_{pq}^{if} A_{pq} \qquad (1)$$

where $D_{pq}^{if}$ is the 1TDM, $A_{pq}$ is the matrix representation of the operator, and both are given for a molecular orbital (MO) basis $\{\phi_p\}$. The 1TDM, in turn, is defined as

$$D_{pq}^{if} = \langle \Psi_i | p^\dagger q | \Psi_f \rangle \qquad (2)$$

where $p^\dagger$ and $q$ are the creation and annihilation operators related to the MOs $\phi_p$ and $\phi_q$. Note that $\Psi_i$ is allowed to be a general, correlated wave function here, meaning that separating these MOs into occupied and virtual orbitals is impossible. Applying the Cauchy-Schwarz inequality to eqn (1),³⁴ we obtain

$$|a_{if}|^2 \le \sum_{pq} |D_{pq}^{if}|^2 \times \sum_{pq} |A_{pq}|^2 = \Omega \times \sum_{pq} |A_{pq}|^2 \qquad (3)$$

where the symbol $\Omega$ denotes the squared Frobenius norm of the 1TDM, that is,

$$\Omega = \sum_{pq} |D_{pq}^{if}|^2 = \|\mathbf{D}^{if}\|^2 \qquad (4)$$

Importantly, we find that the transition property $a_{if}$ necessarily vanishes if $\Omega$ vanishes, that is, if all elements of the 1TDM are zero. Conversely, if any 1TDM element is non-zero, there is at least a conceivable one-electron operator with a non-vanishing transition property. An $\Omega$ value of zero is equivalent to the statement that the state cannot be coupled via a one-electron operator. Thus, a doubly (or higher) excited state exhibits $\Omega = 0$, whereas a purely singly excited state exhibits $\Omega = 1$.

More generally, $\Omega$ can be seen as an effective proportionality factor stating how strongly the transition interacts with one-electron operators. Therefore, a value of $\Omega$ between 0 and 1 can be used to represent a partially doubly excited character.³¹ In practice, the $\Omega$ value is consistent with the fraction of singly excited amplitudes (%T$_1$) alluded to above³⁵ and presents a natural generalization of this concept. An alternative viewpoint, based on ref. 85, is presented in Section S1.† In Section S2 we discuss the possibility of $\Omega$ values larger than 1.†

## 2.2 Definition of doubly excited character via (difference) density matrices

As an alternative to the 1TDM, it is possible to view doubly excited character via the 1-particle state density matrix (1DM) or difference density matrix (1DDM). The 1DM is defined in analogy to eqn (2) as

$$D_{pq}^{ff} = \langle \Psi_f | p^\dagger q | \Psi_f \rangle \qquad (5)$$

Furthermore, the difference density matrix (1DDM) is simply the difference between two state 1DMs

$$\Delta^{if} = \mathbf{D}^{ff} - \mathbf{D}^{ii} \qquad (6)$$

Diagonalization of the 1DDM and separation of the eigenvectors according to their signs provide detachment and attachment densities.³⁵,³⁶ The sum over all positive or negative eigenvalues of the 1DDM—denoted as the promotion number $p$—gives the total number of electrons rearranged during the excitation process. In principle, $p$ could lie between zero and the total number of electrons, but in practice, it usually ranges from 1 to 2. The promotion number appears to be a natural measure for defining a multiply excited character. However, since it is not only affected by the electrons taking part in the primary excitation process but also by secondary orbital relaxation,³⁵,³⁹ $p$ has been considered an unsuitable measure.¹⁴ More generally, doubly excited states usually experience increased $p$ values, but increased $p$ values alone are not a sufficient criterion to assign doubly excited character.

As an alternative measure for double excitation character, the excitation number ($\eta$) was introduced by Barca et al.¹⁴ in the context of the maximum overlap method (MOM). For two single-determinantal wave functions, $\Phi_i$ and $\Phi_f$, $\eta$ is defined as

$$\eta = n - \sum_{jk}^{\text{occ}} |S_{jk}^{if}|^2 \qquad (7)$$

where $S_{jk}^{if}$ is the overlap between the j-th occupied orbitals of $\Phi_i$ and the k-th occupied orbital of $\Phi_f$, $n$ is the total number of electrons in the system. The $\eta$ value ranges from 0 (when $\Phi_i$ is equal to $\Phi_f$) to $n$ (when there is no overlap between any orbitals in $\Phi_i$ and $\Phi_f$). A generalization of $\eta$ to arbitrary wave functions is not trivial, and we discuss this issue in some detail in Section S3.† After several attempts, we suggest using the formula

$$\eta = n_{\text{eff}} - \text{tr}(\mathbf{D}^{ii}\mathbf{D}^{ff})/2 \qquad (8)$$

where $\mathbf{D}^{ii}$ and $\mathbf{D}^{ff}$ are the spin-traced 1DMs. The value of $n_{\text{eff}}$ in this equation is defined as

$$n_{\text{eff}} = \max(\|\mathbf{D}^{ii}\|^2, \|\mathbf{D}^{ff}\|^2)/2 \qquad (9)$$

This expression reduces to eqn (7) for a single-determinantal wave function with doubly occupied spatial orbitals. Furthermore, just like the original expression, it is invariant to a switch between the initial and final states and vanishes if the initial and final states are the same. Moreover, we have verified that this expression yields the expected result in model systems and a variety of realistic computations.

Finally, we want to point out that $p$ and $\eta$ are both based on the 1DMs, and thus, they analyze shifts in the electron density rather than probing the actual wave functions. Crucially, if the initial and final state should have the same 1DMs, then $p$ and $\eta$ would both vanish, formally classifying the state as a “zero-electron transition”. This complicates the assignment in cases
of strong ground-state correlation. For example, states where the HOMO and LUMO are singly occupied (HOMO¹LUMO¹) or states constructed as a linear combination of the configurations with doubly occupied HOMO or doubly occupied LUMO (HOMO²-LUMO², HOMO²+LUMO²), all possess the same 1DMs, and, thus, the *p* and *η* values between them would vanish.

## 2.3 Further descriptors

If the state is predominantly singly excited ($\Omega \approx 1$), then it is meaningful to analyze the 1TDM further and to obtain the natural transition orbitals (NTOs), defined as the singular vectors of the 1TDM.⁴⁰,⁴¹ The number of independent NTO pairs necessary to describe the transition, called the NTO participation ratio,⁴²,⁴³ is defined as

$$PR_{\text{NTO}} = \frac{\left(\sum_{i} \lambda_{i}\right)^{2}}{\sum_{i} \lambda_{i}^{2}} \qquad (10)$$

where $\lambda_i$ are the weights of the NTO pairs. PRNTO can be used to assess the multiconfigurational character of the transition. A value of 1 means that a single configuration state function can express the excitation; higher values imply that this state has a multiconfigurational nature.

Furthermore, we compute the spin-averaged natural orbitals (NOs), defined as the eigenvectors of the 1DM. Their occupation numbers ($n_p$) range from 0 (unoccupied) to 2 (doubly occupied). It is common to characterize the 1DM via $y_0$ and $y_1$, which correspond to the occupation numbers of the lowest unoccupied natural orbitals (LUNO) and LUNO+1, respectively. For example, $(y_0, y_1)$ equals (0,0), (1,0), and (1,1) correspond to a closed shell, pure diradical, and pure tetra-radical characters, respectively.³⁸,⁴⁴ In Section 2.4, we show how these quantities can differentiate between types of doubly excited states. Alternatively, one can also compute the number of effectively unpaired electrons by summing over all NOs of the system. Eqn (12) and (13) show two expressions to obtain the number of unpaired electrons via either

$$n_u = \sum_{p} \min(n_p, 2 - n_p) \qquad (11)$$

or

$$n_{u,nl} = \sum_{p} n_{p}^{2} (2 - n_{p})^{2} \qquad (12)$$

where $n_u$ includes both static and dynamic correlation, while $n_{u,nl}$ suppresses dynamic correlation, thus, focusing on static contributions.³¹,³⁷ In principle, $n_u$ and $n_{u,nl}$ range from zero to the number of electrons. A value of zero represents a closed shell; a value of two a biradical with two open-shell orbitals, four represents a tetra-radical with four open shells, etc.

## 2.4 Classification of singly and doubly excited states

The descriptors presented above provide a toolbox for a comprehensive description of electronic excitation processes. This section shows how they can be combined to give a well-

<figure><img src="image_5.png" /><figcaption>Fig. 2 Limiting cases for singly and doubly excited states using a four-orbital four-electron model. Singly excited states distinguish between single configurational ($S_{SC}$) and multiconfigurational ($S_{MC}$) cases; doubly excited states distinguish between the formal closed-shell ($D_{CS}$) and four open-shells ($D_{OS}$) cases.</figcaption></figure>

defined classification scheme of singly and doubly excited states. For this purpose, we discuss the values of the descriptors for four limiting cases: a single configurational ($S_{SC}$) and a multiconfigurational ($S_{MC}$) singly excited state, and a closed-shell ($D_{CS}$) and an open-shell ($D_{OS}$) doubly excited state (Fig. 2). In addition, we will consider the mixed case ($D_{mix}$) as a case with notable doubly excited character not conforming with any of the four limiting cases.

Herein we use $\Omega$ as the main characteristic to distinguish between singly ($\Omega = 1$) and doubly ($\Omega = 0$) excited states. Alternatively, *p* and *η* can be used where their values correspond to the number of excited electrons. Note, however, that *p* is also strongly affected by orbital relaxation,³⁵,³⁹ and it is unclear how *η* performs in the case of static ground-state correlation.

Within the singly excited states, we distinguish between the single-configurational ($S_{SC}$) and multiconfigurational ($S_{MC}$) limiting cases. In the single-configurational limiting case ($S_{SC}$), the overall excitation can be described as a transition between a single pair of orbitals, e.g., the HOMO → LUMO transition. More generally, we define the $S_{SC}$ case as a state with only a single contributing NTO pair, leading to a value of PRNTO = 1. For the $S_{SC}$ case, $(y_0, y_1)$ is equal to (1, 0) since only one virtual orbital is involved. The multiconfigurational case ($S_{MC}$) is obtained if at least one additional pair of NTOs contributes to the state. In the scheme presented in Fig. 2, $S_{MC}$ is represented by PRNTO = 2 and $(y_0, y_1) = (0.5, 0.5)$. The significance of these differences is discussed in the literature, for example, in the context of excitons and the $L_a/L_b$ states in aromatic molecules.⁴³,⁴⁵⁻⁴⁸

Within the doubly excited states, we distinguish between the closed-shell ($D_{CS}$) and open-shell ($D_{OS}$) limiting cases. In the first case, exemplified by a pure HOMO² → LUMO² transition, two electrons are promoted to the same virtual orbital, and the excited state obtains a closed-shell character ($D_{CS}$) similar to the

<figure><img src="image_6.png" /></figure>
ground state. The second is exemplified by a combined HOMO/ HOMO-1 → LUMO/LUMO+1 transition: two electrons are promoted from two different initial orbitals to two different final orbitals, leaving four open-shell orbitals in total. D<sub>CS</sub> states can only be realized for singlets due to the Pauli principle, whereas D<sub>OS</sub> states can be singlet, triplet, or quintet. Within the presented scheme, the D<sub>CS</sub> and D<sub>OS</sub> states are resolved via NO occupations. The D<sub>CS</sub> limiting case possesses only closed shells and therefore has n<sub>u,nl</sub> = 0 and (y<sub>0</sub>, y<sub>1</sub>) = (0, 0), in analogy to a closed-shell ground state. Here, the LUMO of the ground state becomes a strongly occupied MO of the excited state. The four open-shell orbitals in the D<sub>OS</sub> case, on the other hand, are represented by n<sub>u,nl</sub> = 4 and (y<sub>0</sub>, y<sub>1</sub>) = (1, 1).

The idealized D<sub>CS</sub> state is single-configurational and behaves like a closed-shell ground state. Such a state would be readily described by a single determinant and would be particularly amenable to the maximum overlap method (MOM).<sup>49</sup> Furthermore, a CAS(2,2) active space or a single spin flip from a triplet reference would both suffice to describe such a state. On the other hand, a D<sub>OS</sub> case would always require a more sophisticated treatment, including at least four active orbitals. Similarly, the D<sub>mix</sub> case requires several correlated orbitals to describe its multiconfigurational character (unless the required nondynamic correlation effects can be captured within the exchange-correlation functional employed). We shall explore these issues below in Section 4.2.

## 2.5 Energies of doubly excited states of monomers and dimers

It is instructive to start by presenting the energies of the various states that can be constructed within a two-orbital two-electron model (TOTEM), as shown in Fig. 3 (see also ref. 8, 16 and 50). For simplicity, we consider transitions from HOMO to LUMO. Four spin-adapted wave functions can be constructed within the TOTEM: the ground state (¹G) with a doubly occupied HOMO, the single configurational singly excited states (¹S<sub>SC</sub>/³S<sub>SC</sub>) of singlet and triplet multiplicities, and the closed-shell doubly excited state (¹D<sub>CS</sub>). The relevant energy terms are the one-electron energies of HOMO and LUMO (h<sub>H</sub>, h<sub>L</sub>), the three Coulomb integrals (J<sub>HH</sub>, J<sub>HL</sub>, J<sub>LL</sub>), and the exchange integral (K<sub>HL</sub>). Fig. 3 indicates the different energy terms contributing to the energy. The energy of the ground state (¹G) is determined by the one-electron energy of the HOMO (h<sub>H</sub>) and the Coulomb

<figure><img src="image_7.png" /><figcaption>Fig. 3 Excited-state diagrams constructed within a two-orbital two-electron model. Energy terms are shown in red: one-electron energies of HOMO (h<sub>H</sub>) and LUMO (h<sub>L</sub>), the three Coulomb terms (J<sub>HH</sub>, J<sub>HL</sub>, J<sub>LL</sub>), and the exchange term (K<sub>HL</sub>).</figcaption></figure>

integral between the two electrons located in the HOMO; the energy of the doubly excited state (¹D<sub>CS</sub>) is determined in a completely analogous way only that the HOMO is swapped for the LUMO. Determining the energies of the singlet and triplet S<sub>SC</sub> states is slightly more involved since they are described by two interacting configurations, meaning that an exchange term (K<sub>HL</sub>) also comes into play. In summary, the energies of the states are given as

$$E(^1G) = 2h_H + J_{HH} \quad (13)$$

$$E(^1S_{SC}) = h_H + h_L + J_{HL} + K_{HL} \quad (14)$$

$$E(^3S_{SC}) = h_H + h_L + J_{HL} - K_{HL} \quad (15)$$

$$E(^1D_{CS}) = 2h_L + J_{LL} \quad (16)$$

Before continuing, we note that the one-electron energies used above include the kinetic energy, the nucleus-electron attraction and possibly the interaction with any other electrons present (treated in the sense of a frozen core) but do not consider any terms involving the HOMO and LUMO. An alternative and equally valid viewpoint is provided in ref. 8 by using the orbital energies (ε<sub>H</sub> = h<sub>H</sub> + J<sub>HH</sub>, ε<sub>L</sub> = h<sub>L</sub> + 2J<sub>HL</sub> - K<sub>HL</sub>) that already include interactions within the HOMO and the LUMO.

Using the TOTEM, we can now examine under what circumstances a closed-shell doubly excited state can be of lower energy than a singly excited state. Solving for E(^1D<sub>CS</sub>) < E(^1S<sub>SC</sub>) with the definitions given above, we obtain

$$h_L - h_H < K_{HL} + J_{HL} - J_{LL} \approx K_{HL} \quad (17)$$

where the right-hand side was simplified under the assumption that the J<sub>HL</sub> and J<sub>LL</sub> Coulomb integrals are of similar magnitude.<sup>8</sup> In other words, the doubly excited state becomes favorable if the exchange repulsion is large compared to the difference in one-electron energies between HOMO and LUMO.

Note, however, that eqn (17) and the assumption that the Coulomb integrals are of similar magnitude would also imply that E(^3S<sub>SC</sub>) < E(^1G), i.e. that the triplet lies below the closed shell.<sup>8</sup> K<sub>HL</sub> also couples ¹G and ¹D<sub>CS</sub>, meaning that the singlet ground state would obtain multiconfigurational character if the exchange interaction were, indeed, of the same order of magnitude as the difference in one-electron energies. This discussion highlights that a large exchange interaction relatively favors doubly excited states by pushing ¹S<sub>SC</sub> up in energy. Nevertheless, it also shows that a simple HOMO<sup>2</sup> → LUMO<sup>2</sup> transition cannot be the lowest excited state if the ground state is a closed shell. Indeed, the connection between doubly excited character and static electron correlation in the ground state is emphasized in the literature.<sup>13,14</sup>

Noting that a double excitation from HOMO to LUMO is not feasible for low-energy excited states, we proceed to an alternative type of doubly excited state. This alternative is present in the case of a dimer where two locally excited states can be combined into one doubly excited state of D<sub>OS</sub> type. Such states can be classified according to the spin-multiplicity of the overall
state along with the individual transitions. By combining singlet or triplet states on each monomer and considering all possible spin couplings, one obtains the states ¹(TT), ³(TT),⁵(TT), ³(ST),³(TS), and ¹(SS). Here, the ¹(TT) state is crucial for the singlet fission process, where one high-energy singlet excited state can be converted into two low-energy triplets.³ (See also ref. 52 for a discussion of the ¹(TT) state in the context of spin-exchange internal conversion and ref. 1 for a discussion of intramolecular doubly excited states in singlet fission.) We can use the TOTEM model to evaluate the stability of such a state. If we neglect possible biexciton binding effects, the ¹(TT) state will be the lowest state of singlet multiplicity if twice the excitation energy of ³S<sub>SC</sub> is lower than the excitation energy of ¹S<sub>SC</sub>. More specifically,

$$h_L - h_H < 3K_{HL} + J_{HH} - J_{HL} \approx 3K_{HL} \quad (18)$$

This condition is certainly easier to satisfy than eqn (17). In a dimer, the exchange repulsion associated with ¹S<sub>SC</sub> is avoided, and two exchange integrals are gained due to the two ³S<sub>SC</sub> configurations present.

As discussed below, we found that low-energy D<sub>OS</sub> type states can be readily constructed in dimers. Conversely, we were unable to find any low-lying D<sub>CS</sub>-type states in a variety of investigated molecules. The D<sub>CS</sub> states that were indeed found were of σσ* or nπ* character and trivially lay at about twice the energy of the corresponding singly excited state. Importantly, the low-lying A<sub>g</sub> ππ* states of polyenes and related systems do not fit either the D<sub>OS</sub> or D<sub>CS</sub> limitting case. Therefore, we introduce a third class of excited state, D<sub>mix</sub>, which is characterized by appreciable double excitation character (as determined by Ω and η) but not conforming to either of the limitting cases. Reviewing the TOTEM, we note that it is a simplified model capturing only the D<sub>OS</sub> and D<sub>CS</sub> cases but is unable to account for the energies of D<sub>mix</sub> type states. Indeed, a more involved model combining a triplet-pair state and charge-transfer exciton has been suggested for the latter.⁵¹

### 3. Computational details

Ground state geometry optimizations and vibrational frequencies of ethylene, butadiene, hexatriene, octatetraene, decapentaene, dodecahexatriene, s-tetrazine and the tetracene dimer were obtained at density functional theory (DFT) level with CAM-B3LYP functional,⁵³ cc-pVTZ⁵⁴ basis set, and Grimme's D3 dispersion correction.⁵⁵ DFT calculations were done using Gaussian 16 rev a03.⁵⁶ Vertical excitation energies and wave functions obtained at DFT/MRCI level⁵⁷ employed the def2-TZVP⁵⁸ basis set, except for the tetracene dimer where we use def2-SV(P).⁵⁸ In this approach, the CI expansion is built from Kohn-Sham orbitals using the BH-LYP⁵⁹ functional (as implemented in TURBOMOLE 7.5⁶⁰) and an effective Hamiltonian. Here, we tested two different parametrizations: the original one proposed by Grimme and Waletzke⁵⁷ and R2018 proposed by Marian et al.⁶¹ As shown in Fig. S1 of the ESI,† R2018 does not reproduce the inversion between states ¹B<sub>u</sub> and ²A<sub>g</sub> observed in polyenes with the increase of the number of double bonds in

the system,⁶² while the original parametrization does.⁶³ However, both parametrizations yield similar trends regarding the wave function analysis (Fig. S2†). In this work, we chose to use the original set of parameters. The initial reference space included configurations obtained from single and double excitations of 10 electrons within 10 orbitals.

TDDFT calculations were performed with the BLYP functional⁶⁴,⁶⁵ and 6-31G** basis set⁶⁶ with Q-Chem 5.3.⁶⁷ Multi-reference configuration interaction with single and double excitations (MR-CISD) were carried out using COLUMBUS 7.0.⁶⁸⁻⁷⁰ A complete active reference space (CAS) including all π and π* orbitals was used for polyenes up to four double bonds; polyenes with five and six double bonds were restricted to a CAS(8,8) due to computational cost. A complete active space self-consistent field (CASSCF) with these active spaces was used to construct the orbitals using state-averaging over the first two A<sub>g</sub> and the first B<sub>u</sub> states. MRCI energies are reported using the Pople extensivity correction (+P); 1s orbitals of all carbon atoms were frozen.⁷¹ Vertical excitation energies were also computed at the third-order algebraic diagrammatic construction method (ADC(3)) for the polarization propagator²⁰,⁷² level with the resolution-of-identity approximation and def2-SV(P) basis set, as implemented in Q-Chem.

In ADC(3) calculations, the wave function analysis library (libwfa)³¹,³² was used to obtain Ω-values, participation ratio of the natural transition orbitals (PR<sub>NTO</sub>), occupation of natural orbitals (y<sub>0</sub> and y<sub>1</sub>), number of unpaired electrons (n<sub>u,nl</sub>), and promotion number (p). η at ADC(3) and all descriptors at DFT/MRCI and MRCI levels were obtained externally with TheODORE⁷³ analysis package using a pre-release of version 3.0.

To investigate the cycloaddition of ethylene, we performed a relaxed scan keeping the linear combination of C1–C3 and C2–C4 (Fig. 1) at fixed distances. Ground state geometry optimizations were done at DFT level using B3LYP/cc-pVTZ. Excited-state energies and wave functions were obtained at MRCI+P(8,8)/cc-pVDZ level. All calculations considered D<sub>2h</sub> symmetry. CAS(8,8) wave functions with 6 states in the average (four A<sub>g</sub> and two B<sub>u</sub>) were used as references. MRCI calculations considered 4 frozen core orbitals (belonging to a<sub>g</sub>, b<sub>3u</sub>, b<sub>2u</sub>, and b<sub>1g</sub> representations), and 1 orbital in the active space for each of the 8 irreducible representations of the D<sub>2h</sub> point group.

### 4. Results and discussion

#### 4.1 Formaldehyde dimer

First, we discuss the ten lowest excited states of the formaldehyde dimer at 5 Å intermolecular distance calculated at the ADC(3) level (Table 1). The first six excited states are excitonic combinations of the ³nπ*, ¹nπ*, and ³ππ* states on the monomers. They are predominantly singly excited, as seen by their Ω values near or above 0.9 and η values close to 1.0. Their PR<sub>NTO</sub> values close to 2 (along with appropriate n<sub>u,nl</sub>, y<sub>0</sub>, and y<sub>1</sub> values) allow us to classify them as multiconfigurational (S<sub>MC</sub>) states, according to Fig. 2. The multiconfigurational nature of excitonic states delocalized between two interacting chromophores has been discussed in detail elsewhere,⁴³,⁴⁵ and we shall proceed with the remaining states here. The final four states are
<figure><img src="image_8.png" /></figure>

<table><caption>**Table 1** Excitation energy (ΔE in eV), oscillator strength (f), and wave function descriptors of the lowest excited singlet and triplet states of the formaldehyde dimer at 5 Å intermolecular separation at ADC(3) level</caption><thead><tr><th>State</th><th>Chars</th><th>ΔE (eV)</th><th>f</th><th>Ω</th><th>PR<sub>NTO</sub></th><th>n<sub>u,nl</sub></th><th>y<sub>0</sub></th><th>y<sub>1</sub></th><th>p</th><th>η</th></tr></thead><tbody><tr><td>1<sup>3</sup>A<sub>2</sub></td><td><sup>3</sup>nπ*</td><td>3.67</td><td>—</td><td>0.909</td><td>2.02</td><td>2.46</td><td>0.543</td><td>0.525</td><td>1.15</td><td>1.01</td></tr><tr><td>1<sup>3</sup>B<sub>2</sub></td><td><sup>3</sup>nπ*</td><td>3.67</td><td>—</td><td>0.909</td><td>2.02</td><td>2.46</td><td>0.543</td><td>0.526</td><td>1.15</td><td>1.01</td></tr><tr><td>1<sup>1</sup>A<sub>2</sub></td><td><sup>1</sup>nπ*</td><td>4.12</td><td>0.000</td><td>0.882</td><td>2.01</td><td>2.51</td><td>0.556</td><td>0.537</td><td>1.15</td><td>1.03</td></tr><tr><td>1<sup>1</sup>B<sub>2</sub></td><td><sup>1</sup>nπ*</td><td>4.12</td><td>0.000</td><td>0.882</td><td>2.01</td><td>2.51</td><td>0.555</td><td>0.537</td><td>1.15</td><td>1.03</td></tr><tr><td>1<sup>3</sup>B<sub>1</sub></td><td><sup>3</sup>ππ*</td><td>6.08</td><td>—</td><td>0.945</td><td>2.09</td><td>2.44</td><td>0.526</td><td>0.510</td><td>1.16</td><td>0.95</td></tr><tr><td>1<sup>3</sup>A<sub>1</sub></td><td><sup>3</sup>ππ*</td><td>6.08</td><td>—</td><td>0.945</td><td>2.09</td><td>2.44</td><td>0.523</td><td>0.513</td><td>1.16</td><td>0.95</td></tr><tr><td>2<sup>1</sup>A<sub>1</sub></td><td>nπ*:<sup>1</sup>(TT)</td><td>7.68</td><td>0.000</td><td>0.000</td><td>—</td><td>4.05</td><td>0.973</td><td>0.972</td><td>2.00</td><td>1.95</td></tr><tr><td>2<sup>3</sup>B<sub>1</sub></td><td>nπ*:<sup>3</sup>(TT)</td><td>7.68</td><td>—</td><td>0.000</td><td>—</td><td>4.05</td><td>0.973</td><td>0.972</td><td>2.00</td><td>1.95</td></tr><tr><td>3<sup>3</sup>B<sub>1</sub></td><td>nπ*:<sup>3</sup>(ST)</td><td>8.53</td><td>—</td><td>0.000</td><td>—</td><td>4.05</td><td>0.972</td><td>0.971</td><td>2.00</td><td>1.95</td></tr><tr><td>2<sup>3</sup>A<sub>1</sub></td><td>nπ*:<sup>3</sup>(ST)</td><td>8.53</td><td>—</td><td>0.000</td><td>—</td><td>4.05</td><td>0.972</td><td>0.972</td><td>2.00</td><td>1.95</td></tr></tbody></table>

unambiguously characterized as doubly excited due to their Ω values of exactly 0.000, along with η values near 2. More specifically, we find that these four states comply with the definition of an open-shell doubly excited state (D<sub>OS</sub>) in Fig. 2, and all relevant descriptors are within 0.05 of their idealized values. This discussion highlights that D<sub>OS</sub> states are indeed readily constructed in realistic systems and using a high-level method.

We find that the excitation energy of the 2<sup>1</sup>A<sub>1</sub> state (7.68 eV) is about twice the excitation energy of the singly excited 3nπ* state (3.67 eV), and we, therefore, assign it as the corresponding 1(TT) state. Note, however, that the agreement is not perfect, and the 1(TT) state lies about 0.35 eV higher than expected using the monomer energies. In principle, such a difference could derive from (bi)excitonic interaction effects. However, these are probably negligible at 5 Å, as the states come in pairs of the same energy. Conversely, we ascribe the difference to a lack of internal consistency within ADC(3) in terms of describing singly and doubly excited states at exactly the same level (see Section S6.1 in the ESI†). The next doubly excited state is of triplet multiplicity and almost degenerate with 1(TT); it is assigned as the 3(TT) state. The final two states shown in Table 1 combine singlet and triplet monomer excitations to form a 3(ST) state and are, again, almost degenerate. Their energy (8.53 eV) is significantly higher than the combined 3nπ* and 1nπ* energies (7.79 eV), which we again attribute to a lack of complete internal consistency within ADC(3).

## 4.2 Butadiene and larger polyenes

The photophysics of polyenes is usually discussed in terms of two important states close in energy: the 1<sup>1</sup>B<sub>u</sub> state dominated by the HOMO → LUMO transition and the 2<sup>1</sup>A<sub>g</sub> state with at least partial admixture of the doubly excited HOMO<sup>2</sup> → LUMO<sup>2</sup> configuration.<sup>9,10</sup> Within the following, we study these states in different polyenes using a variety of computational methods. We emphasize that our primary goal is not to obtain an accurate energy gap value between these states but rather to elucidate the nature of the 2<sup>1</sup>A<sub>g</sub> state.

The characterization and excited state ordering of polyenes is still a subject of discussion in the literature. Experimentally, 2<sup>1</sup>A<sub>g</sub> becomes the lowest excited state for polyenes with more than four double bonds.<sup>8,74</sup> Computationally, the energy gap

and order of states are strongly method-dependent.<sup>20,63,75,76</sup> While MS-CASPT2 calculations correctly predict the state inversion, CC3 predicts that 1<sup>1</sup>B<sub>u</sub> remains above 2<sup>1</sup>A<sub>g</sub>, although both methods deliver a small energy gap between those states.<sup>63</sup> ADC(2)-x and ADC(3) always predict 2<sup>1</sup>A<sub>g</sub> state as the lowest excited state, while ADC(2)-s predicts that to be the 1<sup>1</sup>B<sub>u</sub> state for polyenes up to four double bonds.<sup>20</sup> Benchmark studies show that the gap between 2<sup>1</sup>A<sub>g</sub> and 1<sup>1</sup>B<sub>u</sub> in polyenes is reproduced correctly by DFT/MRCI, although the excitation energies are lower than the best theoretical estimate.<sup>61,63,75</sup> Due to these discrepancies, here we use four different computational methods to study polyene excited states: ADC(3), MRCI, DFT/MRCI, and TDDFT/BLYP. For DFT/MRCI we use the original parameterisation because it reproduces the inversion between 1<sup>1</sup>B<sub>u</sub> and 2<sup>1</sup>A<sub>g</sub> expected in larger polyenes; see ESI Section S4.†

In particular, butadiene has become a paradigmatic case for discussing doubly excited states.<sup>8,13,14</sup> Thus, we first focus on this molecule using the ADC(3) method. To obtain a comprehensive picture, we look at several density descriptors (following Fig. 2) along with the percentage of single excitations (%T<sub>1</sub>) within ADC(3). The results are presented in Table 2. We start the discussion with the 1<sup>1</sup>B<sub>u</sub> state, which at this level of theory is the second excited state lying at 6.72 eV. Ω and η values close to 1 unambiguously assign the state as being singly excited. Furthermore, PR<sub>NTO</sub> approximately 1, along with (y<sub>0</sub>, y<sub>1</sub>) = (0.965, 0.039), allows classifying it as an S<sub>SC</sub> single configurational state close to the idealized values presented in Fig. 2.

The assignment of 2<sup>1</sup>A<sub>g</sub> of butadiene, which lies at 6.02 eV, is more involved. An Ω-value of 0.305, a %T<sub>1</sub> value of 31.3%, and an excitation number (η) of 1.537 indicate a partially doubly excited character, which is also supported by a promotion number (p) of 1.672. However, the descriptors are notably different from the idealized case of a doubly excited state (Ω = 0, η = p = 2) shown in Fig. 2. This implies that the admixture of singly excited configurations plays an important role. The dominant contribution to the 2<sup>1</sup>A<sub>g</sub> state is the HOMO<sup>2</sup> → LUMO<sup>2</sup> transition, with a weight of 42%. The HOMO-1 → LUMO (18%) and HOMO → LUMO+1 (13%) transitions come next, followed by many doubly excited configurations, all involving the HOMO-1, HOMO, LUMO, and LUMO+1. Thus, this state is strongly multiconfigurational, with notably different characteristics to the idealized D<sub>CS</sub> case. The NO-based
<table><caption>**Table 2** Excitation energy (ΔE in eV), oscillator strength (f), and wave function descriptors of the lowest excited singlet states of butadiene computed at the ADC(3) level</caption><thead><tr><th>State</th><th>ΔE (eV)</th><th>f</th><th>Ω</th><th>PR<sub>NTO</sub></th><th>n<sub>u,nl</sub></th><th>y<sub>0</sub></th><th>y<sub>1</sub></th><th>p</th><th>η</th><th>%T<sub>1</sub></th></tr></thead><tbody><tr><td>1<sup>1</sup>A<sub>g</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.186</td><td>0.071</td><td>0.053</td><td>—</td><td>—</td><td>—</td></tr><tr><td>2<sup>1</sup>A<sub>g</sub></td><td>6.02</td><td>0.000</td><td>0.305</td><td>1.959</td><td>2.430</td><td>0.890</td><td>0.253</td><td>1.672</td><td>1.537</td><td>31.3</td></tr><tr><td>1<sup>1</sup>B<sub>u</sub></td><td>6.72</td><td>1.739</td><td>0.904</td><td>1.095</td><td>2.073</td><td>0.965</td><td>0.039</td><td>0.998</td><td>0.899</td><td>93.5</td></tr></tbody></table>

characteristics reflect this divergence particularly well. For the idealized D<sub>CS</sub> case, n<sub>u,nl</sub>, y<sub>0</sub>, and y<sub>1</sub> are all zero. However, for the 2<sup>1</sup>A<sub>g</sub> state of butadiene, these values are n<sub>u,nl</sub> = 2.430, and (y<sub>0</sub>, y<sub>1</sub>) = (0.890, 0.253), which are between the limiting cases shown in Fig. 2. For these reasons, we classify this state as D<sub>mix</sub>.

As a next step, we investigate the dependence of the presented results on the electronic structure method and the size of the molecule. For this reason, we computed excitation energies and 1(T)DM descriptors at MRCI, DFT/MRCI, and ADC(3) levels for all-trans-butadiene, all-trans-hexatriene, all-trans-octatetraene, all-trans-decapentaene, and all-trans-dodecahexene (N = 2, 3, 4, 5, 6, respectively, where N is the number of double bonds in the system).

Starting with the excitation energies (Fig. 4A), we find a substantial decrease with increasing N for all methods. However, there is a notable difference in the 2<sup>1</sup>A<sub>g</sub> excitation energies according to the trend ADC(3) < TDDFT < DFT/MRCI<sub>MRCI+P</sub>, spanning a range of up to ≈1 eV. Concerning the Ω values presented in Fig. 4B, all methods aside from TDDFT find substantial double excitation character (Ω < 0.65) for all the molecules. The Ω values generally decrease with increasing N. The only exception is for the largest MRCI computations, which may be affected by size-extensivity problems. Aside from the general trends, there is also a quite notable difference in the obtained Ω values. ADC(3) indicates strong double-excitation character (Ω < 0.30 for N > 2); MRCI delivers intermediate values (0.35 < Ω < 0.5); DFT/MRCI predominantly indicates single excitation character (Ω > 0.5) albeit with notable admixtures of double excitations. TDDFT/BLYP, on the other hand, always predicts a singly excited character (Ω ≈ 1). The low Ω values for 2<sup>1</sup>A<sub>g</sub> contrast with the high Ω values for 1<sup>1</sup>B<sub>u</sub> (Fig. S9†), which are consistently above 0.90 for DFT/MRCI and TDDFT and above 0.85 for ADC(3) and MRCI. Thus, the Ω values clearly distinguish between the wave functions of the 2<sup>1</sup>A<sub>g</sub> and 1<sup>1</sup>B<sub>u</sub> states, assigning a partial doubly excited character to 2<sup>1</sup>A<sub>g</sub>.

As outlined in Section 2.1, Ω has a clear physical meaning by acting as an effective proportionality factor specifying how strongly the states are coupled via one-electron operators. Changes in Ω are related to changes in physically observable transition properties. In the present case, the transition dipole moments of the A<sub>g</sub> states vanish for symmetry reasons, but the differences should show up via enhanced transition quadrupole moments or angular momenta when computed with methods that produce different Ω values. However, this discussion is left to future work.

The excitation numbers (Fig. 5A) agree with the Ω-values in the sense that they always attribute at least partial doubly excited character to the 2<sup>1</sup>A<sub>g</sub> state (η > 1.2). Interestingly, the

MRCI and DFT/MRCI values are both fairly low, with η ≤ 1.3, whereas ADC(3) provides significantly enhanced doubly excited character (η > 1.5). One can understand this discrepancy by

<figure><img src="image_9.png" /><figcaption>**Fig. 4** Excitation energies ΔE (A) and 1TDM norm Ω (B) of the 2<sup>1</sup>Ag state of different polyenes plotted against the number of double bonds (N) calculated at ADC(3), DFT/MRCI (dubbed DFTCI in the insert), MRCI, and TDDFT/BLYP levels.</figcaption></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 5 Excitation number $\eta$ (A) and occupations (B) of the lowest unoccupied natural orbital (LUNO, $y_0$) and LUNO+1 ($y_1$) of $2^1A_g$ state of different polyenes plotted against the number of double bonds (N) calculated at ADC(3), DFT/MRCI (dubbed DFTCI in the insert), MRCI, and TDDFT/BLYP levels.</figcaption></figure>

noticing that these descriptors are differently affected by ground-state correlation (see Section 2.2). Generally speaking, it is not clear how the excitation number, initially developed for comparing single determinantal wave functions, should be interpreted in the case of multiconfigurational states. Nonetheless, it is interesting to discuss the excitation numbers of the $1^1B_u$ excited states (Fig. S10†). For MRCI, these are always below 0.85; for DFT/MRCI and ADC(3), values below 0.95 are obtained. Thus, a clear differentiation between the singly excited $1^1B_u$ and doubly excited $1^1A_g$ states is also present when the excitation numbers are considered. Finally, the $y_0$ values (Fig. 5B) are close

to 1 for all methods, whereas the $y_1$ values are smaller than 0.5. Note that these $y_0/y_1$ values are inconsistent with any limiting cases presented in Fig. 2. Conversely, they illustrate the multi-configurational and partial singly and doubly excited nature of the states, in line with a $D_{mix}$ character. This is, again, markedly different from the singly excited $1^1B_u$ states (Fig. S13†), which for all methods aside from TDDFT, exhibit the idealized values ($y_0 \approx 1, y_1 \approx 0$) expected for $S_{SC}$ states.

In summary, we find at least a partial double excitation character, classified as $D_{mix}$, in the $2^1A_g$ state of all polyenes investigated, revealed by all methods (except TDDFT) and descriptors in agreement with much of the previous literature.⁸⁻¹⁰ Nevertheless, two contradicting viewpoints argue against the double excitation character of butadiene in the literature, which we discuss next.

Shu and Truhlar¹³ have presented butadiene computations at various computational levels to understand the differences between the $2^1A_g$ and $1^1B_u$ states and learn which computational methods are suitable for their description. Crucially, they argued that doubly excited states could not be understood separately, but the multireference character of the ground state promotes the contribution of doubly excited configurations in low-lying excited states. This assessment agrees with the data presented in Table 2 and the discussion in Section 2.5. Nonetheless, we emphasise that the unique properties of polyenes cannot be understood by considering the ground state alone. Otherwise, we would observe similar amounts of double excitations in both $2^1A_g$ and $1^1B_u$. In contrast, the $2^1A_g$ and $1^1B_u$ states possess distinct characters: only $2^1A_g$ obtains double excitations, while $1^1B_u$ retains the singly excited ($S_{SC}$) character. Furthermore, Shu and Truhlar have advocated using local functionals to describe the excited states of polyenes.¹³ They argued that local functionals, such as BLYP, revPBE and M06-L, perform well since they minimize the static correlation error included by the Hartree-Fock exchange in nonlocal functionals.¹³,⁷⁷ Considering only the energies, one finds that TDDFT/BLYP does indeed produce results comparable to the wave-function-based methods (Fig. 4A). However, special care has to be taken due to the differences in $\Omega$-values (Fig. 4B), which would predict significantly altered transition properties.

Subsequently, Barca *et al.*¹⁴ performed DFT/MOM computations on butadiene and analysed them with the help of their original excitation number ($\eta$) definition, as shown in eqn (7). A value of $\eta = 1.022$ was obtained for $2^1A_g$ indicating almost perfect singly excited character. The challenge in interpreting these results is that the correlated computations describe the $2^1A_g$ state as a mixture of the HOMO–1/LUMO, HOMO/LUMO+1, and HOMO²/LUMO² configurations. However, the DFT/MOM method produces only a single open-shell Slater determinant. Thus, it is doubtful whether the MOM method provides a realistic description of this intrinsically multiconfigurational state. Conversely, extending the excitation number to the multi-configurational case (Fig. 5A) highlights the doubly excited character. Barca *et al.* argue that it is not clear a priori whether the doubly substituted determinants in CI wave functions account for electron excitation, electron correlation, or orbital relaxation, not allowing an unambiguous assignment.¹⁴ It is
<figure><img src="image_13.png" /></figure>

precisely for this reason that we have chosen a rigorously defined set of density-matrix-based descriptors to avoid such ambiguities. All the proposed descriptors are invariant to the orbital representation employed. This choice consistently shows the admixture of doubly excited characters in the 2¹Ag state of butadiene. It is important to point out that this contrasts with both Barca et al.'s assignment as a singly excited state, and Shu and Truhlar's argument that the apparent doubly excited character is due to correlation at the ground state.

## 4.3 Cycloaddition of ethylene

Unable to locate the D<sub>CS</sub> case in polyenes, we proceed to a third model, the cycloaddition of ethylene. The dimerization of ethylene to cyclobutane is a typical example illustrating the changes in the wave function character along a chemical reaction. According to the Woodward-Hoffmann rules, the [2 + 2] cycloaddition of ethylene is thermally forbidden due to a change in the ground state electronic configuration, yet it is photochemically allowed.<sup>78,79</sup> The reaction pathway can be explained by the four frontier orbitals shown in Fig. 6. The left represents the case of two isolated ethylene molecules where the frontier orbitals are of π and π* character. As they get closer, the spatial overlap between the orbitals increases until they eventually form the σ and σ* orbitals of cyclobutane shown on the right. Crucially, the HOMO (π<sub>2</sub>) of the isolated ethylene molecules corresponds to the LUMO (σ<sub>2</sub><sup>*</sup>) of cyclobutane and vice versa. Thus, the dimerization requires a change in electronic configuration from π<sub>1</sub><sup>2</sup>π<sub>2</sub><sup>2</sup>π<sub>3</sub><sup>*</sup>0 to σ<sub>1</sub><sup>2</sup>σ<sub>3</sub><sup>2</sup>σ<sub>2</sub><sup>*</sup>0. This process is thermochemically forbidden but can be facilitated via a doubly excited state. In this section, we assess the involvement of this doubly excited state, focusing on how its character changes along the reaction coordinate.

<figure><img src="image_14.png" /><figcaption>Fig. 6 Orbital correlation diagram for two ethylene molecules separated by 4 Å (left) and 1.4 Å (right).</figcaption></figure>

The ethylene dimerization is a practical model illustrating the transition from the D<sub>OS</sub> to D<sub>CS</sub> limits. The closed- or open-shell character is determined by an interplay between LUMO and LUMO+1 energies and the exchange integral, as outlined in Section 2.5. While the exchange term favors D<sub>OS</sub>, D<sub>CS</sub> is favored if LUMO and LUMO+1 are farther apart in energy. At large separations, HOMO (π<sub>2</sub>) and HOMO-1 (π<sub>1</sub>), as well as LUMO (π<sub>3</sub><sup>*</sup>) and LUMO+1 (π<sub>4</sub><sup>*</sup>), are pairwise degenerate. Thus, exchange dominates, and the two electrons are promoted to different unoccupied orbitals in a 1(TT) type state with a π<sub>1</sub><sup>1</sup>π<sub>2</sub><sup>1</sup>π<sub>3</sub><sup>*</sup>1π<sub>4</sub><sup>*</sup>1 configuration. As the distance decreases, the overlap between the ethylene orbitals increases, becoming non-degenerate. Then, the LUMO+1 becomes inaccessible, and a D<sub>CS</sub> state arises. Below, we evaluate the validity of this model using *ab initio* computations.

Generally speaking, three states are relevant in the dimerization process, 1<sup>1</sup>A<sub>g</sub>, 2<sup>1</sup>A<sub>g</sub>, and 1<sup>1</sup>B<sub>u</sub>. Here, the closed shell and doubly excited states—that is π<sub>2</sub><sup>2</sup>π<sub>2</sub><sup>2</sup>π<sub>3</sub><sup>*</sup>0/π<sub>1</sub><sup>2</sup>π<sub>0</sub><sup>0</sup>π<sub>3</sub><sup>*</sup>2/π<sub>1</sub><sup>1</sup>π<sub>2</sub><sup>1</sup>π<sub>3</sub><sup>*</sup>1π<sub>4</sub><sup>*</sup>1 for the dimer and σ<sub>1</sub><sup>2</sup>σ<sub>3</sub><sup>2</sup>σ<sub>2</sub><sup>*</sup>0/σ<sub>1</sub><sup>2</sup>σ<sub>3</sub><sup>0</sup>σ<sub>2</sub><sup>*</sup>2 for cyclobutane—are always of A<sub>g</sub> symmetry. The singly excited states π<sub>2</sub><sup>1</sup>π<sub>2</sub><sup>1</sup>π<sub>3</sub><sup>*</sup>1 and σ<sub>1</sub><sup>2</sup>σ<sub>3</sub><sup>1</sup>σ<sub>2</sub><sup>*</sup>1 are always B<sub>u</sub>. Fig. 7A presents the MRCI+P energies of these states computed along a relaxed scan. The right side shows the case of isolated ethylene molecules, whereas the left side represents the formation of cyclobutane. Starting with the 1<sup>1</sup>A<sub>g</sub> curve, we find that the dimerization is energetically favorable, but a substantial energy barrier of over 2.5 eV is encountered, making the reaction unfeasible in the ground state, as mentioned. Considering the excited states, the singly excited state (1<sup>1</sup>B<sub>u</sub>) is generally below any of the A<sub>g</sub> states, except in the avoided crossing region, where 2<sup>1</sup>A<sub>g</sub> becomes lower in energy. The excited states are fairly flat toward the right, whereas a steep increase of the doubly excited state is seen toward the left once cyclobutane is formed. Indeed, the doubly excited state of interest becomes 4<sup>1</sup>A<sub>g</sub> and reaches an adiabatic energy of 18.5 eV. Fig. 7A suggests a clear mechanism for the photochemical reaction: a photon is initially absorbed by 1<sup>1</sup>B<sub>u</sub>. Subsequently, the two molecules are attracted, forming an excimer. Near the minimum of the excimer, a crossing with 2<sup>1</sup>A<sub>g</sub> is encountered. Furthermore, 2<sup>1</sup>A<sub>g</sub> finally relaxes to the 1<sup>1</sup>A<sub>g</sub> ground state, forming cyclobutane.

To characterize the amount of double excitation character involved, we use the Ω descriptor, as presented in Fig. 7B. The Ω value for 1<sup>1</sup>B<sub>u</sub> is close to 1 throughout the energy profile, highlighting the singly excited nature of this state. Conversely, the Ω value for 2<sup>1</sup>A<sub>g</sub> is close to zero on the right side, highlighting that this state is doubly excited in the limit of the separated dimer. On the left side, the doubly excited character is transferred to the 4<sup>1</sup>A<sub>g</sub> state, as seen by its Ω value close to zero. A large spike in 2<sup>1</sup>A<sub>g</sub> is observed during the first avoided crossing around 2.5 Å. At this point, the 1<sup>1</sup>A<sub>g</sub> and 2<sup>1</sup>A<sub>g</sub> states become multiconfigurational, and a clear definition of doubly excited character becomes more challenging. Note that there is an additional spike in the region around 2.0 Å. It is related to a second avoided crossing involving 2<sup>1</sup>A<sub>g</sub>, 3<sup>1</sup>A<sub>g</sub>, and 4<sup>1</sup>A<sub>g</sub>.

The singly and doubly excited nature of the 1<sup>1</sup>B<sub>u</sub> and 2<sup>1</sup>A<sub>g</sub> states for reactant and product can also be represented by the η and p values (Fig. S16 and S17†). Toward the left and right, they
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 7 Relative energies in eV (A), 1TDM norms $\Omega$ (B), and number of unpaired electrons ($n_{u,nl}$; C) of selected singlet states of the ethylene dimer/cyclobutane system plotted against the intermolecular separation.</figcaption></figure>

represent the singly and doubly excited character, similar to $\Omega$ values. Interestingly, both $\eta$ and $p$ tend towards zero for $1^1B_u$ and $2^1A_g$ at the avoided crossing around 2.5 Å, which can be understood following the discussion in the last paragraph of Section 2.2. The $\pi_1^2\pi_2^2\pi_3^{*0}$ and $\pi_1^2\pi_2^0\pi_3^{*2}$ configurations mix and, as a consequence, all states have the same natural orbitals (with singly occupied $\pi_2$ and $\pi_3^*$ orbitals) and, hence, the same density matrix.

Having verified the overall amount of doubly excited character, we now use $n_{u,nl}$ to obtain a more detailed classification of the states (see Fig. 7C). As expected, $n_{u,nl}$ is close to two for the singly excited $1^1B_u$ state and, aside from the avoided crossing, near zero for the closed-shell $1^1A_g$ state. For the doubly excited $2^1A_g$ state, we find that its $n_{u,nl}$ value is close to 4 for large intermolecular separations, representing the idealized $D_{OS}$ ($1^1TT$, that is, $\pi_1^1\pi_2^1\pi_3^{*1}\pi_4^{*1}$) case. The value steadily decreases as the molecules move together (aside from the avoided crossings). Below the second avoided crossing, the doubly excited state becomes $4^1A_g$; its $n_{u,nl}$ value further decreases until reaching a value of 0.26 for the last point probed. This low $n_{u,nl}$

value, along with the individual ($y_0, y_1$) = (0.16, 0.08) values shown in Fig. S18–S20,† highlights that the state at this geometry does, indeed, closely conform with the idealized $D_{CS}$ state hypothesized in Fig. 2. Note, however, that this doubly excited state occurs at an extremely high vertical excitation energy of 18.5 eV. Indeed, its vertical excitation energy is about twice as high as the singly excited $1^1B_u$ state. Thus, in line with the above discussion (Section 2.5), we can state that a $D_{CS}$-type state can only occur at energies significantly higher than the lowest singly excited state.

## 4.4 Further molecules

Having outlined the different archetypes of doubly excited character in some detail above, we applied the proposed scheme to larger and more complex molecular systems. The selected molecules and a summary of the results are presented in Fig. 8, and more data is presented in the ESI.† We start with a diketopyrrolopyrrole derivative⁷⁶ used in optoelectronics,⁸⁰,⁸¹ OLEDs⁸² and as singlet-fission chromophores⁸³,⁸⁴ and
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 8 Analysis of the lowest singly and doubly excited singlet states of various molecular systems: excitation energies (ΔE in eV), squared 1TDM norms (Ω) and numbers of unpaired electrons (n<sub>u,nl</sub>).</figcaption></figure>

a derivative of its bis-thiophene building block (Fig. 8). A more extended set of derivatives is given in Table S4.† As pointed out previously,⁷⁶ these molecules have an extended polyene backbone with additional functional groups. The lowest doubly excited singlet states (of A<sub>g</sub> symmetry) in these molecules lie at similar energies as the lowest singly excited B<sub>u</sub> singlet states. In both cases shown and all examples in Table S4,† the A<sub>g</sub> states present the clear signature of the D<sub>mix</sub> case (0.3 < Ω < 0.6, 2.5 < n<sub>u,nl</sub> < 2.8). This analysis highlights that the D<sub>mix</sub> case applies to a variety of molecules and is an important model to understand electronic excitations.

Other interesting examples are the tetracene dimer and s-tetrazine. The tetracene dimer at 5 Å intermolecular separation is chosen as a more realistic illustration of dimer excited states relevant to singlet fission. Its lowest singlet excited state at 2.16 eV is a doubly excited ππ* state delocalized over both molecules. With values of Ω = 0.00 and n<sub>u,nl</sub> = 4.05 it almost perfectly aligns with the D<sub>OS</sub> case highlighting the relevance of this case for general dimer excited states. Finally, we investigate s-tetrazine. This molecule possesses a doubly excited nπ* state at 5.11 eV. With values of Ω = 0.00 and n<sub>u,nl</sub> = 0.14 it is a close match to the D<sub>CS</sub> case. Its energy is about twice the energy of S<sub>1</sub>, the singly excited nπ* state. Thus, in agreement with the previous discussions, we find that D<sub>CS</sub> type states are only found well above the first singly excited state.

## 5. Conclusions

This work presents a rigorous and transferable classification scheme for doubly excited states. We propose to define a doubly (or higher) excited state as a state that cannot be coupled to the ground state with any conceivable one-electron operator. This physically meaningful definition can be readily evaluated numerically using the 1TDM norm.

Within the manifold of doubly excited states, we define three cases: the closed-shell (D<sub>CS</sub>) case, where two electrons are promoted together from one orbital to another; the open-shell case (D<sub>OS</sub>), where the excitations occur between two independent orbital pairs; and the multiconfigurational mixed case (D<sub>mix</sub>) possessing only partial doubly excited character and conforming with neither of the above definitions. The underlying energetics are presented, highlighting that the D<sub>OS</sub> limiting case can occur as a low-lying excited state in realistic computations. Conversely, the pure D<sub>CS</sub> case is expected at

significantly higher energies, and only D<sub>mix</sub> is a viable model for low-lying intramolecular doubly excited states. The differentiation between D<sub>OS</sub>, D<sub>CS</sub>, and D<sub>mix</sub> in practical calculations using natural orbital occupations and other readily available density matrix descriptors is outlined. We find that the assignment of the D<sub>OS</sub> and D<sub>CS</sub> cases is generally clear and unambiguous. By contrast, D<sub>mix</sub> type states are sometimes discussed quite controversially in the literature. Other authors¹³,¹⁴ have refrained from using the term “doubly excited” for D<sub>mix</sub> type states reserving it for the D<sub>OS</sub> and D<sub>CS</sub> cases. Ultimately, this is a question of terminology. Nonetheless, this study shows that D<sub>mix</sub> states possess unique wavefunction properties that are clearly differentiated from traditional singly excited states and we, therefore, suggest labelling them as states with partial doubly excited character.

Several practical examples are presented to study the occurrence of different types of doubly excited states and their description with different computational methods. First, we highlight that D<sub>OS</sub>-type states can be readily constructed in dimers where they occur as different combinations of the monomer singlet and triplet states, such as ¹(TT), ³(TT), ³(ST). Computations of the formaldehyde dimer at the ADC(3) level are presented, illustrating that even at this highly correlated level, one obtains 1(T)DM descriptors close to the idealized results.

We proceed to butadiene and larger polyenes to present results on their controversially discussed 2¹A<sub>g</sub> excited states. Computations at the ADC(3), *ab initio* MRCI, and DFT/MRCI levels all agree that admixture of doubly excited character *via* the HOMO² → LUMO² transition plays an important role in the 2¹A<sub>g</sub> excited states of these molecules. At the same time, the description is never close to a D<sub>CS</sub> limiting case, but multi-configurational character and admixture of singly excited configurations play a significant role in line with the D<sub>mix</sub> case.

Unable to locate the D<sub>CS</sub> case in polyenes, we proceed to a third model, the [2 + 2] cycloaddition of ethylene. At large intermolecular separations, a doubly excited D<sub>OS</sub> (¹TT, π<sub>1</sub><sup>1</sup>π<sub>2</sub><sup>1</sup>π<sub>3</sub><sup>*1</sup>π<sub>4</sub><sup>*1</sup>) state is found. Upon dimerization, this state converts into a near-perfect D<sub>CS</sub> (σ<sub>1</sub><sup>2</sup>σ<sub>3</sub><sup>0</sup>σ<sub>2</sub><sup>*2</sup>) state that is strongly dominated by the HOMO² → LUMO² transition. However, this state lies at very high energies (above 18 eV), demonstrating, again, that a low-lying D<sub>CS</sub> state cannot be achieved. Finally, we applied our scheme to an extended set of molecules to highlight that the three archetypes D<sub>CS</sub>, D<sub>OS</sub>, and D<sub>mix</sub> are transferable
<figure><img src="image_19.png" /></figure>

models. Thus, we demonstrate the generality of our scheme for understanding doubly excited states in various systems.

In summary, we present a physically motivated definition of doubly excited character and a classification scheme able to distinguish between its limiting cases, providing a new approach to a long-standing problem. More specifically, we have highlighted challenges in the computational description of doubly excited states of various kinds, outlining the requirements for computational methods to describe them accurately. We emphasize that reproducing excitation energies is not enough to ensure the quality of a computational method to a determined system. A computation should only be deemed accurate if its wave functions and operator expectation values also comply with the reference. Nevertheless, the availability of well-defined and transferable descriptors can provide a solid basis for further discussions of the computational description of the doubly excited character. Furthermore, we hope the presented work can provide new ideas in the science surrounding doubly excited states and that the underlying physics discussed here can provide an improved language to discuss experimental results.

## Data availability

The underlying computational research data (molecular coordinates and input/output files for Q-Chem, Columbus, DFT/MRCI and TheoDORE) are available at Loughborough University's data repository at https://doi.org/10.17028/rd.lboro.22303765.

## Author contributions

Conceptualization, methodology, software: F. P.; data curation, visualization: M. T. do C.; funding acquisition and resources: M. T. do C. and M. B.; investigation, writing (original draft preparation): M. T. do C., F. P.; supervision: F. P., M. B., and J. M. T.; writing (review & editing): M. T. do C., F. P., J. M. T., M. B.

## Conflicts of interest

The authors declare no competing financial interest.

## Acknowledgements

F. P. thanks Ignacio Fernández Galván and Anna I. Krylov for discussions about wave function analysis methods. F. P. has received funding from the UK EPSRC (reference EP/V048686/1) and acknowledges access to Loughborough's high-performance computing system lovelace. M. T. do C., M. B., and J. M. T. have received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No 828753 (BoostCrop). M. B. and J. M. T. thank the European Research Council (ERC) advanced grant SubNano (Grant agreement 832237). M. T. do C. is thankful for the support of the PhD mobility program funded by AMUtech. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

## References

1. M. E. Sandoval-Salinas and D. Casanova, *ChemPhotoChem*, 2021, 5, 282-293.

2. T. Zeng, R. Hoffmann and N. Ananth, *J. Am. Chem. Soc.*, 2014, **136**, 5755-5764.

3. M. B. Smith and J. Michl, *Chem. Rev.*, 2010, **110**, 6891-6936.

4. K. W. Stone, D. B. Turner, K. Gundogdu, S. T. Cundiff and K. A. Nelson, *Acc. Chem. Res.*, 2009, **42**, 1452.

5. B. P. Fingerhut, M. Richter, J. W. Luo, A. Zunger and S. Mukamel, *Phys. Rev. B: Condens. Matter Mater. Phys.*, 2012, **86**, 1-10.

6. C. Heshmatpour, P. Malevich, F. Plasser, M. Menger, C. Lambert, F. Šanda and J. Hauer, *J. Phys. Chem. Lett.*, 2020, **11**, 7776-7781.

7. P. De Silva, *J. Phys. Chem. Lett.*, 2019, **10**, 5674-5679.

8. J. H. Starcke, M. Wormit, J. Schirmer and A. Dreuw, *Chem. Phys.*, 2006, **329**, 39-49.

9. B. Hudson and B. Kohler, *Annu. Rev. Phys. Chem.*, 1974, **25**, 437-460.

10. K. Schulten, I. Ohmine and M. Karplus, *J. Chem. Phys.*, 1976, **64**, 4422-4441.

11. H. A. Frank and R. J. Cogdell, *Photochem. Photobiol.*, 1996, **63**, 257-264.

12. X. Hu, A. Damjanović, T. Ritz and K. Schulten, *Proc. Natl. Acad. Sci. U. S. A.*, 1998, **95**, 5935-5941.

13. Y. Shu and D. G. Truhlar, *J. Am. Chem. Soc.*, 2017, **139**, 13770-13778.

14. G. M. J. Barca, A. T. B. Gilbert and P. M. W. Gill, *J. Chem. Theory Comput.*, 2018, **14**, 9-13.

15. P. F. Loos, M. Boggio-Pasqua, A. Scemama, M. Caffarel and D. Jacquemin, *J. Chem. Theory Comput.*, 2019, **15**, 1939-1956.

16. M. E. Casida and M. Huix-Rotllant, *Annu. Rev. Phys. Chem.*, 2012, **63**, 287-323.

17. A. Dreuw and M. Wormit, *WIREs Computational Molecular Science*, 2015, **5**, 82-95.

18. O. Christiansen, H. Koch and P. Jørgensen, *J. Chem. Phys.*, 1995, **103**, 7429-7441.

19. J. Leitner, A. L. Dempwolff and A. Dreuw, *J. Chem. Phys.*, 2022, **157**, 184101.

20. P. H. P. Harbach, M. Wormit and A. Dreuw, *J. Chem. Phys.*, 2014, **141**, 17.

21. P. F. Loos, A. Scemama, A. Blondel, Y. Garniron, M. Caffarel and D. Jacquemin, *J. Chem. Theory Comput.*, 2018, **14**, 4360-4379.

22. H. Lischka, D. Nachtigallová, A. J. A. Aquino, P. G. Szalay, F. Plasser, F. B. C. MacHado and M. Barbatti, *Chem. Rev.*, 2018, **118**, 7293-7361.

23. D. Casanova and A. I. Krylov, *Phys. Chem. Chem. Phys.*, 2020, **22**, 4326-4342.

24. Z. Rinkevicius, O. Vahtras and H. Ågren, *J. Chem. Phys.*, 2010, **133**, 114104.

25. G. M. J. Barca, A. T. B. Gilbert and P. M. W. Gill, *J. Chem. Theory Comput.*, 2018, **14**, 1501-1509.

26. D. Hait and M. Head-Gordon, *J. Phys. Chem. Lett.*, 2021, **12**, 4517-4529.
<figure><img src="image_20.png" /></figure>

27 C. D. Lin, *Nucl. Instrum. Methods Phys. Res.*, 1985, **240**, 572-576.

28 C. D. Lin, *Adv. Atom. Mol. Phys.*, 1986, **22**, 77-142.

29 N. T. Maitra, *Annu. Rev. Phys. Chem.*, 2022, **73**, 117-140.

30 J. P. Coe and M. J. Paterson, *Mol. Phys.*, 2014, **112**, 733-739.

31 F. Plasser, M. Wormit and A. Dreuw, *J. Chem. Phys.*, 2014, **141**, 024106.

32 F. Plasser, A. I. Krylov and A. Dreuw, *WIREs Computational Molecular Science*, 2022, 1-15.

33 A. I. Krylov, *J. Chem. Phys.*, 2020, **153**, 080901.

34 S. Matsika, X. Feng, A. V. Luzanov and A. I. Krylov, *J. Phys. Chem. A*, 2014, **118**, 11943-11955.

35 F. Plasser, S. A. Bäppler, M. Wormit and A. Dreuw, *J. Chem. Phys.*, 2014, **141**, 24107.

36 M. Head-Gordon, A. M. Grana, D. Maurice and C. A. White, *J. Phys. Chem.*, 1995, **99**, 14261-14270.

37 M. Head-Gordon, *Chem. Phys. Lett.*, 2003, **372**, 508-511.

38 T. Minami and M. Nakano, *J. Phys. Chem. Lett.*, 2012, **3**, 145-150.

39 F. Plasser, S. A. Mewes, A. Dreuw and L. González, *J. Chem. Theory Comput.*, 2017, **13**, 5343-5353.

40 A. V. Luzanov, A. A. Sukhorukov and V. É. Umanskii, *Theor. Exp. Chem.*, 1976, **10**, 354-361.

41 R. L. Martin, *J. Chem. Phys.*, 2003, **118**, 4775-4777.

42 F. Plasser and H. Lischka, *J. Chem. Theory Comput.*, 2012, **8**, 2777-2789.

43 F. Plasser, *J. Chem. Phys.*, 2016, **144**, 194107.

44 M. Nakano, H. Fukui, T. Minami, K. Yoneda, Y. Shigeta, R. Kishi, B. Champagne, E. Botek, T. Kubo, K. Ohta and K. Kamada, *Theor. Chem. Acc.*, 2011, **130**, 711-724.

45 I. Mayer, *Chem. Phys. Lett.*, 2007, **443**, 420-425.

46 F. Bettanin, L. F. A. Ferrão, M. Pinheiro, A. J. A. Aquino, H. Lischka, F. B. C. Machado and D. Nachtigallova, *J. Chem. Theory Comput.*, 2017, **13**, 4297-4306.

47 P. Kimber and F. Plasser, Classification and Analysis of Molecular Excited States, in *Comprehensive Computational Chemistry*, ed. R. Boyd and M. Yanez, 2023 (in print).

48 P. Kimber and F. Plasser, *Phys. Chem. Chem. Phys.*, 2020, **22**, 6058-6080.

49 A. T. B. Gilbert, N. A. Besley and P. M. W. Gill, *J. Phys. Chem. A*, 2008, **112**, 13164-13171.

50 X. K. Chen, D. Kim and J. L. Brédas, *Acc. Chem. Res.*, 2018, **51**, 2215-2224.

51 W. Barford, *Phys. Rev. B*, 2022, **106**, 035201.

52 S. Bai and M. Barbatti, *J. Chem. Theory Comput.*, 2019, **15**, 1503-1513.

53 T. Yanai, D. P. Tew and N. C. Handy, *Chem. Phys. Lett.*, 2004, **393**, 51-57.

54 T. H. Dunning, *J. Chem. Phys.*, 1989, **90**, 1007-1023.

55 S. Grimme, *J. Comput. Chem.*, 2006, **27**, 1787-1799.

56 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson,

D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell Jr, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, *Gaussian 16, Revision A. 03*, Gaussian, Inc., Wallingford, CT, 2016.

57 S. Grimme and M. Waletzke, *J. Chem. Phys.*, 1999, **111**, 5645-5655.

58 F. Weigend and R. Ahlrichs, *Phys. Chem. Chem. Phys.*, 2005, **7**, 3297-3305.

59 A. D. Becke, *J. Chem. Phys.*, 1993, **98**, 1372-1377.

60 S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Rückert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodyński and J. M. Yu, *J. Chem. Phys.*, 2020, **152**, 184107.

61 C. M. Marian, A. Heil and M. Kleinschmidt, *Wiley Interdiscip. Rev.: Comput. Mol. Sci.*, 2019, **9**, 1-31.

62 W. Fuß, Y. Haas and S. Zilberg, *Chem. Phys.*, 2000, **259**, 273-295.

63 M. R. Silva-Junior, M. Schreiber, S. P. A. Sauer and W. Thiel, *J. Chem. Phys.*, 2008, **129**, 104103.

64 A. D. Becke, *Phys. Rev. A*, 1988, **38**, 3098.

65 C. Lee, W. Yang and R. G. Parr, *Phys. Rev. B: Condens. Matter Mater. Phys.*, 1988, **37**, 785.

66 R. Ditchfield, W. J. Hehre and J. A. Pople, *J. Chem. Phys.*, 1971, **54**, 720-723.

67 Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kus, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C. M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. Distasio, H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T. C. Jagau, H. Ji, B. Kaduk, K. Khistaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser,
<figure><img src="image_21.png" /></figure>

P. Manohar, S. F. Manzer, S. P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. Oneill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y. C. Su, A. J. W. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z. Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J. Da Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C. P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill and M. Head-Gordon, *Mol. Phys.*, 2015, **113**, 184–215.

68 H. Lischka, T. Müller, P. G. Szalay, I. Shavitt, R. M. Pitzer and R. Shepard, *Wiley Interdiscip. Rev.: Comput. Mol. Sci.*, 2011, **1**, 191–199.

69 H. Lischka, R. Shepard, R. M. Pitzer, I. Shavitt, M. Dallos, T. Müller, P. G. Szalay, M. Seth, G. S. Kedziora, S. Yabushita and Z. Zhang, *Phys. Chem. Chem. Phys.*, 2001, **3**, 664–673.

70 H. Lischka, R. Shepard, I. Shavitt, R. M. Pitzer, M. Dallos, Th. Müller, P. G. Szalay, F. B. Brown, R. Ahlrichs, H. J. Böhm, A. Chang, D. C. Comeau, R. Gdanitz, H. Dachsel, C. Ehrhardt, M. Ernzerhof, P. Höchtl, S. Irle, G. Kedziora, T. Kovar, V. Parasuk, M. J. M. Pepper, P. Scharf, H. Schiffer, M. Schindler, M. Schüler, M. Seth, E. A. Stahlberg, J.-G. Zhao, S. Yabushita, Z. Zhang, M. Barbatti, S. Matsika, M. Schuurmann, D. R. Yarkony, S. R. Brozell, E. V. Beck, J.-P. Blaudeau, M. Ruckenbauer, B. Sellner, F. Plasser, J. J. Szymczak, R. F. K. Spada and A. Das, *COLUMBUS, an ab initio electronic structure program, release 7.0*, 2019.

71 J. A. Pople, R. Seeger and R. Krishnan, *Int. J. Quantum Chem.*, 1977, **12**, 149–163.

72 J. Schirmer, *Phys. Rev. A*, 1982, **26**, 2395–2416.

73 F. Plasser, *J. Chem. Phys.*, 2020, **152**, 084108.

74 R. J. Cave and E. R. Davidson, *Chem. Phys. Lett.*, 1988, **148**, 190–196.

75 C. M. Marian and N. Gilka, *J. Chem. Theory Comput.*, 2008, **4**, 1501–1515.

76 M. T. do Casal, J. M. Toldo, F. Plasser and M. Barbatti, *Phys. Chem. Chem. Phys.*, 2022, 23279–23288.

77 H. S. Yu, S. L. Li and D. G. Truhlar, *J. Chem. Phys.*, 2016, **145**, 130901.

78 F. Bernardi, M. Olivucci and M. A. Robb, *Acc. Chem. Res.*, 1990, **23**, 405–412.

79 N. D. Epiotis and S. Shaik, *J. Am. Chem. Soc.*, 1978, **100**, 9–17.

80 S. Ghosh, S. Shankar, D. S. Philips and A. Ajayaghosh, *Mater. Today Chem.*, 2020, **16**, 100242.

81 N. Luo, G. Zhang and Z. Liu, *Org. Chem. Front.*, 2021, **8**, 4560–4581.

82 A. Shukla, S. K. M. McGregor, R. Wawrzinek, S. Saggar, E. G. Moore, S. C. Lo and E. B. Namdas, *Adv. Funct. Mater.*, 2021, **31**, 1–11.

83 L. Shen, Z. Tang, X. Wang, H. Liu, Y. Chen and X. Li, *Phys. Chem. Chem. Phys.*, 2018, **20**, 22997–23006.

84 I. Papadopoulos, M. J. Álvaro-Martins, D. Molina, P. M. McCosker, P. A. Keller, T. Clark, Á. Sastre-Santos and D. M. Guldi, *Adv. Energy Mater.*, 2020, **10**, 1–10.

85 S. Klaiman and L. S. Cederbaum, *J. Chem. Phys.*, 2014, **141**, 194102.