Theory of triangulene two-dimensional crystals

Equilateral triangle-shaped graphene nanoislands with a lateral dimension of $n$ benzene rings are known as $[n]$triangulenes. Individual $[n]$triangulenes are open-shell molecules, with single-particle electronic spectra that host $n-1$ half-filled zero modes and a many-body ground state with spin $S=(n-1)/2$. The on-surface synthesis of triangulenes has been demonstrated for $n=3,4,5,7$ and the observation of a Haldane symmetry-protected topological phase has been reported in chains of $[3]$triangulenes. Here, we provide a unified theory for the electronic properties of a family of two-dimensional honeycomb lattices whose unit cell contains a pair of triangulenes with dimensions $n_a,n_b$. Combining density functional theory and tight-binding calculations, we find a wealth of half-filled narrow bands, including a graphene-like spectrum (for $n_a=n_b=2$), spin-1 Dirac electrons (for $n_a=2,n_b=3$), $p_{x,y}$-orbital physics (for $n_a=n_b=3$), as well as a gapped system with flat valence and conduction bands (for $n_a=n_b=4$). All these results are rationalized with a class of effective Hamiltonians acting on the subspace of the zero-energy states that generalize the graphene honeycomb model to the case of fermions with an internal pseudospin degree of freedom with $C_3$ symmetry.


I. INTRODUCTION
The quest for new states of matter that do not occur naturally in conventional materials fuels the study of artificial quantum lattices in a variety of platforms, including cold atoms 1,2 , trapped ions 3 , quantum dots 4,5 , dopants in silicon 6 , functionalized graphene bilayer 7 , moiré heterostructures 8 and adatoms [9][10][11] . These systems may work as quantum simulators of both spin and fermionic Hamiltonians, such as the Hubbard model, and are expected to host strongly correlated electronic phases.
Here, we explore the properties of two-dimensional (2D) artificial crystals with triangulenes as building blocks. Triangulenes are open-shell polycyclic aromatic hydrocarbons. Their low-energy degrees of freedom are provided by electrons that partially occupy zero-energy modes inside a large gap of strongly covalent molecular states formed by π atomic orbitals. In standard conditions, open-shell molecules are very reactive and therefore not suitable for manipulation. However, thanks to the advances in on-surface chemical synthesis [12][13][14] , it has been possible to go around this problem, so that the controlled fabrication of triangulenes of various sizes [15][16][17][18] , as well as triangulene dimers 19 , chains 20 , rings 20,21 and other structures 22 , has been recently demonstrated. On the other hand, the on-surface synthesis of large-area carbon-based crystals has also been established [23][24][25] . Therefore, the synthesis of 2D triangulene crystals seems within reach using state-of-the-art techniques, motivating the present work.
Individual [n]triangulenes have a ground state with spin S = (n − 1)/2, on account of their strong intramolecular exchange 26,27 .
The small spin-orbit coupling of carbon makes their magnetic anisotropy negligible 20,28 .
The intermolecular exchange coupling in [3]triangulene dimers has been determined to be antiferromagnetic [19][20][21]29 . Recently reported 20 experimental results show that chains and rings of [3]triangulenes display the key features of the Haldane phase for antiferromagnetically coupled S = 1 spins, namely a gap in the excitation spectrum and the emergence of fractional spin-1/2 edge states. Here, we study [n a , n b ]triangulene 2D crystals, i.e., honeycomb lattices with a unit cell made of two triangulenes with dimensions n a , n b (see Fig. 1 for the cases n a = n b = 2 and n a = 2, n b = 3).
The rest of this paper is organized as follows. In section II, we review the electronic properties of individual [n]triangulenes. We argue that, based on a nearest-neighbor tight-binding (TB) approximation, [n a , n b ]triangulene 2D crystals should have n a + n b − 2 flat bands at zero energy. In section III, we calculate the energy spectrum using density functional theory (DFT) and we find n a + n b − 2 narrow bands, in disagreement with the expectations based on the nearest-neighbor TB model. We show that adding a third neighbor hopping to the TB model accounts for the DFT results and is essential to understand the weakly dispersive bands. In section IV, we build a minimal TB model that includes a) b) c) FIG. 1: (a) Two variants of [2]triangulene, with A (B) sublattice-type sites colored in red (blue). The top (bottom) triangulene has excess of red (blue) sites. (b) Centrosymmetric triangulene honeycomb lattice with a unit cell formed by a pair of [2]triangulenes with opposite sublattice imbalance. (c) Non-centrosymmetric triangulene 2D crystal with a unit cell formed by an A-type [2]triangulene and a B-type [3]triangulene.
only the zero-energy modes of the triangulenes and the effect of third neighbor hopping, compare it to both DFT and full TB models, and show that it accounts for the energy bands so obtained. These include a graphenelike spectrum for n a = n b = 2, spin-1 Dirac electrons for n a = 2, n b = 3, p x,y -orbital honeycomb physics for n a = n b = 3, and a gapped system with flat valence and conduction bands for n a = n b = 4. In section V, we address the effect of interactions and the spin physics in triangulene 2D crystals. In section VI, we wrap up and present our conclusions.

