Broken-symmetry magnetic phases in two-dimensional triangulene crystals

We provide a comprehensive theory of magnetic phases in two-dimensional triangulene crystals, using both Hubbard model and density functional theory (DFT) calculations. We consider centrosymmetric and non-centrosymmetric triangulene crystals. In all cases, DFT and mean-field Hubbard model predict the emergence of broken-symmetry antiferromagnetic (ferrimagnetic) phases for the centrosymmetric (non-centrosymmetric) crystals. This includes the special case of the [4,4]triangulene crystal, whose non-interacting energy bands feature a gap with flat valence and conduction bands. We show how the lack of contrast between the local density of states of these bands, recently measured via scanning tunneling spectroscopy, is a natural consequence of a broken-symmetry N\'eel state that blocks intermolecular hybridization. Using random phase approximation, we also compute the spin wave spectrum of these crystals, including the recently synthesized [4,4]triangulene crystal. The results are in excellent agreement with the predictions of a Heisenberg spin model derived from multi-configuration calculations for the unit cell. We conclude that experimental results are compatible with an antiferromagnetically ordered phase where each triangulene retains the spin predicted for the isolated species.


I. INTRODUCTION
Triangulenes are graphene fragments with the shape of an equilateral triangle, terminated with zigzag edges and of various sizes, customarily defined in terms of the number n of benzenes in a given edge [1][2][3] .According to single-particle theory, [n]triangulenes host n − 1 nonbonding half-filled zero modes 2 .Coulomb interactions favor the maximal spin configuration, very much like the Hund's first rule in atoms, so that [n]triangulenes are predicted 2,[4][5][6][7] to have a ground state with total spin S = n−1 2 (see Fig. 1a), consistent with Lieb's theorem for the Hubbard model for bipartite lattices at half-filling 8 , and in agreement with Ovchinnikov's rule 9 .
Using inelastic electron tunneling spectroscopy 21 , zerobias Kondo resonances in individual [2]triangulenes 16 , as well as spin excitations in [3]triangulene dimers 17 , rings 18,19 and chains with more than 40 units 18 have been observed.These experiments provide strong evidence that these zero-and one-dimensional supramolecular structures remain open-shell and their low-energy electronic properties can be accounted for by spin Hamiltonians with antiferromagnetic interactions (Fig. 1b).
Spin-restricted density functional theory (DFT) calculations of [n, m]triangulene crystals-i.e., honeycomb 2D crystals whose unit cell is made of a pair of triangulenes with sizes n and m-show the formation of n + m − 2 weakly dispersive energy bands 22 .Using tightbinding models, it has been shown 22 that these bands are made of linear combinations of the in-gap zero modes of the triangulenes, hybridized via third-neighbor hopping.Intermolecular hybridization splits the zero modes into bonding-antibonding pairs, promoting non-magnetic closed-shell electronic configurations.Therefore, in contrast with the case of isolated triangulenes, interactions need to overcome intermolecular hybridization in order to promote open-shell states.This is expected to be harder in the case of the [4, 4]triangulene crystal, for which both spin-restricted DFT 23 and tight-binding calculations predict a narrow-gap insulator, unlike the [2, 2] and [3, 3]  cases, that feature Dirac cones at the Fermi energy.The synthesis of a [4, 4]triangulene 2D lattice has been recently reported 20 , putting this specific system under the spotlight.
In this work, we undertake a systematic study of the electronic properties of triangulene 2D crystals, focusing on the magnetic properties of their ground states.To do so, we go beyond the spin-restricted framework in the case of DFT, and beyond non-interacting tight-binding models.For that matter, we take the natural next step, doing spin-unrestricted DFT calculations and adding a Hubbard term to the tight-binding model used in previous work.The Hubbard model is treated at three levels of approximation: collinear mean-field theory, random phase approximation (RPA) and exact diagonalization of small structures in a restricted space of configurations, FIG. 1.(a) [4]triangulene has a fourfold degenerate ground state with total spin S = 3/2.(b) [4]triangulene dimer is an open-shell singlet with an entangled wave function as a result of antiferromagnetic intermolecular coupling.(c) Two examples of broken-symmetry Néel states for [4,4]triangulene 2D crystals, obtained with a collinear mean-field Hubbard model.The size of the circles represents the magnitude of the local moments, with red/blue colors denoting spin-↑/↓.Total spin is no longer a good quantum number.
the so-called complete active space (CAS) method.
Previous spin-unrestricted DFT calculations predict that [2,2]-(ref. 24) and [3,3]-(ref. 22) triangulene 2D crystals should display antiferromagnetic order, with the two sublattices being polarized in opposite directions.The [4,4]triangulene crystal is different from [2,2] and [3,3] as it features a small band-gap when calculated both with spin-restricted DFT 20,22,23 and with the conventional single-orbital tight-binding model with thirdneighbor hopping 22 .On the basis of this narrow gap, an excitonic insulating phase has been proposed 23 , taking as a reference-state the closed-shell non-magnetic ground state.A major goal of this manuscript is to address whether the [4,4]triangulene crystal also features an antiferromagnetic phase (Fig. 1c), and how this affects the size of the gap and the putative excitonic insulator.
The rest of the paper is organized as follows.In section II we review the different theoretical methods used in this work.In section III we present our results for triangulene dimers within the CAS approach for the Hubbard model.These calculations allow us to derive the effective spin interactions, in the form of polynomials of the Heisenberg coupling, and to estimate the magnitude of the intermolecular exchange couplings.In section IV we present the results of collinear mean-field Hubbard and spin-unrestricted DFT calculations.The results are very similar, validating the Hubbard model, and systematically predict broken-symmetry magnetic phases as the ground state of triangulene 2D crystals.In section V we present our Hubbard model RPA calculations of the spin waves for the [2,3], [3,3] and [4,4] crystals, and compare them with those obtained from the spin models derived from the Hubbard model CAS calculations.In section VI we discuss how the lack of contrast in the local density of states (LDOS) of the conduction and valence bands can be used to identify the emergence of broken-symmetry states, providing an explanation to recent experimental scanning tunneling microscope (STM) spectroscopy results 20 .In section VII we present the conclusions.

II. METHODS
In this section we provide a brief description of the theoretical methods used throughout the paper.

