# Numerical Differentiation using local Chebyshev-Approximation

Stefan H. Reiterer

##  1 Motivation

In applied mathematics, especially in optimization, functions are often only
provided as so called ”Black-Boxes” provided by software packages, or very
complex algorithms, which make automatic differentation very complicated or
even impossible. Hence one seeks the numerical approximation of the
derivative.

Unfortunately numerical differentation is a difficult task in itself, and it
is well known that it is numerical instable. There are many works on this
topic, including the usage of (global) Chebyshev approximations. Chebyshev
approximations have the great property that they converge very fast, if the
function is smooth. Nevertheless those approches have several drawbacks, since
in practice functions are not smooth, and a global approximation needs many
function evalutions.

Nevertheless there is hope. Since functions in real world applications are
most times smooth except for finite points, corners or edges. This motivates
to use a local Chebyshev approach, where the function is only approximated
locally, and hence the Chebyshev approximations still yields a fast
approximation of the desired function. We will study such an approch in this
work, and will provide a numerical example.

Disclaimer: This work was done as private research during my study years and
is not related to my current affiliation at all. No funding was recieved for
this work.

##  2 Chebyshev Polynomials and One Dimensional Functions

First we recall some well known definitions. the following definitions and
properties are taken from [1].

######  Definition 2.1 (Chebyshev-Polynomials).

For $n\in\mathbb{N}_{0}$ the $n$-th Chebyshev polynomial
$T_{n}(x):\,[-1,1]\rightarrow\mathbb{R}$ is given by
$T_{n}(x)=\cos(n\arccos(x))$. From this definition we see that the Chebyshev
polynomial takes it’s extrema at

$x_{k}:=-\cos\left(\frac{k\pi}{n}\right)\text{ for }k=0,\ldots,n.$

The points $(x_{k})_{k=0}^{n}$ are the so called Gauss-Lobatto Grid-points.

The best way to approximate functions via Chebyshev-polynomials lies in the
following

######  Proposition 2.2 (Chebyshev-Series & Co).

Set $\omega(x):=\frac{2}{\pi}(1-x^{2})^{-1/2}$, then

$\left<f,g\right>_{\omega}:=\int\limits_{-1}^{1}f(x)g(x)\omega(x)\,dx$

forms the (weighted) scalar product of the Hilbert-Space

$L_{2,\omega}:=\left\\{f:\,[-1,1]\rightarrow\mathbb{R}\,|\left\|f\right\|_{L_{2,\omega}}=\sqrt{\left<f,f\right>_{\omega}}<\infty\right\\}.$

Then for a function $f\in L_{2,\omega}$ the Chebyshev-Coefficients

$a_{n}:=\left<f,T_{n}\right>_{\omega},$

exists, and the Chebyshev-Series

$\frac{a_{0}}{2}+\sum\limits_{n=1}^{\infty}a_{n}T_{n}(x),$

converges in the $L_{2,\omega}$ sense to $f$.

###### Proof.

These well known facts follow from applying classical Fourier-theory to

$\tilde{f}=f(\cos(t))\in L_{2}(0,\pi)$. ∎

Since Chebyshev-Polynomials are “cosines in disguise” other important
properties from Fourier-theory carry over to Chebyshev-series, like spectral
convergence for smooth functions etc.

Another way to approximate a function is to use the Chebyshev-Interpolation
Polynomial $C_{n}=\sum\limits_{k=0}^{n}=b_{k}T_{k}$, which is uniquely defined
by $f(x_{k})=C_{n}(x_{k})$ for $k=0,\ldots,n$.

Although the Chebyshev-Series yields the best approximation there is a little
known fact between the two approximation types, namely the

######  Proposition 2.3.

Let $f:\,[-1,1]\rightarrow\mathbb{R}$ a function with

$f(x)=\frac{a_{0}}{2}+\sum\limits_{n=1}^{\infty}a_{n}T_{n}(x),$

and

$\frac{\left|a_{0}\right|}{2}+\sum\limits_{n=1}^{\infty}\left|a_{n}\right|<\infty.$

Further let for $f_{N}(x):=\frac{a_{0}}{2}+\sum_{n=1}^{N}a_{n}T_{n}(x)$,

$E_{T}(N):=\sup_{x\in[-1,1]}\left|f(x)-f_{N}(x)\right|,$

be the trunctation error, and Then

$E_{T}(N)\leq\sum\limits_{n=N}^{\infty}\left|a_{n}\right|,$

and for the interpolation error

$\sup_{x\in[-1,1]}\left|f(x)-C_{N}(x)\right|\leq 2E_{T}(N).$

Hence the penalty for using interpolation instead of truncation is at most a
factor of two! Additionally the coefficients $b_{k}$ of the interpolation
polynomial $C_{n}$ are related to the exact coefficents $a_{k}$ by the
identity

$b_{k}=a_{k}+\sum\limits_{j=1}^{\infty}a_{k+2jn}+a_{-k+2jn}.$ (1)

That means the approximated coefficients differ from the exact coefficients by
the aliasing-error, and hence the error vanishes with $\mathcal{O}(a_{n})$

###### Proof.

See [1, Thm. 6& Thm. 21]. ∎

Since with Chebyshev interpolation we are close to the realm of the DFT,
additional methods from signal processing, like denoising could be considered
to handle numerical distortions.

To work with a function $g:\,[a,b]\rightarrow\mathbb{R}$ a generalized
Chebyshev interpolation for $g$ can be achieved by using a linear
transformation $\varphi:\,[a,b]\rightarrow[-1,1]$, and interpolate the
function $f=g\circ\varphi$.

Since computing the derivative of a polynomial can be done exactly we could
use the derivative of $C_{N}$ for some $N$, as the numerical derivative of
$f$. This means

$f^{\prime}\approx C_{N}^{\prime}.$

Then how about the errors? First recall the identity
$T^{\prime}_{n}(x)=nU_{n-1}$. Hence for the truncation error we have

$f^{\prime}(x)-f^{\prime}_{N}(x)=E_{N}^{\prime}(x)=\sum\limits_{n=N}^{\infty}na_{n}U_{n-1}(x).$

