Exchange rules for diradical {\pi}-conjugated hydrocarbons

TL;DR

This paper develops rules to accurately predict zero mode states, spin multiplicity, and shell nature of planar pi-conjugated hydrocarbons, extending beyond traditional theorems using analytical, Hubbard model, and quantum chemistry methods.

Contribution

It introduces a comprehensive set of rules for predicting electronic states and spin configurations in diradical pi-conjugated hydrocarbons, surpassing existing theorems.

Findings

Rules accurately predict zero mode states and spin multiplicity.

Analytical and computational methods validate the predictions.

Extends understanding beyond Lieb's theorem and Ovchinnikov's rule.

Abstract

A variety of planar {\pi}-conjugated hydrocarbons such as heptauthrene, Clar's goblet and, recently synthesized, triangulene have two electrons occupying two degenerate molecular orbitals. The resulting spin of the interacting ground state is often correctly anticipated as S = 1, extending the application of Hund's rules to these systems, but this is not correct in some instances. Here we provide a set of rules to correctly predict the existence of zero mode states, as well as the spin multiplicity of both the ground state and the low-lying excited states, together with their open- or closed-shell nature. This is accomplished using a combination of analytical arguments and configuration interaction calculations with a Hubbard model, both backed by quantum chemistry methods with a larger Gaussian basis set. Our results go beyond the well established Lieb's theorem and Ovchinnikov's rule, as we address the multiplicity and the open-/closed-shell nature of both ground and excited states.

Tables

Table 1. Table 1: Hubbard integrals for the structures from Fig.1

Structure	${𝒰}_1/\overline{𝒰}$	${𝒰}_2/\overline{𝒰}$	$J/\overline{𝒰}$	$t_{12}/\overline{𝒰}$	$\Delta /\overline{𝒰}$	$\overline{𝒰}\left(U\right)$	$\lambda \left(\widehat{ℛ}\right)$
Triangulene	1	1	2	0	0	0.096	$e^{\pm \frac{i2\pi }{3}}$
Clar’s goblet	1	1	0	0	0	0.144	$1$
Heptauthrene	1.220	0.779	1.209	0	0.604	0.106	$\pm 1$

Equations

$${\mathcal{H}}_0\Gamma +\Gamma {\mathcal{H}}_0=0$$

$$⟨{\psi }_n|\Gamma |{\psi }_n⟩=0\left(E_n\ne 0\right)$$

$$Tr\Gamma =\sum _i⟨i|\Gamma |i⟩=N_A-N_B$$

$$Tr\Gamma =\sum _n⟨{\psi }_n|\Gamma |{\psi }_n⟩+\sum _{A}1-\sum _{B}1$$

$${\mathcal{H}}_{\mathcal{U}}$$

$${\mathcal{U}}_{\eta }=U\sum _i|{\varphi }_{\eta }(i){|}^4$$

$$J=2U\sum _i|{\varphi }_1(i){|}^2|{\varphi }_2(i){|}^2$$

$${\mathcal{V}}_{12}=\sum _{\sigma }{n}_{1\sigma }\left(t_{12}{C}_{1\overline{\sigma }}^{†}{C}_{2\overline{\sigma }}+t_{12}^{\ast }{C}_{2\overline{\sigma }}^{†}{C}_{1\overline{\sigma }}\right)$$

$${\mathcal{V}}_{\mathrm{pair}}=\Delta C_{1\uparrow }^{†}{C}_{1\downarrow }^{†}{C}_{2\uparrow }{C}_{2\downarrow }+{\Delta }^{\ast }{C}_{2\uparrow }^{†}{C}_{2\downarrow }^{†}{C}_{1\uparrow }{C}_{1\downarrow }$$

$$t_{12}=U\sum _i|{\varphi }_1(i){|}^2{\varphi }_1(i{)}^{\ast }{\varphi }_2(i)$$

$$\Delta =U\sum _i({\varphi }_1(i{)}^{\ast }{)}^2{\varphi }_2(i{)}^2$$

$$\hat{\mathcal{R}}{\varphi }_{1,2}={\lambda }_{1,2}{\varphi }_{1,2}$$