A. Hubbard model
Following previous work 2, [25][26][27][28] , we use a single-orbital Hubbard model to describe π-magnetism in graphene nanostructures.The Hubbard model 29 can be written as where the indices i, j run over carbon atoms, t i,j stands for the hopping between sites i and j, U is the on-site Hubbard repulsion, c † iσ (c iσ ) denotes the operator that creates (annihilates) an electron in site i with spin projection along a quantization axis σ =↑, ↓, and n iσ = c † iσ c iσ is the corresponding number operator.While the first term in the Hamiltonian describes hopping between different sites, the second deals with the intra-atomic Coulomb repulsion cost associated to having a given site (or, more formally, the corresponding p z -orbital) doubly occupied.
Unless stated otherwise, we consider systems at halffilling and assume that all t i,j are zero except when i and j are first or third neighbors.We denote first and third neighbor hoppings by t and t 3 , respectively.Secondneighbor hopping t 2 introduces charge inhomogeneities that are penalized by the Hartree interaction, so that best agreement with DFT is obtained by assuming t 2 = 0. Throughout this paper we set t = −2.7 eV and t 3 is taken as a free parameter.For triangulene 2D crystals, good agreement with DFT calculations is obtained if we take t 3 0.1t (ref. 22).

B. CAS
Due to the exponential increase in complexity as the size of a quantum system grows, exact diagonalization of many-body problems is only possible for rather small systems.To treat larger systems, approximate solutions have to be introduced, one of them being the configuration interaction method in the CAS approximation.Here, we follow the implementation of the CAS method for the Hubbard model presented in previous works by some of us 7,18,30,31 .First, the single-particle spectrum of a given triangulene structure is obtained.Then, a subset of N MO molecular orbitals (MOs)-containing the zero modes and closest states in energy-is selected, and a complete set of multi-electronic configurations with N e electrons occupying these N MO MOs is considered; the rest of the electrons are assumed to fully occupy the MOs below the active space.The Hubbard Hamiltonian is represented in this restricted basis set and diagonalized numerically.Hereinafter, we shall refer to this procedure as CAS(N MO , N e ).Since there is one electron per π-orbital for triangulenes at charge neutrality, we always consider N e = N MO .

C. Mean-field
In contrast with the CAS method, the mean-field approximation for the Hubbard model makes it possible to include all the single-particle states of molecules and crystals, but interactions are treated approximately.The mean-field theory can be formulated variationally, where the many-body wave function is written in terms of a set of independent electrons that occupy the energy levels of a mean-field Hamiltonian.Here, we impose that the total S z is a good quantum number, thus breaking the spin-rotation invariance present in the original Hubbard model; this is the so-called collinear mean-field approximation, extensively used in the modelling of magnetism in graphene nanostructures 2,25,26,28,[32][33][34][35][36][37][38][39] .In this case, the Hamiltonian takes the form: where the local densities n iσ are computed with the variational wave function.Therefore, the variational wave function and the mean-field Hamiltonian have to be determined in a self-consistent manner.In practice, this is done by iteration, starting from an initial guess for the local densities.In crystals, the local densities are also periodic so that the eigenvalues and eigenvectors of Eq. ( 2) satisfy Bloch's theorem and can be classified in terms of a wave vector k.
In general, we classify the collinear mean-field solutions in two groups: broken-symmetry solutions for which the expectation values of the local spin operators are finite, and non-magnetic solutions, that are isomorphic to the non-interacting case, except from a trivial rigid shift of the energies.Therefore, the mean-field method provides the value of S z that minimizes the energy, the expectation value of the local moments, and a set of energy levels.These three quantities can be compared with DFT.In the case of graphene nano-islands 2 and ribbons 25,28,40 , the predictions of mean-field theory were found to be in good agreement with those of DFT for some values of U .For triangulene 2D crystals we also find a good agreement.Therefore, comparison of DFT and mean-field Hubbard models allows us to obtain an educated guess for U in these systems.
In our mean-field calculations for 2D triangulene crystals, we have considered a 5 × 5 Monkhorst-Pack grid for the k-point sums and a tolerance of 10 −4 for convergence in the local densities.Different initial guesses for the local densities were tested, with the antiferromagnetic guess found to yield the lowest energy solution in all relevant cases.

D. RPA
In order to study spin excitations of 2D triangulene crystals, we use the standard RPA to calculate their transverse spin susceptibility matrix for wave vector Q and frequency Ω (refs. 26,27,41,42), which is the space and time Fourier transform of the spinflip Green function, where i, i are atomic site indices within a unit cell, is the reduced Planck constant, N is the number of unit cells, R, R denote unit cell positions, S + R,i (t) is the time-dependent version (in the Heisenberg picture) of the spin ladder operator S + R,i ≡ c † R,i,↑ c R,i,↓ at time t, S − ≡ (S + ) † , θ(t) is the unit step function, and [•, •] denotes the commutator.The spin-flip Green function depends only on the relative position of unit cells R − R due to the translation symmetry of the crystal.
Within the RPA, we first obtain the mean-field susceptibility χ 0 , which corresponds to taking the average •, • over a self-consistent mean-field state associated with the Hamiltonian given in Eq. ( 2).Then the "interacting" susceptibility can be calculated using the RPA equation, which can be cast in the following matrix form: where [χ] contains the matrix elements χ ii (Q, Ω), and analogously for [χ 0 ].The specific mean-field susceptibilities that are relevant to us are given by Lindhard-like expressions, where ψ k,λ,σ (i) is the wave function coefficient, at site i, of a Bloch eigenstate of band λ with wave vector k and spin σ of the mean-field Hamiltonian.The associated eigenvalues are E k,λ,σ , and f (E) is the Fermi-Dirac distribution function.The sum over k spans the Brillouin zone of the crystal.To calculate χ 0 ii (Q, Ω), we have used 2500 reciprocal space points within the Brillouin zone of the crystal (equivalent to considering N = 2500 unit cells), which guarantees convergence of the k-space sum.All the results are obtained at zero temperature.An empirical broadening of the single-particle states η = 5 meV has been adopted.
The RPA expression for χ also allows to determine the critical value of U , denoted by U c , above which the nonmagnetic solutions are no longer stable.The magnetic instability is signaled by the condition det(1 − U c [χ 0 ]) = 0, with [χ 0 ] calculated at Ω = 0 for the Hamiltonian in Eq. ( 2) with U = 0.The kind of spin arrangement towards which the true self-consistent mean-field solution tends, either ferromagnetic or antiferromagnetic, is indicated by the wave vector Q at which the condition is satisfied for the smallest U c , together with the eigenvector of [χ 0 ] corresponding to its largest eigenvalue 43 .