Considering that $U_{n}$ has it’s extrema at $\pm 1$ with $U_{n-1}(\pm 1)=(\pm
1)^{n}n$, we immediately see that

$\sup_{x\in[-1,1]}\left|f^{\prime}(x)-f^{\prime}_{N}(x)\right|\leq\sum\limits_{n=N}^{\infty}n\left|a_{n}\right|\sup_{x\in[-1,1]}\left|U_{n-1}(x)\right|\leq\sum\limits_{n=N}^{\infty}n^{2}\left|a_{n}\right|=\mathcal{O}(n^{2}a_{n}).$

For the interpolation polynomial $C_{N}$ we get by using the aliasing identity
(1) and rearranging terms like in the proof of [1, Thm. 21] the error
estimation

$\left|f^{\prime}(x)-C^{\prime}_{N}(x)\right|\leq
2\sum\limits_{n=N}^{\infty}n^{2}\left|a_{n}\right|=\mathcal{O}(2n^{2}a_{n}).$

Hence we have the

######  Proposition 2.4.

For $f:\,[-1,1]\rightarrow\mathbb{R}$ smooth we have the error estimation

$\sup_{x\in[-1,1]}\left|f^{\prime}(x)-C^{\prime}_{N}(x)\right|\leq
2E_{T}^{\prime}(N)\mathcal{O}(2n^{2}a_{n}).$

Hence the penalty for using interpolation instead of truncation when
differentiating is at most a factor of two!

##  3 Practical Considerations in the 1D Case and Generalizations

In this section consider a function $f:\,[a,b]\rightarrow\mathbb{R}$, which is
continuous and piecewise smooth, but not differentiable in a set of finite
points $(y_{k})_{k\in I}\subset[a,b]$ (with $I$ finite and allowed to be
empty).

If we interpolate now the function $f$ on the interval $[a,b]$, we know from
Fourier-Theory, that the error will converge only slowly to zero. We also note
that it is easy to see that (from variable transformations) that error
estimation depends on the length of the interval $[a,b]$. with the factor
$(\frac{b-a}{2})^{N}$. Hence if one restricts the function $f$ only on a small
interval, one achieves more precise interpolation results, with lesser degrees
of the interpolation polynomial. If we restrict the function $f$ to some
interval $[c,d]\subset[a,b]$, with

  1. 1.

$d-c<1$ and

  2. 2.

$y_{i}\not\in(c,d)$ for $i\in I$,

then the sequence of the Chebyshev interpolation polynomials
$(C_{N})_{N\in\mathbb{N}}$ converges very fast (with spectral convergence) and
uniformly to the original function in the interval $[c,d]$.

From Proposition 2.4 we know that this then also applies to $C_{N}^{\prime}$,
since the error decays with $\mathcal{O}(n^{2}a_{n})$. Nevertheless we also
see that for practical use it is important to ensure the smoothness on the
observed interval $[c,d]$, because of the term $n^{2}$. This leads to the
following definition

######  Definition 3.1.

Let $f:\,[a,b]\rightarrow\mathbb{R}$ be continuous and piecewise smooth, then
we call the Chebyshev interpolation polynomial
$C_{N,[c,d]}:\,[c,d]\rightarrow\mathbb{R}$ of $f_{|[c,d]}$ the local smooth
Chebyshev polynomial iff

$(y_{k})_{k\in I}\cap(c,d)=\emptyset.$

This motivates the following algortithm to compute derivatives of a function
$f$ (if the points $Y:=(y_{k})_{k\in I}$ are known):

⬇

# input: function f, and point x, estimated lenght h

set of disallowed points Y, number of interpolation points N

def compute_derivative(f,x,h,Y):

if x not in Y:

# classical derivative

c = x-h; d= x+h

# ensure that we are smooth on interval [c,d]

while c < a and b < d and intersection(Y,[c,d]) != {}:

h = make_h_smaller(h)

c = x-h; d= x+h

CN = compute_interpolation(f,c,d,N)

CN_prime = derive_polynomial(CN)

return CN_prime(x)

else: # x in Y

# subgradient

c[0] = x-h; d[0] = x

c[1] = x; d[1] = x+h

while all([c[i] < a and b < d[i] and \

intersection(Y,(c[i],d[i])) != {} for i in [0,1]]):

h = make_h_smaller(h)

c[0] = x; d[0] = x+h

c[1] = x-h; d[1] = x

for i in [0,1]:

CN[i] = compute_interpolation(f,c[i],d[i],N)

CN_prime[i] = derive_polynomial(CN[i])

return [CN_prime[0](x),CN_prime[1](x)]

As one can see we also included more general derivatives based on directional
derivatives like subgradients. It also would be possible to define a weak
derivative witch is motivated by Chebyshev convergence theory defined by
$f^{\prime}(x):=\frac{1}{2}(f^{\prime}(-x)+f^{\prime}(+x))$.

A question which still remains, is what to do when the set $Y$ is unknown.
This is still an open problem, but since one can observe the speed of the
decay of the coefficients $a_{n}$ (or $b_{n}$) it is possible to locate non
differentiable points. There are several works on this topic from spectral
theory.

Now to the $M$ dimensional case: Consider a functional
$f:\Omega\subseteq\mathbb{R}^{M}\rightarrow\mathbb{R}$, which is continuous
and piecewise smooth on the domain $\Omega$. We can compute the directional
derivatives in a direction $h\in\mathbb{R}^{M}$ in a point $x\in\Omega$ by
deriving the one dimensional function $g(t):=f(x+th)$. Hence the methods from
the one dimensional case, can be carried over to the $M$ dimensional setting.

##  4 Examples

as first we define the function $f_{1}$ given by (see also Figure 1):

$f_{1}(x)=\begin{cases}x^{4}\text{ if }x>0,\\\ 0\text{ else.}\end{cases}$

First we compare computing the first derivative with local Chebyshev
approximation with finite differences given by