A. Spectrum and symmetries
In this section, we briefly review the electronic properties of individual [n]triangulenes. We first consider their single-particle properties, as described with the TB Hamiltonian, considering one orbital per site and nearestneighbor hopping, t. Because of the bipartite character of the lattice, the energy spectrum must be composed of two types of states: finite-energy states, with electronhole symmetry, and at least |N A −N B | zero modes, where N A/B is the number of sites of the A/B sublattice. In the case of [n]triangulenes (see Fig. 2a for n = 2, 3, 4), |N A − N B | = n − 1 and we find n − 1 zero modes. We can build triangulenes with excess of either A or B sites (Fig. 1a). We denote their zero mode wave functions by |a and |b , respectively.
The zero modes are separated from the finite-energy states by a large energy splitting, as shown in Fig. 2a. At half filling, relevant for polycyclic aromatic hydrocarbons at charge neutrality, the zero modes are half-filled. Therefore, the low-energy physics of triangulene crystals are expected to be dominated by the electrons occupying the zero-energy states.
We now focus on the wave functions of the zero modes. As we show in Fig. 2b, the zero modes are sublatticepolarized, i.e., they have a non-zero weight exclusively in the majority sublattice 27,30 . Because of the degeneracy of the zero modes for n ≥ 3, their representation is not unique. It is extremely useful to choose |a and |b as eigenstates of the R 2π/3 counterclockwise rotation operator: where e i3ω a,b = 1. We thus have three possible values for the exponents: ω a,b = 0, ±2π/3.

B. Interactions and magnetic properties
The n − 1 zero modes of the [n]triangulenes host n − 1 electrons. Therefore, [n]triangulenes have an open-shell structure. Calculations carried out with DFT 27,31 , quantum chemistry 30 and model Hamiltonians 27,30,32 consistently predict that the many-body ground state maximizes the spin of the electrons in the zero modes, so that 2S = n − 1. When modeled with the Hubbard model, the spin of the ground state can be anticipated using Lieb's theorem 33 The underlying mechanism for the ferromagnetic intratriangulene exchange is a molecular version of the atomic Hund's coupling.  Table I). The zero modes vanish at the minority sublattice, including the corner sites that link triangulenes to their neighbors.
coupled together, Sutherland's theorem 34 warrants a minimal number of |n a − n b | zero modes per unit cell, and the same number of zero-energy (E = 0) flat bands for the corresponding 2D crystal. Specifically, for n a = n b the theorem does not ensure the existence of E = 0 flat bands. However, the binding sites of the unit cell in the triangulene 2D lattice belong to the minority sublattice, whereas the zero modes are hosted in the majority sublattice (Fig. 2b). Therefore, nearest neighbor hopping does not hybridize zero modes of adjacent triangulenes. As a consequence, t does not lift the zero mode degeneracy in a [n a , n b ]triangulene crystal, so that a nearest neighbor TB model predicts n a + n b − 2 flat bands at zero energy. We anticipate that intermolecular hybridization of zero modes in triangulene 2D crystals will be governed by the small third neighbor hopping, leading to weakly dispersive half-filled energy bands.