E. DFT
The DFT calculations have been performed with the local-density approximation, as implemented in Quantum Espresso 44 .We have used norm-conserving pseudopotentials, with a kinetic energy cutoff of 50 Ry and a k-point sampling of 12 × 12 × 1 in a Monkhorst-Pack mesh.To avoid spurious interaction between replicas we have set a vacuum distance of 20 Å.We have set the same experimental lattice parameter and atomic positions for the three cases of magnetic order (non-magnetic, ferromagnetic and antiferromagnetic).The optimized atomic positions have been calculated using the non-magnetic phase and the final structure is planar.

F. LDOS
The LDOS at energy E and position r = (x, y, z) was calculated using the following equation: The delta function was approximated by a Lorentzian of the form where Γ is the half width at half maximum of the Lorentzian function.In our calculations, we took Γ = 8 meV and used a 5 × 5 Monkhorst-Pack grid for the k-point sum.Moreover, we considered a carbon Slater distribution for the 2p z atomic wave function, with r 0 = 0.325 Å (refs. 30,45).For clarity, R denotes a lattice vector and R i is the specific position of site i in the corresponding unit cell.

III. HUBBARD MODEL CAS CALCULATIONS FOR CENTROSYMMETRIC DIMERS
In this section we present the results of CAS calculations for [n]triangulene dimers.The main goal here is to show that, for U/|t| 1, the low-energy spectrum can be mapped to a spin model, providing evidence that the dimers remain open-shell and the triangulenes preserve their magnetic moments.We focus on the cases of n = 4 and n = 3 dimers, as the n = 2 case has been already studied in detail in previous work 31 .
dimers have a null sublattice imbalance, so they may have no zero modes.However, we find 22 that there are 2(n−1) states close to zero energy, on account of the vanishing weight of the zero modes on the intermolecular binding sites.Only third-neighbor hopping leads to a small intermolecular hybridization of the triangulene zero modes 22 .
In Fig. 2 we show the single-particle spectra for [3]and [4]triangulene dimers, obtained by solving the Hamiltonian of Eq. ( 1) with U = 0 and t 3 = t/10.As expected, for the [3]triangulene dimer we find four states close to zero energy, originating from the weak intermolecular hybridization of the two zero modes hosted by each monomer individually, promoted by third neighbor hopping.For the [4]triangulene dimer, a similar result is found, only this time the three zero modes of the monomers hybridize to give six states close to zero energy.We also depict the choice of MOs that will enter in the CAS calculation for each of the molecules.These active spaces were chosen to include an additional pair of orbitals besides the zero modes, as this is crucial to account for the Coulomb-driven superexchange mechanism 31 .
We now discuss our CAS calculations for the [3]-and [4]triangulene dimers.The results for U = |t| and t 3 = t/10 are presented in panels (a) and (c) of Fig. 3.While the [3]-and [4]triangulene monomers are sublattice imbalanced, which according to Lieb's theorem implies a ground state with finite total spin (S = 1 and S = 3/2, respectively), for [n]triangulene dimers the sublattice imbalance vanishes and the ground state has S = 0.For the n = 3 dimer, this ground state is followed by a triplet (S = 1) and a quintet (S = 2).For the [4] triangulene dimer, an additional septet (S = 3) follows the S = 1 and S = 2 manifolds.
In panels (b) and (d) of Fig. 3, we show the CAS results for different values of U , thus allowing to study how the energies of the many-body states are affected by the strength of the on-site Coulomb repulsion.Inspecting this figure, it becomes clear that, for U |t|, the lowenergy excitation order S = 1, S = 2 (and S = 3 for the [4]triangulene dimer) is preserved and, crucially, remains well separated from higher-energy excitations.As U is reduced, however, the low-lying excitations become closer to the high-energy ones, and for a critical value of U a crossover is visible.
In the parameter region where the low-energy manifold is well separated from the higher-energy states, the low-energy spectrum of the triangulene dimers can be modeled by a simple spin Hamiltonian where each triangulene is represented by a spin whose value is that of the ground state of the corresponding monomer.To establish a quantitative comparison, we postulate a non-linear n J (meV) β2 β3 3 27.9 0.12 -4 11.3 0.09 0.007 Heisenberg dimer Hamiltonian, In Appendix A, we derive analytical expressions for the energy levels of this Hamiltonian.By matching these expressions with the results found with CAS for the lowenergy manifolds of the [n]triangulene dimers, we are able to compute J, β 2 , β 3 as a function of U and t 3 .As a reference, in Table I we give their values for U = |t| and t 3 = 0.1t.We see that, for both molecules, the exchange coupling J is in the order of tens of meV, with the n = 3 dimer presenting the larger antiferromagnetic exchange.In both cases, it is found that the biquadratic term (given by β 2 J) is approximately 10% of the bilinear one (J), emphasizing its importance to accurately capture the energy levels with a spin model.As for β 3 , which is only included in the model of the n = 4 dimer (as explained in Appendix A), it is found to be one order of magnitude smaller than β 2 .Thus, we see that the bicubic term only introduces minor corrections to the energy spectrum, which further justifies not accounting for it to describe the n = 3 dimers.
The fact that we can map the low-energy levels of the fermionic CAS calculation to a spin model, together with the fact that, for U |t|, these are well separated from higher-energy excitations provides a strong evidence that the dimers are in the open-shell regime, the triangulenes host local moments, and the singlet ground state arises from the intermolecular antiferromagnetic coupling.This shows that, although intermolecular hybridization is present, the magnetic nature of the triangulenes is preserved and the intermolecular interactions are antiferromagnetic.A comparison of the intermolecular hybridization and the Coulomb energies is provided in Appendix B. As we decrease U , the low-energy excitations and the high-energy ones become closer, and the validity of the model is no longer warranted.The spin model description certainly fails where the crossover between low-and high-energy excitations occurs 47 .