$${\mathcal{U}}_{\eta }$$

$$J$$

$$t_{12}$$

$$\Delta$$

$$E_T$$

$$E_{S1}$$

$$E_{S2}$$

$$E_{S3}$$

Full text

I Introduction

Understanding the so-called exchange interactions in atoms, molecules and crystals is one of the central topics in the study of the electronic properties of matter. In the case of atoms, whether or not a given atom has a magnetic ground state is a long settled issue: only open-shell atoms can have a magnetic moment whose magnitude is established by the so-called Hund’s rules. The nature of inter-atomic exchange is also well understood in a wide class of oxides, since the seminal work of Goodenoughgoodenough63 and Kanamorikanamori59 .

The study of $\pi$-conjugated hydrocarbons with diradical nature (which also includes polycyclic aromatic hydrocarbons (PAHs) or nanographenes) goes back to more than a century agoSchlenk ; Borden77 ; rajcareview ; moritaReview . Experimental research in this area has faced the challenge of the very high chemical reactivity of diradicals. However, recent developments of new synthesis methods, in ultrahigh vacuum surfaces, has made it possible to synthetize highly reactive diradical species, such as triangulenepavlivcek2017 or graphene nanoribbons with zigzag edgeswang2016 as well as other open-shell PAHs NachoNat , and to explore them using scanning probe microscopies.

Diradicals host two degenerate, or almost degenerate, states that lie in the gap between the doubly occupied and empty molecular orbitals that are formed primarily by $\pi$ orbitalsYoneda2009 . Very often they are localized in one of the two interpenetrating triangular sublattices that form the honeycomb bipartite lattice. Naively, one could expect that diradicals always had a ground state with $S=1$, due to the similitude with the case of open-shell atoms. Whereas this is often the case in a variety of systems melle2015 ; Schlenk ; gryn'ova2015 ; Lischka2016 ; Sandoval2019 ; Yoneda2011 , in others there is a violation of the Hund’s ruleYazyev09 ; gryn'ova2015 ; Sheng2013 ; Trinquier2015 ; Malrieu2016 ; Trinquier2018 ; Pogodin1985 . Therefore, the sign of the exchange interaction in this class of diradicals is not always the samerajcareview .

In general, the spin of the ground state in PAHs can be anticipated using the Ovchinnikov’s ruleovchinnikov1978 , that states that the spin $S$ of the ground state is given by $S=\frac{N_{A}-N_{B}}{2}$, where $N_{A,B}$ are the number of C atoms in each of the interpenetrating triangular sublattices that form the honeycomb lattice. Interestingly, Lieb upgraded the Ovchinnikov’s rule into a theoremLieb89 , which states that the exact interacting ground state of the Hubbard model for a bipartite system is given also by $S=\frac{N_{A}-N_{B}}{2}$. Yet another theoremsutherland86 establishes that for a bipartite lattice the number of zero modes is given by $N_{Z}=|N_{A}-N_{B}|$. In consequence, systems with a sublattice imbalance of $N_{A}=N_{B}\pm 2$, such as triangulene, are diradicals with $S=1$.

However, this picture is not complete for several reasons. First, some diradicals, such as the Clar’s goblet or the undistorted (planar) cyclobutadiene have no sublattice imbalance and still have two zero modes and, in agreement with the Ovchinnikov’s rule and Lieb’s theorem, have $S=0$. Second, whereas Lieb’s theoremLieb89 has predictive power on the spin $S$ of the ground state, it has a very poor explaining power: what is the nature of the antiferromagnetic exchange in those diradicals with a GS with $S=0$? Third, the Lieb’s theorem does not provide any information on spin of the excited states, nor on their open-shell vs closed-shell nature.

In this work we address these longstanding research questions in the case of planar conjugated hydrocarbons or nanographenes (as well as its generalization to all kind of diradical $\pi$-conjugated hydrocarbons in the suppl. mat.). First, we provide a unified approach to anticipate if a given bipartite system has zero mode states. This constitutes a prerequisite for open-shell configurations. Second, we provide a set of rules that determine the multiplicity and the open-/closed- shell nature of the lowest in energy many-body states for these nanographene diradicals. We find that the key ingredient that defines the properties in the multi-electronic problem are the transformation properties of the molecular orbitals under the symmetry operations of the point group of the molecule. We therefore establish a set of simple rules that permit to anticipate the relative position and approximate excitation energy of the six lowest energy multi-electron states, a degenerate triplet and three singlets.