III. NON-MAGNETIC ENERGY BANDS: DFT VS TB
We now discuss the energy bands of several [n a , n b ]triangulene 2D crystals, calculated with two methods: spin-unpolarized DFT and TB (see results for n a = n b = 2, 3, 4 in Fig. 3). The spin-unpolarized DFT calculations enforce non-magnetic solutions, so that interactions do not break time-reversal nor sublattice symmetry. Our DFT calculations were carried out with Quantum Espresso 35 , using the Perdew-Burke-Ernzerhof functional 36 . The edge carbon atoms were passivated with hydrogen. The kinetic energy cutoff considered was 30 Ry. The charge density and potential energy cutoffs used were 700 Ry. We employed k-grids of 8 × 8 × 1 for the [3,3]triangulene crystal and of 10×10×1 for the [2,2] and [4,4] cases. No significant deviations from planarity were found.
Below we show that the resulting low-energy bands are associated with the triangulene zero modes and can be accounted for by the TB calculations, as long as third neighbor hopping is included. We start by addressing the DFT bands obtained for the [2,2]-, [3,3]-and [4,4]triangulene crystals ( Fig. 3d-f). The details of these first-principle calculations can be found in Appendix ??. Our results are in agreement with previous work by the group of Feng Liu 37,38 for the [2,2] and [4,4] cases.
The overall picture of the spin-restricted DFT bands for the [n, n]triangulene, with n = 2, 3, 4, is the following. A set of 2(n − 1) weakly dispersive bands is located around the Fermi energy, well separated from higher/lower energy bands. Given that the number of zero modes per unit cell is 2(n − 1), it can be expected that the wave functions of these bands are mostly made of zero modes. In the case of n = 2 ( Fig. 3d), the bands are isomorphic to the π bands of graphene, but with a smaller bandwidth of 605 meV. The n = 3 bands, shown in Fig. 3e, feature two Dirac cones and, in addition, two flat bands at the maximum and minimum of the graphene-like Dirac bands. We note that these bands are identical to those obtained in the p x,y -orbital honeycomb model 39,40 . The spectrum of the [4,4]triangulene crystal ( Fig. 3f, see also a zoom in Fig. 5g) features a small gap of 0.17 eV, with flat valence/conduction bands and a pair of graphene-like bands whose top (bottom) is degenerate with the valence (conduction) flat band at the Γ point. We note that the same set of bands was found in crystalline networks where topological zero-energy states are hosted at the junctions of a network of graphene nanoribbons 41 .
The appearance of non-dispersive bands at finite energy is quite peculiar, as non-bonding states in bipartite lattices usually correspond to sublattice-polarized modes with E = 0. Here, the sublattice-polarized E = 0 states seem to hybridize, forming pairs of bonding-antibonding states that lead to electron-hole symmetric bands, some of which being flat bands with finite energy. The origin of the peculiar dispersion of these low-energy bands, expected to be all flat at E = 0 in the nearest-neighbor TB model, is addressed in the next section.
The dispersion of the low-energy bands obtained with DFT suggests that the zero modes of adjacent triangulenes do hybridize, so that the nearest neighbor TB approximation fails. Second neighbor hopping does not couple states hosted at different sublattices, including the zero modes of adjacent triangulenes, and, in addition, leads to bands that break electron-hole symmetry. Therefore, it is natural to consider a third neighbor hopping, t 3 . Our calculations with t 3 = 0.1t show energy bands in qualitative agreement with those of DFT for the [2, 2]-, [3,3]-and [4,4]triangulene crystals ( Fig. 3g-i). If we take t 3 = 0 (not shown), the n a + n b − 2 low-energy bands become flat, with E = 0. Thus, we conclude that third neighbor hopping accounts for the observed dispersion and is the key energy scale that controls intertriangulene hybridization.

