# Theory of triangulene two-dimensional crystals

R. Ortiz1,2 G. Catarina1,3,4, J. Fernández-Rossier3111On permanent leave from
Departamento de Física Aplicada, Universidad de Alicante, 03690 San Vicente
del Raspeig<EMAIL_ADDRESS>1Departamento de
Física Aplicada, Universidad de Alicante, 03690 San Vicente del Raspeig, Spain
2Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4,
20018 Donostia-San Sebastián, Spain
3Theory of Quantum Nanostructures Group, International Iberian Nanotechnology
Laboratory (INL), Av. Mestre José Veiga, 4715-330 Braga, Portugal
4Centro de Física das Universidades do Minho e do Porto, Universidade do
Minho, Campus de Gualtar, 4710-057 Braga, Portugal

###### Abstract

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 atomsHart et al. (2015); Mazurenko et al. (2017),
trapped ionsBlatt and Roos (2012), quantum dotsHensgens et al. (2017);
Mortemousque et al. (2021), dopants in siliconSalfi et al. (2016),
functionalized graphene bilayerGarcía-Martínez and Fernández-Rossier (2019),
moiré heterostructuresKennes et al. (2021) and adatomsYang et al. (2017);
Khajetoorians et al. (2019); Yang et al. (2021). 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 $\pi$ 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 synthesisCai et al. (2010); Ruffieux et al. (2016); Su et al. (2020),
it has been possible to go around this problem, so that the controlled
fabrication of triangulenes of various sizesPavliček et al. (2017); Mishra et
al. (2019); Su et al. (2019); Mishra et al. (2021a), as well as triangulene
dimersMishra et al. (2020), chainsMishra et al. (2021b), ringsMishra et al.
(2021b); Hieulle et al. (2021) and other structuresCheng et al. (2022), has
been recently demonstrated. On the other hand, the on-surface synthesis of
large-area carbon-based crystals has also been establishedSteiner et al.
(2017); Moreno et al. (2018); Telychko et al. (2021). Therefore, the synthesis
of 2D triangulene crystals seems within reach using state-of-the-art
techniques, motivating the present work.

Figure 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.

Individual $[n]$triangulenes have a ground state with spin $S=(n-1)/2$, on
account of their strong intramolecular exchangeOvchinnikov (1978); Fernández-
Rossier and Palacios (2007). The small spin-orbit coupling of carbon makes
their magnetic anisotropy negligibleLado and Fernández-Rossier (2014); Mishra
et al. (2021b). The intermolecular exchange coupling in $[3]$triangulene
dimers has been determined to be antiferromagneticMishra et al. (2020, 2021b);
Jacob et al. (2021); Hieulle et al. (2021). Recently reportedMishra et al.
(2021b) 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 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 graphene-like 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.

##  II Electronic properties of $[n]$triangulenes

###  II.1 Spectrum and symmetries

Figure 2:  (a) Number of states vs energy of $[n]$triangulene, for $n=2,3,4$,
calculated with a nearest-neighbor TB model, featuring $n-1$ zero-energy
modes. (b) Representation of the zero mode wave functions: the circle size
scales with the modulus squared, the arrows represent the phases. For $n=3,4$,
the zero modes are chosen as eigenvectors of the $C_{3}$ symmetry operator
(see Eq. (1) and Table 1). The zero modes vanish at the minority sublattice,
including the corner sites that link triangulenes to their neighbors.

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 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\rangle$ and $|b\rangle$, 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 sublattice-polarized, i.e., they have a non-zero weight
exclusively in the majority sublatticeFernández-Rossier and Palacios (2007);
Ortiz et al. (2019). Because of the degeneracy of the zero modes for $n\geq
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:

