Are Multicentre Bond Indices and Related Quantities Reliable Predictors of Excited-State Aromaticity?

Abstract

Systematic scrutiny is carried out of the ability of multicentre bond indices and the NOEL-based similarity index d_AB to serve as excited-state aromaticity criteria. These indices were calculated using state-optimized complete active-space self-consistent field wavefunctions for several low-lying singlet and triplet states of the paradigmatic molecules of benzene and square cyclobutadiene and the inorganic ring S₂N₂. The comparison of the excited-state indices with aromaticity trends for individual excited states suggested by the values of magnetic aromaticity criteria show that whereas the indices work well for aromaticity reversals between the ground singlet and first triplet electronic states, addressed by Baird’s rule, there are no straightforward parallels between the two sets of data for singlet excited states. The problems experienced while applying multicentre bond indices and d_AB to singlet excited states are explained by the loss of the information inherently present in wavefunctions and/or pair densities when calculating the first-order density matrix.

1. Introduction

Despite its somewhat vaguely defined qualitative nature, the concept of aromaticity has had huge impacts on organic chemistry, starting with the formulation of the Hückel aromaticity rules [1,2] and encompassing a broad research area including the elucidation of the link between cyclic delocalization and energetic stabilization of conjugated (poly)cyclic hydrocarbons [3,4,5,6,7,8,9,10], the role of cyclic conjugation in inducing the ring currents [11,12,13,14,15,16,17,18,19] responsible for the special magnetic properties of aromatic compounds, and revealing the links between electron counts, orbital topology and selection rules in pericyclic reactions [20,21,22,23]. The fact that the phenomenon of aromaticity can be associated with a very wide range of structural, energetic, and magnetic properties [3,4,5,6,7,8,9,10,11,12,13,14,16,17,18,24,25] has given impetus to numerous attempts to define measures or indices that are intended to characterize the “extent” of aromaticity in quantitative terms [16,17,18,19,24,25,26,27,28,29,30,31,32,33]. However, such efforts have often been plagued by discrepancies between the various types of indices; these have led to the postulation of a multidimensional character for this phenomenon [34,35,36], even implying “orthogonality” between energetic and magnetic measures of aromaticity as manifested by the reported absence of a straightforward link between these two types of aromaticity measure [19,29,37,38,39]. However, it has been demonstrated that such discrepancies are most often observed when trying to juxtapose quantities that are not straightforward to compare. One such example is provided by attempts to correlate the extents of cyclic delocalization in the individual benzene rings in polycyclic aromatic hydrocarbons (PAHs), as given by multicentre bond indices, with the values of nucleus-independent chemical shieldings (NICS) [18]. This fails not least because of the incompatibility between the strictly local character of multicentre indices and the fact that the NICS value for an individual ring is “contaminated” by the interfering contributions of the other rings [37,39]. The agreement between these two types of index is in fact restored when the contaminating contributions are properly taken into account; analogous parallels have been established between multicentre indices, induced ring currents and (when properly accounted for) the energetic effects of cyclic conjugation [19,37,38,39,40,41].

Although initially most studies were focused on aromaticity in ground electronic states, Baird’s pioneering discovery of the reversal of Hückel’s aromaticity rules upon electronic excitation from the singlet ground to the first triplet excited state [42] directed attention to the systematic investigation of excited-state aromaticity [43,44,45,46,47,48,49,50,51]. The importance of such studies for the understanding of the photochemical/physical properties of photoactive materials has prompted the development of experimental and computational tools that are capable of providing reliable estimates of excited-state aromaticity. Amongst the first attempts at theoretical justification of Baird’s discovery of aromaticity reversals in the lowest excited states of cyclic conjugated hydrocarbons is a study by Iljić et al. [43] which looked at the extension of the concept of topological resonance energy (TRE) to low-lying states of cyclic conjugated hydrocarbons. The authors of that study demonstrated that the TRE values for the ground and lowest excited states of conjugated rings reproduce the aromaticity reversal predicted by Baird’s rule. Despite the elegant simplicity of this approach, the calculation of TREs has serious inherent limitations arising from the Hückel molecular orbital (HMO) foundations of the underlying graph-theoretical considerations. Modern quantum chemical calculations are not subject to such limitations and the scope of excited-state aromaticity studies was subsequently extended to formulating Baird-style rules for higher excited states. The most convincing proof of aromaticity and/or antiaromaticity reversals in the first and higher excited states was provided by the results of systematic studies of various magnetic properties with state-specific complete active-space self-consistent field (CASSCF) wavefunctions constructed from gauge-included atomic orbitals (GIAOs) [44,45,46,52]. Given that multicentre bond indices have been applied successfully for the quantitative evaluation of the local aromaticities of individual benzene rings in PAHs [28,29,37,40], it was natural to try to find out whether the same approach could provide a computationally efficient and sufficiently accurate characterization of excited-state aromaticity. The aim of the current work is to carry out a systematic comparative study of the performance and reliability of multicentre bond indices and other first-order density-based quantities for the description and classification of excited-state aromaticity in the paradigmatic molecules benzene and cyclobutadiene, as well as in disulfur dinitride, which has been shown recently to be the first inorganic ring that exhibits changes in aromaticity between different electronic states [52]. As will be shown, it turns out that such quantities have significant difficulties distinguishing properly between singlet diradical and zwitterionic character.

