Theory of triangulene two-dimensional crystals

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 nearest-neighbor 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 electron–hole 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 figure 2(a) 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 (figure 1(a)). We denote their zero mode wave functions by $|a\rangle$ and $|b\rangle$, respectively.

Zoom In Zoom Out Reset image size

Table 1. From left to right: size of the triangulene, labels for the corresponding n − 1 zero-energy modes, phase that each zero mode wave function acquires upon the action of the $R_{2\pi/3}$ operator (see equation (1) and figure 2(b)), and values of each zero mode wave function at the binding sites $j = 1,2$. The label z stands for either a or b.

n	zero mode	ωz	z(1)	z(2)
2	z0	0	$\frac{1}{\sqrt{6}}$	$\frac{-1}{\sqrt{6}}$
3	$z_+$	$+\frac{2\pi}{3}$	$\frac{1}{\sqrt{11}}$	$\frac{-1}{\sqrt{11}}$
3	$z_-$	$-\frac{2\pi}{3}$	$\frac{1}{\sqrt{11}}$	$\frac{-1}{\sqrt{11}}$
4	z0	0	$\frac{1}{\sqrt{12}}$	$\frac{-1}{\sqrt{12}}$
4	$z_+$	$+\frac{2\pi}{3}$	$\frac{1}{\sqrt{21}}$	$\frac{-1}{\sqrt{21}}$
4	$z_-$	$-\frac{2\pi}{3}$	$\frac{1}{\sqrt{21}}$	$\frac{-1}{\sqrt{21}}$

The zero modes are separated from the finite-energy states by a large energy splitting, as shown in figure 2(a). 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 figure 2(b), the zero modes are sublattice-polarized, 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 \geqslant 3$, their representation is not unique. It is extremely useful to choose $|a\rangle$ and $|b\rangle$ as eigenstates of the $R_{2\pi/3}$ counterclockwise rotation operator:

$$\begin{align} R_{2\pi/3}|a\rangle = e^{i\omega_a}|a\rangle, \quad R_{2\pi/3}|b\rangle = e^{i\omega_b}|b\rangle, \end{align} \tag{ 1 }$$

where $e^{i3\omega_{a,b}} = 1$. We thus have three possible values for the exponents: $\omega_{a,b} = 0,\pm 2\pi/3$.

The n − 1 zero modes of $[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], which states that $2S = |N_A-N_B|$. For $[n]$triangulenes, $|N_A-N_B| = n-1$, so that $2S = n-1$. The underlying mechanism for the ferromagnetic intra-triangulene exchange is a molecular version of the atomic Hund's coupling.

Let us now consider the unit cell of a $[n_a,n_b]$triangulene crystal, i.e. a dimer made of two triangulenes of A and B types. Separately, these triangulenes would have $n_{a/b}-1$ zero modes each. When 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 (figure 2(b)). 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.

The minimal Hamiltonian for the $[2,2]$triangulene crystal is a $2\times2$ 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 $\omega_z = 0$. Using that $z(1) = -z(2) = \frac{1}{\sqrt{6}}$ (see table 1), we get $t^{(0)}_{ab} = t^{(1)}_{ab} = t^{(2)}_{ab} = t_3/3$. Thus, we obtain the following effective Bloch Hamiltonian:

$$\begin{align} \mathcal{H}^{[2,2]}_\text{eff}(\phi_1,\phi_2) = \frac{t_3}{3} \begin{pmatrix} 0 & f (\phi_1,\phi_2) \\ f^{\,*} (\phi_1,\phi_2) & 0 \end{pmatrix}, \end{align} \tag{ 8 }$$

with

$$\begin{align} f (\phi_1,\phi_2) = 1 + e^{i\phi_1} + e^{i\phi_2} \end{align} \tag{ 9 }$$

and

$$\begin{align} \phi_{1,2} = \vec{k} \cdot \vec{R}_{1,2}. \end{align} \tag{ 10 }$$

The corresponding energy bands are given by:

$$\begin{align} \begin{array}{*{20}{l}} {\varepsilon _ \pm ^{[2,2]}({\phi _1},{\phi _2})}\\ {\,\,\,\,\,\,\, = \pm \frac{{{t_3}}}{3}\left| {f({\phi _1},{\phi _2})} \right|}\\ {\,\,\,\,\,\,\, = \pm \frac{{{t_3}}}{3}\sqrt {3 + 2\,\cos {\phi _1} + 2\,\cos {\phi _2} + 2\,\cos ({\phi _1} - {\phi _2})} .} \end{array}\end{align} \tag{ 11 }$$

The Hamiltonian obtained in equation (8) is identical to the TB model for the p_z orbitals of graphene, but with an effective hopping given by:

$$\begin{align} t^{[2,2]}_\textrm{eff} = \frac{t_3}{3}, \end{align} \tag{ 12 }$$

which leads to an energy bandwidth of $6t^{[2,2]}_\textrm{eff} = 2t_3$. The bandwidth obtained by our DFT calculations is around 0.6 eV, which implies $t_3~\simeq 0.3$ eV, in line with the values from the literature [43]. A Taylor expansion of equation (8) around the K and K^ʹ points yields Dirac cones with Fermi velocity $\hbar v_\mathrm{F}^{[2,2]} = \frac{5}{2} t_3 d$, where d denotes the carbon–carbon distance.

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₃, instead of by the nearest neighbor hopping.

The $[2,3]$triangulene (see figure 6(a)) is the smallest non-centrosymmetric triangulene 2D crystal. 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 1, the effective Bloch Hamiltonian can be written as:

$$\begin{align} & \mathcal{H}^{[2,3]}_\text{eff} (\phi_1,\phi_2)\nonumber\\ & \quad\!\!\!\!\! = \begin{pmatrix} 0 & \mathcal{F}^{[2,3]}_{a_0,b_+} (\phi_1,\phi_2) & \mathcal{F}^{[2,3]}_{a_0,b_-} (\phi_1,\phi_2) \\ {\mathcal{F}^{[2,3]}_{a_0,b_+} }^* (\phi_1,\phi_2) & 0 & 0 \\ {\mathcal{F}^{[2,3]}_{a_0,b_-} }^* (\phi_1,\phi_2) & 0 & 0 \end{pmatrix}, \end{align} \tag{ 13 }$$

with

$$\begin{align} \mathcal{F}^{[2,3]}_{a_0,b_\pm} (\phi_1,\phi_2) = \frac{2~t_3}{\sqrt{6 \times 11}} f(\phi_1 \pm \theta, \phi_2 \mp \theta) \end{align} \tag{ 14 }$$

and $\theta = 2\pi/3$. The corresponding eigenvalues read as:

$$\begin{align} \epsilon^{[2,3]}_0 (\phi_1,\phi_2) = 0 \end{align} \tag{ 15 }$$

and

$$\begin{align} &\epsilon^{[2,3]}_\pm (\phi_1,\phi_2)\nonumber \\ & \quad = \pm \sqrt{\left| \mathcal{F}^{[2,3]}_{a_0,b_+} (\phi_1,\phi_2) \right|^2 + \left| \mathcal{F}^{[2,3]}_{a_0,b_-} (\phi_1,\phi_2) \right|^2} \nonumber \\ & \quad = \pm \frac{2~t_3}{\sqrt{33}} \sqrt{3 - \cos \phi_1 - \cos \phi_2 - \cos(\phi_1 - \phi_2) }. \end{align} \tag{ 16 }$$

In figure 6(d), we plot these energy bands, which are verified to yield an excellent agreement with the low-energy bands of the full TB model (figures 6(b) and (c)).

Zoom In Zoom Out Reset image size

At the Γ point, we have $\phi_1 = \phi_2 = 0$ and therefore $\mathcal{F}^{[2,3]}_{a_0,b_\pm} (0,0) = 0$, which accounts for the triple degeneracy of the bands at that point. A Taylor expansion of equation (13) around Γ leads to the following expression for the dispersive bands:

$$\begin{align} \epsilon^{[2,3]}_\pm (\vec{k}) \simeq \pm 6~\sqrt{\frac{3}{11}} t_3 d |\vec{k}|, \end{align} \tag{ 17 }$$