IV. 2D CRYSTALS: DFT AND MEAN-FIELD HUBBARD MODEL CALCULATIONS
In this section we undertake the study of magnetic properties in 2D triangulene crystals.For that matter, we compare DFT-based calculations, both spin-restricted and spin-polarized, with mean-field Hubbard model results.We consider ferromagnetic (FM) and antiferromagnetic (AF) broken-symmetry solutions, as well as nonmagnetic (NM) states.In all cases considered, we find that the lowest energy configuration corresponds to the AF solution.
A. DFT for the [4, 4]triangulene crystal We now discuss the electronic properties of the [4,4]triangulene crystal, as described with DFT-based calculations.We note that both the spin-restricted and the AF cases of the [2,2]-and [3,3]triangulene crystals were addressed in previous works 22,24 .In both systems, it was found that the NM solution is an excited state and describes a zero-gap semiconductor with two Dirac cones and a narrow bandwidth.The spin-polarized AF solution opens up a large gap and is the ground state.
Previous work 20,22,23 has shown that the spinunpolarized [4,4]triangulene crystal is a narrow-gap semiconductor with flat valence and conduction bands.Here, we go beyond the NM framework and study two magnetic phases, AF and FM.We find that the AF phase has smaller energy than both the FM (E F M − E AF = 0.171 eV) and the NM (E N M − E AF = 0.457 eV).It is thus apparent that DFT calculations confirm the open-shell nature of the [4]triangulenes when covalently bonded to form a 2D honeycomb crystal.If we model the energy difference between the AF and FM phases with a classical Heisenberg model on a honeycomb lattice, we get 6JS 2 = 0.171 eV.Using S = 3/2, we pull out J = 12.7 meV.For the [3,3]triangulene crystal, a similar analysis 22 AF = 0.159 eV and J [3,3] = 26.5 meV.
In Fig. 4, we show the energy bands for the three configurations (NM, FM, AF) of the [4,4]triangulene crystal, together with the distribution of the magnetic moments in the FM and AF solutions.The three solutions are gapped, but the size of the gap increases in the magnetic phases, specially in the AF case.The FM bands have a similar line shape than the NM bands, except for the top of the conduction band.The AF bands are much narrower than the NM bands.This relates to the quenching of intermolecular hybridization due to the opposite-sign spin splitting of the zero modes of adjacent molecules.
We note that, whereas the magnetic moments lie predominantly in the majority sublattice of each triangulene, there is a smaller magnetization with opposite sign in the minority sublattice, coming presumably from electrons in non-zero modes.Moreover, we find that the magnetization per triangulene shares a similar pattern for both FM and AF solutions, and the values obtained are compatible with the predictions for individual triangulenes.

B. Mean-field Hubbard model results
We now present our results for the [2,2]-, [2,3]-, [3,3]-, and [4,4]triangulene 2D crystals, obtained using the collinear mean-field approximation to the Hubbard model at half-filling.For the centrosymmetric [n, n]trianguelene crystals, we find that, for U above an n-dependent critical value U c (n) below which the ground state solutions are NM (see Subsection IV D), the lowest energy solutions are AF, in agreement with DFT calculations.As for the non-centrosymmetric [2,3] case, the ground state obtained is always ferrimagnetic.
In Fig. 5, we show the energy bands for both NM and ground state (magnetic) configurations, obtained with U = 0 and U = |t|, respectively.Two features are immediately apparent.First, the dispersion of the AF bands is narrower compared to the NM case.This is a consequence of suppressed intermolecular hybridization, on account of the opposite-sign spin splitting in the two triangulenes of the unit cell.Second, the separation between conduction and valence bands increases in the magnetic phases.Thus, the [2,2]-, [2,3]-, and [3,3]triangulene crystals, gapless for U = 0, become gapped when magnetic order appears.In the case of the [4,4]triangulene crystal, gapped for U = 0, the interactions increase the gap by more than a factor of 3. The gap of the magnetically ordered phases reflects the fact that every triangulene is full-shell in a spin-channel, so that the addition of a new electron is only possible in the minority spin channel, that became spin-split.
In Fig. 5, we also show the local magnetic moments corresponding to the ground state mean-field solutions.The magnetization pattern is such that carbon sites in different sublattices are magnetized with opposite sign.For U |t|, the magnetic moments per triangulene are close to the values expected from Lieb's theorem for individual triangulenes, and in qualitative agreement with those of DFT.We note that mean-field theory is not constrained by Lieb's theorem, that applies to exact solutions.
For the non-centrosymmetric [2,3]triangulene crystal, the magnetic order appears for arbitrarily small values of U .This is expected on account of the flat band at the Fermi energy.For small values of U , magnetic moments are only present in the larger unit, that hosts the flat-band states.As U is ramped up, the magnitude of the magnetic moments in both units increases towards values close to those of the isolated triangulenes, and a ferrimagnetic ground state is obtained.

