Exploitation of Baird aromaticity and Clar’s rule for tuning the triplet energies of polycyclic aromatic hydrocarbons

: Polycyclic aromatic hydrocarbons (PAH) are a prominent substance class with a variety of applications in molecular materials science. Their electronic properties crucially depend on the bond topology in ways that are often highly non-intuitive. Here, we study, using density functional theory, the triplet states of four PAHs based on the biphenylene motif ﬁnding dramatically different triplet excitation energies for closely related isomeric structures. These differences are rationalised using a qualitative description of Clar sextets and Baird quartets, quantiﬁed in terms of nucleus independent chemical shifts, and represented graphically through a recently developed method for visualising chemical shielding tensors (VIST). These results are further interpreted in terms of a 2D rigid rotor model of aromaticity and through an analysis of the natural transition orbitals involved in the triplet excited states showing good consistency between the different viewpoints. We believe that this work constitutes an important step in consolidating these varying viewpoints of electronically excited states.

The properties of PAHs can be widely tuned using the bond topology [17][18][19] as well as by doping with heteroatoms [9,[20][21][22][23]. To do so effectively, it is particularly helpful to establish clear structure-property relationships bridging between the molecular structure and the observable properties. As one option, this can be achieved by considering the biradical character and unpaired electrons [22,[24][25][26][27]. A particularly attractive alternative option is given via aromaticity noting that, following Hückel's [28] and Clar's [29] rules, qualitative predictions of molecular properties can already be made following simple rules based on the bond topolgy and the number of electrons [18,19,22]. Even more Baird's rule [30] of excited state aromaticity allows to effectively predict triplet energies [31][32][33][34].
However, it is not trivial to move from the qualitative description of aromaticity based on bond topology to quantifiable aromaticity descriptors and a number of methods with different merits and shortcomings have been put forward for this task. On the one hand, geometric criteria are used, such as the harmonic oscillator model of aromaticity [35]. On the other hand, magnetic properties such as induced current densities [36,37] and nucleus independent chemical shifts (NICS) [38] are used as direct probes of the electronic structure. However, a special challenge in the application of these methods is that the underlying quantities (the current density susceptibility and the chemical shielding) are both represented by tensor fields making it difficult to view all of the relevant information together. To meet this challenge, we have recently developed the VIST (visualisation of chemical shielding tensors) method [39] providing graphical information about the full chemical shielding tensor, thus, providing detailed information on (anti)aromaticity along with its local variations and anisotropy. The method is exemplified in Fig. 1: The principal axes of the shielding tensor at a given point in space are, first, determined via an eigenvalue decomposition and subsequently visualised via dumb-bells whose size depends on the associated eigenvalue.
Using the VIST method, it is the purpose of this work to illustrate the concept of Baird aromaticity in a class of PAHs based on the biphenylene motif (see Fig. 2). Recently, the biphenylene unit [40] has been considered for molecular conductance [8] and singlet fission [33]. Computational studies on the biphenylene unit integrated into larger PAHs have highlighted its effect on Baird aromaticity [41] and biradical character [42]. Here, we will investigate how minor structural variations in biphenylene derivatives can have a dramatic effect on triplet energies and analyse this effect in terms of bond topology, shielding tensors, and molecular orbitals. For this purpose, we construct three derivatives of biphenylene by formally fusing two benzene rings to the molecule. These molecules (3)(4)(5) along with biphenylene (2) and cyclobutadiene (1) are shown in Fig. 2. Resonance structures of these molecules are shown with Clar sextets highlighted in blue and Baird quartets in red. Following Ref. 41, we hypothesise that the structures with simultaneous Clar sextets and Baird quartets (i.e. 4 and 5) will exhibit particularly pronounced Baird aromaticity leading to lowered triplet energies. Within this work, we start the discussion with a theory section illustrating the concept of Baird aromaticity within the MO picture, discussing its effects on the chemical shielding, and providing the basics of the VIST method used in the graphical representation. Subsequently, the computational results are presented starting with a basic consideration of energies, proceeding to a discussion of the shielding tensors in the CBD and benzene building blocks, and continuing with a detailed discussion of the properties of molecules 2-5 considering shielding tensors as well as molecular orbitals. Before concluding, we present a discussion of the relation between Baird aromaticity and other ways of describing excited-state electronic structure.

Aromaticity -The rigid rotor viewpoint
In light of the discussion to follow, it is expedient to consider the phenomenon of aromaticity in the context of a simple 2D rigid rotor (see Ref. 43 for a related discussion based on Hückel theory). The eigenstates of a 2D rigid rotor [44] are schematically represented in Fig. 3 (b). The first state is non-degenerate whereas all other levels come in pairs of two degenerate states sorted according to an increasing value of the angular momentum quantum number m l . The lowest state can be identified with the s-orbital (m l = 0) of the hydrogen atom. The next levels have the same in-plane angular wavefunctions as the p x /p y (|m l | = 1), d xy /d x 2 −y 2 (|m l | = 2) pairs etc increasing in energy and angular momentum.
Inspection of the molecular orbitals (MOs) of planar conjugated cyclic molecules shows that the MOs have similar symmetry properties as these rigid rotor eigenstates, i.e. there is one MO with no nodal plane, two MOs with one nodal plane, two MOs with two nodal planes etc. These MOs are shown are shown in Fig. 3 (a) for benzene and (c,d) cyclobutadiene (CBD) highlighting their similar structure among each other and with the idealised 2D rigid rotor. Since the lowest level is non-degenerate and all others are doubly degenerate, it follows that one needs an odd number of electron pairs (4n + 2 electrons) to fill up any given level of degenerate orbitals, as shown in Fig. 3 (a). On the other hand, if an even number of electron pairs (4n electrons) are present, then the highest level will only be half-filled -usually accompanied by symmetry-breaking -as shown for CBD in Fig. 3 (c). In line with general chemical knowledge, filled levels provide enhanced chemical stability and, hence, cyclic systems with 4n + 2 electrons are particularly stable and systems with 4n electrons are not. This is the essence of Hückel's rule.  In Figure 3 (d), the case for the lowest triplet state in CBD is shown. An electron was promoted from the HOMO to the LUMO leading to two singly occupied orbitals, whose degeneracy is restored as CBD obtains four-fold symmetry. The occupation patterns now look similar to benzene in the sense that both quasidegenerate orbitals are evenly occupied. Indeed, CBD in its triplet state is considered aromatic according to Baird's rule [30], which states a reversal between aromaticity and antiaromaticity in the triplet state. Due to the low HOMO/LUMO gap between the two quasidegenerate orbitals in singlet CBD one expects a low triplet excitation energy. However, due to the high degree of spatial overlap between HOMO and LUMO inducing enhanced exchange repulsion [45][46][47], one expects singlet excitation energies to be appreciably higher. The combination of low triplet energies with high exchange splitting provides a possibility of realising the energetics required for the singlet fission process [3] and, indeed, Baird aromaticity has been considered as a promising route toward new molecules capable of singlet fission [32][33][34]48].

Definition of the magnetic shielding tensors
The magnetic shielding tensor σ(R) at a given point in space R is defined by the relation where B ext and B ind are vectors describing the applied external and induced field, respectively [49]. Thus, roughly speaking, the shielding is the negative of the ratio between external and induced field. In chemical systems, the induced field tends to be about a million times smaller than the external field and, hence, the shielding is on the order of a few ppm. More specifically, the shielding is given as a tensor field represented by a 3 × 3 matrix where all matrix elements depend on the position R.
The nucleus independent chemical shift (NICS), which serves as a prominent aromaticity criterion [38,50], is -in its basic form -defined as the negative of the spherically averaged magnetic shielding In practical computations, the shielding is usually evaluated as the mixed second derivative of the energy with respect to an external magnetic field B β and the nuclear magnetic moment µ γ [51][52][53] where β and γ are Cartesian coordinates (x, y, z). Note that in Eq. (4) and all following equations atomic units are employed. Practically evaluating this expression leads to the following form [51,54,55] where α is the fine-structure constant, Ψ 0 /Ψ I are the ground/excited state wavefunctions, E 0 /E I are their energies,L β /L γ are angular momentum operators with respect to the gauge origin and the point probed, and r = (x, y, z) is the distance vector from the gauge origin. The total shielding is usually considered a sum of two terms σ dia γβ (diamagnetic) and σ para γβ (paramagnetic) noting, however, that the distinction is not unique and their relative magnitude depends on the gauge origin chosen [53].
The first term σ dia γβ , denoted diamagnetic, only depends on the distribution of the electrons in space. To get a qualitative understanding of this term, we may set the origin into the point probed (R = 0) and evaluate one of its diagonal components where ρ(r) is the electron density, r = |r|, and θ is the angle from the z-axis. This term, related to the Lamb formula, represents the textbook view of NMR spectroscopy -the chemical shielding is proportional to the electron density ρ(r) in the vicinity of the nucleus of interest [44,55].
The second term σ para γβ , denoted paramagnetic, has a more complex form measuring how a perturbation of the wavefunction by an external magnetic field induces an angular momentum. It is usually given in a perturbative expansion in terms of the wavefunctions Ψ I and energies E I of the excited states of the system [52,54] as shown in Eq. (5). Viewing, again, only one diagonal component at the gauge origin, we obtain Whereas, the diamagnetic shielding given by Eq. (6) is a non-specific term present whenever there is electron density, one finds that the paramagnetic shielding [Eq. (7)] depends on the precise shapes and energies of the orbitals involved.
It is possible to construct a semi-quantitative model for the shielding contributions in (anti)aromatic systems based on the 2D rigid rotor discussed in Sec. 2.1. For this purpose, we assume that the electron density is distributed in a ring of radius R 0 around the gauge origin (r = R 0 ). In addition we assume this ring to be in the xy-plane, which means sin 2 θ = 1. Inserting these values into Eq. (6), then leads to a diamagnetic shielding of where N el = ρ(r)dr is the number of electrons contributing to the ring current.
In the case of an antiaromatic system, we can provide an estimate of the paramagnetic shielding by considering the HOMO/LUMO contribution. As shown in Fig. 3 (c), the HOMO and LUMO of an antiaromatic system both relate to the same angular momentum quantum number m l (see also Refs 43,50). The corresponding real orbitals are given as π −0.5 cos(m l φ) and π −0.5 sin(m l φ) where φ is the azimuthal angle in the xy-plane (cf. Ref 44). In addition, we can estimate the excitation energy of the first state via the HOMO/LUMO gap, i.e., E 1 -E 0 ≈ E HL (thus ignoring any post-MO contributions to the excitation energy [47]). Using these approximations, one obtains the following paramagnetic shielding.
The analogous expression vanishes for aromatic systems as the matrix element of the angular momentum operator between functions with different m l vanishes under the assumptions stated. Note, that this is also true for a Baird aromatic system, as shown in Fig. 3 (d), considering that due to the Pauli principle any excitation into a higher lying triplet state must involve orbitals of higher m l .
The results are summarised in Tab. 1. For an aromatic system only the diamagnetic contribution plays a role and, hence, enhanced shielding going along with negative NICS values is expected at the centre of an aromatic ring. For an antiaromatic system diamagnetic shielding and paramagnetic deshielding both play a role and if the latter predominates, then one finds the characteristic deshielding found at the centre of antiaromatic rings. Table 1. Diamagnetic and paramagnetic shielding contributions computed at the centre of a 2D rigid rotor in terms of the number of electrons (N el ), the radius (R 0 ), the angular momentum quantum number (m l ) of the HOMO and LUMO, and the HOMO/LUMO gap (E HL ).

Diamagnetic shielding
Finally, it is interesting to discuss how these terms change with increasing ring size R 0 . If one assumes that with increasing ring size also the number of electrons N el increases, and specifically that N el /R 0 is constant, then one finds that σ dia zz should not have a strong dependence on the ring size. In addition, we can assume that the m l quantum number of the HOMO and LUMO of an aromatic system is proportional to the number of electrons and, thus, the ring size. Hence, also the N el m l 2 /R 0 3 term should be approximately constant meaning that the paramagnetic contribution is also insensitive to the ring size. These considerations illustrate the suitability of NICS values to probe rings of various sizes [38].

Visualisation of shielding tensors (VIST)
Eq. (2) shows that the chemical shielding is a tensor field going over all space (cf. Ref. 52), i.e. it is given by a 3 × 3 tensor containing 9 independent values at every point in space. Whereas scalar fields can be represented in 3D space via isosurfaces and vector fields via arrows, one has to construct a more involved representation for a tensor field. To do so, we have developed the VIST (visualisation of chemical shielding tensors) method [39]. VIST proceeds by computing the principal axes q (1) , q (2) , q (3) of the shielding tensor via an eigenvalue decomposition For visualisation, we construct a local coordinate system oriented according to the eigenvectors q (i) and visualise the three components as dumb-bells whose size and length depends on the absolute value of the associated eigenvalue |t (i) | and whose color depends on the sign (blue or red), see Fig. 1. Through encoding the eigenvectors and eigenvalues in the graphical representation, we are able to represent the full information given in the 3 × 3 tensor graphically. Specifically, we draw the length L of the axis and the radius R of the sphere as where t (i) is given in ppm and L and R are given in Å. To compare these results to the NICS values, it is worth noting that in analogy to Eq. (3) the NICS value is a third of the sum of the three eigenvalues according to Reviewing Eq. (1), we find that the principal axes of the shielding tensor are the directions in space where the induced magnetic field is parallel to the external magnetic field. If we set B ext = B 0 q (i) and insert this into Eq. (1), we find by applying Eq. (10) that
Chemical shielding tensors for S 0 and T 1 states were computed using RKS and UKS, respectively, at the respective optimised geometries using TPSSh/def2-TZVP as implemented in Gaussian 09 [67] using gauge-including atomic orbitals [51] and applying tight SCF convergence criteria. The shielding tensors were analysed using the VIST (visualisation of chemical shielding tensors) method [39] employing a pre-release version of TheoDORE 2.4 [68]. The visual molecular dynamics (VMD) [69] program was used as a graphical backend for creating the tensor representations in connection with the molecular structures and NTO pictures.
The underlying research data are available via a separate repository [70]: Molecular coordinates; input/output files for Gaussian, TheoDORE, and Q-Chem.

Energetics
We start with a discussion of the basic energetics of the molecules presented in Fig. 2 considering their geometries optimised in their lowest singlet (S 0 ) and triplet (T 1 ) states (Tab. 2). Starting with the adiabatic triplet excitation energy (∆E(S 0 /T 1 )), we find that it almost vanishes for cyclobutadiene (CBD, 1) in line with CBD's expected antiaromaticity. Moving to biphenylene a substantial rise in triplet energy to 1.92 eV is seen. Interestingly, increasing the π-system to 3 further increases the triplet energy to 2.64 eV. Moving to the isomeric molecules 4 and 5 yields a dramatic drop in the excitation energy to just above 1 eV. Table 2. Adiabatic excitation energies to the first triplet state ∆E(S 0 /T 1 ) and total energies E(S 0 ) and E(T 1 ) for the molecules studied (all energies given in eV). a

Molecule
∆E(S 0 /T 1 ) E(S 0 ) E(T 1 ) Due to the fact that molecules 3, 4, and 5 are isomeric, we can analyse their energetics in some more detail. In Tab. 2 the total energies of their singlet E(S 0 ) and triplet E(T 1 ) states, all referenced against the S 0 energy of 3, are listed. This shows that the formal interconversion from 3 to 4 or 5 is strongly endothermic (≈0.5 eV) in the singlet state. However, in the triplet state the energetic ordering is reversed meaning that 4 and 5 are stabilised by more than 1 eV with respect to 3.
The above-mentioned trends clearly run counter to the standard picture stating the excitation energies should decrease with an extension of the π-system as confinement effects are reduced. Conversely, the concept of ground state antiaromaticity in connection with triplet state Baird aromaticity provides an attractive option for explanation. Reviewing the energetics in Tab. 2, we hypothesise that moving from 1 to 3 the antiaromaticity and Baird aromaticity of CBD are both "blurred" as the π-system is extended leading to increased triplet energies. Conversely, reviewing Fig. 2, we find that out of all structures considered only 4 and 5 possess resonance structures with a Baird quartet along with two Clar sextets. Following Ref. 41 we hypothesise that 4 and 5, therefore, possess enhanced ground-state antiaromaticity as well as triplet state Baird aromaticity when compared to 3 explaining both the S 0 and T 1 energies in Tab. 2.
Simply looking at the energetics, it is difficult to assess whether (anti)aromaticity is really the underlying cause for the dramatic differences in energy or whether this is just a coincidence. A different viewpoint and more detailed insight is required to further substantiate the hypothesis. We will endeavour to provide this viewpoint via the computation of nucleus independent chemical shifts [38] and specifically our newly developed method for the visualisation of the underlying chemical shielding tensors [39].

Shielding tensors for the building blocks
We start the discussion of shielding tensors with the building blocks, benzene and cyclobutadiene. In Fig. 4 shielding tensors, computed at the TPSSh/def2-TZVP level, are shown 1 Å above the the centre of the ring. The tensor for benzene Fig. 4 (a) has a particularly simple shape with just one dominant out-of-plane component (-29.5 ppm) and vanishing in-plane components (<1 ppm). This shape reflects the aromaticity of benzene. The dominant component agrees well with NICS(1) zz literature values computed at the B3LYP (-28.5 ppm) [49], PBE0 (-29.5 ppm) [39], MP2 (-30.4 ppm) [71], and CASSCF (-28.8 ppm) [71] levels highlighting that NICS values are robust descriptors of the electronic structure.
The shape of the shielding tensor is in qualitative agreement with the model of Eq. (8) showing that a circular charge distribution should yield a strong out-of-plane contribution of the shielding tensor and vanishing in-plane components. Even more, we can approximate the numerical value. If we consider an effective radius of R 0 = 2.65 a.u. (1.4 Å) and assume that half the π-electrons contribute to the shielding effectively at the point probed (N el = 3), then we can use Eq. (8) to estimate the shielding as This value is indeed in good agreement with the value of -29.5 ppm given in Fig. 4 (a) showing that the qualitative model developed is applicable to understand the physics involved. Note, however, that the extremely good match of less than 1 ppm deviation is certainly just coincidence.
Proceeding to CBD (Fig. 4), we find a strongly deshielded (red, +54.1 ppm) contribution representing the strong antiaromaticity in this system. CASSCF computations [71] place the corresponding NICS(1) zz value at a similar and only slightly lower value of +39.3 ppm highlighting that the method chosen performs reasonably well despite the expected multireference character. In parallel to the molecular plane, one slightly shielded (blue) and one slightly deshielded (red) contribution is found. This anisotropy of the shielding reflects the broken symmetry of CBD. The shielded component points toward the CC double bonds (bond length of 1.33 Å) and the deshielded component toward the single bonds (1.57 Å) highlighting the different properties of the two bonds.
Proceeding to the qualitative model of Tab. 2, we again consider half the electrons as being involved (N el = 2). Furthermore we set m l = 1, R 0 = 1.89 a.u., and E HL = 0.076 a.u. and we find by application of Eqs (6) and (7) the following value for the shielding (15) This is, again, in qualitative agreement with the value of +54.1 ppm obtained in the computation highlighting that the viewpoint adopted here is a valid approximation of the rather involved shielding expressions.
Proceeding to the triplet state of CBD [ Fig. 4) (c)], we find that the dominant component of the shielding becomes negative (-17.6 ppm) indicating Baird aromaticity for the triplet state. This phenomenon can be understood by considering Fig. 3 (d): By promoting an electron from the HOMO to the quasidegenerate LUMO, the system now occupies a filled shell meaning enhanced stability. Reviewing Eq. (7), we find that all orbitals up to m l = 1 are filled and, therefore, no excited states are present that could be coupled to the ground state viaL z meaning that the paramagnetic shielding at the ring centre vanishes. Thus, only the diamagnetic shielding σ dia zz is left predicting a value of -28.2 ppm according to Eq. (15), which again at least produces the same order of magnitude.