2. Computational Methodology

2.1. Electronic Structure Calculations

The aromaticity of the low-lying electronic states of benzene, square cyclobutadiene and disulfur dinitride has been analysed using a range of magnetic criteria including NICS calculated with CASSCF-GIAO wavefunctions at fixed ground electronic state geometries [27,44,45,46,52]. To enable direct comparisons, we use the same levels of theory and the same geometries in the current work. All excited electronic state properties discussed in this work correspond to vertical excitations because we chose to use identical ground-state geometries for all electronic states of a given molecule. All CASSCF calculations on benzene and cyclobutadiene reported in this paper were carried out within the 6-311++G(2d,2p) basis, whereas use was made of the cc-pVTZ basis for S₂N₂.

It is important in this work to focus on vertical excitations not least because the electronic wavefunction changes much more rapidly than the molecular geometry. By examining the excited-state wavefunction at the ground-state geometry we can establish whether a given vertically excited state is intrinsically aromatic, antiaromatic or non-aromatic. If it turns out to be aromatic, then it will of course tend to retain a geometry that is similar to that of the ground state. If, on the other hand, it is antiaromatic then it is likely to experience a geometry distortion that leads to a lower energy, less antiaromatic and closer to non-aromatic geometry. The same does of course apply for systems which are antiaromatic in their electronic ground states, such as the ground state of square (D_4h) cyclobutadiene; the relaxation of the geometry of that system to rectangular (D_2h) decreases the antiaromaticity, not only making this state much more non-aromatic but also rendering it significantly less interesting to study as an example of an antiaromatic molecule.

The S₀ (1 ¹A_1g), T₁ (1 ³B_1u), S₁ (1 ¹B_2u) and S₂ (1 ¹B_1u) electronic states of benzene were described using state-optimized π-space CASSCF(6,6) wavefunctions (with “6 electrons in 6 orbitals”). We used the experimental D_6h gas-phase ground-state geometry with C−C and C−H bond lengths of 1.3964 Å and 1.0831 Å, respectively, which was obtained through analysis of the ν₄ vibration-rotation bands of C₆H₆ and C₆D₆ [53].

The calculations for the S₀ (1 ¹B_1g), T₁ (1 ³A_2g), S₁ (1 ¹A_1g) and S₂ (1 ¹B_2g) electronic states of square (D_4h) cyclobutadiene employed state-optimized π-space CASSCF(4,4) wavefunctions (with “4 electrons in 4 orbitals”). We used C−C and C−H bond lengths of 1.447 Å and 1.076 Å, respectively, that were optimized with the cc-pVTZ basis through a multireference averaged quadratic coupled cluster (MR-AQCC) approach, taking orbitals from corresponding state-averaged π-space CASSCF(4,4) wavefunctions that included the ground state, lowest triplet state and two lowest singlet excited states (SA-4-CASSCF) [54].

The calculations on the S₀ (1 ¹A_g), T₁ (1 ³B_3u) and S₁ (1 ¹A_u) electronic states of S₂N₂ were carried out using state-optimized CASSCF(22,16) wavefunctions (with “22 electrons in 16 orbitals”). For this purpose we used the D_2h semi-experimental equilibrium geometry established by Perrin et al. [55], with R(SN) = 1.64182 Å and ∠(NSN) = 91.0716°, in a coordinate system that places N at positions (±1.171748 Å, 0.0, 0.0), and S atoms at positions (0.0, ±1.150035 Å, 0.0), respectively.

All of the CASSCF calculations required for the present work were primarily carried out using Gaussian 03 [56] but, purely for our convenience, the same wavefunctions were also obtained using MOLPRO [57,58]. For reasons that we have explained, it was important to use the ground-state geometries for all of these calculations. We note in passing that accurate excited-state geometry optimizations of antiaromatic states would require methods such as CASPT2, given that those based on a closed-shell reference do not describe correctly the biradical character. Such studies are outside the scope of the present work but may be considered when CASPT2 analytical gradients and Hessians become widely available, making the optimization and characterization of excited-state local minima and saddle points very much faster and more reliable.

2.2. Multicentre Bond Indices

Such indices were originally introduced [59,60] as mono-, bi- tri- and generally k-centre contributions resulting at the closed-shell SCF or Kohn–Sham level of theory from the identity (1)