C. Comparison between mean-field and DFT
In this section, we briefly compare the results of the mean-field Hubbard models with those of DFT, for the [2, 2], [3, 3] and [4, 4] crystals.The DFT results for the [2, 2] crystals are taken from ref. 24 .As for the [3,3] crystals, DFT calculations were reported in ref. 22 by two of us.Since the comparison of the NM phases was already established in previous work 22 , we focus on the magnetic phases.
Qualitatively, both levels of theory are in agreement.They both predict AF solutions as the ground state, with magnetic moments close to those predicted for isolated triangulenes.Moreover, both in mean-field and DFT the band-gap of the magnetic solutions is much larger than the NM cases, and the band dispersion is narrower.
Given the uncertainty over the best value of U , we make no attempt to find the value of U for which this agreement is better, and we take U = |t| as a reasonable guess.It is apparent that the mean-field bands obtained with U = |t| are in good agreement with the DFT calculations.The same is also verified for the magnetization patterns (see figures (4)d,e and lower panels (5).A quantitative comparison between the mean-field theory for U = |t| and DFT is provided in Table II.We find a fairly good agreement that justifies the use of Hubbard models for this type of system.Specifically, for the [4, 4] case, we obtain a good agreement in: (i) the energy difference between the different magnetic phases (with the NM configuration featuring the highest energy of the three); (ii) the band gaps of the NM and AF solutions, with both levels of theory predicting an increase of the band gap by a similar factor in the AF phase; (iii) the S z per triangulene (discussed above); (iv) the absolute value of the magnetization, defined by |M tot | = gµ B i |S z (i)|, where g = 2 is the electron g-factor and µ B stands for the Bohr magneton.

D. Critical value of U
We now discuss the minimal value of U that makes the NM solution unstable within the mean-field Hubbard approximation.This can be obtained in two ways.First, by comparing the NM and the magnetic solutions of a mean-field calculation as a function of U and finding the critical value U c above which the disordered phase becomes an excited state.Second, a faster approach, discussed in Section II D and adopted here, where we look for the value U = U c for which the non-interacting RPA susceptibility diverges.The results are shown in Fig. 6 for the [n, n] crystals with n = 2, 3, 4. We note that, for a System Quantity DFT Mean-field [2, 2]  EF M − EAF (eV) 0.11 a 0.109 [2, 2]  ENM − EAF (eV) 0.12 a 0.097 [3, 3]  EF M − EAF (eV) 0.159 b 0.137 [4, 4]  EF M − EAF (eV) 0.171 0.133 [4, 4]  ENM − EAF (eV) 0.457 0.508 NM [4, 4]  Gap (eV) 0.185 0.148 AF [4, 4]  Gap (eV) 0.716 0.625 AF [4, 4]  |Mtot| (µB) 8.89 9.01 AF [4, 4] Sz per triangulene 1. 53  1.49For the [2, 2]triangulene crystal, whose low-energy single-particle Hamiltonian maps exactly to that of a honeycomb model 22 , the effective first neighbor hopping is given by t = |t 3 |/3 and the effective Hubbard interaction is Ũ = U/6 (ref. 22).Therefore, by renormalizing the U c = 2.2|τ | equation, we can estimate a critical value U c = 4.4|t 3 | for the [2,2] crystal, in good agreement with the U c numerically obtained (see Fig. 6).
Given that the low-energy spectrum of the [3, 3]triangulene crystal also features graphene-like bands in the neighborhood of the Fermi energy, we can also compare the numerically obtained U c with that of the honeycomb crystal.For the [3, 3] crystal, the effective first neighbor hopping is given by t = 2|t 3 |/11 and the effective Hubbard is Ũ U/11 (ref. 22).Therefore, we also estimate a critical value U c 4.4|t 3 |.The fact that our numerical estimates for U c (t 3 ) are slightly different (see Fig. 6) reflects the fact that U c is also influenced by the flat bands away from the Fermi energy.
The critical values of U for t 3 = 0.1t are in the range of U c 0.45|t| < 1.2 eV.Estimates of atomic U for carbon are higher than this, in the range of 3.5 eV (ref. 6).Mean-field theories are known to underestimate U c .For instance, Quantum Monte Carlo methods 48 predict U c = 4.5|τ | for the Hubbard model on the honeycomb lattice.Even if U c is twice as large as the values predicted by mean-field, magnetic order should appear in the triangulene crystals.
Interestingly, the value of U c is very similar for the [3, 3] and [4, 4] crystals.This result further supports the picture that, once moderately large interactions are included, the fact that the NM bands of the [4, 4] crystal have a band-gap does not seem to have a dramatic effect on its electronic properties, and the [4, 4]triangulene crystal is (antiferro) magnetic.

V. COLLECTIVE SPIN EXCITATIONS IN 2D TRIANGULENE CRYSTALS
A. RPA for the Hubbard model The choice of a "ground state" with broken spin rotation symmetry implies the existence of gapless Goldstone modes, the spin waves.Here we obtain the spin wave spectra of 2D triangulene lattices by computing the transverse spin susceptibility χ(Q, Ω) of the Hubbard Hamiltonian (Eq.( 1)) in the RPA, as discussed in Section II D. The spin wave frequencies are associated with the poles of χ(Q, Ω).For a given wave vector, two poles occur at energies ± Ω(Q), due to the opposite directions of the spins in the two magnetic sublattices. 49rom those we can build a spin wave dispersion relation, shown in Figs.7a,b for the [3, 3] and [4, 4] crystals, and in Fig. 7c for the [2, 3] crystal.We note that, for the centrosymmetric cases ( [3, 3], [4, 4]), the two modes are degenerate, in contrast with the [2, 3] for which we find an acoustic and an optical branch of spin waves.In these figures, the symbols represent the locations of the poles of χ(Q, Ω) for a few wave vectors along two highsymmetry directions in the honeycomb Brillouin zone.It is apparent that the bandwidth of the magnon spectrum is larger for the [3, 3] crystal, in agreement with the larger values of intermolecular exchange obtained with the CAS calculations for the dimers.

B. Comparison with spin models
We now compare the RPA results with those of a Heisenberg spin model with first-neighbor exchange J, calculated in the linear spin-wave approximation 50,51 .The calculation (not shown) is standard 51 .The spin operators are expressed in terms of Holstein-Primakoff (HP) bosons 50 , taking the quantization axis parallel to the classical ground state (AF for the [3, 3], [4, 4], ferrimagnet for the [2, 3]), where the classical magnetization of each [n]triangulene is 2S = n − 1.The resulting bosonic Hamiltonian is truncated so that only terms bilinear in the HP bosons are kept.This bilinear Hamiltonian can be solved exactly, by means of a paraunitary canonical transformation.For a lattice with two spins per unit cell, such as the honeycomb, two spin-wave branches are obtained, given by: where , S A and S B denote the spin of the triangulenes in sublattice A and B, 3φ k = 1+e ia1•k +e ia2•k , a 1,2 are the lattice vectors of the honeycomb lattice and J is the intermolecular exchange.
It is apparent that in the AF case we have S A = S B and the two branches become degenerate.It is also apparent that, for k = (0, 0), the lower energy branch vanishes, complying with the Goldstone theorem.
In Figs.7a-c, we compare the magnon dispersion calculated from the fermionic RPA theory with the spin wave dispersion of Eq. ( 12) .Taking S from the meanfield calculation, we determine the value of intermolecular exchange J that provides the best fitting to the RPA calculation within the fermionic model.We find that the RPA curves lies exactly on top of the spin-wave curves, providing additional support to the notion that the lowenergy exciations of 2D triangulene crystals can be described with spin model Hamiltonians, very much like the one-dimensional triangulene spin chain.
We can determine the dependence of J on U and t 3 by repeating this procedure for different values of those parameters.In Fig. 8 we plot J, so obtained, as a function of U/|t| with t 3 = t/10 for the [3, 3] and [4, 4]  cases.The general behavior is qualitatively very similar to the results from CAS calculations for dimers, as seen in Figs.10a,c.In fact, even the actual values of J given by RPA and CAS are reasonably similar for 0.5 U/|t| 1.5.This qualitative good agreement backs-up the robustness of the main underlying picture of this work: despite the intermolecular hybridization between triangulenes in the 2D crystals considered here, they retain their magnetic moment.

VI. PREDICTIONS FOR STM SPECTROSCOPY
We now discuss experimental consequences of the magnetic order discussed in the previous sections.Given that, so far, triangulenes structures are studied with STM 15,19,20 , we focus properties that can be probed with this technique.STM dI/dV can reveal two different properties of the surface 21 : LDOS and inelastic excitations.In the case of nanographenes, LDOS features are revealed as prominent peaks at large voltages, in the range of hundreds of meV, corresponding to resonant tunneling accross specific energy levels of the molecules.

A. Probing LDOS
Here we discuss the LDOS at the energy of the valence and conduction bands, that can be measured by means of STM spectroscopy.The LDOS is sensitive to the interatomic coherence: by virtue of Eq. ( 8), the LDOS is proportional to the square of the MO wave function, that in turn is a linear combination of atomic orbitals.[4,4] crystal, c) [2,3] crystal, for which the two spin wave branches (corresponding to different polarities) are non-degenerate due to the ferrimagnetic nature of the ground state.In contrast, for the [3, 3] and [4, 4] (antiferromagnetic) crystals, the two polarities are degenerate, thus a single dispersion is shown for each.In a), b) and c) panels, the dots have been extracted from magnon spectral densities (the imaginary part of the transverse spin susceptibility) and the solid curves are fits to nearest-neighbor Heisenberg models.d) Integrated magnon density of states for the nearest-neighbor antiferromagnetic Heisenberg model on the honeycomb lattice, obtained from Eq. ( 12) with SA = SB = S. J (meV) [3,3], RPA [4,4], RPA [4,4], CAS [3,3], CAS FIG. 8. Effective exchange as a function of U , obtained by fitting the spin wave dispersion relation to a nearestneighbor Heisenberg model in the linear spin wave approximation (dots), and from CAS calculations on single dimers (solid curves).
Therefore, LDOS is sensitive to the relative phases of the weights of the MO at different atoms.Specifically, for the non-interacting bands of a bipartite lattice, electron-hole symmetric states, such as valence and conduction bands have opposite relative phases between adjacent atoms.More formally, let us denote the wave function of a conduction band MO as where ψ k,λ,σ (A/B) encodes the MO weight on all the atoms in the unit cell that belong to the A and B sublattices.Then, for a bipartite lattice, the wave function of the electron-hole conjugate state λ in the valence band is given by 36 Thus, electron-hole conjugate MO wave functions have the same probability amplitudes but opposite phases at one sublattice.As a result, LDOS will have a enhancement/depletion at the regions connecting atoms with different sublattices.Specifically, at the bonding region between any pair of atoms, we can truncate Eq. ( 10) keeping only contribution of the two closest atoms, a and b, that, by definition, belong to different sublattices.This leads to We can now compute the difference of the LDOS computed in the bonding region between two atoms, for which eq. ( 15) holds, evaluated at energies +E λ and −E λ .The contributions to the difference of LDOS at these two energies from states with the same wave vector k in a pair of electron-hole conjugate bands λ, λ will be given by From Eq. ( 16) it is apparent that the LDOS contrast between electron-hole bands is controlled by the weights (and crucially the respective phases) of the wave functions of the MOs of adjacent atoms (which belong to different sublattices).
In the case of conduction and valence bands of triangulenes, the weight of the wave functions inside each triangulene is all on a single sublattice.Therefore, the LDOS contrast is only seen at the inter-triangulene binding sites.These are shown for the [4, 4] triangulene crystal in Fig. 9.For the non-interacting case we find a depletion of the LDOS at the intermolecular binding sites at the conduction band energy and a corresponding enhancement of the valence band (see Fig. 9a,b).
We now discuss how interactions, described at the mean-field level, change this picture.The brokensymmetry Néel states result in the presence of a staggered exchange potential.As a result, for a given spin direction, the on-site energy of two adjacent triangulenes is no longer the same.Consequently, the wave functions of valence and conduction bands no longer have the same weight on both sublattices 36 , i.e.Eqs. ( 13) and ( 14) relating the wave functions of valence and conduction bands no longer hold.In the interacting cases the MOs become sublattice biased, and in the very strong coupling limit, completely sublattice polarized.This ultimately reduces the amplitude of the bonding-anti-bonding interference effect, as shown in Fig. 9c,d .This is in agreement with the experimental observations of Delgado et al. 20 .
We note here that the reduction of the LDOS contrast between valence and conduction bands in the interacting cases relates to the reduced bandwidth of the interacting bands.The spin-dependent staggered potential creates an energy barrier for intermolecular hybridization.
In the work of Delgado et al., 20 , the observed reduced contrast of LDOS at the binding site is attributed to an excitonic insulator state that arises on account of the small gap obtained from the spin-unpolarized DFT calculations.As both our DFT and mean-field results show, the spin polarized solution has lower energy, a larger gap (that makes the excitonic insulator state less likely) and, more important, already accounts for the reduced LDOS contrast in terms of the sublattice symmetry breaking of the AF solution.