Our results are based on an analysis using theory at three levels of complexity. First, an analytical description of the Hubbard model for these compounds when the active space is restricted to the two in-gap states only. Second, numerical calculations, still within the Hubbard model, in a larger active space that includes occupied and virtual molecular orbitals. Third, ab initio quantum chemistry methodscusinato2018 , including complete active space (CAS) configuration interaction (CI) calculations, carried using a gaussian orbital basis, and further corrected by second-order N-electron Valence Perturbation Theory (NEVPT2).

II Zero modes (non-bonding states).

We first revisit the problemlonguet50 ; Borden77 of predicting the existence of non-bonding zero modes in a PAH. Our discussion applies when the system is described with a tight-binding model with one $p_{z}$ orbital per atom in the first neighbor hopping approximation. This defines a tight-binding model in a bipartite graph and permits to provide a unified and compact picture that accounts for a number of results derived over the yearssutherland86 ; inui1994 ; Fajtlowicz2005OnMM ; Fernandez2007 ; weikprb2016 . Bipartite graphs can be drawn as the superposition of two interpenetrating sublattices (see Fig.1). We define the $\Gamma$ matrix that assigns a $+1$ to the carbon sites of the $A$ sublattice and a $-1$ to the carbon sites of the $B$ sublattice.

The Hückel (tight-binding) Hamiltonian ${\cal H}_{0}$ that describes the $\pi$ orbitals has the so-called chiral symmetryRyu02 ; delplace2011 , given by the equation:

[TABLE]

This anti-commutation relation, different from usual commutators associated to symmetries, entails several consequences relevant for the ensuing discussion. First, let us consider the eigenstates ${\cal H}_{0}|\psi_{n}\rangle=E_{n}|\psi_{n}\rangle$. Eq.(1) implies that $|\psi^{\prime}_{n}\rangle\equiv\Gamma|\psi_{n}\rangle$ is also an eigenstate of ${\cal H}_{0}$ with eigenvalue $-E_{n}$. If $E_{n}\neq 0$, $|\psi_{n}\rangle$ and $|\psi^{\prime}_{n}\rangle$ have to be orthogonal, as they are eigenstates with different eigenvalues. Their orthogonality can be written up as:

[TABLE]

which implies an equal weight on the two sublattices.

Hamiltonians ${\cal H}_{0}$ with chiral symmetry can also host states with $E_{n}=0$, that we label as $\phi_{\eta}$. It can be easily demonstrated (see suppl. mat.) that these zero modes are sublattice polarized and, therefore, are eigenstates of $\Gamma$. As a result, they satisfy $\langle\phi_{\eta}|\Gamma|\phi_{\eta}\rangle=\pm 1$, depending on whether they are localized in the $A$ or $B$ sublattice. Thus, Eq.(1) classifies the eigenstates of ${\cal H}_{0}$ in three groups, according to the expectation value of $\Gamma$, that can be [math], for bonding and antibonding states, and $+1$ or $-1$ for sublattice polarized non-bonding states.

With this background, we now provide a straightforward demonstration of the well known resultLonguetHiggings ; Borden77 ; sutherland86 ; inui1994 that for a given bipartite graph with $N_{A}-N_{B}\neq 0$ there are, at least, $N_{Z}=|N_{A}-N_{B}|$ zero modes. For that matter, we first compute the trace of the $\Gamma$ operator, calculated in the atomic orbital basis:

[TABLE]

where $i$ labels the atomic sites. We now write down the trace in the basis of eigenstates of ${\cal H}_{0}$. We break it down in terms of three types, namely, states with $E_{n}\neq 0$ that do not contribute to the trace, and $A$ or $B$ zero modes, that contribute with $\pm 1$:

[TABLE]