$f^{\prime}(x)\approx\frac{f(x+h)-f(x-h)}{2h}=:f_{h}(x),$

by computing the values at $x=0$ and $x=0.5$. See Table 1 for results. As we
can observe the errors of $f_{h}$ and $C_{3,[x\pm h]}$ are equal. This comes
without surprise, since the derivative of the $3$rd Chebyshev polynomial is
exactly the central difference quotient.

$x$ | $h$ | $f_{h}$ | $C_{3,[x\pm h]}$ | $C_{5,[x\pm h]}$
---|---|---|---|---
$0.5$ | $1e-3$ | $2e-6$ | $2e-6$ | $2e-14$
| $1e-4$ | $1.9e-8$ | $1.9e-8$ | $2.6e-13$
| $1e-5$ | $1.9e-10$ | $1.9e-10$ | $1.8e-12$
$0.0$ | $1e-3$ | $5e-10$ | $5e-10$ | $1.4e-10$
| $1e-4$ | $5e-13$ | $5e-13$ | $2.6e-13$
| $1e-5$ | $5e-16$ | $5e-16$ | $1.5e-16$
Table 1: Errors for $f_{1}$

Next we will modify $f_{1}$ in the following way: For some $\varepsilon>0$,
and a randomly standard normal distributed variable $X$ we define $f_{2}$ by

$f_{2}(x)=f_{1}(x)+\varepsilon X$

and then try to differentiate it at $0.5$. the results can be seen at Table 2.

$x$ | $h$ | $f_{h}$ | $C_{3,[x\pm h]}$ | $C_{5,[x\pm h]}$ | $C_{7,[x\pm h]}$
---|---|---|---|---|---
$0.5$ | $1e-1$ | $2e-2$ | $2e-2$ | $5.5e-9$ | $4.2e-8$
| $1e-2$ | $2.4e-6$ | $2e-4$ | $4e-7$ | $8e-7$
| $1e-3$ | $2.4e-6$ | $2e-6$ | $1.8e-6$ | $7.3e-6$
| $1e-4$ | $1e-6$ | $4e-6$ | $2e-5$ | $2e-5$
| $1e-5$ | $2.8e-5$ | $9e-5$ | $2e-5$ | $7e-4$
| $1e-6$ | $7e-5$ | $1e-3$ | $4e-5$ | $2e-3$
Table 2: Errors for $f_{2}$

The errors are maximal values after taking some samples (There is some space
for improvement here, e.g. better statistics, but we get the picture…) We
Observe the following behavior: While $f_{h}$ and $C_{3,[x\pm h]}$ behave like
expected, we see that making $h$ smaller, does not make the approximation
quality of $C_{5,[x\pm h]}$ and $C_{7,[x\pm h]}$ better, but even worse. This
can be explained by the fact, that for smooth $f$, we know that the function
behaves like it’s Taylor polynomial of lower order, while the higher order
terms are neglectable. This also means that computing a higher order Chebyshev
polynomial does not make sense, and we even get into more trouble. hence when
choosing a suitable $h$ for the Chebyshev approximation should be not too
small, and yet not too big.

Figure 1: Function $f_{1}$

Finally we try the Local Chebyshev Method on the 2D Rosenbrock function

$R_{a,b}(x)=(a-x_{0})^{2}+b(x_{1}-x_{0}^{2})^{2},$

with parameters $a=1$ and $b=100$, and try to find the minimum with help of
gradient methods. The argmin is $x_{0}=(1,1)$ Additionally we add some
disturbances $\delta(x)$ and $\varepsilon(x)$ in the interval $[-1e-6,1e-6]$.
While $\delta$ is some randomly normal distributed variable, $\varepsilon$ is
a jump function which changes it’s sign (which is fatal for finite
differences). As optimization method we used Steepest Descent with Armijo-rule
and $\left\|\nabla f\right\|<1e-3$ as terminating criterion. We used $h=1e-6$
for the finite difference method, and order $5$ polynomials with $h=1e-4$ for
the Local Chebyshev Method. and The number of iterations and the computed
results (rounded up to 3 digits) can be seen in Table 3. The results suggest,
that computing the gradient with the Local Chebyshev Method makes the
algorithm more robust. Also it was observed that the stepsize $h$ should be
not too small.

Function | Method | Iteration Numbers | Result
---|---|---|---
$R_{a,b}$ | Exact Gradient | $1427$ | $(0.999,0.998)$
| Finite Differences | $1427$ | $(0.999,0.998)$
| Local-Chebyshev | $1428$ | $(0.999,0.998)$
$R_{a,b}+\delta$ | Finite Differences |  $19999$ (interrupted) | $(0.986,0.973)$
| Local-Chebyshev | $1090$ | $(0.996,0.992)$
$R_{a,b}+\varepsilon$ | Finite Differences |  $19999$ (interrupted) | $(0.486,0.234)$
| Local-Chebyshev | $1565$ | $(0.9990.998)$
Table 3: Iteration numbers for optimization

## References

  * [1] Boyd, John P. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.


# Supplementary materials
Superexchange Mechanism in Coupled Triangulenes Forming spin-1 Chains

Yasser Saleem  Institut für Physikalische Chemie, Universität Hamburg,
Grindelallee 117, D-20146 Hamburg, Germany  Torben Steenbock  Institut für
Physikalische Chemie, Universität Hamburg, Grindelallee 117, D-20146 Hamburg,
Germany  Emha Riyadhul Jinan Alhadi  Institute of Physics, Faculty of Physics,
Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5,
87-100 Toruń, Poland  Weronika Pasek  Institute of Physics, Faculty of
Physics, Astronomy and Informatics, Nicolaus Copernicus University,
Grudziadzka 5, 87-100 Toruń, Poland  Gabriel Bester  Institut für
Physikalische Chemie, Universität Hamburg, Grindelallee 117, D-20146 Hamburg,
Germany  Pawel Potasz  Institute of Physics, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Toruń,
Poland

(12 March, 2024)

††preprint: APS/123-QED

##  I Density functional theory