that accounts for the single Dirac cone. The corresponding Fermi velocity, $\hbar v_\mathrm {F}^{[2,3]} = 6~\sqrt{\frac{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 [44].

For $n_a = n_b = 3$, there are two zero modes per triangulene, that we label as $|a_\pm \rangle$ and $|b_\pm \rangle$. The pairs of zero modes $|z_+ \rangle$ and $|z_- \rangle$ can be distinguished by the following property: $R_{2\pi/3} |z_\pm \rangle = e^{\pm i 2~\pi / 3} |z_\pm \rangle$. Thus, we can think of the phases $\omega_z = \pm 2\pi/3$ as a pseudospin degree of freedom.

Using the minimal model, the $2\times2$ hopping matrix of the effective Bloch Hamiltonian reads as:

$$\begin{align} \tau^{[3,3]} (\phi_1, \phi_2) = \begin{pmatrix} \mathcal{F}^{[3,3]}_{a_+,b_+} (\phi_1,\phi_2) & \mathcal{F}^{[3,3]}_{a_+,b_-} (\phi_1,\phi_2) \\ \mathcal{F}^{[3,3]}_{a_-,b_+} (\phi_1,\phi_2) & \mathcal{F}^{[3,3]}_{a_-,b_-} (\phi_1,\phi_2) \end{pmatrix}, \end{align} \tag{ 18 }$$

where

$$\begin{align} \mathcal{F}^{[3,3]}_{a_+,b_+} (\phi_1,\phi_2) = \mathcal{F}^{[3,3]}_{a_-,b_-} (\phi_1,\phi_2) = \frac{2t_3}{11} f (\phi_1,\phi_2), \end{align} \tag{ 19 }$$

$$\begin{align} \mathcal{F}^{[3,3]}_{a_+,b_-} (\phi_1,\phi_2) = \frac{2t_3}{11} f (\phi_1 +\theta,\phi_2 - \theta) \end{align} \tag{ 20 }$$

and

$$\begin{align} \mathcal{F}^{[3,3]}_{a_-,b_+} (\phi_1,\phi_2) = \frac{2t_3}{11} f (\phi_1 -\theta,\phi_2 + \theta). \end{align} \tag{ 21 }$$

Thus, we see that the pseudospin conserving terms lead to the same type of matrix elements than in the conventional honeycomb graphene model, whereas the pseudospin flip terms contain additional $\theta = 2\pi/3$ phases.

Diagonalizing the minimal Hamiltonian, we obtain

$$\begin{align} \epsilon_{1,\pm}^{[3,3]} (\phi_1,\phi_2) = \pm 3~t_\mathrm{eff}^{[3,3]}, \end{align} \tag{ 22 }$$

that accounts for the flat bands at finite energy, and

$$\begin{align} \epsilon_{2,\pm}^{[3,3]} (\phi_1,\phi_2) = \pm t_\mathrm{eff}^{[3,3]} | f (\phi_1,\phi_2)|, \end{align} \tag{ 23 }$$

that accounts for the graphene-like bands, with an effective hopping $t_\mathrm{eff}^{[3,3]} = 2~t_3 / 11$. The dispersive graphene-like bands feature Dirac cones at the K and K^ʹ points ($\phi_1 = -\phi_2 = \pm 2~\pi /3$), with Fermi velocity $\hbar v_F ^{[3,3]} = \frac{21}{11} t_3 d$, which is smaller than those obtained for the $[2,2]$ and $[2,3]$ cases.

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–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 electron–hole partner. Using numerical diagonalization for finite size clusters, we have verified that, indeed, every hexagonal ring hosts two localized states with energies $\pm 3~t_\mathrm{eff}^{[3,3]}$. These states are bonding, in the sense that both sublattices participate, and are thus different from E = 0 flat bands, which are sublattice-polarized.

In contrast to the previous triangulene 2D crystals, the energy bands of the $[4,4]$triangulene feature a gap at half-filling. 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 figures 5(g)–(i). The three topmost valence bands, and their conduction band electron–hole partners, are identical to the bands of a Kagome lattice: a flat valence/conduction band that is degenerate at the Γ point with the top/bottom of graphene-like bands. The lowest (highest) energy conduction (valence) band is thus dispersionless, which is of particular interest for the field of flat band physics.

Given that for the $[4,4]$triangulene crystal the number of zero modes per triangulene (3) is commensurate with the coordination number of the honeycomb lattice (3), it is natural to associate the band-gap to the formation of an orbital valence bond solid (VBS). This idea is further reinforced by the fact that a gapped band structure is also obtained for the $[7,7]$triangulene crystal, which has 6 zero modes per triangulene. In order to test this idea, and following previous work [41], we build a basis set of zero modes that are localized around one corner. To do so, we take advantage of the isomorphism of the states of a N = 3 ring TB model and the three zero modes chosen as eigenstates of the C₃ operator, $|z_0\rangle, |z_\pm \rangle$, that for the sake of this discussion we relabel as $| \omega \rangle$ with $\omega = 0,\pm\frac{2\pi}{3}$, respectively. In the N = 3 ring we have:

$$\begin{align} |\omega \rangle = \frac{1}{\sqrt{3}}\sum_{\alpha = 1}^{3} e^{i \omega \alpha }|c_\alpha \rangle. \end{align} \tag{ 24 }$$

We now invert this equation to obtain the states $|c_\alpha\rangle$ as linear combinations of the C₃ eigenstates. This results in three zero modes peaked around one corner (see figure 7). In the basis of corner states $|c_\alpha\rangle$, with $\alpha = 1,2,3$, the resulting inter-triangulene Hamiltonians $\tau_\alpha$ across the $\alpha = 1,2,3$ corners can be approximately written as:

$$\begin{align} \tau_1 = \tilde{t} \left(\begin{array}{ccc} 1 & \eta & \eta\\ \eta & 0 &0\\ \eta & 0 & 0 \end{array} \right), & \ \tau_2 = \tilde{t} \left(\begin{array}{ccc} 0 & \eta & 0\\ \eta & 1 &\eta\\ 0 & \eta & 0 \end{array} \right), \nonumber\\ & \ \tau_3 = \tilde{t} \left(\begin{array}{ccc} 0 & 0& \eta\\ 0 & 0 &\eta\\ \eta & \eta & 1 \end{array} \right),\end{align} \tag{ 25 }$$

with $\tilde{t} = 0.035{\text{t}}$, η = 0.097. The terms dropped in these matrices are all order $10^{-2} \tilde{t}$. The dominant $\tilde{t}$ matrix elements correspond to face to face corner states, whereas the $\eta\tilde{t}$ terms correspond to hopping from 'second neighbour corners'. If we had η = 0, then we would have a dispersionless VBS, where every electron is paired to just one electron. The finite value of η brings about the dispersion.

Zoom In Zoom Out Reset image size

Analytical expressions for the energy bands can be obtained [41], yielding:

$$\begin{align} \epsilon^{[4,4]}(\phi_1,\phi_2) = \pm \tilde{t}(1-2\eta), \ \pm \tilde{t} \left( 1+ \eta \pm \eta \left| f(\phi_1,\phi_2) \right| \right). \end{align} \tag{ 26 }$$

These analytical expressions are in very good agreement with the bands computed by numerical diagonalization of the Hamiltonian without dropping small terms in the $\tau_\alpha$ matrices above. We have also verified that the states in the flat bands are also ring states, very much like those of the $[3,3]$triangulene crystal.

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, $\ell$. For atoms, the orbital degeneracies are given by $2\ell+1$, with $\ell = 0,1,2,\ldots$, 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,\ldots$, and their symmetry properties are determined by the discrete $R_{2\pi/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 [45].

Finally, we note that even though the dispersive low-energy bands are narrow, interactions can turn these systems into Mott insulators at half-filling, as we discuss in the next section.

We now discuss qualitatively the effect of interactions. Unless otherwise stated, we focus on the half filling case. At the single $[n]$triangulene level, interactions have a huge impact on the electrons at the zero modes, determining the strong intramolecular ferromagnetic exchange that leads to $2S = n-1$. In the case of $[n_a,n_b]$triangulenes, Coulomb repulsion competes with the inter-triangulene hybridization. We quantify this competition in terms of the ratio between two energy scales, that we can define for each zero mode. On the one hand, we have the effective Coulomb interaction,

$$\begin{align} \tilde U_a = U \sum_i |a(i )|^4. \end{align} \tag{ 27 }$$

This effective interaction is thus given by the atomic Hubbard U times the inverse participation ratio (IPR). On the other hand, the average inter-triangulene hopping energy,

$$\begin{align} K_a = \frac{1}{n_b-1}\sum_{b} t^{(0)}_{ab}, \end{align} \tag{ 28 }$$

where $t^{(0)}_{ab}$ is the intracell hopping, given in equation (5), and $n_b-1$ is the number of zero modes of the n_b triangulene. The average of the ratio $\tilde U_a/K_a$ over the a zero modes provides a single figure of merit to quantify the strength of the interactions. The ratios depend both on the values of atomic energy scales U and t₃ and on the properties of the molecular orbitals, $t^{(0)}_{ab}, \tilde U_a$.

For instance, in the case of the $[2,2]$triangulene crystal, with just a single zero mode per triangulene, the IPR of the zero mode is $1/6$, so that the ratio yields $r\equiv \frac{\tilde{U}}{t^{(0)}_{ab}} = \frac{U/6}{t_3/3} = \frac{U}{2t_3}$. For U = 1.5 t and $t_3 = 0.1{\text{t}}$, we have $r\simeq 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 [46]. 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

$$\begin{align} J^{[2]} = 4\frac{\left(t^{[2,2]}_\textrm{eff}\right)^2}{\tilde{U}} = 6\times \frac{4}{9} \frac{t_3^2}{U} \simeq 48\,\textrm{meV}, \end{align} \tag{ 29 }$$

for t = 2.7 eV, U = 1.5 t and $t_3 = 0.1{\text{t}}$. As we shall discuss in a forthcoming work [47], in addition to the kinetic exchange, there are other contributions to the antiferromagnetic exchange in this system.

For $[n,n]$triangulenes with $n = 3,4,5$ we find the average of ratios $\tilde U_a/K_a$ to be 3.9, 2.9 and 3.1, respectively, if we take typical values [19, 20] of U = 1.5 t and $t_3 = 0.1{\text{t}}$. Thus, interactions cannot be neglected and will very likely drive these systems into a Mott insulating phase.

We now discuss the results of our spin-polarized DFT calculations on the $[3,3]$triangulene crystal. 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 one-dimensional (1D) chains of triangulene crystals [20]. The results are shown in figure 8 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 figure 8(b), 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.

Zoom In Zoom Out Reset image size

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 $6 J^{[3]} S^2 = E_{\mathrm {FM}}-E_{\mathrm {AFM}} = 159.2$ meV. Using S = 1 we infer $J^{[3]} = 26.5$ meV. This analysis, that overlooks additional contributions such as biquadratic exchange, is to be compared with the result of J = 14 meV, recently inferred for the S = 1 triangulene dimer, in the same approximation [19].

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 filling case. These show a transition to a Neel ordered Mott insulating phase for U > 4.2 t. The existence of a spin liquid phase has been studied extensively using Heisenberg spin models [48]. 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 $S = |n_a-n_b|/2$ 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 calculations in finite size clusters [49]. 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 Affleck–Kennedy–Lieb–Tasaki (AKLT) model for $S = 3/2$, that contains Heisenberg (bilinear) terms as well as biquadratic and bicubic spin couplings [50]. Importantly, the AKLT model can be solved exactly and features a spin VBS. The ground state of the AKLT model is a universal resource for measurement based quantum computing (MBQC) [51]. Importantly, the AKLT model has a gap in the excitation spectrum [52–54], 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, and even bicubic [50] terms in the honeycomb case. In one dimension, the 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 diagonalization of the spin model [55] 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 one dimension 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 [56]. This matter deserves further study.