Comparison of Eq.(3) and (4) establishes $N_{zA}-N_{zB}=N_{A}-N_{B}$, where $N_{zA/zB}$ is the number of zero modes localized in a given sublattice. Importantly, Eq.(4) gives the minimum number of zero modes, as we can certainly have $N_{zA}+N_{zB}>N_{A}-N_{B}$. Thus, in this work we consider systems, such as the so-called triangulene and heptauthrene (Fig.1a, c), that have sublattice imbalance of $2$ and two zero modes in the majority sublattice (Fig.2 a, c, d, f) and also systems with two zero modes and $N_{A}=N_{B}$ (Fig.2b, e).

We have found two different classes of situations in which systems with $N_{A}=N_{B}$ have two zero modes:

• •

When there is a symmetry operation $\hat{\cal R}$ that conmutes both with ${\cal H}_{0}$ and $\Gamma$, such that the different irreducible representations of $\hat{\cal R}$ have a partial trace of $\Gamma$ different from zeroKoshino . This is the case of $4n$ anulenes, such as cyclobutadiene. Importantly, Jahn-Teller mechanism is operative here, and a structural distortion of the PAH can lift the degeneracy of the zero modes.

• •

When the structure is formed by fragments that, when considered separately, have zero modes whose wavefunction have no weight on the junction region Borden77 . This is the case of Clar’s gobletclar72 (Fig.2b, e). Here the zero modes are robust with respect to Jahn-Teller distorsions.

The fact that, for $N_{A}=N_{B}$, the two zero modes are hosted in different sublattices has important consequences in the many-body Hamiltonian, as we discuss below.

III PAH diradicals: a Hubbard model description.

III.1 Hubbard model in a minimal Active Space.

The discussion of the effect of electron-electron interactions starts with the description of the Hubbard model that, in the context of nanographenes, has been shownFernandez2007 ; Fernandez08 to give results very similar to those of Density Functional Theory calculations. We start with the minimal Hilbert space for two electrons in two molecular orbitals, $\phi_{1}$ and $\phi_{2}$, that are linear combination of the $\pi$ atomic orbitals of the PAH’s C atoms. These molecular orbitals are obtained by diagonalization of a single-particle Hamiltonian, ${\cal H}_{0}$, that describes first neighbor single orbital tight-binding model. The resulting spectra and molecular orbitals for a variety of nanographenes are shown in Fig.2. All of them have two zero modes that are occupied with two electrons at half filling. When interactions are ignored, the two-fold orbital degeneracy of the in-gap states results in a six-fold degeneracy of the multi-electronic ground state, that is lifted by Coulomb interactions.

When treated in the Hubbard approximation, the interaction term is simply given by $U\sum_{i}n_{i\uparrow}n_{i\downarrow}$, where $n_{i\sigma}$ stands for the occupation of the atomic $\pi$ orbital with spin $\sigma$ at carbon $i$. In contrast, it can be seen that the same expression, projected over the two in-gap molecular orbitals results in a Hamiltonian with four different terms (see suppl. mat.):

[TABLE]

where $n_{\eta\sigma}=C^{\dagger}_{\eta\sigma}C_{\eta\sigma}$ is the occupation operator of the molecular orbital $\phi_{\eta}$, $\vec{S}_{\eta}=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}C^{\dagger}_{\eta\sigma}\vec{\tau}_{\sigma,\sigma^{\prime}}C_{\eta\sigma^{\prime}}$ is the spin operator associated to that orbital (see suppl. mat.) and $C^{\dagger}_{\eta\sigma}$ is the operator that creates and electron with spin $\sigma$ in the molecular orbital $\phi_{\eta}$.

Hamiltonian (5) has four types of interactions. First, Hubbard-like terms, with energy ${\cal U}_{\eta}-\frac{J}{4}$, that describe the Coulomb energy penalty of double occupation of states $\phi_{1}$ and $\phi_{2}$. The energy scales ${\cal U}_{\eta}$ are given by:

[TABLE]

where $\sum_{i}|\phi_{\eta}(i)|^{4}$ is the so-called inverse participation ratio and is a metric of the extension of the orbitalOrtiz2018 .

The exchange integral is given by:

[TABLE]

