Superexchange Mechanism in Coupled Triangulenes Forming Spin-1 Chains

We show that the origin of the antiferromagnetic coupling in spin-1 triangulene chains, which were recently synthesized and measured by Mishra et al. (Nature2021, 598, 287−29234645998 ), originates from a superexchange mechanism. This process, mediated by intertriangulene states, opens the possibility to control parameters in the effective bilinear-biquadratic spin model. We start from the derivation of an effective tight-binding model for triangulene chains using a combination of tight-binding and Hartree–Fock methods fitted to hybrid density functional theory results. Next, correlation effects are investigated within the configuration interaction method. Our low-energy many-body spectrum for NTr = 2 and NTr = 4 triangulene chains agree well with the bilinear-biquadratic spin-1 chain antiferromagnetic model when indirect coupling processes and superexchange coupling between triangulene spins are taken into account.

2][3][4] Furthermore, in finite systems, characteristic spin-S = 1 2 edge excitations emerge 5 .Affleck, Kennedy, Lieb, and Tasaki (1987) proposed the AKLT model to elucidate the origin of the Haldane gap, introducing an isotropic one-dimensional spin-1 model with a unique ground state that also confirms the presence of an energy gap to excited states 3 .This model, characterized by an additional biquadratic exchange term, allows for a representation of the ground state in a valence-bond basis.The AFM Heisenberg spin-1 chain, along with models like AKLT, belongs to a broader class of bilinear-biquadratic (BLBQ) spin-1 chain models, which display diverse magnetic orderings dependent on the energetic ratio between linear and biquadratic exchange terms (see Refs. 6-8 and references therein).
Spin models are effective models of the low-energy subspaces of a Fermionic system [9][10][11] .
They can be derived from other simplified models such as the Hubbard model, which is an approximation to the full interacting Hamiltonian containing all two-body Coulomb interactions.AFM coupling in spin models can be explained by the kinetic exchange mechanisms, a coupling between singly and doubly occupied sites by effective hopping.This process favours AFM ordering of electron spins, since electrons with opposite spins are able to freely hop from one site to the next without being repelled due to Pauli exclusion 12 .Another exchange mechanism, known as superexchange has been proposed to explain AFM order in transition metal oxides, [12][13][14][15] .In that case, two cation orbitals are separated by an anion, as such, direct hopping between cations is quenched, while hopping between localized d-orbitals (located on the cation) occurs through an intermediate p-state (located on the anion) 12,14,[16][17][18][19] .
Recently, in 2021 Mishra et al. 20 observed signatures of fractional edge states, and gapped spin excitations in spin-1 triangulene chains.An AFM coupling between neighboring spins was described as a direct process in an effective Hubbard model with an enhanced nextnext-nearest-neighbor (NNNN) hopping commonly labelled t 3 in literature 10,20,21 .While that model captures major physical properties of the system, we show that the microscopic mechanism of the AFM coupling is an indirect superexchange process that includes states localized on the connection of triangulene molecules.This opens up opportunities to control paramaters of the spin model in a variety of ways.
In our work, we start from a combination of tight-binding (TB) and Hartree Fock (HF) methods and fit the results to density functional theory (DFT), in order to derive an effective single-particle model of triangulene chains that includes the presence of valence band electrons.The obtained spatial variation of TB hopping integrals takes into account geometric effects.We include correlations within the complete active space (CAS) using the configuration interaction (CI) method.We analyze the low-energy many-body spectrum as a function of CAS size and compare it with experimental results from Ref. 20.We find that the direct kinetic exchange contribution is too small to explain the spin-1 AFM model spectrum and the exchange mechanism that dominates the value of the coupling constant J is due to superexchange.Finally, we fit the many-body spectrum to the BLBQ model for N Tr = 2 and N Tr = 4 triangulene chains, confirming that at low-energies, the system behaves as a spin-1 chain, and allows us to extract effective spin parameters.
Triangular graphene quantum dots (TGQDs) with zigzag edges exhibit promising potential as building blocks for spin chains.This is because the single-electron spectrum for these quantum dots collapses to a degenerate shell at the Fermi level, a phenomenon resulting from the sublattice imbalance.Consequently, it is predicted that these quantum dots possess a magnetic spin-polarized ground state (GS) at half-filling [22][23][24][25][26][27][28][29][30][31] , a property consistent with Lieb's theorem 32 .The degeneracy of this shell is proportional to the length of a single edge of the triangular structure 28 .Triangulene, which is the second smallest TGQD with zigzag edges, containing 22 carbon atoms (also referred to as Clar's hydrocarbon), is of particular interest in this work and is shown in Fig. 1(a).Triangulene has two degenerate states at the Fermi level and consequently possesses a spin-1 ground state.It was predicted to exist and attempted to be synthesized in 1953 by Clar and Stewart 33 .Triangulene faced synthesis challenges due to the molecule's high reactivity arising from the two unpaired electrons.Overcoming this challenge more than 60 years later, Pavliček et al. 34 successfully synthesized triangulene in 2017.Following this milestone, noteworthy experimental work focused on the synthesis of triangulene as well as larger-sized TGQDs emerged [35][36][37] .
The first step in constructing a synthetic spin-1 chain using triangulene involved combining two triangulene molecules to form a triangulene dimer (TD) as illustrated in Fig. 1