Shielding tensors for biphenylene derivatives
Having discussed the individual benzene and CBD molecules, it is of interest to observe what happens when these are fused together obtaining the biphenylene molecule (2). Starting with the VIST plots for the singlet ground state [ Fig. 5 (a)], we find that the antiaromaticity of the central CBD ring and the aromaticity of the outer benzene rings are both reduced in magnitude when compared to the building blocks. The reduction in the primary shielding components is modest for CBD (from +54.1 to +38.2 ppm) whereas a reduction by about two thirds occurs for the benzene components (-29.5 to -11.1 ppm). Next, we proceed to the T 1 state of biphenylene. TDDFT indicates that this state is dominated by the HOMO/LUMO transition with respect to the ground state. To obtain a rigorous graphical representation of the transition, we use the natural transition orbitals (NTOs) of the S 0 → T 1 transition, shown in Fig. 5 (c). The NTOs clearly reflect the shape of the HOMO and LUMO of CBD (see Fig. 3).
This indicates that the transition occurs largely between the two quasi-degenerate orbitals on CBD, which fits with the model of Baird aromaticity.
The T 1 shielding tensors, shown in Fig. 5 (b), reflect at first glance the profound changes the system undergoes upon excitation. The central CBD ring becomes strongly shielded with a value of -35.2 ppm that is higher than, both, isolated CBD in the triplet and benzene in the singlet state. Also the shielding at the benzene rings is enhanced compared to the singlet. These changes highlight the Baird aromaticity of the triplet state in line with the NTOs presented.
Having established the properties of biphenylene, it is of interest to see how the local aromaticity changes upon addition of two benzene rings on the sides obtaining 3 as shown in Fig. 6. Viewing the singlet ground state properties [ Fig. 6 (a)], we find that the central biphenylene unit bears close resemblance to isolated biphenylene [ Fig. 5 (a)] with the exception that the antiaromaticity at the central CBD unit is slightly reduced (to +27.2 ppm). The outer two benzene rings resemble isolated benzene showing strong shielding (-24.4 ppm). The shielding tensors for the lowest triplet state of 3 are shown in Fig. 6 (b). These have a strikingly different appearance to the analogous values for biphenylene [ Fig. 5 (b)]. The central CBD and outer benzene rings have almost vanishing out-of-plane values (below 5 ppm in magnitude) whereas the intermediate benzene ring has a notably deshielded value of +9.7 ppm. To understand this behaviour, we turn to the associated NTOs as shown in Fig. 6 (c). These NTOs differ in two important ways when compared to biphenylene [ Fig. 5 (c)]. First, we find that two NTO transitions contribute significantly to this state, whereas for biphenylene there was only one dominant transition. Second, we find that the dominant NTO pair, making up 57% of the overall transition has an entirely different shape when compared biphenylene. There is only little involvement of CBD and, interestingly its orbitals contribute in the opposite way as expected. The acceptor orbital possesses some contribution of the fully bonding CBD π-orbital, shown at the bottom in Fig. 3 (c). By contrast, the donor possesses some contribution of the fully anti-bonding CBD π * -orbital (not shown in Fig. 3). The dominant NTO pair, thus, highlights the strong mixing of the CBD MOs with the other MOs in the system. This consideration shows that this transition does not agree with the Baird aromaticity model laid out in Fig. 3 (d). On the contrary, the NTOs on the two sides bear resemblance to the HOMO and LUMO of naphthalene (see, e.g. Ref. 72). The emerging antiaromaticity can, therefore, be understood in the same sense as an individual napthalene molecule in its triplet state becomes antiaromatic [73]. Going back to Fig. 6 (c), one finds that only the second NTO pair (contributing 28%) is of the expected shape analogous to biphenylene.
In summary, we find that the difference in shielding tensors between 2 and 3 is well-reflected by the difference in the underlying NTOs highlighting the consistency between aromaticity descriptors and the MO picture. From a more methodological point of view we want to emphasise that the involvement of two NTO pairs in the state implies multiconfigurational character of the triplet state [72,74] meaning that the single determinant description provided by UKS is not able to cover all relevant aspects of the wavefunction. Nonetheless, Fig. 6 clearly shows the differences in electronic structure between 2 and 3.
Proceeding to 4, shown in Fig. 7, we find a notably different behaviour when compared to the isomeric 3. Viewing the shielding tensors in the ground state, Fig. 7 (a), we find enhanced antiaromaticity at the CBD ring showing deshielding of +61.1 ppm, which is even higher than for isolated CBD [ Fig. 4 (b)]. This enhanced antiaromaticity agrees well with the higher relative energy as shown in Tab. 2. Both benzene rings come out as clearly shielded, albeit with lower values than isolated benzene. Proceeding to the triplet state [ Fig. 7 (b)], we find a strongly aromatic CBD ring (-28.6 ppm) similarly to biphenylene. This is accompanied by a weakly shielded intermediate benzene ring and a strongly shielded (-25.1 ppm) outer benzene ring. Reviewing the shielding tensors, these can be seen to reflect the sextet-quartet-sextet structure indicated in Fig. 2. It is precisely this feature that has been associated with enhanced stability for triplet excited states [41]. Viewing the NTOs of the S 0 → T 1 transition, we find that, as opposed to 3, there is only one dominant NTO pair already making up 93% of the overall transition. Moreover, this NTO pair closely resembles the one in biphenylene. In particular there is a strong contribution between the quasi-degenerate CBD frontier orbitals, which is linked with the Baird aromaticity of this compound, as discussed above.
Proceeding to 5, shown in Fig. 8, we find a largely similar structure in terms of, both, the shielding tensors and NTOs as for 4 and, again, entirely different properties to 3. This highlights that the qualitative considerations in terms of quartets and sextets, as presented in Fig. 2, are a powerful way to estimate the properties of these molecules and that other details in the bonding pattern do not affect the electronic structure in an important way.