and is a metric of the overlap of the two zero modes. It is apparent that $J=0$ for disjoint zero modes.

Second, a ferromagnetic exchange term between the spins in orbitals (notice that $J\geq 0$). The last two terms are the density assisted hopping terms,

[TABLE]

and the pair hopping term:

[TABLE]

The Hubbard matrix elements that control these Coulomb assisted hoppings are:

[TABLE]

and

[TABLE]

Both $t_{12}$ and $\Delta$ can be complex numbers.

At this point, we make two crucial observations. First, since $\phi_{1}$ and $\phi_{2}$ are degenerate, there is not a unique representation for them and the values of the Hubbard integrals ${\cal U}_{\eta}$, $J$, $\Delta$ and $t_{12}$ depend on our choice of representation. In the following we choose $\phi_{1}$ and $\phi_{2}$ so that they diagonalize a symmetry operator of the point group of the molecule. This permits to obtain closed expressions for the spectrum of Hamiltonian (5). To do so, the first step is to find a symmetry operator $\hat{\cal R}$ that conmutes with the single particle Hamiltonian, $[\hat{\cal R},{\cal H}_{0}]=0$ and with the sublattice operator $\Gamma$. We choose $\phi_{1}$ and $\phi_{2}$ such that:

[TABLE]

In table I the relevant symmetry operators and the eigenvalues $\lambda_{1,2}$ are listed for the three systems of interest. Symmetries include a 120∘ rotation in the case of triangulene, and reflection around the mirror symmetry axis in the case of Clar’s goblet and heptauthrene.

Our second central observation is the fact that Hubbard integrals, defined in eqs.(6, 7, 10, 11) have to remain invariant under the symmetry operation. Hence, if we replace $\phi_{1}$ and $\phi_{2}$ by ${\cal\hat{R}}\phi_{1}=\lambda_{1}\phi_{1}$ and ${\cal\hat{R}}\phi_{2}=\lambda_{2}\phi_{2}$, we have:

[TABLE]

Thus, the prefactor in the right hand side of Eq.(13) have to be identical to one, otherwise the corresponding Hubbard integral vanishes. As it can be inferred from table I, the Hubbard integral $t_{12}=0$ for $C_{3}$ and heptauthrene systems, out of symmetry. Obviously, $t_{12}$ also vanishes for all diradicals with disjoint zero modes. When $t_{12}=0$, we obtain analytical expressions for the six eigenvalues of Eq.(5):

[TABLE]

where $\overline{{\cal U}}\equiv\frac{1}{2}\left({\cal U}_{1}+{\cal U}_{2}\right)$. The first line corresponds to the triplet. The energy $E_{S1}$ corresponds to the open-shell singlet, and $E_{S2}$ and $E_{S3}$ to closed-shell singlets. These equations constitute our first important result.

III.2 Hubbard model Reduced active space: theory

In the following we apply this theory to three different types of diradicals, with and without sublattice imbalance and with different point group symmetries. In all cases, we determine the open-/closed-shell of the multielectronic states by the overlap of the multielectronic wave function with the configurations with a well defined occupation of the single particle states, always taking the zero modes as eigenstates of a symmetry operator.

III.2.1 $C_{3}$ diradicals: Hubbard model.

We consider first the case of systems with $N_{A}-N_{B}=2$ and $C_{3}$ symmetry, such as trianguleneFernandez2007 and trimethylenemethane (see suppl. mat.). In this case, the relevant symmetry is the in-plane 120*∘* rotation, $\hat{\cal R}(\frac{2\pi}{3})$. The diagonal representation of $\hat{\cal R}(\frac{2\pi}{3})$ in the subspace of zero modes has two eigenvalues $e^{\pm i\epsilon}$ with $\epsilon=\frac{2\pi}{3}$. Therefore, given that the rotation matrix is real, we can easily see that $\phi_{1}(i)=\left(\phi_{2}(i)\right)^{*}$, i.e., overlap between zero modes is maximal. From here, we obtain ${\cal U}_{1}={\cal U}_{2}=\frac{J}{2}$, which implies that exchange is the dominant energy scale and $\Delta=0$ and $t_{12}=0$.