IV. EFFECTIVE LOW-ENERGY THEORY
The TB results of the previous section illustrate that [n a , n b ]triangulene 2D crystals provide a platform to generate a wide class of non-trivial energy bands, including graphene-like Dirac dispersion, p x,y -orbital honeycomb physics and a gapped system with flat valence and conduction bands. The emergence of these peculiar results calls for a deeper comprehension.
In order to ascertain the ultimate origin of the peculiar properties of the low-energy bands of [n a , n b ]triangulenes, we define an effective Hamiltonian in terms of the n a + n b − 2 zero modes in the unit cell, calculated with the TB model with t 3 = 0 57 . These zero modes are split into two groups consisting of n a − 1 |a and n b − 1 |b states localized in the A and B sublattices of adjacent triangu- lenes, respectively. Third neighbor hopping connects |a and |b at three pairs of atomic dimers at the corners of the triangulenes (Fig. 4). On the other hand, t 3 does not couple zero modes of the same triangulene, as they belong to the same sublattice. The resulting effective model defines a honeycomb lattice with n a + n b − 2 states per unit cell, whose Bloch Hamiltonian can be written as: where k denotes the 2D crystal momentum, 0 [a/b] are square matrices of dimension n a/b − 1 with all entries equal to zero, and τ is an (n a − 1) × (n b − 1) rectangular matrix with entries given by: In the previous expression, R 0 = (0, 0) is the null vector associated to the intracell hopping terms, and R 1,2 are the primitive vectors that define the hexagonal lattice of [n a , n b ]triangulenes: where a 1,2 are the primitive vectors of the graphene hexagonal lattice. The ( k-independent) hopping matrix elements t (m) ab are defined below. As usual, the hopping between A and B states of the honeycomb lattice has one intracell (the m = 0 term of the sum in Eq. (3)) and two intercell (the m = 1, 2 terms of the sum in Eq. (3)) contributions.
The intracell hopping matrix elements can be expressed as: (1) and Fig. 2b), and values of each zero mode wave function at the binding sites j = 1, 2.
The label z stands for either a or b.
where a(j) = j|a and b(j) = j|b denote the components of the zero mode wave functions |a and |b at the atomic sites j = 1, 2, following the labeling depicted in Fig. 4. Importantly, we can always impose that a (1) In the following, we shall assume this gauge fixing. For that reason, we will use the notation z for cases where we want to label either a or b.
Using the point group symmetry of the triangulenes, the intercell hopping matrix elements can be expressed in terms of the intracell ones. Given the values of z(1) and z(2), together with the phases ω z that each zero mode wave function acquires when applying the R 2π/3 operator (see Table I for n = 2, 3, 4), the remaining values of z(j ) and z(j ) can determined by making use of the C 3 symmetry. Specifically, we use that z(j) = e iωz z(j ) = e −iωz z(j ) to obtain: and t (2) which fully clarifies how to build the effective Hamiltonian of Eq. (2).
The effective Hamiltonian described above is the minimal TB model that is expected to capture the n a + n b −2 narrow bands of [n a , n b ]triangulenes. This minimal model describes a honeycomb lattice whose unit cell has n a − 1 and n b − 1 states in the A and B sublattice-type sites, respectively. We can think of the different states of the same sublattice as a pseudospin degree of freedom with C 3 symmetry, thereby generalizing the honeycomb model of graphene to the case of fermions with an additional internal structure. The magnitude of the effective hopping is controlled by the weight of the zero modes at the intracell binding sites j = 1, 2. The intercell hopping matrix elementns depend on the phase differences ω a −ω b . When nonzero, we can interpret the phases ω z as Peierls phases associated to a pseudospin-dependent magnetic field. Importantly, the model must still be time-reversal invariant.
In Fig.5c,f,i, we show the energy bands obtained with the minimal TB model for n a = n b = 2, 3, 4. taking t 3 = 0.1t. We find an excellent agreement with the full TB model (Fig.5b,e,h), which in turn agrees with the spinunpolarized DFT results (Fig.5a,d,g). Below, we take an in-depth look into the effective Hamiltonian of Eq. (2), which allow us to have a better understanding of the band dispersion for n a = n b = 2, 3, 4, as well as for the non-centrosymmetric case n a = 2, n b = 3 (not addressed yet).