Discussion
Before concluding, we want to discuss the relation between the chemical shielding tensors used above and other common ways of describing the electronic structure of PAHs. In the above presentation, we have shown that Baird aromaticity provides a powerful way of rationalising variations in the triplet energies of the molecules studied. In addition, we have highlighted that the energies can also be explained within the MO picture following Fig. 3 and by realising that the orbitals involved do indeed possess the shapes expected from the rigid rotor model. Using a different viewpoint, Ryerson et al. [27] have pointed out that also a biradical model is highly suitable for discussing the energetics in related systems. It should, therefore, be stressed that Baird aromaticity is just one possible language of describing the phenomena seen. Nonetheless, it is the strength of the aromaticity rules that they allow qualitative predictions based on simply counting the number of electrons and analysing bonding patterns as exemplified in the striking differences between molecules 3 vs 4/5.
From a more formal viewpoint it is worth pointing out that there is a specific hierarchy between frontier orbitals, biradical character/NTOs, and the shielding tensors. Frontier orbitals are simply intermediates in approximate theories and it is well-known that their shapes and energies change dramatically with the level of theory chosen [47,75]. Measures for biradical character and the NTOs can be defined based on the wavefunctions alone [24,66,76] and are, thus, well-defined without any explicit reference to a computational level of theory. However, out of this list, only the chemical shielding tensors are a physically observable quantity. The chemical shielding at the atomic positions is of course routinely measured in 1 H and 13 C NMR experiments. Even more, it is possible to approximate the NICS via NMR experiments using complexated 6 Li [77] and, more generally, the shielding defined as a ratio of two magnetic fields [Eq. (1)] is always at least in principle an observable.
More generally, we suggest using aromaticity as one viewpoint in the elucidation of excited-state electronic structure in a similar sense as exciton theory and valence bond theory can be helpful viewpoints in elucidating otherwise hidden properties of excited states [74,78]. However, it should be understood that these are just different ways of describing the same underlying physics and that Baird aromaticity is not somehow detached from other ways of describing excited-state electronic structure.
On a different note, it should be pointed out that the practice of computing only one NICS value per ring in multi-ring systems, as done above, has been contested [50]. We believe that for the purposes of the current study, i.e. comparing the electronic structure of closely related systems and illustrating the excitation processes, this practice provides an intuitive and representative account of the underlying electronic structure. However, for a more detailed and transferable representation, it may be necessary to compute more fine-grained NICS-XY-scans [79].

Conclusions
Within this work a class of PAHs based on the biphenylene motif was studied showing that the bond topology had dramatic effects on their triplet excitation energies. Strikingly, a difference of more than 1.5 eV in adiabatic excitation energies was found for 3 when compared to isomeric 4 and 5. These results were rationalised in the context of Baird aromaticity using a qualitative description of Clar sextets and Baird quartets, quantified in terms of NICS values, and graphically represented using a recently developed method for visualising chemical shielding tensors (VIST). A comparison with natural transition orbitals showed a consistent picture highlighting that Baird aromaticity should not be seen as a detached phenomenon but as one possible way to categorise excited states. More generally, we believe that excited-state aromaticity constitutes a fascinating and important phenomenon and hope that this work will stimulate future studies by showing how it can be readily visualised and how it integrates into the framework of excited-state electronic structure theory.
Funding: This research was funded by the EPSRC, grant number EP/V048686/1.

Conflicts of Interest:
The author declares no conflict of interest.