$$Tr{(PS)}^k=2^{\left(k-1\right)}N={\displaystyle \sum }_A{\mathsf{\Delta }}_A^{\left(k\right)}+{\displaystyle \sum }_{A<B}{\mathsf{\Delta }}_{AB}^{\left(k\right)}+{\displaystyle \sum }_{A<B<C}{\mathsf{\Delta }}_{ABC}^{\left(k\right)}+\dots {\displaystyle \sum }_{A<B<C\dots K}{\mathsf{\Delta }}_{ABC\dots K}^{\left(k\right)}$$

(1)

in which P and S denote the charge-density and overlap matrices, respectively, and N is the number of electrons. The usefulness of these indices for structural elucidations arises from the interesting nontrivial finding that their values mimic sensitively the presence and/or absence of bonding interactions between individual atoms in a molecule. Thus, for example, in the case of molecules that are well described by the familiar classical Lewis model of localized two-centre two-electron bonds, the corresponding 2-centre bond indices, which coincide in this case with the well-known Wiberg–Mayer indices [61,62], attain non-negligible values only between classically bonded atoms while the corresponding values of the indices for pairs of classically nonbonded atoms are negligible. Such indices are also very useful for molecules whose descriptions transcends the classical Lewis model by involving instead bonding interactions that are delocalized over more than two atomic centres. The indices retain the ability in such cases to detect and reveal just those atomic fragments engaged in the delocalized bonding, whereas the values for the remaining fragments are very small.

In various earlier papers the definition was modified [29,63,64] and the indices were instead normalized to N, as shown in Equation (2).

$$\frac{1}{2^{\left(k-1\right)}}Tr{(PS)}^k=N={\displaystyle \sum }_A{\mathsf{\Delta }}_A^{\left(k\right)}+{\displaystyle \sum }_{A<B}{\mathsf{\Delta }}_{AB}^{\left(k\right)}+{\displaystyle \sum }_{A<B<C}{\mathsf{\Delta }}_{ABC}^{\left(k\right)}+\dots {\displaystyle \sum }_{A<B<C\dots K}{\mathsf{\Delta }}_{ABC\dots K}^{\left(k\right)}$$

(2)

A slight disadvantage of this alternative definition is that the resulting values of the indices tend to decrease rather dramatically with increasing k. We have thus chosen for the present work to use the original definition (Equation (1)) in which the normalization sum is $Tr{(PS)}^k$ (which is equal to $2^{k-1}N$ at the closed-shell SCF level). Otherwise, even relevant indices would be rather small for k > 3. An obvious alternative to the different values returned by Equations (1) and (2) would be to quote instead the proportion of the quantity being “partitioned”.

In the case of k-centre indices the general definition (Equation (1)) leads to Equation (3)

$$MC{I}_{ABC\dots K}={\mathsf{\Delta }}_{ABC\dots K}^{\left(k\right)}={\displaystyle \sum }_i{\Gamma }_i\left[{\displaystyle \sum }_{\mu }^A{\displaystyle \sum }_{\nu }^B{\displaystyle \sum }_{\lambda }^C\dots {\displaystyle \sum }_{\tau }^K{(PS)}_{\mu \nu }{(PS)}_{\nu \lambda }{(PS)}_{\lambda \sigma }\dots {(PS)}_{\tau \mu }\right]$$

(3)

where the permutation operator ${\Gamma }_i$ ensures that the index includes all of the terms that correspond to different permutations of atomic labels.

The above general formula can also be straightforwardly extended beyond the scope of the Hartree-Fock approximation. The formula remains formally the same, except of course that the idempotent charge-density matrix is replaced by the corresponding correlated first-order density matrix [65,66]. The normalization sum $Tr{(PS)}^k$ is of course no longer straightforwardly linked to the total number of electrons, as in the case of the closed-shell SCF approximation. The above definitions of multicentre indices that are based on a Mulliken-like partitioning can easily be generalized to the framework of QTAIM [67] analysis [30,68], such that Equation (4) is transformed to:

$$MC{I}_{ABC\dots K}={\displaystyle \sum }_i{\Gamma }_i\left[{\displaystyle \sum }_{\alpha }{\displaystyle \sum }_{\beta }{\displaystyle \sum }_{\gamma }\dots \dots {\displaystyle \sum }_{\tau }{\eta }_{\alpha }{\eta }_{\beta }{\eta }_{\gamma }\dots {\eta }_{\tau }\langle \alpha |\beta {\rangle }_A\langle \beta |\gamma {\rangle }_B\langle \gamma |\delta {\rangle }_C\dots \langle \tau |\alpha {\rangle }_K\right]$$

(4)

Here the symbol $\langle \lambda |\sigma {\rangle }_X$ denotes the domain-condensed overlap of natural orbitals $\lambda$ and $\sigma$ (i.e., integration over the QTAIM atomic domain of atom $X$), ${\eta }_{\lambda }$ denotes the occupation number of $\lambda$, and the summations again run over all permutations of atomic labels.