$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,$ (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$.

###  II.2 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
DFTFernández-Rossier and Palacios (2007); Wang et al. (2009), quantum
chemistryOrtiz et al. (2019) and model HamiltoniansFernández-Rossier and
Palacios (2007); Güçlü et al. (2009); Ortiz et al. (2019) 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 theoremLieb (1989), 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.

###  II.3 Intermolecular hybridization of zero modes

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 theoremSutherland (1986) 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

Figure 3:  (a-c) Atomic structure of (a) $[2,2]$-, (b) $[3,3]$\- and (c)
$[4,4]$triangulene crystals. Red and blue colors denote the two sublattices.
(d-f) Spin-unpolarized DFT bands of (d) $[2,2]$-, (e) $[3,3]$\- and (f)
$[4,4]$triangulene crystals. (g-i) Energy bands obtained with the TB model,
including third neighbor hopping $t_{3}=0.1t$, for (g) $[2,2]$-, (h) $[3,3]$\-
and (i) $[4,4]$triangulene crystals.

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 EspressoGiannozzi et al. (2009), using the Perdew-Burke-Ernzerhof
functionalPerdew et al. (1996). 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 $\vec{k}$-grids of
$8\times 8\times 1$ for the $[3,3]$triangulene crystal and of $10\times
10\times 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 LABEL:appendix:DFT_details. Our results are in agreement with
previous work by the group of Feng LiuZhou and Liu (2020); Sethi et al. (2021)
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 $\pi$ 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 modelWu et al. (2007); Wu and Sarma (2008). 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 $\Gamma$ 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 nanoribbonsTamaki et al.
(2020).

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 inter-
triangulene 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.

Figure 4:  Scheme illustrating how third neighbor hopping $t_{3}$ leads to
hybridization of zero-energy states in adjacent triangulenes via three pairs
of atomic dimers. The zero mode wave functions at atomic sites labeled with
the same numerals ($j=1,2$) but with different number of primes are related by
$C_{3}$ symmetry.

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$333See Ref.Catarina and Fernández-
Rossier (2022) for a similar derivation applied to a toy model for 1D crystals
of [3]triangulenes.. These zero modes are split into two groups consisting of
$n_{a}-1$ $|a\rangle$ and $n_{b}-1$ $|b\rangle$ states localized in the $A$
and $B$ sublattices of adjacent triangulenes, respectively. Third neighbor
hopping connects $|a\rangle$ and $|b\rangle$ 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:

$\mathcal{H}_{\text{eff}}(\vec{k})=\begin{pmatrix}0^{[a]}&\tau(\vec{k})\\\
\tau^{\dagger}(\vec{k})&0^{[b]}\end{pmatrix},$ (2)

where $\vec{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 $\tau$
is an $(n_{a}-1)\times(n_{b}-1)$ rectangular matrix with entries given by:

$\tau_{ab}(\vec{k})=\sum_{m=0,1,2}e^{i\vec{k}\cdot\vec{R}_{m}}t^{(m)}_{ab}.$
(3)

In the previous expression, $\vec{R}_{0}=(0,0)$ is the null vector associated
to the intracell hopping terms, and $\vec{R}_{1,2}$ are the primitive vectors
that define the hexagonal lattice of $[n_{a},n_{b}]$triangulenes:

$\vec{R}_{1,2}=(n_{a}+n_{b}+1)\vec{a}_{1,2},$ (4)

where $\vec{a}_{1,2}$ are the primitive vectors of the graphene hexagonal
lattice. The ($\vec{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:

$t^{(0)}_{ab}=t_{3}\sum_{j=1,2}a(j)b^{*}(j),$ (5)

where $a(j)=\langle j|a\rangle$ and $b(j)=\langle j|b\rangle$ denote the
components of the zero mode wave functions $|a\rangle$ and $|b\rangle$ at the
atomic sites $j=1,2$, following the labeling depicted in Fig. 4. Importantly,
we can always impose that $a(1)=b(1)\in\mathbb{R}$, from which it follows that
$a(2)=b(2)\in\mathbb{R}$. 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$.

$n$ | zero mode | $\omega_{z}$ | $z(1)$ | $z(2)$
---|---|---|---|---
2 | $z_{0}$ | 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 | $z_{0}$ | 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}}$
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}$ rotation operator (see Eq. (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$.

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 $\omega_{z}$ that each
zero mode wave function acquires when applying the $R_{2\pi/3}$ operator (see
Table 1 for $n=2,3,4$), the remaining values of $z(j^{\prime})$ and
$z(j^{\prime\prime})$ can determined by making use of the $C_{3}$ symmetry.
Specifically, we use that
$z(j)=e^{i\omega_{z}}z(j^{\prime})=e^{-i\omega_{z}}z(j^{\prime\prime})$ to
obtain:

$t^{(1)}_{ab}=t_{3}\sum_{j^{\prime}=1^{\prime},2^{\prime}}a(j^{\prime})b^{*}(j^{\prime})=e^{i(\omega_{b}-\omega_{a})}t^{(0)}_{ab}$
(6)

and

$t^{(2)}_{ab}=t_{3}\sum_{j^{\prime\prime}=1^{\prime\prime},2^{\prime\prime}}a(j^{\prime\prime})b^{*}(j^{\prime\prime})=e^{i(\omega_{a}-\omega_{b})}t^{(0)}_{ab},$
(7)

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 $\omega_{a}-\omega_{b}$. When nonzero, we can interpret the phases
$\omega_{z}$ as Peierls phases associated to a pseudospin-dependent magnetic
field. Importantly, the model must still be time-reversal invariant.

Figure 5:  Energy bands of $[n,n]$triangulene crystals, for (a-c) $n=2$, (d-f)
$n=3$ and (g-i) $n=4$, obtained with (a,d,g) spin-unpolarized DFT, (b,e,g) the
full TB model and (c,f,i) the minimal TB model defined by Eq. (2). The TB
results were obtained with $t_{3}=0.1t$.

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 spin-unpolarized
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).

###  IV.1 $[2,2]$triangulene crystal

The minimal Hamiltonian for the $[2,2]$triangulene crystal is a $2\times 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
$\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:

$\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},$ (8)

with

$f(\phi_{1},\phi_{2})=1+e^{i\phi_{1}}+e^{i\phi_{2}}$ (9)

and

$\phi_{1,2}=\vec{k}\cdot\vec{R}_{1,2}.$ (10)

The corresponding energy bands are given by:

$\displaystyle\epsilon^{[2,2]}_{\pm}(\phi_{1},\phi_{2})$
$\displaystyle=\pm\frac{t_{3}}{3}\left|f(\phi_{1},\phi_{2})\right|$ (11)
$\displaystyle=\pm\frac{t_{3}}{3}\sqrt{3+2\cos\phi_{1}+2\cos\phi_{2}+2\cos(\phi_{1}-\phi_{2})}.$
(12)

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:

$t^{[2,2]}_{\textrm{eff}}=\frac{t_{3}}{3},$ (13)

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 Tran et al.
(2017). A Taylor expansion of Eq. (8) around the $K$ and $K^{\prime}$ points
yields Dirac cones with Fermi velocity $\hbar
v_{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_{3}$, instead of by the nearest neighbor hopping.

###  IV.2 $[2,3]$triangulene crystal

The $[2,3]$triangulene is the smallest non-centrosymmetric 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 1, the effective Bloch
Hamiltonian can be written as:

$\mathcal{H}^{[2,3]}_{\text{eff}}(\phi_{1},\phi_{2})=\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},$
(14)

with

$\mathcal{F}^{[2,3]}_{a_{0},b_{\pm}}(\phi_{1},\phi_{2})=\frac{2t_{3}}{\sqrt{6\times
11}}f(\phi_{1}\pm\theta,\phi_{2}\mp\theta)$ (15)

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

$\epsilon^{[2,3]}_{0}(\phi_{1},\phi_{2})=0$ (16)

and

$\displaystyle\epsilon^{[2,3]}_{\pm}(\phi_{1},\phi_{2})$
$\displaystyle=\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}}$
(17)
$\displaystyle=\pm\frac{2t_{3}}{\sqrt{33}}\sqrt{3-\cos\phi_{1}-\cos\phi_{2}-\cos(\phi_{1}-\phi_{2})}.$
(18)

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 $\Gamma$ 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 Eq. (14) around
$\Gamma$ leads to the following expression for the dispersive bands:

$\epsilon^{[2,3]}_{\pm}(\vec{k})\simeq\pm
6\sqrt{\frac{3}{11}}t_{3}d|\vec{k}|,$ (19)

that accounts for the single Dirac cone. The corresponding Fermi velocity,
$\hbar v_{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 $\Gamma$, the minimal Hamiltonian can be mapped into a spin-1
Dirac modelMizoguchi et al. (2021).

Figure 6:  (a) Atomic structure of the $[2,3]$triangulene crystal. (b) Energy
bands of the $[2,3]$triangulene crystal, obtained with the full TB model,
using $t_{3}=0.1t$. (c) Zoom of panel (b) onto the three lowest-energy bands.
(d) Energy bands of the $[2,3]$triangulene crystal, obtained with the minimal
TB model, taking $t_{3}=0.1t$.

###  IV.3 $[3,3]$triangulene crystal

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 i2\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\times 2$ hopping matrix of the effective Bloch
Hamiltonian reads as:

$\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},$
(20)

where

$\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}),$
(21)
$\mathcal{F}^{[3,3]}_{a_{+},b_{-}}(\phi_{1},\phi_{2})=\frac{2t_{3}}{11}f(\phi_{1}+\theta,\phi_{2}-\theta)$
(22)

and

$\mathcal{F}^{[3,3]}_{a_{-},b_{+}}(\phi_{1},\phi_{2})=\frac{2t_{3}}{11}f(\phi_{1}-\theta,\phi_{2}+\theta).$
(23)

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