A. [2, 2]triangulene crystal
The minimal Hamiltonian for the [2,2]triangulene crystal is a 2 × 2 matrix that is isomorphic to the TB model for a monoatomic honeycomb lattice with one orbital per site and nearest neighbor hopping. The reason for that is the fact that there is only one zero mode per triangulene, with ω z = 0. Using that z(1) = −z(2) = 1 √ 6 (see Table I), we get t (0) ab = t 3 /3. Thus, we obtain the following effective Bloch Hamiltonian: with and The corresponding energy bands are given by: = ± t3 3 3 + 2 cos φ1 + 2 cos φ2 + 2 cos(φ1 − φ2).
The Hamiltonian obtained in Eq. (8) is identical to the TB model for the p z orbitals of graphene, but with an effective hopping given by: which leads to an energy bandwidth of 6t In short, we see that third neighbor hopping in [2,2]triangulene crystals produces an intermolecular hybridization of the zero modes that leads to the formation of an artificial graphene lattice, with a narrower bandwidth governed by t 3 , instead of by the nearest neighbor hopping.

B. [2, 3]triangulene crystal
The [2,3]triangulene is the smallest noncentrosymmetric crystal (Fig. 6a). It has a sublattice imbalance of |N A − N B | = 1 that accounts for the presence of one flat band at E = 0. Within the minimal model, it has three states per unit cell. Using the results of Table I, the effective Bloch Hamiltonian can be written as: with F [2,3] a0,b± (φ 1 , φ 2 ) = and θ = 2π/3. The corresponding eigenvalues read as: In Fig. 6d, we plot these energy bands, which are verified to yield an excellent agreement with the low-energy bands of the full TB model (Fig. 6b,c). At the Γ point, we have φ 1 = φ 2 = 0 and therefore F [2,3] a0,b± (0, 0) = 0, which accounts for the triple degeneracy of the bands at that point. A Taylor expansion of Eq. (14) around Γ leads to the following expression for the dispersive bands: that accounts for the single Dirac cone. The corresponding Fermi velocity, v [2,3] F = 6 3 11 t 3 d, is larger than that obtained in the [2, 2]triangulene crystal. It must also be noted that, in the neighborhood of Γ, the minimal Hamiltonian can be mapped into a spin-1 Dirac model 43 .

C. [3, 3]triangulene crystal
For n a = n b = 3, there are two zero modes per triangulene, that we label as |a ± and |b ± . The pairs of zero modes |z + and |z − can be distinguished by the following property: R 2π/3 |z ± = e ±i2π/3 |z ± . Thus, we can think of the phases ω z = ±2π/3 as a pseudospin degree of freedom.
The emergence of flat bands at finite energy is remarkable. These bands are associated to states that are localized around any given hexagonal ring in the honeycomb lattice [39][40][41] . Every unit cell of the [3,3]triangulene crystal must contribute with one state to each of the two flat bands. Thus, this makes half a state per triangulene and band. Every triangulene participates in three hexagonal rings in a honeycomb lattice. Therefore, every triangulene contributes with 1/6 of state per band to a given hexagonal ring. Thus, the six triangulenes that form a ring give rise to one localized state, plus its electronhole partner. Using numerical diagonalization for finite size clusters, we have verified that, indeed, every hexagonal ring hosts two localized states with energies ±3t [3,3] eff . These states are bonding, in the sense that both sublattices participate, and are thus different from E = 0 flat bands, which are sublattice-polarized.