B. Probing magnons
In contrast to the large-bias LDOS measurements discussed above, STM spectroscopy can reveal inelastic excitations as bias-symmetric steps, at bias voltages below 100 meV, whose height is dramatically smaller than the resonant peaks.The underlying mechanism for these inelastic steps is inelastic cotunneling of electrons 21,52 .In spin systems, inelastic electron tunneling spectroscopy (IETS) can probe spin transitions between the ground state and excited states that satisfy the rule for the change of total spin ∆S = 0, ±1.Therefore, IETS is optimal to probe magnons 53,54 .We expect that d 2 I dV 2 will have a line shape that reflects the density of states of magnon excitations.In Fig. 7d, we show the magnon LDOS associated to the dispersion energy from Eq. ( 12) for S A = S B = S, relevant for [n, n]triangulene crystals.It is apparent that the magnon LDOS features an outstanding Van Hove singularity, at the energy √ 8JS, corresponding to the M points in the Brillouin zone.Therefore, we anticipate the presence of a prominent feature at eV = ± √ 8JS energy in the d 2 I dV 2 spectra.For the [4, 4]triangulene crystal, taking J 9 meV (see Table III), and S = 3/2, the steps are expected at ±38 meV.

VII. SUMMARY AND CONCLUSIONS
The main goal of this paper is to describe the consequences of electron-electron interactions in triangulene 2D crystals.Specifically, we address the question of whether triangulenes retain their magnetic moments when forming 2D crystals that entail intermolecular hybridization.This is particularly relevant in the case of [4, 4] triangulene crystals, for which the single-particle model 22 predicts an insulating state that, naively, may quench the emergence of magnetism.
We employ both spin-unrestricted DFT calculations and Hubbard model with three different approximations: 1. Multiconfigurational calculations of triangulene dimers.
2. Mean-field approximation of the 2D crystals, whose results are in qualitative agreement with DFT calculations.
3. RPA calculations of the spin excitations.
Importantly, these different methods permit us to perform cross-validations.For instance, both the spinpolarized DFT and mean-field Hubbard model yield very similar results for all the key quantities (see table II).The Hubbard-model RPA calculations, building on top of mean-field solutions, predict an excitation spectra that can be fitted very well to a Heisenberg model, with just a single fitting parameter, the effective exchange (see figure (7)a,b,c).In turn, this exchange is in qualitative agreement with the one obtained from CAS calculations, for a range of values of U (see figure (8)).
Our main conclusions are: 1. Triangulenes retain their magnetic moment when forming in two-dimensional crystals, according to both DFT and Hubbard model calculations.
2. Triangulene crystals are insulating, on account of the electron-electron interactions.This is supported both by our mean-field Hubbard model and our spin-polarized DFT calculations.In the case of the [4, 4] crystal, the size of the gap, calculated with DFT, comes out 3.9 times larger than the spin-unpolarized gap, which calls for a revision of the predictions 20,23 of an excitonic insulator state based on the smaller gap of the non-magnetic ground state.
3. Two-dimensional triangulene crystals are magnetically ordered, either antiferrromagnetically, in the centrosymmetric case, or ferrimagnetically, for noncentrosymmetric crystals.This statement is based on both DFT and mean-field Hubbard calculations.
4. The value U = |t| gives a very good agreement between mean-field Hubbard model and DFT for several quantities, such as the intermolecular exchange, the magnetic moment and the band-gap.
We have not tried to fine-tune U/|t| to improve that agreement, but we can be sure that U |t| is a good ball-park reference for this important ratio, in agreement with previous work 15 .
5. The low energy spin excitations obtained from the RPA fermionic calculations are very well described with Heisenberg Hamiltonians (see Fig. 7).The exchange interactions so obtained are in qualitative agreement with those obtained from CAS, meanfield and DFT (see table III).
6. Intermolecular exchange features non-linear interactions, beyond Heisenberg.This is found by comparing CAS calculations for the Hubbard model with spin models.For U |t|, the values of the non-linear terms are in qualitative agreement with previous work for S = 1 15 .For S = 3/2 triangulenes, the value of the non-linear interactions are very small, so that it is very unlikely that the system realizes the AKLT model for the honeycomb lattice 55 , that would be for relevant measurementbased quantum computing 56 .7. Magnetically ordered states reduce the intermolecular hybridization in triangulene 2D crystals.This has two consequences: • First, the bandwith of magnetically ordered triangulene crystals are narrower than those of the non-interacting case.• Second, a reduction in the difference of the LDOS measured at the valence and conduction band energies at the intermolecular binding sites.This reduction has been observed experimentally in recent work 20 in [4, 4] triangulene crystals.
We now briefly discuss the robustness of the predicted antiferromagnetic states.By construction, both the mean-field and DFT calculations predict brokensymmetry solutions or non-magnetic solutions.Both quantum and thermal fluctuations can destroy long range order in two dimensions.
At T = 0, broken symmetry states are robust in this class of system.Using Quantum Monte Carlo, it was shown that the Hubbard model in the honeycomb lattice, at half-filling, features antiferromagnetic long-range order 48 .Since the effective model for the [2, 2] crystal is a Hubbard Hamiltonian that maps into a S = 1/2 Heisenberg model, given that quantum fluctuations scale with 1/S, and given that larger triangulenes have larger S, we expect that at T = 0 the ground state of the centrosymmetric triangulene crystals also feature Néel long-range order.
In contrast, thermal fluctuations are expected to destroy long-range order, on account of Mermin-Wagner theorem 57 .However, the spin correlation length may be larger than system size for small crystals reported experimentally 20 .Therefore, the broken-symmetry solutions remain a good approximation for these systems, as in the case of one-dimensional edge magnetism in graphene ribbons 58 .
Our results, together with previous experimental work 18,20 , should pave the way for the design of other nanographene molecular crystals 22,59,60  Magnons in magnetic 2D materials for a novel electronics (2D MAGNONICS).This study forms part of the Advanced Materials programme and was supported by MCIN with funding from European Union NextGenera-tionEU (PRTR-C17.I1) and by Generalitat Valenciana, project SPINO2D, reference MFA/2022/009.