The resulting spectrum, shown in Fig.3a, has a ground state triplet, in agreement with Lieb’s theorem, with energy $E_{T}=-\frac{J}{4}$. The three excited singlets are arranged in a low energy closed-shell doublet, with energies $E_{S2}=E_{S3}=\frac{J}{4}$, and a high energy open-shell singlet with energy $E_{S1}=3\frac{J}{4}$. The spectrum of excited states of the triangulene is quite peculiar for two reasons. First, the open-shell singlet is the highest energy excitation, within the manifold of states considered in this approximation. This reflects the large value of exchange that arises from the maximal overlap of the zero modes. Second, the closed-shell singlets are degenerate.

III.2.2 Diradicals with $N_{A}-N_{B}=2$ and reflection symmetry.

We now discuss the case of PAHs with a different point group. We consider the heptauthrene diradical, that has a reflection symmetry that preserves sublattice. The eigenvalues of this symmetry operator are $\pm 1$; thus, its application leaves all the Hubbard integrals unaffected except for $t_{12}$ and $t_{21}$ that change sign and, therefore, must vanish identically. Since the eigenvalues of the symmetry operator are real, the states $\phi_{1}$ and $\phi_{2}$ are also real, which automatically entails $J=2\Delta$. However, we have ${\cal U}_{1}\neq{\cal U}_{2}$. The resulting many-body spectrum, within the reduced active space, is shown in Fig.3c. The ground state has $S=1$, complying with the Lieb’s theorem. The excited states follow a more conventional arrangement, with three non-degenerate singlets, with the open-shell singlet in between the two closed-shell singlets.

III.2.3 $N_{A}=N_{B}$ diradicals.

We now apply Hamiltonian (5) to the case of diradicals with $N_{A}=N_{B}$ such as the Clar’s goblet clar72 ; Yazyev09 ; Fajtlowicz2005OnMM and the cyclobutadieneKoshino ; SchumannR (see suppl. mat.). The single-particle zero modes of diradicals with $N_{A}=N_{B}$ are disjoint, i.e., they are located in different atomsKenneth . As a result, the only non-zero energy scales of the restricted many-body Hamiltonian (5) are ${\cal U}_{1,2}$. The resulting interacting spectrum (Fig.3b), in the minimal active space, presents a quartet ground state, formed by the open-shell $S=0$ and $S=1$ states, and two closed-shell excited $S=0$ states. Extension of the Hilbert space is thus necessary to resolve the singlet-triplet energy difference, that must favour a $S=0$ ground state to comply with Lieb’s theoremLieb89 .

IV Hubbard model calculations with an extended Hilbert space.

The description of the diradicals using a restricted Hilbert space sometimes gives spectra with peculiar degeneracies, such as the first excited state of the $C_{3}$ diradicals or the ground state degeneracy of the $N_{A}=N_{B}$ diradicals. We now explore whether these degeneracies are artefacts of the truncated Hilbert space, or they are byproducts of symmetry. For that matter, we have carried out exact diagonalizations using a configuration interaction method for the Hubbard model, extending the active space to include valence and conduction states. We denote a complete active space (CAS) with $\cal N$ electrons and $\cal M$ molecular orbitals, each with a two-fold spin degeneracy, as CAS($\cal N$, $\cal M$). For ${\cal N}={\cal M}=4$ we have a total of 70 many-body states.

Our results for CAS(4,4) for the three systems of interest are shown in the second row of Fig.(3). In the case of triangulene, the main peculiar properties of the spectrum obtained in the minimal model are preserved here. This is also the case of trimethylenemethane (see suppl. mat.), for which the CAS(4,4) calculation covers the complete Hilbert space. We thus conclude that the peculiar symmetry obtained in the analytical CAS(2,2) model is not a specific feature of a truncated active space. Our results for heptauthrene, using CAS(4,4) are also in line with the analytical model.