D. [4, 4]triangulene crystal
In contrast to the previous triangulene 2D crystals, the energy bands of the [4,4]triangulene feature a gap at halffilling. As in the previous cases, the spin-unpolarized DFT bands close to E = 0 are well described both by the full TB Hamiltonian and by the minimal low-energy model, as shown in Fig. 5g-i. The three topmost valence bands, and their conduction band electron-hole partners, are identical to the bands of a Kagome lattice: a flat band that is degenerate at the Γ point with the top/bottom of two graphene-like bands. The lowest (highest) energy conduction (valence) band is thus dispersionless, which is of particular interest for the emergent field of flat band physics.
Given that, for the [4,4]triangulene crystal, the number of zero modes per triangulene (three) is a multiple of the coordination number of the honeycomb lattice (three), it is natural to associate its band gap to the formation of an orbital valence bond solid. This idea is further reinforced by the fact that a gapped band structure is also obtained for the [7,7]triangulene crystal, which has six zero modes per triangulene.
E. Overall picture of the single-particle states The low-energy bands of [n a , n b ]triangulene crystals are remarkable for several reasons. First, they only disperse because of a third neighbor hopping, that would normally play a residual role, yet it is dominant here. Second, they feature flat bands at finite energy, which are interesting on their own right. Third, they provide a generalization of the honeycomb TB model, with a single orbital per site, to a more general class of honeycomb models where fermions have an internal pseudospin degree of freedom.
At this point, we can draw an analogy between this pseudospin degree of freedom and the orbital angular momentum of atoms, . For atoms, the orbital degeneracies are given by 2 + 1, with = 0, 1, 2, ..., and their symmetry properties are governed by the spherical harmonics. For [n]triangulenes, the orbital degeneracies are given by n − 1, with n = 2, 3, 4, ..., and their symmetry properties are determined by the discrete R 2π/3 operator. Very much like the angular momentum wave function of the valence electrons in atoms shapes their electronic properties, the pseudospin of the fermions in [n a , n b ]triangulene crystals dictates the resulting band structures. In particular, the [2,2] case leads to artificial graphene, [2,3] to the S = 1 Dirac Hamiltonian, and [3,3] to the p x,y -orbital honeycomb model, all of which featuring Dirac cones at the Fermi energy. The [4,4]triangulene, in contrast, is a band insulator, with flat valence and conduction bands. It must be noted that, even though we have been always considering triangulene 2D crystals at half-filling, other options could be possible if carriers are injected into them, for instance via the application of a gate voltage or through chemical doping 44 .
Finally, we note that even though the dispersive lowenergy bands are narrow, interactions will turn these systems into Mott insulators at half-filling, as we discuss in the next section.

A. General considerations
We now discuss qualitatively the effect of interactions. Unless otherwise stated, we focus on the half-filling case. At the single triangulene level, interactions have a huge impact on the electrons at the zero modes, and they determine the strong intramolecular ferromagnetic exchange that leads to 2S = n − 1, for [n]triangulenes. In the case of triangulene crystals, Coulomb repulsion competes with inter-triangulene hybridization. We quantify this competition in terms of the ratio of two energy scales, that we can define for each zero mode. On the one hand we have the effective Coulomb interaction, On the other hand, the average inter-triangulene hopping energy, where τ ab (0) is the expression given in Eq. for a given unit cell, we can then average the ratioŨ a /K a over the zero modes. The results depend both on the values of atomic energy scales U and t 3 and on properties of the molecular orbitals, such as their inverse participation ration (IPR), given by i |a(i)| 4 , and the hopping matrix elements. We find that, for all triangulene systems considered here, this ratio is systematically much larger than 1 if we take typical values 19,20 of U = 1.5t and t 3 = 0.1t. Thus, we expect that, at half filling, triangulene 2D crystals will be Mott insulators. Perhaps the only exception could be the [4,4]triangulene that has a gap before interactions are included.
We illustrate our point with the case of the [2,2]triangulene dimer, for which the IPR of the zero mode is 1/6, so that the ratio r ≡Ũ t [2,2] eff = U/6 t3/3 = U 2t3 . For U = 1.5t and t 3 = 0.1t, we have r 7.5, clearly in the insulating side of the metal-insulator transition predicted for r > 4.5 for the Hubbard model at half filling in the honeycomb lattice 45 . In the Mott insulating regime, with the charge fluctuations frozen, we expect the [2,2]triangulene crystal to behave like a S = 1/2 honeycomb lattice with antiferromagnetic exchange for t = 2.7eV , U = 1.5t and t 3 = 0.1t. As we shall discuss in a forthcoming work 46 , in addition to the kinetic exchange, there are other contributions to the antiferromagnetic exchange in this system.