Appendix A: Derivation of spin Hamiltonian parameters
In this Appendix, we derive the analytical expressions for the energy levels of the spin model dimer Hamiltonian of Eq. (11).Then, by matching these expressions with the numerical results obtained with CAS for [n]triangulene dimers, we study how the parameters of the spin model, i.e.J, β 2 and β 3 , depend on U and t 3 .
In terms of the total spin S and the spin of the triangulenes s, the eigenvalues of the spin model of Eq. ( 11) are given by where S can take the values S = 0, 1, ..., 2s.Thus, for the n = 3 (n = 4) case, S can take values up to S = 2 (S = 3).For the [4]triangulene dimer, the spectrum of the spin model has four multiplets, with S = 0, 1, 2, 3.
For the [3]triangulene dimer, we assume β 3 = 0 since the S = 1 dimer model can take the values S = 0, 1, 2, and we can only fit two energy parameters out of three multiplets.As we found in previous work 18 , the model with β 3 = 0 can account for experimental observations of a large number of structures.The excitation energies for the n = 3 case are related to J and β 2 as follows: These equations can be easily inverted to obtain J and β 2 for a given fermionic calculation.
In the case of the [4]triangulene dimer, with s = 3/2, we obtain the following equations: As before, the system of equations can be inverted in order to obtain expressions for J, β 2 and β 3 in terms of the excitations energies; this then allows us to match the spin model with the fermionic calculation and obtain the dependence of the parameters with U and t 3 .
In Fig. 10 we present the values of J, β 2 and β 3 obtained for the two considered molecules, as a function of U , for possible values of t 3 .In each panel, the data is only presented up to the critical value of U for which the crossover between low-and high-energy excitations occurs; for smaller U the extraction of the spin model parameters is not valid.From these figures, one clearly sees that, for both dimers, the intermolecular antiferromagnetic exchange J is in the order of a few tens of meV, with the n = 3 dimer presenting a stronger intermolecular exchange than the n = 4.In both molecules we find that for U |t|, the parameter β 2 , which quantifies the weight of the quadratic term relative to the linear one in the model Hamiltonian, takes values up to approximately 1/5, emphasizing its importance to accurately describe these molecules with a spin model.The value of β 3 , describing the strength of the cubic term relative to the leading one, is found to be much smaller than β 2 , indicating that it only introduces a small correction in the spectrum of the [4]triangulene dimer.Moreover, the fact that β 3 β 2 further justifies our choice of setting β 3 = 0 for the n = 3 dimer.
At last, we note that, for the n = 3 dimer, the value of β 2 approaches 1/3 asymptotically as U decreases.In that limit, the spin model-which corresponds to the wellknown AKLT model 55  model.We note that the singlet-triplet gap cannot vanish for the Hubbard model, as Lieb's theorem 8 states that the ground state is unique.The low-energy MOs are bonding and antibonding linear combinations of zero modes.Therefore, the intermolecular hybridization is proportional to the splitting between electron-hole symmetric single-particle energies, ±E m , given by δ m = 2|E m |.The addition energy associated to the double occupancy of a single-triangulene zero mode is given by the product of the atomic Hubbard repulsion parameter U with the inverse participation ratio (IPR) 7 , where the zero mode wave functions can be obtained from the MOs |m ± , associated to the states with energies ±E m , through the equation: We note that the zero modes so obtained are not necessarily identical to those obtained from the solution of the individual triangulene problem, on account of the degeneracy of the zero mode manifold, that allows one to define different zero mode bases.The values of U depend on that choice.It is found that for centrosymmetric triangulene dimers, U In Table IV we show the values of r m for the n = 3 and n = 4 triangulene dimers, obtained with U = |t| and t 3 = 0.1t.All ratios are smaller than 1, in some cases much smaller.For comparison, the ratios for the closest energy MOs not formed with zero modes are ∼ 10 for n = 3 and ∼ 11 for n = 4.This analysis backs up the open-shell nature of triangulene dimers.