Instead of using orbitals it is of course possible to calculate indices separately for α and β spin-orbitals. Such an approach was reported in earlier extensions of multicentre indices to open-shell systems [69]. The total index is of course then the sum of the corresponding α and β contributions. We note that a recent study dealing with the application of multicentre bond indices to the excited-state aromaticity of benzene and cyclobutadiene [47] used natural spin-orbitals (NSOs) instead of natural orbitals (NOs), even for singlet states. The multicentre indices (MCIs) reported in the present study were calculated with our own programs using the QTAIM approach [67], with the required domain-condensed overlaps generated using the AIMAll program [70].

Using NSO rather than NO expansions is of course straightforward for systems with nonzero spin because one may use combinations of the charge-density and spin-density matrices to generate the different NSO expansions for the α and β one-electron densities. On the other hand, the spin-density matrix is null for singlet states and so the NSO expansions of the α and β one-electron densities must coincide. Given that the NO occupations are split equally between the α and β NSOs we may refer to this as the “half-electron scheme”. Because the α and β NSOs are the same, and coincide with the NOs, the total α and β k-centre multicentre indices calculated using this scheme are straightforwardly related to those calculated from the NO occupations (Equation (4)) via the trivial scaling factor ${\left(\frac{1}{2}\right)}^{k-1}$.

A potential problem for singlet states can easily be appreciated by noticing that α and β NSO occupations do not distinguish between the combination of “singlet diradical” determinants $\phi \overline{\psi }$ and $\overline{\phi }\psi$ and the combination of “closed shell” determinants $\phi \overline{\phi }$ and $\psi \overline{\psi }$. As a consequence, the resulting multicentre indices do not take explicit account of which states have a high degree of diradical character and which of them do not. The degree of diradical character could of course be an important feature for considerations of aromaticity. As such the “half-electron” approach, although formally the correct one, might appear to be slightly questionable when considering, for example, an inherently diradical species such as the singlet ground state of square cyclobutadiene. As was demonstrated in the seminal study by Salem and Rowland [71], the diradicals represented in the simplest model by two degenerate orbitals occupied by two electrons form four electronic states, namely two biradical states (singlet and triplet) and two zwitterionic singlet states:

$$\begin{gathered} S_1,T_1=\frac{1}{\sqrt{2}}\left(\left|\varphi \overline{\psi }\right|\pm \left|\psi \overline{\phi }\right|\right) \\ {Z}_1\pm Z_2=\frac{1}{\sqrt{2}}\left(\left|\phi \overline{\phi }\right|\pm \left|\psi \overline{\psi }\right|\right) \end{gathered}$$

(5)

Although all of these states differ at the level of the pair density, and thus also of the energy, all of the spinless one-electron densities coincide:

$$\rho \left(1\right)={\phi }^2\left(1\right)+{\psi }^2\left(1\right)$$

(6)

Such observations made it seem attractive to consider taking explicit account of singlet diradical character by means of artificial modifications of the α and β NSO occupations. In the simple case of the singlet ground state of square cyclobutadiene, which features two singly occupied orbitals, we could for example consider that the first of them is pure α and the other one is pure β spin. We use the label “diradical scheme” for this somewhat artificial approach in which the α and β densities are now allowed to be different, albeit they still add to the correct total. In actual practice we did unfortunately find that manipulations of this sort were far from satisfactory. There were particular complications and uncertainties for cases such as states of benzene which feature two pairs of degenerate orbitals (each corresponding of course to one of the E irreducible representations in D_6h symmetry). We were also concerned that some invariances to orbital rotations might be lost and we noticed that using analogous manipulations for triplet states resulted in “artificial” NSO occupations that bear no obvious resemblance to the actual ones. As a consequence, we reluctantly mostly abandoned this diradical scheme and so we focus here on our results that were obtained for the singlet and triplet states with the actual NSO occupation numbers. Nonetheless, because of this inability of the first-order density matrix to reflect important features that are present in the wavefunction and the pair density, we considered it useful to take into account in our considerations also the eventual manifestation in the wavefunction of diradical character, given that it could be very important in the evaluation of the degree of aromaticity. For this purpose and in order to provide additional insights into the nature of the individual excited states of the molecules studied we also quantify the contributions to the occupation numbers that arise from diradical character. (Note that, instead of using the actual NSO occupations for the singlet states, we could have used Equation (4) with the NO occupations, rescaling the resulting 4- and 6-centre indices by $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.$ and $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$32$}\right.$, respectively).

2.3. Similarity of Excited States

As an auxiliary tool to assess the aromaticity and/or antiaromaticity of individual excited states we also used a molecular similarity index that is based on the number of overlapping electrons (NOEL), as was introduced some time ago by Cioslowski [72]. In essence, the index of similarity between two molecules A and B is defined in terms of the first-order density matrices of the corresponding molecules as