(b).
This was accomplished experimentally 20,38 and investigated theoretically 9,10,39,40 in recent years.To understand the single-particle properties of the TD, we employ a p z orbital TB model, akin to previous studies 10, [41][42][43] .One assumption in these TB models, is that that the hopping parameters between two sites depends solely on the distance between them.However, when graphene is cut into a finite flake, electron hopping can be dramatically modified especially along the edge.We correct the hoppings by solving the HF equation given by 29,30,42  where c † i,σ /c i,σ creates/annihilates a p z electron on site i with spin σ. t ilσ are the twodimensional (2D) infinite system hopping elements.We take the following standard accepted values for t ilσ : t = −2.8eV as the nearest-neighbor (NN) hopping, t ′ = −0.1 eV as the nextnearest-neighbor (NNN) hopping and t 3 = −0.07eV is the NNNN hopping 44 .ρ 0 jkσ ′ are the density matrix elements for a 2D infinite graphene sheet obtained previously 29,42 .ρ jkσ ′ are the density matrix elements for the finite system computed with respect to the HF ground state.The Coulomb matrix elements ⟨ij|V |kl⟩ are screened by a dielectric constant κ = 3 which gives good agreement with DFT calculations, as seen in Fig. 1, see details of calculations in Supplemental Materials (SM) 45 .
For the TD shown in Fig. 1(b), there are N s = 44 atomic sites, and in the singleorbital model here, this corresponds to N e = 44 electrons for the charge neutral system.
It is convenient for the CI calculations that will follow, to remove four electrons from the TD, similar to what was done in previous work 27,29,30 and add them back when performing CI calculations.Therefore, we self-consistently solve the HF equation, as represented in Eq. 1, for N ↑ = 20 spin-up and N ↓ = 20 spin-down p z electrons, which is justified by noticing the presence of the energy gap between the lowest N st = 20 energy states and four degenerate shell states, seen in the tight-binding spectrum 10 .We obtain HF quasiparticle orbitals and compare them with the results from DFT calculations for a +4 cation TD.Fig. 1(c) show agreement between DFT and HF energy levels, and Fig. 2 shows agreement between their corresponding wavefunctions.We use one fitting parameter between the two calculations, which is a static screening constant κ = 3, which reduces the strength of Coulomb interaction by 1/κ.Examination of Fig. 1(c) reveals the presence of four distinct states in close proximity to the Fermi level, and a substantial energy gap separates these states from others.The wavefunctions of states No. 5 (the first state below the degenerate shell, Fig. 2(b) and (g)) and No. 10 (the first state above the degenerate shell, Fig. 2(f) and (l)), are localized at the connection between the two triangulenes.The wavefunctions of four degenerate shell states significantly differs from the wavefunctions of two degenerate states from an isolated triangulene as illustrated in Fig. 2(b-e) for the HF calculations and (h-k) for DFT calculations (compared with the Fig. 1 inset in SM 45 ).In particular, after HF self-consistent calculations, the degenerate shell states No. 6 (Fig. 2

(b) and (h)) and
No. 7 (Fig. 2(c) and (i)) strongly hybridize with the inter-triangulene states.This implies that in order to analyze coupling between the spins of isolated trianglulenes, we must take into account at least the six states shown in Fig. 2.
As seen, in the last line of Eq. 1, the HF equation reduces to a single effective hopping parameter tilσ which accounts for finite-sized modifications of the TB hoppings.It is given It is apparent in Eq. 2 that when ρ jkσ ′ = ρ 0 jkσ ′ , the effective hopping parameter reduces to that of the extended system hopping parameters t ilσ .Furthermore, deviations of the hopping elements will largely depend on the density matrix elements, as such, to obtain effective hopping paramaters that more closely resemble those that should be used in a TB model, we recompute the density matrix ρ jkσ ′ at charge neutrality.Thus, all states below the dotted line in Fig. 1(c) are occupied.Fig. 3 shows how the extended system hopping parameters are modified in HF for the TD.The line connecting two atoms represents the hopping between the two sites while the color represents the magnitude of the hopping.We see in Fig. 3(a) that in the center of the structure the effective NN hopping is almost equal to the 2D infinite The units for each color plot are measured with respect to their 2D infinite system hopping values.
system value t.Away from the center is where the largest modifications occur.There is an enhancement of the NN hopping element of about 10% near the bottoms of each triangulene molecule, and a reduction of about 25% on the connection between the two triangulenes.
Similarly, NNN t ′ hopping element in the center of each individual triangulene resembles that of the 2D infinite system, but a 300% enhancement of the NNN hopping element along the bottoms of each triangulene molecule is observed, as seen in Fig. 3(b).Finally in Fig. 3(c), even though, we observe an overall enhancement of the NNNN t 3 hopping, we observe a reduction of it on the connection between the two triangulene molecules.This contradicts a strong enhancement of t 3 in the effective Hubbard model considered in previous works 10,20,46 .
We notice that the renormalization of the hopping elements is dominantly attributed to lack of the charge uniformity brought upon by the edge effects of the TD.
In order to include the effects of correlations and obtain excited states of the TD we combine CI with the HF method.By rotating the many-body Hamiltonian to the basis of HF states one obtains 30,47 where ϵ HF pσ are the energies of the HF quasiparticle orbitals obtained by solving Eq. 1. b † p,σ /b p,σ creates/annihilates an electron on the HF orbital p with spin σ. τ pqσ are the double counting corrections which arises to avoid double counting of the interactions that were already contained at the mean-field level 30,47 .This effectively lowers the contribution of the last term in Eq. 3.
To solve the problem exactly for N Tr = 2, one would need to include N e = 44 electrons on N st = 44 states, which is currently impossible to solve, due to an unmanageably large many-body Hilbert space.Thus, we solve the CI Hamiltonian given by Eq. 3 in the basis of configurations by considering CAS(4,4), CAS (6,6), CAS (10,10), and CAS (14,14) states, as shown in Fig. 4(a), where we use the notation CAS(N st ,N e ) with N st /2 taken above and N st /2 taken below the dotted line shown in Fig. 1(c).
In the CAS(4,4) results, the ground state is ferromagnetic (FM) with S = 2, and corresponds to a Hund's rule like filling of electrons on the degenerate shell consisting of four states.This occurs when the Coulomb exchange amongst electrons on the shell of four states dominates compared to the effective AFM exchange coupling J between the spin-1 quasiparticles.When we include, into our active Hilbert space, the inter-triangulene HF states, N st = 6, shown in Fig. 2 meV, respectively) and could be improved by using a more advanced model for the dielectric screening.
We fit our CAS (14,14) spectrum to the BLBQ model given by where J is the bilinear exchange coupling constant between neighboring spin sites, and β is it's biquadratic partner.By fitting to the BLBQ model, we extract the effective exchange terms J = 8.21 meV and β = 0.016.Finally, in Fig. 4(b) we compare many-body spectra for N Tr = 4 (we take all eight degenerate shell states and six inter-triangulene states in the calculations) and the corresponding energy spectrum of the spin Hamiltonian given by Eq. 4.
We find good agreement between energy spectra in both models for each total spin eigenstate proving the validity of spin-1 model description (with almost the same parameters as for TD case, the effective exchange terms are J = 8.29 meV and β = 0.016).Small differences in the excited state energies can be attributed to a finite Hilbert space.We note that calculations that include more HF states is beyond our numerical capabilities.Fig. 4(c) shows the manybody spin density for the lowest energy triplet state, which resemble characteristic edge states in the Haldane phase.It is worth mentioning, that the excess of this spin density is similar to the one from the isolated triangulene (see SM 45 ).
The AFM superexchange mechanism is usually derived using higher order perturbation theory.We notice that after self-consistent HF calculations, we already are in a basis of strongly hybridized states of the degenerate shell and the inter-triangulene states.This hybridization comes from interactions between electrons populating all filled states and virtual electrons on a degenerate shell, and thus can not be uncoupled by any rotation within the effective low-dimensional Hilbert space; one needs to include all the single-particle energy states, not only the four degenerate shell states (for should also be responsible for AFM coupling in other triangulene systems, in particular, in recently proposed two-dimensional triangulene crystals 40,48 .Following this path, such newly created AFM spin lattice systems should reveal similarities to antiferromagnetism in transition metal oxides 16,17,19 , where superexchange mechanism plays a dominant role.
Superexchange interaction in anion-bridged dinuclear transition-metal complexes can be described by the Goodenough-Kanamori rules [49][50][51] .In analogy to these systems, it might be possible to establish similar rules for understanding the spin coupling in triangulene spin chains.However, the orbitals originating from the degenerate shell of triangulene are significanlty different than the d-orbitals of the metal ions.Thus, determining if similar rules apply here requires more detailed studies.
All density functional theory (DFT) calculations were carried out with the Turbomole 7.5 program package 1 employing the B3LYP global hybrid exchangecorrelation functional with 20% of the exact Hartree-Fock exchange 2-7 , Ahlrich's triple-zeta split-valence basis set with polarization functions on all atoms, def2-TZVP 8 , and the empirical dispersion correction of Grimme in the third generation 9 .Additionally, we employ the multipole-accelerated resolution of identity approximation for Coulomb integrals (MARIJ) 8,10-12 to speed-up the calculation of the Coulomb integrals in the selfconsistent field (SCF) algorithm.
First, we perform structure optimizations on the spin ground states approximated by Broken-Symmetry (BS) determinants 13 , where the S = 1 spins of adjacent triangulenes are aligned anti-parallel.In the structure optimizations convergence criteria of 10 −7 Hartree for the energy and 10 −4 Hartree/Bohr for the gradient were employed.On top of the optimized structures, we perform single-point calculations for non-spin polarized cations that serve as the reference in the full configuration interaction (CI) calculations.The cations are obtained by removing all unpaired electrons, two for each triangulene subunit, from the systems.In these single-point calculations,the same convergence criterion for the energy of 10 −7 Hartree is employed.

II. CALCULATION OF COULOMB MATRIX ELEMENTS IN THE SITE BASIS
The Coulomb matrix elements ⟨ij|V |kl⟩ are given explicitly as 14 where, ⃗ r 1 , ⃗ r 2 are the coordinates of electron 1 and electron 2, κ is the dielectric constant, and R y is the Rydberg constant.ϕ i (⃗ r) are p z Slater orbitals centered on atom i, they are given as where ⃗ r i is the position of carbon atom i, and ξ = 3.25 15 .The lengths are in units of Bohrs.The integrals are solved efficiently in real space using the VEGAS integration alogorithim contained in the GNU scientific library 16 .Table I shows all quantum mechanical scattering elements considered and their magnitude for κ = 1.When electrons are far enough away, we take a classical limit of the Coulomb elements, mainly we also take the elements:

III. CONFIGURATION INTERACTION BASIS
The many-body Hilbert space can be divided into smaller subspaces with total spin S and azimuthal spin S z .We construct the basis, in the occupation number representation, distributing particles among singleparticle states labeled with spin S z .The total number of possible configurations N cf for N e particles distributed on N st single particle states with a given spin N ↓ or N ↑ , arXiv:2403.07774v1 [cond-mat.mtrl-sci]12 Mar 2024 where N e = N ↓ + N ↑ is determined by a product of binomial coefficients, N cf = Nst N ↑ • Nst N ↑ .We do not rotate the Hamiltonian matrix to a S basis as this is an additional computational cost, and instead determine the ground state from calculations of expectation value of total spin S for each energy eigenstate.

IV. TRIANGULENE
We analyze a single triangulene molecule, a triangular graphene quantum dot with zigzag edges shown in Fig. 1(a) of the main text.This quantum dot has broken sublattice symmetry as seen by counting the number of red balls (carbon atoms belonging to sublattice A) and the number of gray balls (carbon atoms belonging to sublattice B).This is also a bipartite lattice, and as such, Lieb's theorem applies 17 .
We start by performing DFT calculations, and find the ground state to be S = 1 in agreement with Lieb's theorem and previous experimental and theoretical work [18][19][20] .We then perform HF calculations by solving the HF equation.We take κ = 3 giving the results shown in Fig. S1.At the top of the VB there are two degenerate states that are spin up, and a large gap that separates the spin down states.The splitting between the spin up and spin down states arises from a net spin-polarization.These states are localized at the edge, and the spin density of triangulene (shown as an inset in Fig. S1) shows the localized spin-1 quasi particle tends to localize at the edge of the triangle.

V. COMPARISON BETWEEN DIFFERENT MANY-BODY INTERACTING FERMIONIC MODELS -RESULTS FOR TWO TRIANGLES
We compare results obtained within different manybody Hamiltonians by restricting Coulomb matrix elements in real space to the dominant ones.Within the Hubbard model, one gets an effective mean-field Hamiltonian given as where the Hubbard parameter U = ⟨ii|V |ii⟩, ni = c † iσ c iσ , and ⟨n i,−σ ⟩ = ρ ii,σ ′ , are diagonal elements of density matrix with σ ′ = −σ.For the extended Hubbard model one has where V ij = ⟨ij|V |ji⟩.The Coulomb matrix elements in the basis of mean-field energies ⟨pq|V |rs⟩ are obtained The inset shows the spin density of the HF ground state.
using a basis rotation where b pσ = i B ip c iσ , and B ip are eigenvectors obtained by solving the mean-field Hamiltonian.The form of the many-body Hamiltonian is the same for Hubbard and extended Hubbard, and full interacting models given by Eq. 3 in the main article, but the models differ in the parameters used defined by ϵ HF pσ , τ pqσ and Coulomb elements given by Eq. 5. A real space spin density distribution for a given many-body state D, used when computing the spin density shown in Fig. 4 in the main article, is calculated using the formula where |D⟩ = Fig.S2 shows the many-body spectrum for the Hubbard model, extended Hubbard model and the fully interacting model for N Tr = 2 structure.Calculations in each case were done for CAS (6,6).We find that the Hubbard model tends to overestimate the spin gaps, with the extended Hubbard model capturing most of the quantitative features of the full model.Notice that we are showing only the three lowest energy states, as the excited states are separated from these states by a large gap on the order of hundreds of meV.

VI. EXTENDED HUBBARD MODEL RESULTS
FOR NTr = 4 TRIANGULENES Using the extended Hubbard model, we analyze the N Tr = 4 structure and discuss properties of longer triangulene chains.In Fig. S3, we show the Hartree-Fock energy spectrum indicating the energy gaps between the inter-triangulene states (sectors B) and the degenerate shell states (sector A and C), ∆ AB and ∆ BC .Within this model, the degenerate shell states split into a set of four double degenerate states.There are three intertriangulene states below and three above the degenerate shell.The energy gap ∆ AB is around twice smaller than ∆ BC , and this relation is approximately true also for longer chains, see the inset.These gaps determine the role of excitations and need to be related to scattering Coulomb matrix elements.For κ = 3 used here, we find some of scattering elements between sector A and B, and between B and C as large as 0.2 eV (e.g.⟨AA|v|BB⟩ and ⟨AB|v|BC⟩).Comparing them to energy gaps ∆ AB ∼ 1.8 eV, ∆ BC ∼ 0.7 eV, one can clude that the states from sector C are as important as the states from sector A, see also Section VIII, where we discuss the superexchange mechanism.
Fig. S4 shows charge densities in these three energy sectors; the densities are normalized by the number of states in a given sector.Sector A has localized charge density in the center of the chain, the central intertriangulene connection and on the first and last triangulenes.On can notice a higher charge density on carbon atoms from the same sublattice as the degenerate shell states, which confirms hybridization with states from sector B (within a tight-binding model, charge density from sector A is similar to that from sector C).In sector B charge density is mainly localized on the edges of triangulenes, and only one sublattice.Sector C has charge density on the three connections between the four triangulenes.
Fig. S5 shows the many-body spectra obtained using the Hamiltonian given by Eq. 3 in the main article with states from only sector B (8 states) and with all three sectors A, B, and C (14 states).Similar to the N Tr = 2 case, inclusion of the inter-triangulene states increases the splitting between singlet and triplet, and triplet and quintuplet.One can also see the order of total spin states (we show only the ten lowest energy states for the 8 state calculation) agrees with the order of states in BLBQ model for the lowest five states, compared with Fig. S6.The sixth state in the extended Hubbard model calculations with 14 states has S = 2 (seventh has S = 0).We attribute these differences between this extended Hubbard model and the BLBQ model to a CAS that is not large enough to converge excited states in this extended Hubbard model calculation.

VII. SPIN MODEL RESULTS
The antiferromagnetic spin-1 Heisenberg model with the added biquadratic operator can effectively describe the chain of triangulane molecules.The effective Hamiltonian is: where J is the coupling constant, β is the biquadratic term amplitude, and Si is a spin vector of the ith molecule in the chain.
Using the definition of the ladder operators Ŝ± i = Ŝx i ± i Ŝy i , one can rewrite the Hamiltonian as: The single spin-1 site has three possible Ŝz which gives us the matrix: Similarly, Hamiltonian matrices for other S z subspaces are: The full Hamiltonian matrix is block diagonal, and its eigenvalues for the singlet, triplet, and quintuplets are After shifting eigenvalues by J(1 + β), we see that only the singlet state depends on the β: E(S) = −3J(1 − β), E(T ) = −2J, and E(Q) = 0.For β < 1 3 the ground state is a singlet and for β > 1 3 it is a triplet state.For β = 1 3 , the singlet and triplet states are degenerate, as shown in Fig. S6(a).
Fig. S6(b) shows numerical results for a chain of N = 4 spin-1 sites.The transition between the singlet and triplet ground state again occurs at β = 1 3 , but the gap between triplet and quintuplet now depends on β.

VIII. SUPEREXCHANGE MECHANISM ANALYSIS
Isolated triangulene has a triplet ground state and can be represented by an effective spin-1 site.We analyze the processes responsible for coupling between neighboring spin-1 states on the example of two triangulenes.Our methodology for treating Coulomb interaction relies on a two-step process.First, we include interactions at the Hartree-Fock level for a closed shell system, and next, we populate unoccupied states up to charge neutrality and diagonalize the many-body Hamiltonian within a restricted subspace.For the N Tr = 2 system, after the HF basis rotation, the degenerate shell states from two triangles form symmetric and antisymmetric linear combinations of states from each triangle, as seen in the wavefunctions of Fig. 3  states and valence and conduction band inter-triangulene states.
This hinders the perturbative analysis.The two isolated triangulenes have four perfectly degenerate edge states, two for each triangulene.The perfect strategy would be to rotate the HF degenerate shell states back to the isolated triangulene basis states.In that case, one could show that coupling between triangulene spin-1 states are through the inter-triangulene states -an indirect superexchange mechanism, one of the main results of this work.However, after HF self-consistent calculations, the degenerate shell states are too strongly hybridized with the inter-triangulene states.
We analyze the superexchange mechanism at the many-body level.We take six HF states to construct our many-body Hilbert space, the highest valence band state, four degenerate shell states and the lowest conduction band state.These states are labelled 5-10 in Fig. 2(c) of the main text, in this analysis, for simplicity we shift the labels to 1-6.The charge neutral system has N el = 6 electrons.For S z = 0, one can construct in total, N cf = 400 configurations.The low-energy subspace, which we call later a single occupation subspace, corresponds to double occupation of the valence band state and single occupation of four degenerate shell states.One can construct six such low-energy configurations.In the occupation representation, these states can be written as where state 1 is the valence band state lying just below the four degenerate shell states, states 2 − 5 are the four degenerate shell states, and state 6 is the state just above the four degenerate states (unoccupied here).The six configurations defined in Eq. 12 can then be rotated into two total spin S = 0 singlets, three S = 1 triplets and one S = 2 quintuplet.We obtain the lowest energy state E 0 (S) within each total spin single occupation subspace by diagonalizing a 2 × 2 matrix within S = 0 subspace and a 3 × 3 matrix within S = 1 subspace, at the same time appropriately rotating the full Hilbert space within each total spin sector.This procedure is related to the strong coupling between configurations given by Eq. 12, due to the hybridization between the degenerate shell states and inter-triangulene states.Next, the obtained lowest energy states within each total spin subspace, are corrected by second order perturbation contributions due to the coupling to the rest of states E i (S) from the manybody Hilbert space (which have been rotated to the total spin subspaces), beyond the single occupation subspace (S = 0 subspace contains in total N cf = 175 configurations, S = 1 contains N cf = 189 configurations, and S = 2 contains N cf = 35 configurations).The perturbative Hamiltonian is written as where Ψ 0 and Ψ i are the many-body wave function corresponding to energy E 0 and E i , respectively, and H is the many-body Hamiltonian given by Eq. (3) in the main article.Although this procedure is second order in perturbation theory, if the hybridization between the edge states of the individual triangulenes and the intertriangulene states were weak, then perturbative treatment of the lowest energy total spin states (corresponding to singlet, triplet and quintuplet) would be appropriate.
In that case, fourth order perturbation theory would be required in order to describe coupling of the degenerate states of two triangulenes (and thus the spin-1 quasiparticles) through the inter-triangulene states.
The energy spectra after diagonalization within the single occupation subspace, and after second order perturbation correction is compared to the exact many-body spectrum in Fig. S7.The order of total spin states after basis rotation within a single occupation subspace is the same as within truncated N st = 4 states with the S = 2 quintuplet as the ground state, see Fig. 4(a) in the main article.Coupling of the single occupation subspace to higher energy configurations leads to a change of the order of states that now agrees with predictions within the two spin-1 Heisenberg Hamiltonian, with the singlet as the ground state.Furthermore, results after this coupling, agrees as well as with the experiments 21 .The energies are close to the exact energies obtained after diagonalization of the full many-body Hamiltonian.We notice that configurations mainly contributing to the perturbation includes all six HF states, thus both the state below and the state above the four degenerate shell states are important, which confirms the existence of an indirect AFM superexchange mechanism.

FIG. 1 .
FIG. 1. Geometrical representation of (a) triangulene and (b) the triangulene dimer.The balls denote location of carbon atoms with the red balls corresponding to sublattice A and the gray balls corresponding to sublattice B. The sticks connecting carbon atoms refer to the alternating single and double bonds.In (c) we show the single-particle spectrum of the triangulene dimer for states near the Fermi level.The blue horizontal lines are HF results, while the red crosses correspond to DFT results obtained with the B3-LYP exchange-correlation functional.The dotted line corresponds to the Fermi level at half-filling

FIG. 2 .
FIG. 2. Hartree-Fock (a-f) and DFT (g-l) wavefunctions for the six states located closest the the Fermi level at half-filling.(a-f) and (g-l) corresponds to the states 5-10 shown in Fig. 1(c).

FIG. 3 .
FIG. 3. Charge neutral effective HF hopping elements for (a) nearest neighbor t (b) next nearest neighbor t ′ (c) next-next nearest neighbor t 3 .The black dots correspond to the location of carbon atoms, and the line connecting them corresponds to the effective hopping between carbon atoms.

FIG. 4 .
FIG. 4. (a) Many-body spectrum for N Tr = 2 triangulenes as a function CAS size.The x-axis denotes the number of states taken relative to the dotted line shown in Fig. 1 (c), i.e. six states means taking 3 states above and 3 states below the dotted line.The "x" data points overlaying the data for 14 states corresponds to the spin model results for J = 8.21 meV and β = 0.016.b) Comparison of the many-body spectrum for N Tr = 4 triangulenes with corresponding BLBQ spin model results.In (c) we show the many-body spin density for the lowest energy triplet state for N Tr = 4 indicated by a circle in (b).