$\epsilon_{1,\pm}^{[3,3]}(\phi_{1},\phi_{2})=\pm 3t_{\mathrm{eff}}^{[3,3]},$
(24)

that accounts for the flat bands at finite energy, and

$\epsilon_{2,\pm}^{[3,3]}(\phi_{1},\phi_{2})=\pm
t_{\mathrm{eff}}^{[3,3]}|f(\phi_{1},\phi_{2})|,$ (25)

that accounts for the graphene-like bands, with an effective hopping
$t_{\mathrm{eff}}^{[3,3]}=2t_{3}/11$. The dispersive graphene-like bands
feature Dirac cones at the $K$ and $K^{\prime}$ 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 latticeWu et al. (2007); Wu and Sarma (2008); Tamaki et al. (2020).
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
3t_{\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.

###  IV.4 $[4,4]$triangulene crystal

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 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 $\Gamma$ 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.

###  IV.5 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, $\ell$. For atoms, the
orbital degeneracies are given by $2\ell+1$, with $\ell=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\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 dopingWang et al. (2022).

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

##  V Effect of interactions

###  V.1 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,

$\tilde{U}_{a}=U\sum_{i}|a(i)|^{4}.$ (26)

On the other hand, the average inter-triangulene hopping energy,

$K_{a}=\frac{1}{n_{b}}\sum_{b}\tau_{ab}(0),$ (27)

where $\tau_{ab}(0)$ is the expression given in Eq. (3), evaluated at the
$\Gamma$ point. In order to obtain a single figure of merit for a given unit
cell, we can then average the ratio $\tilde{U}_{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 $\sum_{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
valuesMishra et al. (2020, 2021b) 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\equiv\frac{\tilde{U}}{t^{[2,2]}_{\rm
eff}}=\frac{U/6}{t_{3}/3}=\frac{U}{2t_{3}}$. For $U=1.5t$ and $t_{3}=0.1t$, 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 latticeSorella and Tosatti (1992). 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

$J^{[2]}=4\frac{(t^{[2,2]}_{\rm
eff})^{2}}{\tilde{U}}=6\times\frac{4}{9}\frac{t_{3}^{2}}{U}\simeq 48{\rm
meV},$ (28)

for $t=2.7eV$, $U=1.5t$ and $t_{3}=0.1t$. As we shall discuss in a forthcoming
workJacob et al. , in addition to the kinetic exchange, there are other
contributions to the antiferromagnetic exchange in this system.

###  V.2 Spin-polarized DFT calculations for the $[3,3]$triangulene crystal

Figure 7:  a) Energy bands, computed with spin-polarized DFT, for the
$[3,3]$triangulene crystal. b) Computed magnetic moment density, showing
antiferromagnetic order.

We now discuss the results of our spin-polarized DFT calculations on the
$[3,3]$triangulene crystal (see details in Appendix
LABEL:appendix:DFT_details). 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 crystalsMishra et al. (2021b). 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 LiuZhou and Liu (2020), 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_{FM}-E_{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
approximationMishra et al. (2020).

###  V.3 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 modelsGong et al. (2013). In the case of $[2,3]$ at half
filling, Lieb theoremLieb (1989) 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
clustersGong et al. (2015). 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 couplingsAffleck et al. (1987). 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
computingWei et al. (2011) (MBQC). Importantly, the AKLT model has a gap in
the excitation spectrumGarcia-Saez et al. (2013); Pomata and Wei (2020); Lemm
et al. (2020), 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 bicubicAffleck et al. (1987) 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 modelGanesh et al. (2011) 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 chainsMishra et al. (2021b).

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 grapheneCao
et al. (2018). 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. 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. 2.

The zero modes vanish identically at the inter-triangulene binding sites. As a
result intermolecular hybridization is governed by third neighbor hopping
$t_{3}$.

  3. 3.

Centrosymmetryc $[n,n]$triangulene crystals have $2(n-1)$ narrow bands. For
$n\geq 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. 4.

Non centrosymmetric $[n_{a},n_{b}]$triangulene crystals, with $n_{a}\neq
n_{b}$, feature $|n_{a}-n_{b}|$ flat bands at $E=0$, because of the bipartite
character of the latticeSutherland (1986); Ortiz et al. (2019)

The results above are validated comparing spin-unpolarized DFT calculations
with tight-binding Hamiltonians at two levels: a full-lattice model that
includes first and third neighbor hopping and a minimal tight-binding 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 half-filling these
crystals are expected to be Mott insulators, with their low energy physics
governed by the spin degrees of freedom. Application of Lieb theoremLieb
(1989) 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}\neq
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 structuresMishra et al. (2021b) 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.

Acknowledgements We acknowledge financial support from the Ministry of Science
and Innovation of Spain (grant No. PID2019-109539GB-41), from Generalitat
Valenciana (grant No. Prometeo2021/017) and from Fundacao Para a Ciencia e a
Tecnologia, Portugal (grant No. PTDC/FIS-MAC/2045/2021). R. Ortiz acknowledges
funding fromGeneralitat Valenciana and Fondo Social Europeo (Grant No.
ACIF/2018/175). G. Catarina acknowledges financial support from Fundação para
a Ciência e a Tecnologia (Grant No. SFRH/BD/138806/2018).

## References

  * Hart et al. (2015) R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse, and R. G. Hulet, Nature 519, 211 (2015), ISSN 1476-4687.
  * Mazurenko et al. (2017) A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017), ISSN 1476-4687.
  * Blatt and Roos (2012) R. Blatt and C. F. Roos, Nature Physics 8, 277 (2012).
  * Hensgens et al. (2017) T. Hensgens, T. Fujita, L. Janssen, X. Li, C. J. Van Diepen, C. Reichl, W. Wegscheider, S. Das Sarma, and L. M. K. Vandersypen, Nature 548, 70 (2017), ISSN 1476-4687.
  * Mortemousque et al. (2021) P.-A. Mortemousque, E. Chanrion, B. Jadot, H. Flentje, A. Ludwig, A. D. Wieck, M. Urdampilleta, C. Bäuerle, and T. Meunier, Nature Nanotechnology 16, 296 (2021), ISSN 1748-3395.
  * Salfi et al. (2016) J. Salfi, J. A. Mol, R. Rahman, G. Klimeck, M. Y. Simmons, L. C. L. Hollenberg, and S. Rogge, Nature Communications 7, 11342 (2016), ISSN 2041-1723.
  * García-Martínez and Fernández-Rossier (2019) N. A. García-Martínez and J. Fernández-Rossier, Phys. Rev. Research 1, 033173 (2019).
  * Kennes et al. (2021) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021).
  * Yang et al. (2017) K. Yang, Y. Bae, W. Paul, F. D. Natterer, P. Willke, J. L. Lado, A. Ferrón, T. Choi, J. Fernández-Rossier, A. J. Heinrich, et al., Phys. Rev. Lett. 119, 227206 (2017).
  * Khajetoorians et al. (2019) A. A. Khajetoorians, D. Wegner, A. F. Otte, and I. Swart, Nature Reviews Physics 1, 703 (2019).
  * Yang et al. (2021) K. Yang, S.-H. Phark, Y. Bae, T. Esat, P. Willke, A. Ardavan, A. J. Heinrich, and C. P. Lutz, Nature Communications 12, 993 (2021), ISSN 2041-1723.
  * Cai et al. (2010) J. Cai, P. Ruffieux, R. Jaafar, M. Bieri, T. Braun, S. Blankenburg, M. Muoth, A. P. Seitsonen, M. Saleh, X. Feng, et al., Nature 466, 470 (2010).
  * Ruffieux et al. (2016) P. Ruffieux, S. Wang, B. Yang, C. Sánchez-Sánchez, J. Liu, T. Dienel, L. Talirz, P. Shinde, C. A. Pignedoli, D. Passerone, et al., Nature 531, 489 (2016).
  * Su et al. (2020) J. Su, M. Telychko, S. Song, and J. Lu, Angewandte Chemie 132, 7730 (2020).
  * Pavliček et al. (2017) N. Pavliček, A. Mistry, Z. Majzik, N. Moll, G. Meyer, D. J. Fox, and L. Gross, Nature Nanotechnology 12, 308 (2017).
  * Mishra et al. (2019) S. Mishra, D. Beyer, K. Eimre, J. Liu, R. Berger, O. Groning, C. A. Pignedoli, K. Müllen, R. Fasel, X. Feng, et al., Journal of the American Chemical Society 141, 10621 (2019).
  * Su et al. (2019) J. Su, M. Telychko, P. Hu, G. Macam, P. Mutombo, H. Zhang, Y. Bao, F. Cheng, Z.-Q. Huang, Z. Qiu, et al., Science advances 5, eaav7717 (2019).
  * Mishra et al. (2021a) S. Mishra, K. Xu, K. Eimre, H. Komber, J. Ma, C. A. Pignedoli, R. Fasel, X. Feng, and P. Ruffieux, Nanoscale 13, 1624 (2021a).
  * Mishra et al. (2020) S. Mishra, D. Beyer, K. Eimre, R. Ortiz, J. Fernández-Rossier, R. Berger, O. Gröning, C. A. Pignedoli, R. Fasel, X. Feng, et al., Angewandte Chemie International Edition (2020).
  * Mishra et al. (2021b) S. Mishra, G. Catarina, F. Wu, R. Ortiz, D. Jacob, K. Eimre, J. Ma, C. A. Pignedoli, X. Feng, P. Ruffieux, et al., Nature 598, 287 (2021b), ISSN 1476-4687.
  * Hieulle et al. (2021) J. Hieulle, S. Castro, N. Friedrich, A. Vegliante, F. R. Lara, S. Sanz, D. Rey, M. Corso, T. Frederiksen, J. I. Pascual, et al., Angewandte Chemie International Edition 60, 25224 (2021).
  * Cheng et al. (2022) S. Cheng, Z. Xue, C. Li, Y. Liu, L. Xiang, Y. Ke, K. Yan, S. Wang, and P. Yu, Nature Communications 13, 1705 (2022), ISSN 2041-1723.
  * Steiner et al. (2017) C. Steiner, J. Gebhardt, M. Ammon, Z. Yang, A. Heidenreich, N. Hammer, A. Görling, M. Kivala, and S. Maier, Nature communications 8, 1 (2017).
  * Moreno et al. (2018) C. Moreno, M. Vilas-Varela, B. Kretz, A. Garcia-Lekue, M. V. Costache, M. Paradinas, M. Panighel, G. Ceballos, S. O. Valenzuela, D. Peña, et al., Science 360, 199 (2018).
  * Telychko et al. (2021) M. Telychko, G. Li, P. Mutombo, D. Soler-Polo, X. Peng, J. Su, S. Song, M. J. Koh, M. Edmonds, P. Jelínek, et al., Science advances 7, eabf0269 (2021).
  * Ovchinnikov (1978) A. A. Ovchinnikov, Theoretica Chimica Acta 47, 297 (1978).
  * Fernández-Rossier and Palacios (2007) J. Fernández-Rossier and J. J. Palacios, Physical Review Letters 99, 177204 (2007).
  * Lado and Fernández-Rossier (2014) J. L. Lado and J. Fernández-Rossier, PRL 113, 027203 (2014).
  * Jacob et al. (2021) D. Jacob, R. Ortiz, and J. Fernández-Rossier, Phys. Rev. B 104, 075404 (2021).
  * Ortiz et al. (2019) R. Ortiz, R. Á. Boto, N. García-Martínez, J. C. Sancho-García, M. Melle-Franco, and J. Fernández-Rossier, Nano Lett. 19, 5991 (2019).
  * Wang et al. (2009) W. L. Wang, O. V. Yazyev, S. Meng, and E. Kaxiras, Physical review letters 102, 157201 (2009).
  * Güçlü et al. (2009) A. D. Güçlü, P. Potasz, O. Voznyy, M. Korkusinski, and P. Hawrylak, Phys. Rev. Lett. 103, 246805 (2009).
  * Lieb (1989) E. H. Lieb, Physical Review Letters 62, 1201 (1989).
  * Sutherland (1986) B. Sutherland, Phys. Rev. B 34, 5208 (1986).
  * Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Zhou and Liu (2020) Y. Zhou and F. Liu, Nano Letters 21, 230 (2020).
  * Sethi et al. (2021) G. Sethi, Y. Zhou, L. Zhu, L. Yang, and F. Liu, Physical Review Letters 126, 196403 (2021).
  * Wu et al. (2007) C. Wu, D. Bergman, L. Balents, and S. D. Sarma, Physical review letters 99, 070401 (2007).
  * Wu and Sarma (2008) C. Wu and S. D. Sarma, Physical Review B 77, 235107 (2008).
  * Tamaki et al. (2020) G. Tamaki, T. Kawakami, and M. Koshino, Physical Review B 101, 205311 (2020).
  * Tran et al. (2017) V.-T. Tran, J. Saint-Martin, P. Dollfus, and S. Volz, AIP Advances 7, 075212 (2017).
  * Mizoguchi et al. (2021) T. Mizoguchi, Y. Kuno, and Y. Hatsugai, Phys. Rev. B 104, 035161 (2021).
  * Wang et al. (2022) T. Wang, A. Berdonces-Layunta, N. Friedrich, M. Vilas-Varela, J. P. Calupitan, J. I. Pascual, D. Peña, D. Casanova, M. Corso, and D. G. de Oteyza, Journal of the American Chemical Society 144, 4522 (2022), pMID: 35254059, eprint https://doi.org/10.1021/jacs.1c12618, URL https://doi.org/10.1021/jacs.1c12618.
  * Sorella and Tosatti (1992) S. Sorella and E. Tosatti, EPL (Europhysics Letters) 19, 699 (1992).
  * (46) D. Jacob, R. Ortiz, G. Catarina, and J. Fernández-Rossier, _Theory of intermolecular exchange in nanographenes_.
  * Gong et al. (2013) S.-S. Gong, D. Sheng, O. I. Motrunich, and M. P. Fisher, Physical Review B 88, 165138 (2013).
  * Gong et al. (2015) S.-S. Gong, W. Zhu, and D. Sheng, Physical Review B 92, 195110 (2015).
  * Affleck et al. (1987) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
  * Wei et al. (2011) T.-C. Wei, I. Affleck, and R. Raussendorf, Phys. Rev. Lett. 106, 070501 (2011).
  * Garcia-Saez et al. (2013) A. Garcia-Saez, V. Murg, and T.-C. Wei, Phys. Rev. B 88, 245118 (2013).
  * Pomata and Wei (2020) N. Pomata and T.-C. Wei, Phys. Rev. Lett. 124, 177203 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.177203.
  * Lemm et al. (2020) M. Lemm, A. W. Sandvik, and L. Wang, Phys. Rev. Lett. 124, 177204 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.177204.
  * Ganesh et al. (2011) R. Ganesh, D. Sheng, Y.-J. Kim, and A. Paramekanti, Physical Review B 83, 144414 (2011).
  * Cao et al. (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Catarina and Fernández-Rossier (2022) G. Catarina and J. Fernández-Rossier, Physical Review B 105, L081116 (2022).