$$d_{AB}={{\displaystyle \int }}^{ }{\left|{\Gamma }^A\left(x,x^{\prime }\right)-{\Gamma }^B\left(x,x^{\prime }\right)\right|}^2 dx d{x}^{\prime }$$

(7)

where the density matrices ${\Gamma }^X$ are conveniently represented by NSO expansions. The smaller the values of $d_{AB}$, the more similar are the first-order densities of the molecules A and B.

NOEL-based comparisons of systems with different geometries would involve also the optimization of the mutual positions of the two molecules, so as to maximize the similarity. In the present work, however, the comparisons of the different states of a given molecules are much more straightforward because of our decision to consider vertical excitations, i.e., fixed geometries. The above general formula then reduces to

$$\begin{array}{cc}{d}_{AB}=d_{AB}^{\alpha }+d_{AB}^{\beta }=& \\ & Tr\left(P_A^{\alpha }S\right)\left(P_A^{\alpha }S\right)+Tr\left(P_B^{\alpha }S\right)\left(P_B^{\alpha }S\right)-2Tr\left(P_A^{\alpha }S\right)\left(P_B^{\alpha }S\right)+\hfill \\ & Tr\left(P_A^{\beta }S\right)\left(P_A^{\beta }S\right)+Tr\left(P_B^{\beta }S\right)\left(P_B^{\alpha \beta }S\right)-2Tr\left(P_A^{\beta }S\right)\left(P_B^{\beta }S\right)\hfill \end{array}$$

(8)

where $P_X^{\alpha }$ and $P_X^{\beta }$ denote the α and β one-electron density matrices, respectively, for a particular electronic state of a given molecule and S is the overlap matrix.

3. Results and Discussion

As our primary tool for the evaluation of excited-state aromaticity we used the multicentre bond indices whose calculation requires knowledge of the first-order density matrix provided in a quantum chemical calculation via natural orbitals and their occupation numbers. In view of the potential problems mentioned in the previous section, the use of quantities based on the first-order density matrix might not always be a completely satisfactory approach: This matrix is not able to reflect all of the features of a more complicated wavefunction and, in certain cases, the features not carried over could be of crucial importance. The relevance of this concern can be illustrated using simple considerations applied to wavefunctions exhibiting diradical character which are often encountered when describing excited electronic states.

It is of course entirely straightforward to construct an expansion of a CASSCF wavefunction in terms of determinants built from NOs so as to reproduce the already known NO occupation numbers. Then we can determine also the net contributions arising from determinants in which a particular NO is singly occupied. Especially for singlet states, the results provide a useful quantitative measure of the extent of diradical character. In most cases, sufficient qualitative information can be obtained just by examining the compositions of the most important determinants in the expansion and, as shown below, doing so is essential when evaluating the reliability of multicentre indices and NOEL-based similarity values as aromaticity criteria.

The need for a more detailed analysis of the nature of each individual excited state is highlighted by the observation that states of very different character can have fairly similar patterns of NO occupation numbers, as can be seen in

Table 1. Active-space NO occupation numbers for the S₀, S₁, S₂ and T₁ states of C₆H₆ and active-space NSO occupation numbers for the T₁ state (all in descending order). Values in brackets show the net contributions from all determinants in which a particular NO is singly occupied.

	S0	S1	S2	T1	T1
					α	β
η1	1.961 (0.028)	1.863 (0.107)	1.952 (0.039)	1.910 (0.078)	0.986	0.924
η2	1.902 (0.046)	1.445 (0.500)	1.494 (0.051)	1.464 (0.487)	0.966	0.498
η3	1.902 (0.046)	1.445 (0.500)	1.494 (0.051)	1.464 (0.487)	0.966	0.498
η4	0.100 (0.046)	0.569 (0.500)	0.524 (0.051)	0.536 (0.487)	0.503	0.034
η5	0.100 (0.046)	0.569 (0.500)	0.524 (0.051)	0.536 (0.487)	0.503	0.034
η6	0.036 (0.028)	0.109 (0.087)	0.012 (0.007)	0.090 (0.079)	0.077	0.013

,

Table 2. Active-space NO occupation numbers for the S₀, S₁, S₂ and T₁ states of square C₄H₄ and active-space NSO occupation numbers for the T₁ state (all in descending order). Values in brackets show the net contributions from all determinants in which a particular NO is singly occupied.

	S0	S1	S2	T1	T1
					α	β
η1	1.905 (0.073)	1.835 (0.031)	1.994 (0.001)	1.914 (0.072)	0.993	0.921
η2	1.000 (1.000)	1.005 (0.000)	1.000 (0.000)	1.000 (1.000)	0.964	0.036
η3	1.000 (1.000)	1.005 (0.000)	1.000 (0.000)	1.000 (1.000)	0.964	0.036
η4	0.095 (0.073)	0.155 (0.031)	0.006 (0.001)	0.086 (0.072)	0.079	0.007