Tr = 2 )
and inter-triangulene valence and conduction band states.This makes it difficult to explicitly identify the superexchange processes responsible for the antiferromagnetic coupling between two S = 1 spins within perturbation theory.The degenerate states from isolated triangulene molecules are significantly different than the four HF degenerate shell states in the TD.We extend the discussion about identification of superexchange mechanism in the SM 45 .The presence of indirect coupling between spins in triangulene spin-1 chains opens a new possibility to control parameters in BLBQ model given by Eq. 4. While it might be difficult, or even impossible, to go beyond the Haldane phase on the phase diagram 8 , slight variations of β might be possible by tuning the energy gap between the inter-triangulene states and the degenerate shell states, which translates into control of the singlet-triplet splitting.The superexchange mechanism discussed in our work for triangulene spin-1 chains FIG. S1.Single-particle spectrum of triangulene NTr = 1 at half-filling.The horizontal bars are HF results, while the crosses correspond to DFT results.Red color corresponds to spin up levels while blue color corresponds to spin down levels.The inset shows the spin density of the HF ground state.

cA
FIG. S2.Many-body spectrum obtained for the Hubbard model, extended Hubbard model, and fully interacting model for CAS(6,6) calculations.
FIG. S3.Hartree-Fock spectrum from the extended Hubbard model for NTr = 4. Three sets of states, A, B, and C are indicated, with energy gaps between them, ∆AB and ∆BC.The inset shows the scaling of the energy gaps with the number of triangulenes NTr in a chain.


FIG. S4.Charge densities of Hartree-Fock spectrum from extended Hubbard model for NTr = 4 in three energy sectors indicated in Fig. S3.

3 FIG
FIG. S6.The BLBQ model results as a function of β parameter for spin chains with (a) N = 2 and (b) N = 4 sites.The blue dashed line indicates the spectra from the Fermionic model in the main text.
i projections: |↑⟩ = |1⟩, |↓⟩ = |−1⟩, and |0⟩.To obtain analytical results, we construct the basis for N=2 chain sites, composed of all combinations of Ŝz i values.The Ŝz = 0 subspace is spanned by three basis vectors: |0, 0⟩, |−1, 1⟩, |1, −1⟩.Acting with the Hamiltonian operator on this subspace results in: FIG. S7.Perturbation analysis of coupling between spins.0 indicates the lowest energy total spin states obtained from diagonalization of many-body Hamiltonian matrix within a single occupation subspace of the degenerate shell states.Perturbation indicates energies of these states after energy corrections included within 2nd order perturbation theory.Exact are energies obtained after diagonalization full many-body Hamiltonian.