All density functional theory (DFT) calculations were carried out with the
Turbomole 7.5 program package TURBOMOLE  employing the B3LYP global hybrid
exchange–correlation functional with 20% of the exact Hartree-Fock exchange
Dirac1929 ; Slater1951 ; Vosko1980 ; Becke1988 ; Lee1988 ; Becke1992 ,
Ahlrich’s triple-zeta split-valence basis set with polarization functions on
all atoms, def2-TZVP Weigend2006 , and the empirical dispersion correction of
Grimme in the third generation Grimme2010 . Additionally, we employ the
multipole-accelerated resolution of identity approximation for Coulomb
integrals (MARIJ) Eichkorn1995 ; Eichkorn1997 ; Sierka2003 ; Weigend2006  to
speed-up the calculation of the Coulomb integrals in the self-consistent field
(SCF) algorithm.

First, we perform structure optimizations on the spin ground states
approximated by Broken-Symmetry (BS) determinants Noodleman1981 , where the
$S=1$ spins of adjacent triangulenes are aligned anti-parallel. In the
structure optimizations convergence criteria of $10^{-7}$ Hartree for the
energy and $10^{-4}$ Hartree/Bohr for the gradient were employed. On top of
the optimized structures, we perform single-point calculations for non-spin
polarized cations that serve as the reference in the full configuration
interaction (CI) calculations. The cations are obtained by removing all
unpaired electrons, two for each triangulene subunit, from the systems. In
these single-point calculations,the same convergence criterion for the energy
of $10^{-7}$ Hartree is employed.

##  II Calculation of Coulomb Matrix Elements in the SITE BASIS

The Coulomb matrix elements $\braket{ij}{V}{kl}$ are given explicitly as
graphenebook

$\braket{ij}{V}{kl}=\int\int
d\vec{r}_{1}d\vec{r}_{2}\phi^{*}_{i}(\vec{r}_{1})\phi^{*}_{j}(\vec{r}_{2})\frac{2R_{y}}{\kappa|\vec{r}_{1}-\vec{r}_{2}|}\phi_{k}(\vec{r}_{2})\phi_{l}(\vec{r}_{1}),$
(1)

where, $\vec{r}_{1}$, $\vec{r}_{2}$ are the coordinates of electron 1 and
electron 2, $\kappa$ is the dielectric constant, and $R_{y}$ is the Rydberg
constant. $\phi_{i}(\vec{r})$ are $p_{z}$ Slater orbitals centered on atom
$i$, they are given as

$\phi_{i}(\vec{r})=\left(\frac{\xi^{5}}{32\pi}\right)^{\frac{1}{2}}ze^{-\frac{\xi}{2}\left|\vec{r}-\vec{r}_{i}\right|},$
(2)

where $\vec{r}_{i}$ is the position of carbon atom $i$, and $\xi=3.25$ xipaper
. The lengths are in units of Bohrs. The integrals are solved efficiently in
real space using the VEGAS integration alogorithim contained in the GNU
scientific library gsl . Table 1 shows all quantum mechanical scattering
elements considered and their magnitude for $\kappa=1$. When electrons are far
enough away, we take a classical limit of the Coulomb elements, mainly we also
take the elements:
$\braket{ij}{V}{ji}=\frac{2R_{y}}{|\vec{R}_{i}-\vec{R}_{j}|}$ for atoms
$i$,$j$ beyond NNN.

Element | Value (eV)
---|---
$\langle 1,1|V|1,1\rangle$ | $17.311$
$\langle 1,2|V|2,1\rangle$ | $8.942$
$\langle 1,3|V|3,1\rangle$ | $5.583$
$\langle 1,1|V|1,2\rangle$ | $3.027$
$\langle 1,2|V|3,1\rangle$ | $1.664$
$\langle 1,2|V|1,2\rangle$ | $0.774$
$\langle 1,1|V|2,2\rangle$ | $0.774$
$\langle 2,2|V|1,3\rangle$ | $0.535$
$\langle 1,2|V|2,3\rangle$ | $0.356$
$\langle 1,1|V|1,3\rangle$ | $0.306$
$\langle 1,1|V|1,4\rangle$ | $0.113$
$\langle 1,2|V|2,4\rangle$ | $0.113$
$\langle 1,3|V|3,4\rangle$ | $0.122$
Table 1: Table of scattering Coulomb matrix elements of $p_{z}$ electrons for
$\kappa=1$. The numbers 1,2,3, and 4 correspond to atoms where the difference
in the numbers corresponds to the neighbor distance i.e. a difference of 1
corresponds to NN, a difference of 2 between numbers corresponds to NNN, and a
difference of 3 corresponds to NNNN.

##  III Configuration Interaction basis

The many-body Hilbert space can be divided into smaller subspaces with total
spin ${\bf S}$ and azimuthal spin $S^{z}$. We construct the basis, in the
occupation number representation, distributing particles among single-particle
states labeled with spin $S^{z}$. The total number of possible configurations
$N_{\rm cf}$ for $N_{\rm e}$ particles distributed on $N_{\rm st}$ single
particle states with a given spin $N_{\downarrow}$ or $N_{\uparrow}$, where
$N_{\rm e}=N_{\downarrow}+N_{\uparrow}$ is determined by a product of binomial
coefficients, $N_{\rm cf}=\binom{N_{\rm st}}{N_{\uparrow}}\cdot\binom{N_{\rm
st}}{N_{\uparrow}}$. We do not rotate the Hamiltonian matrix to a ${\bf S}$
basis as this is an additional computational cost, and instead determine the
ground state from calculations of expectation value of total spin $S$ for each
energy eigenstate.

##  IV Triangulene

We analyze a single triangulene molecule, a triangular graphene quantum dot
with zigzag edges shown in Fig. 1(a) of the main text. This quantum dot has
broken sublattice symmetry as seen by counting the number of red balls (carbon
atoms belonging to sublattice A) and the number of gray balls (carbon atoms
belonging to sublattice B). This is also a bipartite lattice, and as such,
Lieb’s theorem appliesLiebTheorems .