FIG. 3 .
FIG. 3. Results obtained with CAS for [3]triangulene (a,b) and [4]triangulene (c,d) dimers.Panels (a) and (c) show the energy of the many-body states, obtained with U = |t| and t3 = t/10.Panels (b) and (d) show the energy difference between the ground state and the first few excited states, as a function of U , for t3 = t/10.
) where S A , S B are the vectors of the spin operators for the individual [n]triangulenes, taken to be S A = S B ≡ s = 1 and S A = S B ≡ s = 3/2 for n = 3 and n = 4, respectively.

FIG. 4 .
FIG. 4.Electronic structure, obtained with DFT, of[4,4]triangulene 2D crystals for the following cases: (a) nonmagnetic, (b) ferromagnetic and (c) antiferromagnetic.The corresponding magnetizations of the FM and AF cases are shown in panels (d) and (e), respectively; the value of Sz per triangulene is also indicated.Red/blue colors correspond to spin-up/down.In (c), spin-up and spin-down bands are degenerate.

FIG. 5 .
FIG.5.Electronic properties of (a)[2,2]-, (b)[2,3]-, (c)[3,3]-and (d)[4,4]triangulene 2D crystals.Top panels show the NM energy bands, obtained using a tight-binding model with t3 = 0.1t; horizontal black lines denote the Fermi energy.Middle panels show the energy bands of the ground state solution of a collinear mean-field Hubbard model with U = |t|; red/blue colors denote spin-up/down.The corresponding magnetizations are shown in the bottom panels, where the size of the circles represents the magnitude of the local moments.The value of Sz per triangulene is also indicated.

FIG. 7 .
FIG.7.Spin wave dispersion relations within the RPA for [n, m]triangulene crystals.a)[3,3] crystal, b)[4,4] crystal, c)[2,3] crystal, for which the two spin wave branches (corresponding to different polarities) are non-degenerate due to the ferrimagnetic nature of the ground state.In contrast, for the[3, 3] and[4, 4] (antiferromagnetic) crystals, the two polarities are degenerate, thus a single dispersion is shown for each.In a), b) and c) panels, the dots have been extracted from magnon spectral densities (the imaginary part of the transverse spin susceptibility) and the solid curves are fits to nearest-neighbor Heisenberg models.d) Integrated magnon density of states for the nearest-neighbor antiferromagnetic Heisenberg model on the honeycomb lattice, obtained from Eq. (12) with SA = SB = S.

FIG. 10 .
FIG.10.Parameters of the spin model, obtained by matching the spin dimer energy levels to the CAS calculations for triangulene dimers, as a function of U , for t3 = 0.06t, t3 = 0.1t and t3 = 0.14t.Panels (a) and (b) show the results for the[3]triangulene dimer, and panels (c)-(e) show the results for the[4]triangulene dimer.
Appendix B: Comparison of energy scales controlling open-shell nature of [n]triangulene dimers A preliminary estimate of the open-shell nature of [n]triangulene dimers can be obtained by analyzing the ratio between intermolecular hybridization energy of the zero modes and the effective addition energy.
m ≡ U m .Thus, considering only these two energies, every electron-hole symmetric pair maps into an effective Hubbard model dimer at halffilling, with effective hopping τ m = δ m /2 and Hubbard repulsion U m .Depending on the ratio between these two energies61 , model can describe a closed-shell configuration, for large r m , where two electrons occupy the bonding state with energy −|E m |, or an open-shell system where double occupancy of the zero modes is inhibited.We also note that the representation of the many-body Hamiltonian on the zero mode basis contains other interacting terms that couple the effective Hubbard dimers.

TABLE I .
Exchange coupling parameters (J, β2, β3) obtained by equating the eigenvalues of the non-linear Heisenberg dimer Hamiltonian of Eq. (11) to the CAS results obtained for [n]triangulene dimers with U = |t| and t3 = 0.1t.

TABLE II
48: Ref.24; b : Ref.22.honeycombHubbardmodel at half-filling, the mean-field critical value for the NM to AF transition is U c = 2.2|τ | (ref.48),where τ is the first neighbor hopping of the honeycomb lattice.

TABLE IV .
-has a vanishing singlet-triplet gap, but it is not a faithful description of the fermion Energy scales for [n]triangulene dimers, obtained with t3 = 0.1t.MO index m refers to the rank of a given molecular orbital when these are ordered in increasing energy order; the indices in parenthesis refer to the electron-hole symmetric partners.δ is the energy splitting between electronhole symmetric orbitals.IPR is the inverse participation ratio defined in Eq. (B1) of the sublattice mode of Eq. (B2).r is defined in Eq. (B3) and its value, obtained assuming U = |t|, indicates the closed-or open-shell nature of a given pair of MOs.