and

Table 3. Active-space NO occupation numbers for the S₀, S₁ and T₁ states of S₂N₂ and active-space NSO occupation numbers for the T₁ state (all in descending order). Values in brackets show the net contributions from all determinants in which a particular NO is singly occupied.

	S0	S1	T1	T1
				α	β
η1	1.997 (0.002)	1.997 (0.002)	1.998 (0.002)	0.999	0.999
η2	1.997 (0.002)	1.996 (0.003)	1.997 (0.002)	0.999	0.998
η3	1.996 (0.003)	1.994 (0.005)	1.995 (0.004)	0.999	0.998
η4	1.991 (0.005)	1.990 (0.007)	1.994 (0.006)	0.998	0.996
η5	1.986 (0.012)	1.984 (0.013)	1.993 (0.005)	0.998	0.994
η6	1.981 (0.016)	1.968 (0.025)	1.983 (0.016)	0.997	0.991
η7	1.980 (0.010)	1.965 (0.028)	1.982 (0.015)	0.997	0.984
η8	1.964 (0.028)	1.960 (0.031)	1.967 (0.027)	0.996	0.983
η9	1.962 (0.030)	1.960 (0.035)	1.965 (0.024)	0.991	0.983
η10	1.960 (0.030)	1.911 (0.034)	1.965 (0.025)	0.984	0.982
η11	1.908 (0.030)	1.009 (0.979)	1.014 (0.984)	0.983	0.017
η12	0.126 (0.047)	0.995 (0.970)	1.008 (0.986)	0.983	0.017
η13	0.041 (0.031)	0.126 (0.060)	0.037 (0.026)	0.019	0.016
η14	0.040 (0.029)	0.058 (0.045)	0.037 (0.026)	0.019	0.016
η15	0.038 (0.030)	0.044 (0.032)	0.035 (0.026)	0.019	0.014
η16	0.033 (0.024)	0.043 (0.031)	0.031 (0.022)	0.017	0.012

. Such similarities are displayed, for example, by the S₁ and S₂ states of benzene, as well as by all three singlet states of square cyclobutadiene that we examined. The absence of pronounced differences between the patterns of NO occupation numbers is a cause for concern because it is not clear how the multicentre indices, as well as the NOEL-based $d_{AB}$ values, will be able to distinguish properly between such electronic states unless the shapes of the NOs change sufficiently between states.

The net contributions from all determinants in which a particular NO is singly occupied to the wavefunctions for the various electronic states of benzene, square cyclobutadiene and disulfur dinitride are shown in Table 1, Table 2 and Table 3. Clearly, both the S₁ state and, of course, the T₁ state of benzene exhibit significant levels of diradical character, unlike the S₀ and S₂ states (see Table 1). In the case of square cyclobutadiene (Table 2) the states with significant levels of diradical character are S₀ and again, of course, T₁, whereas there only minor traces of such character in the S₁ and S₂ states which appear to be zwitterionic [71]. Moving on to S₂N₂, we can see from Table 3 that that it is only the S₀ ground state that has slight diradical character, whereas there is strong diradical character in both of the excited states. Inspection of the symmetries of the NOs for the electronic states of this molecule reveals that whereas the dominant diradical character in T₁ comes from two unpaired π electrons, much as in the corresponding states of C₆H₆ and C₄H₄, that in the S₁ state is associated with the coupling of an unpaired σ electron to an unpaired π electron. Therefore, we can expect that the NOEL-based similarity index for S₂N₂ will show significant differences between the valence σ system of the S₁ state and those of the S₀ and T₁ states.

The values of multicentre QTAIM indices for C₆H₆, C₄H₄ and S₂N₂ calculated using NSOs are shown in

Table 4. 6-centre MCIs for the S₀, S₁, S₂ and T₁ states of C₆H₆, 4-centre MCIs for the S₀, S₁, S₂ and T₁ states of square C₄H₄ and 4-centre MCIs for the S₀, S₁ and T₁ states of S₂N₂. QTAIM/6-311++G(2d,2p) for C₆H₆ and C₄H₄; QTAIM/cc-pVTZ for S₂N₂. Values in brackets show the π-only (C₆H₆ and C₄H₄) and π-only valence (S₂N₂) components of the total index.

	C6H6	C4H4	S2N2
S0	0.0160 (0.0159)	0.0360 (0.0336)	0.0525 (0.0504)
S1	0.0016 (0.0014)	0.0263 (0.0240)	0.0128 (0.0164)
S2	0.0029 (0.0027)	0.0477 (0.0454)
T1	0.0030 (0.0029)	0.0774 (0.0749)	−0.0045 (−0.0069)