We start by performing DFT calculations, and find the ground state to be $S=1$
in agreement with Lieb’s theorem and previous experimental and theoretical
work pavlivcek2017synthesis ; mishra2019synthesis ; ortiz2022theory . We then
perform HF calculations by solving the HF equation. We take $\kappa=3$ giving
the results shown in Fig. S1. At the top of the VB there are two degenerate
states that are spin up, and a large gap that separates the spin down states.
The splitting between the spin up and spin down states arises from a net spin-
polarization. These states are localized at the edge, and the spin density of
triangulene (shown as an inset in Fig. S1) shows the localized spin-1 quasi
particle tends to localize at the edge of the triangle.

Figure S1:  Single-particle spectrum of triangulene $N_{\rm Tr}=1$ at half-
filling. The horizontal bars are HF results, while the crosses correspond to
DFT results. Red color corresponds to spin up levels while blue color
corresponds to spin down levels. The inset shows the spin density of the HF
ground state.

##  V Comparison between different many-body interacting fermionic models -
results for two triangles

We compare results obtained within different many-body Hamiltonians by
restricting Coulomb matrix elements in real space to the dominant ones. Within
the Hubbard model, one gets an effective mean-field Hamiltonian given as

$\begin{split}H^{H}_{MF}&=\sum\limits_{i,l,\sigma}t_{il\sigma}c_{i\sigma}^{\dagger}c_{l\sigma}+U\sum\limits_{i,\sigma}\left(\langle\hat{n}_{i,-\sigma}\rangle-\frac{1}{2}\right)\hat{n}_{i\sigma},\end{split}$
(3)

where the Hubbard parameter $U=\braket{ii}{V}{ii}$, $\hat{n}_{\rm
i}=c^{\dagger}_{i\sigma}c_{i\sigma}$, and
$\langle\hat{n}_{i,-\sigma}\rangle=\rho_{ii,\sigma^{\prime}}$, are diagonal
elements of density matrix with $\sigma^{\prime}=-\sigma$. For the extended
Hubbard model one has

$\begin{split}H^{EH}_{MF}&=\sum\limits_{i,l,\sigma}t_{il\sigma}c_{i\sigma}^{\dagger}c_{l\sigma}+U\sum\limits_{i,\sigma}\left(\langle\hat{n}_{i,-\sigma}\rangle-\frac{1}{2}\right)\hat{n}_{i\sigma}\\\
&+\sum\limits_{i\neq j,\sigma,\sigma^{\prime}}V_{\rm
ij}\left(\rho_{jj\sigma^{\prime}}-\frac{1}{2}\right)\hat{n}_{\rm
i\sigma},\end{split}$ (4)

where $V_{\rm ij}=\braket{ij}{V}{ji}$. The Coulomb matrix elements in the
basis of mean-field energies $\braket{pq}{V}{rs}$ are obtained using a basis
rotation

$\braket{pq}{V}{rs}=\sum\limits_{i,j}V_{\rm
ij}B^{*}_{ip}B^{*}_{jq}B_{jr}B_{is},$ (5)

where $b_{p\sigma}=\sum\limits_{i}B_{ip}c_{i\sigma}$, and $B_{ip}$ are
eigenvectors obtained by solving the mean-field Hamiltonian. The form of the
many-body Hamiltonian is the same for Hubbard and extended Hubbard, and full
interacting models given by Eq. 3 in the main article, but the models differ
in the parameters used defined by $\epsilon^{HF}_{p\sigma}$, $\tau_{pq\sigma}$
and Coulomb elements given by Eq. 5. A real space spin density distribution
for a given many-body state $D$, used when computing the spin density shown in
Fig. 4 in the main article, is calculated using the formula

$\begin{split}\braket{D}{\hat{n}_{i\sigma}}{D}&=\sum\limits_{p,q}B^{\prime*}_{pi}B^{\prime}_{qi}\braket{D}{b^{\dagger}_{p\sigma}b_{q\sigma}}{D}\\\
&=\sum\limits_{p,q}B^{\prime*}_{pi}B^{\prime}_{qi}\sum\limits_{c,d}A^{*}_{c}(D)A_{d}(D)\braket{c}{b^{\dagger}_{p\sigma}b_{q\sigma}}{d},\end{split}$
(6)

where $|D\rangle=\sum\limits_{c}A_{c}|c\rangle$ and $A_{c}$ are expansion
coefficients, $|c\rangle=\prod_{p,\sigma}b^{\dagger}_{p\sigma}|0\rangle$ are
occupation configuration states and $|0\rangle$ is our vacuum state
corresponding to HF state with all valence band states filled (all states
below CAS). $B^{\prime}_{pi}$ are coefficients from inverse transformation,
$c_{i\sigma}=\sum\limits_{p}B^{\prime}_{pi}b_{p\sigma}$, where
$B^{\prime}_{pi}=B^{*}_{ip}$.

Fig. S2 shows the many-body spectrum for the Hubbard model, extended Hubbard
model and the fully interacting model for $N_{\rm Tr}=2$ structure.
Calculations in each case were done for CAS(6,6). We find that the Hubbard
model tends to overestimate the spin gaps, with the extended Hubbard model
capturing most of the quantitative features of the full model. Notice that we
are showing only the three lowest energy states, as the excited states are
separated from these states by a large gap on the order of hundreds of meV.

Figure S2: Many-body spectrum obtained for the Hubbard model, extended Hubbard
model, and fully interacting model for CAS(6,6) calculations.

##  VI Extended Hubbard model results for $N_{\rm Tr}=4$ triangulenes