Unlike the analytical model, our CAS(4,4) calculations give non-degenerate $S=0$ ground state for the diradicals with $N_{A}=N_{B}$, such as the Clar’s goblet and the cyclobutadiene. This complies with Lieb’s theoremLieb89 . The lowest excited state is a triplet. We have verified (see suppl. mat.) that the singlet-triplet splitting, scales as $E_{S=1}-E_{S=0}\propto U^{2}/t>0$. This scaling indicates that the antiferromagnetic exchange is driven by virtual excursions to excited states, with electrons in the conduction states and/or holes in the valence states, driven by the Hubbard term. It must be noted that this is different from the usual kinetic exchangeAnderson59 scaling $t^{2}/U$, that has also been discussed for the case of in-gap zero modes hybridized by hopping in nanographenesOrtiz2018 .

V Ab initio calculations.

At this point we further explore to which point the predictions of the Hubbard model are different from more sophisticated yet costly ab initio calculations including nuclear coordinates relaxation, more atomic orbitals per atom, H atoms passivating edge C atoms, and long-range Coulomb interactions. The structures and molecular orbital basis set used to carry out the multi-configurational approaches are obtained using a Density-Functional Theory calculation based with the hybrid exchange-correlation functional B3LYP (see suppl. mat.).

We have carried out ab initio calculations for three different diradicals: triangulene, Clar’s Goblet and heptauthrene (see Fig.3). In all three cases, the ground state and first excited state have the same spin than anticipated in the previous sections. In the case of triangulene, the symmetry of the spectrum also comes in the 3-2-1 sequence obtained using the Hubbard model. The relative value of the energy difference between triplet and open-shell singlet is also in line with the Hubbard model: largest for triangulene, and much smaller for the Clar’s goblet.

In the case of the Clar’s goblet, the ab initio calculations predict a second excited triplet with energy smaller than that of the two singlets, in contrast with the Hubbard model predictions. This discrepancy is probably related to the underestimation of the energy overhead associated to the double occupancy of a molecular orbital in the closed-shell configuration when electronic repulsion is treated in the Hubbard approximation. As a result, the Hubbard model underestimates the energy of the two closed-shell singlets, compared to the second open-shell triplet. Barring this discrepancy in the higher order excited states of the Clar’s goblet, the main features of the Hubbard model are confirmed with the more sophisticated, but much more computationally expensive, ab initio methods.

VI Exchange rules.

We are now in position to conclude with a number of rules for the exchange interactions in PAHs diradicals, based on the analysis of the results obtained with various methods:

Lieb’s rule: the spin of the ground state, $S$, is determined by the sublattice imbalance, $|N_{A}-N_{B}|=2S$. 2. 2.

The spin and orbital degeneracies of the six lowest energy states is predicted, qualitatively, by Hamiltonian (5). The eigenvalues energy scales are governed by ${\cal U}_{1,2}$, $J$, $\Delta$ that, in turn, are strongly conditioned by the point symmetry group of the PAH diradical. 3. 3.

Ferromagnetic exchange $J$ is maximal for $C_{3}$ diradicals, on account of the maximal overlap of the zero modes. The $C_{3}$ symmetry also imposes that the lowest energy excited state turns out to be given by a degenerate pair of closed-shell states. 4. 4.

Ferromagnetic exchange $J$ is minimal for $N_{A}=N_{B}$ diradicals, on account of the disjoint nature of their zero modes.

We acknowledge J. L. Lado for fruitful discussions. J. F.-R. and R. O. acknowledge financial support from MINECO-Spain (Grant No. MAT2016-78625-C2) and from the Portuguese “Fundação para a Ciência e a Tecnologia” (FCT) for the project P2020-PTDC/FIS-NAN/4662/2014. J.F.-R., M. M.-F. and N. G.-M. acknowledge support from the P2020-PTDC/FIS-NAN/3668/2014. J. F.-R. acknowledges support from UTAPEXPL/NTec/0046/2017 projects, as well as Generalitat Valenciana funding (Prometeo2017/139). R. O. and J. C. S.-G. acknowledge ACIF/2018/175 (Generalitat Valenciana and Fondo Social Europeo). M. M.-F. and R. B. would like to acknowledge the Portuguese “Fundação para a Ciência e a Tecnologia” (FCT) for the project IF/00894/2015 and FCT Ref. UID/CTM/50011/2019 for CICECO - Aveiro Institute of Materials. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 664878.