. For singlet states we have used the actual NSO occupations, i.e., the “half-electron” scheme. While inspecting the numbers in this table, it is useful to adopt the largest total 6-centre MCI value of 0.016 for the archetypal example of an aromatic system, the ground state of benzene, as a yardstick for assessing the values of this index for the other states of this molecule. Whereas the value of the index for the first excited singlet state is smaller than its S₀ counterpart by an order of magnitude, those for the S₂ and T₁ states which turn out to very close are, instead, smaller by a factor of about five. It should be mentioned that, in keeping with expectations, the total values of the indices are dominated by their π-only components for all of the states included in Table 4. The situation in the case of S₂N₂ is similar to that in C₆H₆: Once again, the largest value of the 4-centre MCI is that for the ground state, whereas the excited-state indices are considerably smaller. The 4-centre MCIs for the electronic states of square C₄H₄ follow a different pattern: The largest value corresponds to the lowest triplet state T₁, whereas the indices for the three singlet states are considerably smaller.

The changes in the values of the MCIs in Table 4 between S₀ and T₁ states are fully consistent with Baird’s original rule [42] which addresses aromaticity reversals involving the singlet ground and lowest triplet electronic states only. However, when singlet states come into play, there are some notable discrepancies from the behaviour expected on the basis of the results from previous studies which discuss in detail ground- and excited-state magnetic properties, including several types of NICS [44,45,52]. Let us start with benzene. The S₁ state of C₆H₆ is correctly classified as antiaromatic, with S₁ being more antiaromatic than T₁ which is in line with predictions based on magnetic properties [44,45] and several multicentre and delocalization indices [47]. However, S₂ in C₆H₆ is predicted to be just as antiaromatic as T₁ whereas the magnetic properties of this state, calculated at the same geometry and level of theory, strongly suggest that it is even more aromatic than the ground state S₀ [44]. Incidentally, S₂ in C₆H₆ was classified as more antiaromatic than S₁ in [47] but this was due to analysing the doubly-degenerate S₄ rather than S₂ (for details, see [44]). We continue our analysis with square C₄H₄. According to the 4-centre MCI values, S₁ is the most antiaromatic state of this molecule whereas the isotropic shielding isosurface for this state and other magnetic properties, calculated at the same geometry and level of theory, show clearly that it is aromatic [44]. S₂ is predicted to be less antiaromatic than S₀ which is in agreement with the results of magnetic property calculations [44,45]. Finally, according to the respective 4-centre MCI values, in S₂N₂ S₁ is less antiaromatic than T₁ whereas magnetic property calculations suggest the opposite ordering for these two states [52]. One further observation is that the values of the 6-centre MCI for benzene are, for most states, of smaller magnitude than the corresponding 4-centre MCIs for square cyclobutadiene and disulfur dinitride, which precludes comparisons between MCIs for different rings.

The discrepancies between the current MCI-based assessments of the aromaticities of singlet excited states and those coming from calculations of various magnetic second-order response properties [44,45] underline our concerns as to whether first-order density-based indices are capable of reflecting more subtle effects whose description relies on the more detailed information inherent to wavefunctions or pair densities for particular electronic states.

Similar concerns are associated with comparisons utilizing the NOEL-based quantity d_AB. This quantity was originally designed [72] as a quantitative measure of the similarity between the electron densities of different molecules but, in this study, we use it to investigate the similarity between the electron densities of different excited states of one molecule. The values of the NOEL-based similarity index $d_{AB}$ for the low-lying electronic states of benzene, square cyclobutadiene and disulfur dinitride, calculated using NSOs, are summarized in

Table 5. Similarity indices $d_{AB}$ calculated at the CASSCF(6,6)/6-311++G(2d,2p) level for low-lying electronic states of C₆H₆. Values in brackets are π-only contributions to the total index.

State	S1	S2	T1
S0	0.4402 (0.4396)	0.3639 (0.3612)	0.8323 (0.8319)
S1		0.0453 (0.0433)	0.4465 (0.4464)
S2			0.4866 (0.4842)

,

Table 6. Similarity indices $d_{AB}$ calculated at the CASSCF(4,4)/6-311++G(2d,2p) level for low-lying electronic states of square C₄H₄. Values in brackets are π-only contributions to the total index.

State	S1	S2	T1
S0	0.0121 (0.0117)	0.0557 (0.0532)	0.8672 (0.8672)
S1		0.0417 (0.0408)	0.8754 (0.8752)
S2			0.9086 (0.9067)

and

Table 7. Similarity indices $d_{AB}$ calculated at the CASSCF(22,16)/cc-pVTZ level for low-lying electronic states of S₂N₂. Values in brackets are (π-only, σ-only) valence contributions to the total index.

State	S1	T1
S0	0.9596 (0.4796, 0.4799)	1.8061 (1.7896, 0.0166)
S1		2.8024 (2.2926, 0.5098)

. These results demonstrate clearly, in keeping with expectations, that the total $d_{AB}$ indices for C₆H₆ and C₄H₄ are, in all cases, dominated by their π-only components. In the case of S₂N₂, $d_{AB}$ indices involving the S₁ state show large σ contributions due to the composition of the wavefunction for this state (see above).