Using the extended Hubbard model, we analyze the $N_{\rm Tr}=4$ structure and
discuss properties of longer triangulene chains. In Fig. S3, we show the
Hartree-Fock energy spectrum indicating the energy gaps between the inter-
triangulene states (sectors $B$) and the degenerate shell states (sector $A$
and $C$), $\Delta_{\rm AB}$ and $\Delta_{\rm BC}$. Within this model, the
degenerate shell states split into a set of four double degenerate states.
There are three inter-triangulene states below and three above the degenerate
shell. The energy gap $\Delta_{\rm AB}$ is around twice smaller than
$\Delta_{\rm BC}$, and this relation is approximately true also for longer
chains, see the inset. These gaps determine the role of excitations and need
to be related to scattering Coulomb matrix elements. For $\kappa=3$ used here,
we find some of scattering elements between sector $A$ and $B$, and between
$B$ and $C$ as large as $0.2$ eV (e.g. $\braket{AA}{v}{BB}$ and
$\braket{AB}{v}{BC}$). Comparing them to energy gaps $\Delta_{\rm AB}\sim 1.8$
eV, $\Delta_{\rm BC}\sim 0.7$ eV, one can conclude that the states from sector
$C$ are as important as the states from sector $A$, see also Section VIII,
where we discuss the superexchange mechanism.

Figure S3: Hartree-Fock spectrum from the extended Hubbard model for $N_{\rm
Tr}=4$. Three sets of states, $A$, $B$, and $C$ are indicated, with energy
gaps between them, $\Delta_{\rm AB}$ and $\Delta_{\rm BC}$. The inset shows
the scaling of the energy gaps with the number of triangulenes $N_{\rm Tr}$ in
a chain.

Fig. S4 shows charge densities in these three energy sectors; the densities
are normalized by the number of states in a given sector. Sector $A$ has
localized charge density in the center of the chain, the central inter-
triangulene connection and on the first and last triangulenes. On can notice a
higher charge density on carbon atoms from the same sublattice as the
degenerate shell states, which confirms hybridization with states from sector
$B$ (within a tight-binding model, charge density from sector $A$ is similar
to that from sector $C$). In sector $B$ charge density is mainly localized on
the edges of triangulenes, and only one sublattice. Sector $C$ has charge
density on the three connections between the four triangulenes.

Figure S4: Charge densities of Hartree-Fock spectrum from extended Hubbard
model for $N_{\rm Tr}=4$ in three energy sectors indicated in Fig. S3.

Fig. S5 shows the many-body spectra obtained using the Hamiltonian given by
Eq. 3 in the main article with states from only sector $B$ (8 states) and with
all three sectors $A$, $B$, and $C$ (14 states). Similar to the $N_{\rm Tr}=2$
case, inclusion of the inter-triangulene states increases the splitting
between singlet and triplet, and triplet and quintuplet. One can also see the
order of total spin states (we show only the ten lowest energy states for the
8 state calculation) agrees with the order of states in BLBQ model for the
lowest five states, compared with Fig. S6. The sixth state in the extended
Hubbard model calculations with 14 states has $S=2$ (seventh has $S=0$). We
attribute these differences between this extended Hubbard model and the BLBQ
model to a CAS that is not large enough to converge excited states in this
extended Hubbard model calculation.

Figure S5: Many-body spectrum from extended Hubbard model for $N_{\rm Tr}=4$
for CAS(8,8) and CAS(14,14).

##  VII Spin model results

Figure S6: The BLBQ model results as a function of $\beta$ parameter for spin
chains with (a) $N=2$ and (b) $N=4$ sites. The blue dashed line indicates the
spectra from the Fermionic model in the main text.

The antiferromagnetic spin-1 Heisenberg model with the added biquadratic
operator can effectively describe the chain of triangulane molecules. The
effective Hamiltonian is:

$H=J\sum_{i}[(\bar{S}_{i}\cdot\bar{S}_{i+1})+\beta(\bar{S}_{i}\cdot\bar{S}_{i+1})^{2}],$
(7)

where J is the coupling constant, $\beta$ is the biquadratic term amplitude,
and $\bar{S}_{i}$ is a spin vector of the ith molecule in the chain.

Using the definition of the ladder operators
$\hat{S}_{i}^{\pm}=\hat{S}_{i}^{x}\pm i\hat{S}_{i}^{y}$, one can rewrite the
Hamiltonian as:

$\begin{split}H&=J\sum_{i}[\frac{1}{2}(\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}+\hat{S}_{i}^{-}\hat{S}_{i+1}^{+})+\hat{S}_{i}^{z}\hat{S}_{i+1}^{z}\\\
&+\frac{\beta}{4}(\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}+\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}\\\
&+\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}+\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}\hat{S}_{i}^{+}\hat{S}_{i+1}^{-})\\\
&+\frac{\beta}{2}(\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}\hat{S}_{i}^{z}\hat{S}_{i+1}^{z}+\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}\hat{S}_{i}^{z}\hat{S}_{i+1}^{z}\\\
&+\hat{S}_{i}^{z}\hat{S}_{i+1}^{z}\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}+\hat{S}_{i}^{z}\hat{S}_{i+1}^{z}\hat{S}_{i}^{-}\hat{S}_{i+1}^{+})\\\
&+\beta(\hat{S}_{i}^{z}\hat{S}_{i+1}^{z})^{2}].\end{split}$ (8)

The single spin-1 site has three possible $\hat{S}_{i}^{z}$ projections:
$\ket{\uparrow}=\ket{1}$, $\ket{\downarrow}=\ket{-1}$, and $\ket{0}$. To
obtain analytical results, we construct the basis for N=2 chain sites,
composed of all combinations of $\hat{S}_{i}^{z}$ values. The $\hat{S}^{z}=0$
subspace is spanned by three basis vectors: $\ket{0,0}$, $\ket{-1,1}$,
$\ket{1,-1}$. Acting with the Hamiltonian operator on this subspace results
in:

$\begin{split}&\hat{H}\ket{-1,1}=J(1-\beta)\ket{0,0}+J\beta\ket{1,-1}+J(2\beta-1)\ket{-1,1},\\\
&\hat{H}\ket{0,0}=2J\beta\ket{0,0}+J(1-\beta)\ket{1,-1}+J(1-\beta)\ket{-1,1},\\\
&\hat{H}\ket{1,-1}=J(1-\beta)\ket{0,0}+J(2\beta-1)\ket{1,-1}+J\beta\ket{-1,1},\\\
\end{split}$ (9)

which gives us the matrix:

$H_{S^{z}=0}=J\begin{bmatrix}(2\beta-1)&(1-\beta)&\beta\\\
(1-\beta)&2\beta&(1-\beta)\\\ \beta&(1-\beta)&(2\beta-1)\end{bmatrix}.$ (10)

Similarly, Hamiltonian matrices for other $S^{z}$ subspaces are:

$\begin{split}H_{|S^{z}|=1}&=J\begin{bmatrix}\beta&1\\\
1&\beta\end{bmatrix}\\\
H_{|S^{z}|=2}&=J\begin{bmatrix}(1+\beta)\end{bmatrix}.\end{split}$ (11)

The full Hamiltonian matrix is block diagonal, and its eigenvalues for the
singlet, triplet, and quintuplets are $E(S)=-2J(1-2\beta)$,
$E(T)=-J(1-\beta)$, and $E(Q)=J(1+\beta)$. After shifting eigenvalues by
$J(1+\beta)$, we see that only the singlet state depends on the $\beta$:
$E(S)=-3J(1-\beta)$, $E(T)=-2J$, and $E(Q)=0$. For $\beta<\frac{1}{3}$ the
ground state is a singlet and for $\beta>\frac{1}{3}$ it is a triplet state.
For $\beta=\frac{1}{3}$, the singlet and triplet states are degenerate, as
shown in Fig. S6(a).

Fig. S6(b) shows numerical results for a chain of $N=4$ spin-1 sites. The
transition between the singlet and triplet ground state again occurs at
$\beta=\frac{1}{3}$, but the gap between triplet and quintuplet now depends on
$\beta$.

##  VIII Superexchange mechanism analysis

Figure S7: Perturbation analysis of coupling between spins. $0$ indicates the
lowest energy total spin states obtained from diagonalization of many-body
Hamiltonian matrix within a single occupation subspace of the degenerate shell
states. Perturbation indicates energies of these states after energy
corrections included within 2nd order perturbation theory. Exact are energies
obtained after diagonalization full many-body Hamiltonian.

Isolated triangulene has a triplet ground state and can be represented by an
effective spin-1 site. We analyze the processes responsible for coupling
between neighboring spin-1 states on the example of two triangulenes. Our
methodology for treating Coulomb interaction relies on a two-step process.
First, we include interactions at the Hartree-Fock level for a closed shell
system, and next, we populate unoccupied states up to charge neutrality and
diagonalize the many-body Hamiltonian within a restricted subspace. For the
$N_{\rm Tr}=2$ system, after the HF basis rotation, the degenerate shell
states from two triangles form symmetric and antisymmetric linear combinations
of states from each triangle, as seen in the wavefunctions of Fig. 3 of the
main article. Moreover, we observe strong mixing between the four degenerate
shell states and valence and conduction band inter-triangulene states.

This hinders the perturbative analysis. The two isolated triangulenes have
four perfectly degenerate edge states, two for each triangulene. The perfect
strategy would be to rotate the HF degenerate shell states back to the
isolated triangulene basis states. In that case, one could show that coupling
between triangulene spin-1 states are through the inter-triangulene states -
an indirect superexchange mechanism, one of the main results of this work.
However, after HF self-consistent calculations, the degenerate shell states
are too strongly hybridized with the inter-triangulene states.

We analyze the superexchange mechanism at the many-body level. We take six HF
states to construct our many-body Hilbert space, the highest valence band
state, four degenerate shell states and the lowest conduction band state.
These states are labelled 5-10 in Fig. 2(c) of the main text, in this
analysis, for simplicity we shift the labels to 1-6. The charge neutral system
has $N_{\rm el}=6$ electrons. For $S_{\rm z}=0$, one can construct in total,
$N_{\rm cf}=400$ configurations. The low-energy subspace, which we call later
a single occupation subspace, corresponds to double occupation of the valence
band state and single occupation of four degenerate shell states. One can
construct six such low-energy configurations. In the occupation
representation, these states can be written as

$\begin{split}|I\rangle&=|1\downarrow,1\uparrow,2\downarrow,3\downarrow,4\uparrow,5\uparrow\rangle\\\
|II\rangle&=|1\downarrow,1\uparrow,2\uparrow,3\uparrow,4\downarrow,5\downarrow\rangle\\\
|III\rangle&=|1\downarrow,1\uparrow,2\downarrow,3\uparrow,4\downarrow,5\uparrow\rangle\\\
|IV\rangle&=|1\downarrow,1\uparrow,2\uparrow,3\downarrow,4\uparrow,5\downarrow\rangle\\\
|V\rangle&=|1\downarrow,1\uparrow,2\downarrow,3\uparrow,4\uparrow,5\downarrow\rangle\\\
|VI\rangle&=|1\downarrow,1\uparrow,2\uparrow,3\downarrow,4\downarrow,5\uparrow\rangle,\end{split}$
(12)

where state $1$ is the valence band state lying just below the four degenerate
shell states, states $2-5$ are the four degenerate shell states, and state $6$
is the state just above the four degenerate states (unoccupied here). The six
configurations defined in Eq. 12 can then be rotated into two total spin $S=0$
singlets, three $S=1$ triplets and one $S=2$ quintuplet. We obtain the lowest
energy state $E_{0}(S)$ within each total spin single occupation subspace by
diagonalizing a $2\times 2$ matrix within $S=0$ subspace and a $3\times 3$
matrix within $S=1$ subspace, at the same time appropriately rotating the full
Hilbert space within each total spin sector. This procedure is related to the
strong coupling between configurations given by Eq. 12, due to the
hybridization between the degenerate shell states and inter-triangulene
states. Next, the obtained lowest energy states within each total spin
subspace, are corrected by second order perturbation contributions due to the
coupling to the rest of states $E_{i}(S)$ from the many-body Hilbert space
(which have been rotated to the total spin subspaces), beyond the single
occupation subspace ($S=0$ subspace contains in total $N_{\rm cf}=175$
configurations, $S=1$ contains $N_{\rm cf}=189$ configurations, and $S=2$
contains $N_{\rm cf}=35$ configurations). The perturbative Hamiltonian is
written as