B. Spin-polarized DFT calculations for the [3, 3]triangulene crystal
We now discuss the results of our spin-polarized DFT calculations on the [3,3]triangulene crystal (see details in Appendix ??). They provide a realistic assessment of the emergence of local moments and their exchange interactions but can only describe broken symmetry magnetic states, and therefore cannot account for quantum disordered spin phases that are known to occur for 1D chains of triangulene crystals 20 . The results are shown in Fig. 7 for the lowest energy magnetic configuration, where the triangulenes have local moments with opposite magnetizations. The four lowest energy bands clearly show that a very large gap has opened where the Dirac bands once stood. Additionally, from the magnetic moments distribution across the molecules, shown in Fig. 7b, it is apparent that local moments are hosted by the zero modes. We note that spin-polarized DFT results for the [2,2]triangulene crystal obtained by Zhou and Liu 37 , also feature an antiferromagnetic insulator. Therefore, the available spin-polarized DFT results support the picture that [n, n]triangulenes are Mott insulators with antiferromagnetically coupled local moments.
From our calculations we can infer the value of the intermolecular exchange J [3] . For that matter, we compare the energy difference between the ferromagnetic and the antiferromagnetic solutions obtained with DFT, with the result of a classical Heisenberg model. The energy difference per unit cell is given by 6J [3] S(S + 1) = E F M − E AF M = 159.2meV. Using S = 1 we infer J [3] = 26.5meV. This analysis, that overlooks additional contributions such as biquadratic exchange, is to be compared with the result of J = 14meV, recently inferred for the S = 1 triangulene dimer, in the same approximation 19 .

C. Spin phases in honeycomb crystals
We now briefly discuss the type of interacting ground states that can be expected to arise for the triangulene crystals at half filling. In the case of [2,2]-crystals, the relevant reference are the quantum Monte Carlo calculations for the Hubbard model in the half-limit case. These show a transition to a Neel ordered Mott insulating phase for U > 4.2t. The existence, or else, of a spin liquid phase has been studied extensively using Heisenberg spin models 47 . In the case of [2,3] at half filling, Lieb theorem 33 ensures a ground state with S = 1/2 per unit cell. According to Lieb theorem the non-centrosymmetryc crystals [n a , n b ] will have a spin 2S = |n a − n b | per unit cell, and will be either ferromagnetic or ferrimagnetic.
The S = 1 Heisenberg model in the honeycomb lattice has also been studied using density matrix renormalization group (DMRG) calculations in finite size clusters 48 . These calculations predict a Neel phase with long-range order, unless second-neighbor antiferromagnetic interactions, that promote frustration, are significant. For bipartite lattices, second-neighbor exchange coupling would be ferromagnetic. Therefore, a Neel phase is likely to occur for the [3,3]triangulene crystal at half filling.
In the case of the [4,4] triangulene crystal, the relevant spin model is the S = 3/2 Heisenberg Hamiltonian in the honeycomb lattice with antiferromagnetic interactions. Interestingly, this Hamiltonian is not far from the AKLT model for S = 3/2, that contains Heiseberg (bilinear) terms as well as biquadratic and bicubic spin couplings 49 . Importantly, the AKLT model can be solved exactly and features a spin valence bond-solid (VBS). The ground state of the AKLT model is a universal resource for measurement based quantum computing 50 (MBQC). Importantly, the AKLT model has a gap in the excitation spectrum [51][52][53] , so that cooling down a system described with the AKLT model at low enough temperatures would be enough to prepare a universal resource for MBQC.
The AKLT point is a singular point in a class of models, with bilinear and biquadratic couplings (BLBQ) and even bicubic 49 models in the honeycomb case. For 1D, VBS phase is robust in a region of parameters containing both AKLT and the pure Heisenberg limit. For the 2D honeycomb AKLT model, finite-size diagonalisation of the spin model 54 found a VBS to Neel transition when AKLT is linearly transformed into a Heisenberg model, but the precise location of the transition is hard to establish due to finite-size effects. In any event, the interesting question of whether the [4,4]triangulene crystal provides a realization of the VBS in two dimensions remains open. We note that the VBS phase predicted for the S = 1 AKLT in 1D, and it has been observed for [3]triangulene chains 20 .
Deviations from half-filling, introducing carriers in the Mott insulating phases, will very likely bring very interesting electronic phases, as it happens in the case of cuprates and in the case of twisted bilayer graphene 55 . This matter deserves further study.