The data in Table 5 suggest some similarity between the S₀ and S₂ states of benzene, in line with the expected aromaticity of S₂ [44], as well as very little similarity between the S₀ and T₁ states, in agreement with Baird’s rule. However, the surprisingly high level of similarity between the S₁ and S₂ states which have been classified as antiaromatic and aromatic, respectively [44], is very much out of line with the rather different magnetic properties of these states. Somewhat surprising are also the comparable levels of similarity between the S₁ and T₁ states, both of which are supposed to be antiaromatic, and the S₂ and T₁ states, which are supposed to be aromatic and antiaromatic, respectively [44].

Other similarity assessments of questionable utility can be found amongst the data for square cyclobutadiene that are presented in Table 6, starting with the high level of similarity between the S₀ and S₁ states which is unexpected, in view of the predicted aromaticity reversal between these states [44]. Both S₁ and T₁ are expected to be aromatic [44], but the level of similarity between these states is comparable to that between S₂ and T₁, which are expected to be antiaromatic and aromatic, respectively.

The $d_{AB}$ indices are doing a better job in the case of S₂N₂ (see Table 7): S₀ and S₁ are quite dissimilar, and so are S₀ and T₁, as expected for comparisons between wavefunctions corresponding to aromatic and antiaromatic states. S₁ and T₁ come out as very dissimilar which is not unrealistic, as these states have been predicted to show very different levels of antiaromaticity [52]. As has been mentioned, the σ-only valence contributions to $d_{AB}$ are large in all comparisons involving the first singlet excited state, due to the composition of the S₁ wavefunction (see above).

We have shown that the multicentre indices (MCIs) examined in this work perform well for aromaticity reversals involving the singlet ground and lowest triplet electronic states which are covered by Baird’s original rule [42]. Our attempts to apply these indices to aromaticity reversals involving singlet excited states were less satisfactory. While this may seem disappointing, since aromaticity/antiaromaticity switching can be predicted even using simple topological resonance energies [43], it should be emphasized that TRE-based studies do not distinguish between singlet and triplet excited states, and all MCI problems arise when dealing with singlet excited states. When dealing with singlet ground states, 6-centre indices have been found to correlate very well with the energetic stabilization resulting from cyclic delocalization in individual benzene rings in polycyclic aromatic hydrocarbons [29,40]; multicentre indices have also been reported as a reliable measure of aromaticity in all-metal clusters [73].

One potential source of the problems experienced when trying to apply multicentre indices to singlet excited states can be associated with the reasons behind the very good performance of MCIs for the ground states of polycyclic aromatic hydrocarbons. The correlation between MCIs and energetic stabilization stems from Coulson’s integral formula [74] which describes quantitatively the extent of energetic stabilization/destabilization associated with cyclic conjugation. However, this formula can be applied only to the ground states of conjugated hydrocarbons and, as there is no equivalent formula for excited states, there is also no straightforward way of measuring the energetic effects resulting from cyclic conjugation in other than ground states (TRE is applicable only to the lowest excited state). The absence of an energy-based justification of excited-state MCIs may have adverse impact on their performance in comparison to ground-state MCIs. On the other hand, cyclic conjugation in excited states can be thought to induce excited-state ring currents and the nature of these currents (paratropic vs diatropic) is decisive for excited aromaticity. These ring currents can be integrated using the Bio-Savart law (as shown, for example, in [75]), producing excited-state magnetic shielding tensors such as those calculated and analysed in [44,45,46,52] which explains why the magnetic properties of excited states provide reliable measures of excited-state aromaticity.

In addition to these somewhat qualitative arguments, more detailed theoretical considerations can be used to identify additional factors affecting the performance of the multicentre indices for singlet excited states. For this purpose, it is useful to refer again to the paper by Salem and Rowland [71] dealing with the electronic structure of diradicals. As we mentioned previously, although all four biradical and zwitterionic states for the simple two-orbital model (see Equation (5)) differ in energy (and, consequently, in wavefunction and in pair density), their one-electron densities are exactly the same. One straightforward implication is that the first-order density matrices for different electronic states of real systems could omit important details, the absence of which would result in multicentre indices giving misleading information about the extent of similarity between these states. Although the discussion in [71] is focused on inherently diradical species, similar problems, arising from details not available within the first-order density matrix, are apparently more general since, as demonstrated in this study, partial diradical character is evident even in the excited-state wavefunctions of a paradigmatic molecule such as benzene. On the other hand, the undeniable usefulness of multicentre indices as a measure of ground-state aromaticity [29,30,32,40,73] can be attributed to the fact that, at the closed-shell SCF and Kohn–Sham levels of theory, the first-order density matrix determines all higher-order densities so that energy-related quantities are described correctly.