$E^{(2)}(S)=\sum_{i=1}^{N_{\rm
cf}(S)}\frac{|\langle\Psi_{0}(S)|H|\Psi_{i}(S)\rangle|^{2}}{E_{0}(S)-E_{\rm
i}(S)},$ (13)

where $\Psi_{0}$ and $\Psi_{i}$ are the many-body wave function corresponding
to energy $E_{0}$ and $E_{i}$, respectively, and $H$ is the many-body
Hamiltonian given by Eq. (3) in the main article. Although this procedure is
second order in perturbation theory, if the hybridization between the edge
states of the individual triangulenes and the inter-triangulene states were
weak, then perturbative treatment of the lowest energy total spin states
(corresponding to singlet, triplet and quintuplet) would be appropriate. In
that case, fourth order perturbation theory would be required in order to
describe coupling of the degenerate states of two triangulenes (and thus the
spin-1 quasiparticles) through the inter-triangulene states.

The energy spectra after diagonalization within the single occupation
subspace, and after second order perturbation correction is compared to the
exact many-body spectrum in Fig. S7. The order of total spin states after
basis rotation within a single occupation subspace is the same as within
truncated $N_{\rm st}=4$ states with the $S=2$ quintuplet as the ground state,
see Fig. 4(a) in the main article. Coupling of the single occupation subspace
to higher energy configurations leads to a change of the order of states that
now agrees with predictions within the two spin-1 Heisenberg Hamiltonian, with
the singlet as the ground state. Furthermore, results after this coupling,
agrees as well as with the experiments FaselNaturespinchain . The energies are
close to the exact energies obtained after diagonalization of the full many-
body Hamiltonian. We notice that configurations mainly contributing to the
perturbation includes all six HF states, thus both the state below and the
state above the four degenerate shell states are important, which confirms the
existence of an indirect AFM superexchange mechanism.

## References

  * (1) Ahlrichs, R., Bär, M., Häser, M., Horn, H. & Kölmel, C.  Electronic structure calculations on workstation computers: The program system turbomole.  _Chem. Phys. Lett._ 162, 165–169 (1989).
  * (2) Dirac, P. A. M.  Quantum mechanics of many-electron systems.  _Proc. Royal Soc. (London) A_ 123, 714–733 (1929).
  * (3) Slater, J. C.  A simplification of the hartree-fock method.  _Phys. Rev._ 81, 385 (1951).
  * (4) Vosko, S. H., Wilk, L. & Nusair, M.  Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis.  _Can. J. Phys._ 58, 1200–1211 (1980).
  * (5) Becke, A. D.  Density-functional exchange-energy approximation with correct asymptotic behavior.  _Phys. Rev. A_ 38, 3098 (1988).
  * (6) Lee, C., Yang, W. & Parr, R. G.  Development of the colle-salvetti correlation-energy formula into a functional of the electron density.  _Phys. Rev. B_ 37, 785 (1988).
  * (7) Becke, A. D.  Density-functional thermochemistry. i. the effect of the exchange-only gradient correction.  _J. Chem. Phys._ 96, 2155–2160 (1992).
  * (8) Weigend, F.  Accurate coulomb-fitting basis sets for h to rn.  _Phys. Chem. Chem. Phys._ 8, 1057–1065 (2006).
  * (9) Grimme, S., Antony, J., Ehrlich, S. & Krieg, H.  A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu.  _J. Chem. Phys._ 132 (2010).
  * (10) Eichkorn, K., Treutler, O., Öhm, H., Häser, M. & Ahlrichs, R.  Auxiliary basis sets to approximate coulomb potentials.  _Chem. Phys. Lett._ 240, 283–290 (1995).
  * (11) Eichkorn, K., Weigend, F., Treutler, O. & Ahlrichs, R.  Auxiliary basis sets for main row atoms and transition metals and their use to approximate coulomb potentials.  _Theor. Chem. Acc._ 97, 119–124 (1997).
  * (12) Sierka, M., Hogekamp, A. & Ahlrichs, R.  Fast evaluation of the coulomb potential for electron densities using multipole accelerated resolution of identity approximation.  _J. Chem. Phys._ 118, 9136–9148 (2003).
  * (13) Noodleman, L.  Valence bond description of antiferromagnetic coupling in transition metal dimers.  _J. Chem. Phys._ 74, 5737–5743 (1981).
  * (14) Güçlü, A. D., Potasz, P., Korkusinski, M., Hawrylak, P. _et al._ _Graphene quantum dots_ (Springer, 2014).
  * (15) Ransil, B. J.  Studies in molecular structure. ii. lcao-mo-scf wave functions for selected first-row diatomic molecules.  _Rev. Mod. Phys._ 32, 245–254 (1960).  URL https://link.aps.org/doi/10.1103/RevModPhys.32.245.
  * (16) Galassi, M. _et al._ _GNU scientific library_ (Network Theory Limited Godalming, 2002).
  * (17) Lieb, E. H.  Two theorems on the hubbard model.  _Phys. Rev. Lett._ 62, 1201–1204 (1989).  URL https://link.aps.org/doi/10.1103/PhysRevLett.62.1201.
  * (18) Pavliček, N. _et al._ Synthesis and characterization of triangulene.  _Nature Nanotechnology_ 12, 308–311 (2017).  URL https://doi.org/10.1038/nnano.2016.305.
  * (19) Mishra, S. _et al._ Synthesis and characterization of $\pi$-extended triangulene.  _Journal of the American Chemical Society_ 141, 10621–10625 (2019).  URL https://doi.org/10.1021/jacs.9b05319.
  * (20) Ortiz, R., Catarina, G. & Fernández-Rossier, J.  Theory of triangulene two-dimensional crystals.  _2D Materials_ 10, 015015 (2022).
  * (21) Mishra, S. _et al._ Observation of fractional edge excitations in nanographene spin chains.  _Nature_ 598, 287–292 (2021).  URL https://doi.org/10.1038/s41586-021-03842-3.