VI. SUMMARY AND CONCLUSIONS
We have shown that graphene triangulenes are ideal building blocks for bottom-up design of carbon-based 2D honeycomb crystals with both narrow and flat bands at the Fermi energy. Our study permits to set the rules for the rational design of the bands, including the number of narrow bands, and the presence of flat bands either at the Fermi energy or nearby. The design principles for honeycomb crystals with a unit cell made of two triangulenes, of dimensions n a , n b (see Figure 1) are: 1. Isolated [n]triangulenes host n − 1 zero modes that are localized in one sublattice (either a or b, depending on the orientation) (see Figure 1) 2. The zero modes vanish identically at the intertriangulene binding sites. As a result intermolecular hybridization is governed by third neighbor hopping t 3 .
3. Centrosymmetryc [n, n]triangulene crystals have 2(n − 1) narrow bands. For n ≥ 3 they always feature pairs of electron-hole conjugate low energy flat bands, corresponding to states localized in supramolecular hexagonal rings. These flat bands are very different from E = 0 bands, as they are not sub-lattice polarized and they feature intermolecular hybridization.
4. Non centrosymmetric [n a , n b ]triangulene crystals, with n a = n b , feature |n a − n b | flat bands at E = 0, because of the bipartite character of the lattice 30,34 The results above are validated comparing spinunpolarized DFT calculations with tight-binding Hamiltonians at two levels: a full-lattice model that includes first and third neighbor hopping and a minimal tightbinding model where only the zero modes of the triangulenes are included. Using the latter approach we are able to derive analytical expressions for the energy bands of the [2,2], [2,3], and [3,3] triangulene crystals.
We have also addressed the effect of electron-electron interactions on these crystals. Both analytical estimates of the relevant energy scales and spin-polarized DFT calculations permit to anticipate that interactions will dominate the electronic properties of these crystals. At halffilling these crystals are expected to be Mott insulators, with their low energy physics governed by the spin degrees of freedom. Application of Lieb theorem 33 valid for the Hubbard model in bipartite lattices at half filling permits to anticipate that compensated ([n, n]) triangulene crystals have a ground state with S = 0, and therefore antiferromagnetic intermolecular exchange. We have estimated the magnitude of this super-exchange for the [2,2] and [3,3] crystals, and found very large values, in the range of tens of meV. For non-centrosymmetric [n a , n b ] crystals with n a = n b Lieb theorem predicts that the ground state will have 2S = n a − n b per unit cell. Thus, these crystals will be ferromagnets or ferrimagnets.
Doping these Mott insulators away from half-filling seems an almost certain recipe to discover non-trivial correlated electronic phases. We hope our work, together with recent work showing the great interest of triangulene one dimensional structures 20 will motivate experimental work to figure out synthetic routes that permit to create triangulene 2D crystals, and the exploration of their electronic properties in the lab.