# Supplemental Material for “Multi-orbital Kondo screening in strongly
correlated polyradical nanographenes”

Aitor Calvo-Fernández  Departamento de Física, Universidad del País Vasco UPV-
EHU, 48080 Leioa, Spain  Donostia International Physics Center (DIPC), 20018
Donostia-San Sebastián, Spain  Diego Soler-Polo  Institute of Physics, Academy
of Sciences of the Czech Republic, Cukrovarnicka 10, Prague 6, CZ 16200, Czech
Republic  Andrés Pinar Solé  Institute of Physics, Academy of Sciences of the
Czech Republic, Cukrovarnicka 10, Prague 6, CZ 16200, Czech Republic  Shaotang
Song  Department of Chemistry, National University of Singapore, Singapore
117543, Singapore  Oleksander Stetsovych  Institute of Physics, Academy of
Sciences of the Czech Republic, Cukrovarnicka 10, Prague 6, CZ 16200, Czech
Republic  Manish Kumar  Institute of Physics, Academy of Sciences of the Czech
Republic, Cukrovarnicka 10, Prague 6, CZ 16200, Czech Republic  Guangwu Li
Department of Chemistry, National University of Singapore, Singapore 117543,
Singapore  Jishan Wu  Department of Chemistry, National University of
Singapore, Singapore 117543, Singapore  Jiong Lu  Department of Chemistry,
National University of Singapore, Singapore 117543, Singapore  Institute for
Functional Intelligent Materials, National University of Singapore, 117544,
Singapore  Asier Eiguren  Departamento de Física, Universidad del País Vasco
UPV-EHU, 48080 Leioa, Spain  Donostia International Physics Center (DIPC),
20018 Donostia-San Sebastián, Spain  EHU Quantum Center, Universidad del País
Vasco UPV-EHU, 48080 Leioa, Spain  María Blanco-Rey  Departamento de Polímeros
y Materiales Avanzados: Física, Química y Tecnología, Universidad del País
Vasco UPV-EHU, 20018 Donostia-San Sebastián, Spain  Donostia International
Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain  Centro de Física
de Materiales CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain
Pavel Jelínek  Institute of Physics, Academy of Sciences of the Czech
Republic, Prague, Czech Republic  Regional Centre of Advanced Technologies and
Materials, Czech Advanced Technology and Research Institute (CATRIN), Palacký
University Olomouc,Olomouc 78371, Czech Republic

##  S1 Experiments

The low-temperature SPM experiments were performed using a commercial low-
temperature (1.2 K) STM/nc-AFM microscope Specs-JT with a Kolibri sensor (f
$\approx$ 1 MHz) with a Standford Lock-in amplifier. High-resolution AFM and
dI/dV maps were performed using a functionalized tip with a single carbon
monooxide (CO) molecule. CO was dosed to the Au(111) at 4K. Typical CO
coverages of two-three molecules per 20 $\times$ 20 nm2 sample scan area were
used. CO was picked up from the Au(111) by placing the tip on top of the
molecule and increasing the sample bias in constant height STM mode until the
CO functionalizes the tip ($\approx$ 200 pA at 130mV). The high-resolution AFM
imaging was performed at constant height mode with bias voltage of 1 mV with
oscillation amplitude of 50 pm. Lock-in parameters during dI/dV measurements:
frequency: 723 Hz, modulation amplitude: 0.5 mV, time constant: 30 ms.

Nickelocene (NiCp2) was dosed from a tantalum pocket kept at room temperature
onto the prepared samples on Au(111) with submonolayer coverage at T $<$ 4K.
To pick up the NiCp2 on the STM tip, the sample was scanned at low current (10
pA STM setpoint and below 5 mV sample bias) until the nickelocene
functionalizes the tip. Only stable NiCp2 tips were used for IETS
measurements. Additionally, the NiCp2 tips were tested prior measurements by
performing d2I/dV2(z) on Au(111) to observe the characteristic IETS peaks at
$\approx$4 mV of the unperturbed NiCp2 molecule.

The fitting of the Kondo peaks was done with a Python script [1] that includes
fitting routines with a Frota function [2] and a thermally broadened Fano
function [3], where an operator is defined for the convolution of the Fano
line-shape with a derivative of Fermi-Dirac distribution. The fitting was
performed with Model class in LMFIT, using Levenberg-Marquardt least-squares
method.

###  S1.1 Magnetic origin of zero-bias peaks

To demonstrate the magnetic origin of zero energy peaks observed in STS
spectra, shown in Figure 1 in the main text, we carried out STS measurements
with nickelocene-functionalized tip. Figure S1 displays the STS spectra
measured along tip approach towards the sample, which clearly show a strong
renormalization of the spin-excitation of the nickelocene tip at 5 meV along
the tip approach due to exchange interaction between the spin of nickelocene
tip and products 1 and 2, respectively.

Figure S1:  $d^{2}I/dV^{2}$ IETS spectra of 1 and 2 on Au(111) acquired with
NiCp2-functionalized tip. (a, e) Chemical schemes of 1 and 2 molecules.
Constant height STM images of 1 at 3.4 mV (b) and 2 at 1 mV (f), image size:
2.5 $\times$ 2.5 nm2. (c, g) Experimental $d^{2}I/dV^{2}(z)$ IETS taken with
NiCp2-functionalized tip in the positions marked by blue crosses in top and
bottom panels of (b, f) respectively. (d, h) $d^{2}I/dV^{2}$ curves for
selected heights in (c, g) marked with the horizontal lines. Red/blue colored
lines in (c, g) correspond to red/blue spectra in (d, h). While the NiCp2
$d^{2}I/dV^{2}(z)$ IETS remains unaltered when approaching to a non magnetic
site of the molecules (b top, f top, it deforms on the egde providing a
fingerprint of the local magnetic signal originated from the Kondo effect[4]
(b bottom, f bottom). Figure S2: Spin-excitation signal from the doublet
ground state to the first excited quartet state of molecule 1. a) scheme of
the molecule, blue lines highlight the location of the defect. Simultaneously
recorded high resolution CO tip AFM (b) and ${\it dI/dV}$ map (c) images of
the molecule 1 acquired at 24 meV, size: 3 $\times$ 3 nm2. It should be noted
that the signal on the left edge in the experimental dI/dV map is caused by
the presence of the Kondo resonance located at the Fermi energy. Thus, the
appearance of the signal on the left edge has nothing to do with spin
excitation from the ground state to the first excited state. d) simulated spin
excitation map obtained from the Natural Transition Orbitals for the spin flip
operator between the ground state doublet and the first excited quartet state,
see Figure S5.

##  S2 Calculations

###  S2.1 Natural transition orbitals

To calculate Natural Transition Orbitals (NTO) orbitals between two states,
$|\Psi_{-}\rangle$ and $|\Psi_{+}\rangle$ we construct for each $j$ and $k$
the matrix element:

$T_{j,k}=\langle\Psi_{-}|\hat{f}^{\dagger}_{j\downarrow}\hat{f}_{k\uparrow}|\Psi_{+}\rangle,$

where, as above, $j$ is an index labeling molecular orbitals and
$\hat{f}^{\dagger}_{j\sigma}$ are the creation operators for electrons in
molecular spin-orbital $j\sigma.$ We then proceed, following Martin [5], to
diagonalize the matrix $TT^{\dagger}$ to obtain the NTOs. Finally, we use the
PP-STM code [6] to calculate corresponding dI/dV maps using s-like tip.

Figure S3: Single electron Hückel orbitals used for the Hubbard-CAS
calculations in product 1: 1 (HOMO-1), 2(HOMO), 3(SOMO), 4(LUMO), 5(LUMO+1).
The numbering of orbitals is according to the Figure 2 in the main text.
Figure S4: Single electron Hückel orbitals used for the Hubbard-CAS
calculations in product 2. First row, from left to right displays as follows :
1 (HOMO-2), 2(HOMO-1) and 3(HOMO). Second row, from left to right: 4 (LUMO),
5(LUMO+1), 6(LUMO+2). The numbering of orbitals is according to the Figure 2
in the main text. Figure S5:  Calculated NTO and dI/dV maps of the doublet to
quartet spin excitation of molecule 1 a) Three dominant (from left to right)
NTO orbitals for the doublet-quartet spin excitation of product 1, each with
their corresponding weights, b) resulting simulated dI/dV map obtained from
these NTOs. Figure S6:  Calculated NTO and dI/dV maps of the doublet to
doublet spin flip of molecule 1. a)Three dominant (from left to right) NTO
orbitals for the spin flip between two degenerate ground state doublets that
yields the Kondo effect in product 1 with their respective weights, b)
simulated dI/dV map obtained from these orbitals  Figure S7: Comparison of
dI/dV maps obtained from NTO and the Kondo channel for molecule 1. a) NTO
orbitals with their weigth and b) corresponding simulated dI/dV map. c) the
Kondo channel obtained from the diagonalization of the $J$ matrix, and d) the
corresponding simulated dI/dV map. Figure S8: Comparison of dI/dV maps
obtained from NTO and two Kondo channels for molecule 2. a) NTO orbitals with
their weight and b) corresponding simulated dI/dV map, c) two Kondo channels
obtained from the diagonalization of the $J$ matrix and d) the corresponding
simulated dI/dV map.  Figure S9:  Calculated NTO and dI/dV maps of the doublet
to doublet spin flip of molecule 2. a)Two dominant (from left to right) NTO
orbitals for the spin flip between two degenerate ground state doublets that
yields the Kondo effect in product 2 with their respective weights, b)
simulated dI/dV map obtained from these two NTOs.

###  S2.2 Spin-spin correlation of unpaired electrons

We evaluated the spin correlations between the unpaired electrons in the
degenerate ground states of products 1 and 2 as follows:

$\langle\Psi_{0}|\vec{\hat{S_{j}}}\cdot\vec{\hat{S_{k}}}|\Psi_{0}\rangle-\langle\Psi_{0}|\vec{\hat{S_{j}}}|\Psi_{0}\rangle\cdot\langle\Psi_{0}|\vec{\hat{S_{k}}}|\Psi_{0}\rangle,$

where the indices $j$ and $k$ design the molecular orbitals. Here negative
correlation represents antiferromagnetic coupling and positive ferromagnetic
coupling between unpaired spins. We should note that the strong entanglement
between three unpaired electrons is independent of the choice of the MO basis
set.

The correlations between the unpaired spins are summarized in the Tables S1
and S2, where the molecular orbital ”3” is the single-occupied zero-mode that
survives on the butterfly after introducing the defect.

| 2 | 3 | 4
---|---|---|---
2 | 0.050 | -0.060 | 0.014
3 | -0.060 | 0.127 | -0.062
4 | 0.014 | -0.062 | 0.052
Table S1: Spin-spin correlations in 1 for the doublet ground state in basis of Hückel molecular orbitals of the active space shown in Fig.S3.  | 3 | 4
---|---|---
3 | 0.25 | -0.24
4 | -0.24 | 0.25
Table S2: Spin-spin correlations in the $\bf 2$ for the doublet ground state
in basis of Hückel molecular orbitals of the active space shown in Fig.S4.

###  S2.3 Electronic structure of product 1: CASCI calculations

Here we describe the wave-function of the ground and excited states of 1
obtained from CASCI calculations for neutral ($N$) and charged ($N\pm 1$)
states.

Neutral molecule:

The 1 wave-function in its doublet ground state in the $S_{z}=1/2$ subspace is
given by:

$\displaystyle|M_{0}\rangle\approx 0.70|2\uparrow 0\rangle+0.47|0\uparrow
2\rangle+0.38|\uparrow\downarrow\uparrow\rangle+0.21|\uparrow\uparrow\downarrow\rangle,$

whereas the $S_{z}=-1/2$ is obtained by applying the spin lowering operator.

The quadruplet excited state 1, located at 25 meV of the doublet ground state,
has a somewhat simpler structure in the subspace $S_{z}=1/2$:

$|M_{1},S_{z}=1/2\rangle\approx
0.57|\uparrow\uparrow\downarrow\rangle+0.57|\uparrow\downarrow\uparrow\rangle+0.57|\downarrow\uparrow\uparrow\rangle.$

In the two equations above, orbitals 1 and 5 are dismissed, since their
occupations are 2 and 0, and do not change.

Molecule with one extra electron: The two main multiplets in the subspace
$S_{z}=0$ are

$|M^{+}_{0},S_{z}=0\rangle\approx 0.64|22\uparrow\downarrow
0\rangle+0.64|22\downarrow\uparrow 0\rangle+0.13|2\uparrow 2\downarrow
0\rangle+0.13|2\downarrow 2\uparrow
0\rangle-0.14|2\uparrow\downarrow\uparrow\downarrow\rangle+0.14|2\downarrow\uparrow\downarrow\uparrow\rangle,$
$|M^{+}_{1},S_{z}=0\rangle\approx-0.72|22200\rangle+0.38|22020\rangle,$

whereas in the $S_{z}=1$ subspace we have

$|M^{+}_{0},S_{z}=1\rangle\approx-0.91|22\uparrow\uparrow
0\rangle-0.18|2\uparrow 2\uparrow
0\rangle-0.19|2\uparrow\uparrow\uparrow\downarrow\rangle-0.11|22\uparrow
0\uparrow\rangle-0.11|2\uparrow\uparrow 20\rangle.$

Their energies with respect to the ground state of the neutral molecule
(CAS(5,5) for the Hubbard hamiltonian) are respectively: 487 meV and 491 meV.

Molecule with one less electron ($N-1$): The two main multiplets are

$|M^{-}_{0},S_{z}=0\rangle\approx 0.70|2\uparrow\downarrow
00\rangle+0.70|2\downarrow\uparrow 00\rangle,$
$|M^{-}_{0},S_{z}=0\rangle\approx
0.87|22000\rangle+0.32|20200\rangle+0.30|20020\rangle.$

Their energies with respect to the ground state of the neutral molecule
(CAS(5,5) for the Hubbard hamiltonian) are respectively: -80 meV and -78 meV.

###  S2.4 Electronic structure of product 2: CASCI calculations

The triplet and singlet many-body wavefunctions extracted from a CAS(6,6) are
(here we only show occupancies of the orbitals 3 (HOMO) and 4(LUMO) because
they are the most important ones within 96% accuracy for the triplet ground
state and the first excited singlet state):

Neutral molecule:

$|M_{0},S_{z}=0\rangle\approx
0.7|\uparrow\downarrow\rangle+0.7|\downarrow\uparrow\rangle.$

In the subspace $S_{z}=1$ it is simply to 96% of the Slater determinant:

$|M_{0},S_{z}=1\rangle\approx|\uparrow\uparrow\rangle.$

The excited singlet $|M_{1}\rangle$, on the other hand, is simply:

$|M_{1}\rangle\approx 0.7|20\rangle+0.7|02\rangle$

(here the occupation numbers always refer to the orbitals 3(HOMO) and
4(LUMO)).

Molecule with one extra electrons: The two main multiplets in the subspace are

$|M_{0}^{+},S_{z}=1/2\rangle\approx 0.78|2\uparrow\rangle+0.52|\uparrow
2\rangle,$
$|M_{1}^{+},S_{z}=1/2\rangle\approx-0.54|2\uparrow\rangle+0.81|\uparrow
2\rangle.$

The energies of these multiplets with respect to the neutral Ground State are:
2.37eV and 2.41 eV

Molecule with one less electron: The two main multiplets are

$|M_{0}^{-},S_{z}=1/2\rangle\approx 0.78|0\uparrow\rangle+0.52|\uparrow
0\rangle,$
$|M_{1}^{-},S_{z}=1/2\rangle\approx-0.54|0\uparrow\rangle+0.81|\uparrow
0\rangle.$

The energies of these multiplets with respect to the neutral Ground State are:
-1.91eV and -1.89eV.

##  S3 Ionic Anderson Hamiltonian in the channel-energy representation

In the discrete momentum space, the hybridization and conduction terms of the
ionic Anderson Hamiltonian are

$\displaystyle\hat{H}_{\text{con}}=\sum_{\mathbf{k},\sigma}\epsilon_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k}\sigma}\hat{c}_{\mathbf{k}\sigma},$
(1)
$\displaystyle\hat{H}_{\text{hyb}}=\sum_{M_{i}m_{i}}\sum_{M_{f}m_{f}}\sum_{\alpha,\mathbf{k},\sigma}V_{\alpha\mathbf{k}}\bra{M_{f}m_{f}}f^{\dagger}_{\alpha\sigma}\ket{M_{i}m_{i}}\ket{M_{f}m_{f}}\bra{M_{i}m_{i}}c_{\mathbf{k}\sigma}+\text{h.c.},$
(2)

where $\epsilon_{\mathbf{k}}$ is the energy of the conduction electrons,
$\hat{c}_{\mathbf{k}\sigma}$ ($\hat{c}^{\dagger}_{\mathbf{k}\sigma}$)
annihilates (creates) a conduction electron with momentum $\mathbf{k}$ and
spin $\sigma$, and $V_{\alpha\mathbf{k}}$ is the amplitude of an electron
tunneling from the MO $\alpha$ to a conduction state with momentum
$\mathbf{k}$. We begin by re-writing the sums in the continuous momentum
space,

$\displaystyle\sum_{\mathbf{k}}\epsilon_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k}\sigma}\hat{c}_{\mathbf{k}\sigma}$
$\displaystyle\rightarrow\int
d^{3}k\epsilon_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{k}\sigma}\hat{a}_{\mathbf{k}\sigma},$
(3) $\displaystyle\sum_{\mathbf{k}}V_{\alpha k}\hat{a}_{\mathbf{k}\sigma}$
$\displaystyle\rightarrow\nu^{\frac{1}{2}}\int
d^{3}kV_{\alpha\mathbf{k}}\hat{a}_{\mathbf{k}\sigma}.$ (4)

where $a^{(\dagger)}_{\mathbf{k}\sigma}$ is the annihilation (creation)
operator in the continuous $\mathbf{k}$-space, normalized as
$\\{a^{\dagger}_{\mathbf{k}\sigma},a_{\mathbf{k}^{\prime}\sigma}\\}=\delta(\mathbf{k}-\mathbf{k}^{\prime})$,
and $\nu$ is the density of $\mathbf{k}$ points. In order to change to the
energy representation, we define the operators
$\hat{c}^{(\dagger)}_{\epsilon\hat{k}\sigma}=|\nabla_{\mathbf{k}}\epsilon_{\mathbf{k}}|^{-\frac{1}{2}}a^{(\dagger)}_{\mathbf{k}\sigma}$,
where $\hat{k}=\mathbf{k}/|\mathbf{k}|$ and the right hand side is evaluated
at $\epsilon_{\mathbf{k}}=\epsilon$. This allows us to write the equations
above as

$\displaystyle\int
d^{3}k\epsilon_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{k}\sigma}\hat{a}_{\mathbf{k}\sigma}=\int
d\epsilon\int_{\epsilon}ds_{\mathbf{k}}\epsilon\hat{c}^{\dagger}_{\epsilon\hat{k}\sigma}\hat{c}_{\epsilon\hat{k}\sigma}$
(5) $\displaystyle\nu^{\frac{1}{2}}\int
d^{3}kV_{\alpha\mathbf{k}}\hat{a}_{\mathbf{k}\sigma}=\int
d\epsilon\int_{\epsilon}\frac{ds_{\mathbf{k}}}{\sqrt{v(\mathbf{k})}}V_{\alpha\mathbf{k}}\hat{c}_{\epsilon\hat{k}\sigma},$
(6)

where $\int_{\epsilon}ds_{\mathbf{k}}$ represents an integral over the
momentum space surface of $\mathbf{k}$ points for which
$\epsilon_{\mathbf{k}}=\epsilon$, and
$v^{-1}(\mathbf{k}):=v(\epsilon,\hat{\mathbf{k}})=\nu|\nabla_{\mathbf{k}}\epsilon|^{-1}_{\epsilon=\epsilon_{\mathbf{k}}}$.
Based on Eq. 6, we define the channel-energy representation operators as

$\hat{c}_{\epsilon\alpha\sigma}=\frac{1}{[\rho(\epsilon)]^{\frac{1}{2}}V_{\alpha}(\epsilon)}\int_{\epsilon}\frac{ds_{\mathbf{k}}}{\sqrt{v(\mathbf{k})}}V_{\alpha\mathbf{k}}\hat{c}_{\epsilon\hat{k}\sigma},$
(7)

where

$V_{\alpha}(\epsilon)=\left[\frac{\int_{\epsilon}\frac{ds_{\mathbf{k}}}{v(\mathbf{k})}|V_{\alpha\mathbf{k}}|^{2}}{\int
ds_{\mathbf{k}}}\right]^{\frac{1}{2}}$ (8)

is the average of $|V_{\alpha\mathbf{k}}|^{2}$ over the
$\epsilon_{\mathbf{k}}=\epsilon$ surface and

$\rho(\epsilon)=\frac{\int_{\epsilon}\frac{ds_{\mathbf{k}}}{v(\mathbf{k})}}{\int_{\epsilon}ds_{\mathbf{k}}}$
(9)

is the density of states per unit of energy. In general, these operators do
not provide an orthonormal basis, since

$\\{\hat{c}^{\dagger}_{\epsilon\alpha\sigma},\hat{c}_{\epsilon^{\prime}\alpha^{\prime}\sigma^{\prime}}\\}=\frac{\delta(\epsilon-\epsilon^{\prime})\delta_{\sigma\sigma^{\prime}}}{\rho(\epsilon)V_{\alpha}(\epsilon)V_{\alpha^{\prime}}(\epsilon)}\int_{\epsilon_{\mathbf{k}}=\epsilon}\frac{ds_{\mathbf{k}}}{v(\mathbf{k})}V^{*}_{\alpha\mathbf{k}}V_{\alpha^{\prime}\mathbf{k}}.$
(10)

Nonetheless, we know from experimental measurements that that the
hybridization between the MOs $\alpha$ is small, so we assume that

$\int_{\epsilon}\frac{ds_{\mathbf{k}}}{v(\mathbf{k})}V_{\alpha\mathbf{k}}V_{\alpha^{\prime}\mathbf{k}}\approx\rho(\epsilon)V_{\alpha}(\epsilon)^{2}\delta_{\alpha\alpha^{\prime}},$
(11)

which leaves us with an orthonormal basis that fulfils

$\\{c^{\dagger}_{\epsilon\alpha\sigma},c_{\epsilon^{\prime}\alpha^{\prime}\sigma}\\}=\delta(\epsilon-\epsilon^{\prime})\delta_{\alpha\alpha^{\prime}}\delta_{\sigma\sigma^{\prime}}.$
(12)

Thus, in the channel-energy representation the full ionic Anderson Hamiltonian
is

$\displaystyle\hat{\mathcal{H}}=\hat{H}_{\text{mol}}+\hat{H}_{\text{con}}+\hat{H}_{\text{hyb}},$
(13)
$\displaystyle\hat{H}_{\text{mol}}=\sum_{M}\sum_{m}\epsilon_{M}\ket{Mm}\bra{Mm},$
$\displaystyle
H_{\text{con}}=\sum_{\alpha\sigma}\int\epsilon\hat{c}^{\dagger}_{\epsilon\alpha\sigma}\hat{c}_{\epsilon\alpha\sigma}d\epsilon$
$\displaystyle
H_{\text{hyb}}=\sum_{M_{i}m_{f}}\sum_{M_{f}m_{i}}\sum_{\alpha\sigma}\int
d\epsilon[\rho(\epsilon)]^{\frac{1}{2}}V_{\alpha}(\epsilon)\bra{M_{f}m_{f}}f^{\dagger}_{\alpha\sigma}\ket{M_{i}m_{i}}\ket{M_{f}m_{f}}\bra{M_{i}m_{i}}c_{\epsilon\alpha\sigma}+\text{h.c.},$

where we have introduced the $\hat{H}_{\text{mol}}$ containing the energies
$\epsilon_{M}$ of each molecular multiplet $\ket{Mm}$ The Hamiltonian used in
the main paper is obtained from considering the hybridization and the density
of states as constants by performing the changes
$V_{\alpha}(\epsilon)\rightarrow V$ and $\rho(\epsilon)\rightarrow\rho$.

##  S4 Effective Kondo coupling

###  S4.1 Schrieffer-Wolff transformation

The effective Kondo interaction for the ionic Anderson Hamiltonian to second
order in perturbation theory takes the form

$\displaystyle\hat{H}_{K}$
$\displaystyle=\sum_{m_{0},m_{0}^{\prime}}\int_{-D}^{D}d\epsilon\int_{-D}^{D}d\epsilon^{\prime}T(m_{0},\epsilon\alpha
i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})\ket{M_{0},m_{0}}\bra{M_{0},m_{0}^{\prime}}c^{\dagger}_{\epsilon\alpha
i}c_{\epsilon^{\prime}\alpha^{\prime}i^{\prime}}$ (14)

where $T(m_{0},\alpha\epsilon i\leftarrow
m_{0}^{\prime},\alpha^{\prime}\epsilon^{\prime}i^{\prime})$ is the second
order scattering amplitude corresponding to virtual processes where an
incoming electron with quantum numbers $\epsilon^{\prime}$ (energy),
$\alpha^{\prime}$ (channel) and $i^{\prime}$ (spin) is scattered off the
molecule into a state with quantum numbers $\epsilon$, $\alpha$ and $\sigma$,
and the molecule spin changes its spin from $m_{0}^{\prime}$ to $m_{0}$. This
amplitude can be written as a sum of contributions from the intermediate
excited multiplets,

$\displaystyle T(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$
$\displaystyle=\sum_{M_{c}}T^{N_{c}}_{S_{c}}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$ (15)

where $T^{N_{c}}_{S_{c}}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$ is the amplitude
for processes where the intermediate molecular state belongs to the $M_{c}$
multiplet, which has occupation and spin quantum numbers $N_{c},S_{c}$. The
Feynman diagrams for such processes are given schematically in Fig. S10.

$M_{0}$$\epsilon\alpha
i$$M_{0}$$\epsilon^{\prime}\alpha^{\prime}i^{\prime}$$M_{c}$

$M_{0}$$\epsilon\alpha
i$$M_{0}$$\epsilon^{\prime}\alpha^{\prime}i^{\prime}$$M_{c}$

Figure S10: Schematic Feynman diagrams for virtual scattering process through
the intermediate molecular multiplet $M_{c}$ with $N_{c}=N_{0}+1$ (left) and
$N_{c}=N_{0}-1$ (right).

The general expressions for scattering through an intermediate state with one
more electron ($N_{c}=N_{0}+1$) and one less electron ($N_{c}=N_{0}-1$) are

$\displaystyle T^{N_{0}+1}_{S_{c}}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})=$
$\displaystyle-\rho
V^{2}\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}}{\epsilon_{M_{c}}-\epsilon_{M_{0}}}$
(16)
$\displaystyle\times(\tfrac{1}{2},i;S_{0},m_{0}+i^{\prime}-i|S_{c},m_{0}+i^{\prime})(\tfrac{1}{2},i^{\prime};S_{0},m_{0}|S_{c},m_{0}+i^{\prime}),$
$\displaystyle T^{N_{0}-1}_{S_{c}}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})=$
$\displaystyle\rho
V^{2}\frac{\bra{M_{0}}|f^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f_{\alpha^{\prime}}|\ket{M_{0}}}{\epsilon_{M_{c}}-\epsilon_{M_{0}}}$
(17)
$\displaystyle\times(\tfrac{1}{2},i;S_{c},m_{0}-i^{\prime}|S_{0},m_{0}-i^{\prime}+i)(\tfrac{1}{2},i^{\prime};S_{c},m_{0}-i^{\prime}|S_{0},m_{0}),$

where we have used the Wigner-Eckart theorem to decompose the excitation
matrix elements as

$\bra{M_{1},m_{1}}f^{\dagger}_{\alpha
i}\ket{M_{2},m_{2}}=\bra{M_{1}}|f^{\dagger}_{\alpha}|\ket{M_{2}}(\frac{1}{2},i;S_{2},m_{2}|S_{1},m_{1})^{*},$
(18)

where $\bra{M_{1}}|f^{\dagger}_{\alpha}|\ket{M_{2}}$ is the reduced matrix
element with respect to spin symmetry. The Clebsch-Gordan coefficients needed
to evaluate the terms in Eqs. 16 and 17 are

$\displaystyle(\tfrac{1}{2},\pm\tfrac{1}{2};S,M|S+\tfrac{1}{2},M\pm\tfrac{1}{2})$
$\displaystyle=\sqrt{\frac{1}{2}\left(1\pm\frac{M\pm\frac{1}{2}}{S+\tfrac{1}{2}}\right)}$
(19)
$\displaystyle(\tfrac{1}{2},\pm\tfrac{1}{2};S,M|S-\tfrac{1}{2},M\pm\tfrac{1}{2})$
$\displaystyle=\mp\sqrt{\frac{1}{2}\left(1\mp\frac{M\pm\tfrac{1}{2}}{S+\tfrac{1}{2}}\right)}.$
(20)

We begin by evaluating $T^{N_{0}+1}_{S_{0}-\frac{1}{2}}(m_{0},\epsilon\alpha
i\leftarrow m_{0}^{\prime},\epsilon\alpha^{\prime}i^{\prime})$, which yields

$\displaystyle T^{N_{0}+1}_{S_{0}-\frac{1}{2}}(m_{0},\epsilon\alpha i$
$\displaystyle\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})=-\rho
V^{2}\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha^{\prime}}\ket{M_{0}}}{\epsilon_{M_{c}}-\epsilon_{M_{0}}}$
(21)
$\displaystyle\times\begin{cases}&-\frac{1}{S_{0}+\frac{1}{2}}\frac{1}{2}m_{0},\text{
if }i=i^{\prime}=\frac{1}{2},m_{0}=m_{0}^{\prime},\\\
&+\frac{1}{S_{0}+\frac{1}{2}}\frac{1}{2}m_{0},\text{ if
}i=i^{\prime}=-\frac{1}{2},m_{0}=m_{0}^{\prime},\\\
&-\frac{1}{S_{0}+\frac{1}{2}}\frac{1}{2}\sqrt{S_{0}(S_{0}+1)-m_{0}(m_{0}-1)},\text{
if }i=\frac{1}{2},i^{\prime}=-\frac{1}{2},m_{0}=m_{0}^{\prime}-1,\\\
&-\frac{1}{S_{0}+\frac{1}{2}}\frac{1}{2}\sqrt{S_{0}(S_{0}+1)-m_{0}(m_{0}+1)},\text{
if }i=-\frac{1}{2},i^{\prime}=\frac{1}{2},m_{0}=m_{0}^{\prime}+1,\end{cases}$
$\displaystyle=\rho
V^{2}\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha^{\prime}}\ket{M_{0}}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})(S_{0}+\frac{1}{2})}\frac{1}{2}S^{a}_{m_{0}m_{0}^{\prime}}\sigma^{a}_{ii^{\prime}}$

where we have ignored the contributions to the scattering potential. In order
to obtain the final expression we have made use of the relations

$\displaystyle\hat{S}^{a}\sigma^{a}_{ii^{\prime}}=\frac{1}{2}(\hat{S}^{+}\delta_{i,i+1}+\hat{S}^{-}\delta_{i,i-1})+\hat{S}^{z}\sigma^{z}_{ii^{\prime}},$
(22) $\displaystyle
S^{+}_{m_{0}m_{0}^{\prime}}=\sqrt{S_{0}(S_{0}+1)-m_{0}(m_{0}+1)}\delta_{m_{0},m_{0}^{\prime}+1}$
(23) $\displaystyle
S^{-}_{m_{0}^{\prime}m_{0}}=\sqrt{S_{0}(S_{0}+1)-m_{0}(m_{0}-1)}\delta_{m_{0},m_{0}^{\prime}-1}$
(24) $\displaystyle
S^{z}_{m_{0}m_{0}^{\prime}}=m_{0}\delta_{m_{0},m_{0}^{\prime}},$ (25)

where $\sigma$ are the Pauli matrices and the index $a$ labels the components
of the vectors, which are to be summed over when repeated. Following the same
procedure for the remaining $N_{c},S_{c}$ combinations, we obtain

$\displaystyle T^{N_{0}-1}_{S_{0}-\frac{1}{2}}(m_{0},\epsilon\alpha
i\leftarrow m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i)$
$\displaystyle=\rho
V^{2}\frac{\bra{M_{0}}|f^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|f_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})S_{0}}\frac{1}{2}S^{a}_{m_{0}m_{0}^{\prime}}\sigma^{a}_{ii^{\prime}},$
(26) $\displaystyle T^{N_{0}+1}_{S_{0}+\frac{1}{2}}(m_{0},\epsilon\alpha
i\leftarrow m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i)$
$\displaystyle=-\rho
V^{2}\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})(S_{0}+\frac{1}{2})}\frac{1}{2}S^{a}_{m_{0}m_{0}^{\prime}}\sigma^{a}_{ii^{\prime}}$
(27) $\displaystyle T^{N_{0}-1}_{S_{0}+\frac{1}{2}}(m_{0},\epsilon\alpha
i\leftarrow m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i)$
$\displaystyle=-\rho
V^{2}\frac{\bra{M_{0}}|f^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|f_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})(S_{0}+1)}\frac{1}{2}S^{a}_{m_{0}m_{0}^{\prime}}\sigma^{a}_{ii^{\prime}}$
(28)

Substituting these scattering amplitudes into Eq. 14 and introducing the
operator

$\hat{\psi}_{\alpha i}=\int d\epsilon\rho^{\frac{1}{2}}\hat{c}_{\epsilon\alpha
i},$ (29)

we can write the Kondo interaction as

$\hat{H}_{\text{K}}=\frac{J_{\alpha\alpha^{\prime}}}{2}\hat{S}^{a}\cdot\hat{\psi}^{\dagger}_{\alpha
i}\sigma^{a}_{ii^{\prime}}\hat{\psi}_{\alpha^{\prime}i^{\prime}},$ (30)

where

$\displaystyle\hat{S}^{a}=\sum_{m_{0},m_{0}^{\prime}}S^{a}_{m_{0},m_{0}^{\prime}}\ket{M_{0}m_{0}}\bra{M_{0}m_{0}^{\prime}},$
(31) $\displaystyle
J_{\alpha\alpha^{\prime}}=\sum_{M_{c}}\frac{V^{2}}{\epsilon_{M_{c}}-\epsilon_{M_{0}}}C_{\alpha\alpha^{\prime}}^{M_{c}}$
(32) $\displaystyle
C_{\alpha\alpha^{\prime}}^{M_{c}}=\begin{cases}\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha}|\ket{M_{0}}}{(S_{0}+\frac{1}{2})},&\text{if
}N_{c}=N_{0}+1,S_{c}=S_{0}-\frac{1}{2},\\\
\frac{\bra{M_{0}}|f^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|f_{\alpha}|\ket{M_{0}}^{*}}{S_{0}},&\text{if
}N_{c}=N_{0}-1,S_{c}=S_{0}-\frac{1}{2},\\\
-\frac{\bra{M_{0}}|f_{\alpha}|\ket{M_{c}}\bra{M_{c}}|f^{\dagger}_{\alpha}|\ket{M_{0}}}{(S_{0}+\frac{1}{2})},&\text{if
}N_{c}=N_{0}+1,S_{c}=S_{0}+\frac{1}{2},\\\
-\frac{\bra{M_{0}}|f^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|f_{\alpha}|\ket{M_{0}}^{*}}{(S_{0}+1)},&\text{if
}N_{c}=N_{0}-1,S_{c}=S_{0}+\frac{1}{2}.\end{cases}$ (33)

###  S4.2 Scaling equations for the coupling constants

In order to obtain the scaling equations for the coupling constants
$J_{\alpha\alpha^{\prime}}$ and $j_{a}$, we apply the poor man’s scaling
method [7] to the Kondo interaction derived in the previous section. By
integrating over second-order scattering processes where the intermediate
electron lies in the energy range $|\epsilon|\in[D-|\delta D|]$, we obtain the
differential scaling $\delta H_{\text{K}}$ of the Kondo interaction,

$\delta\hat{H}_{K}=\sum_{m_{0},m_{0}^{\prime}}\int_{-D}^{D}d\epsilon\int_{-D}^{D}d\epsilon^{\prime}\delta
T(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})\ket{M_{0},m_{0}}\bra{M_{0},m_{0}^{\prime}}\hat{c}^{\dagger}_{\epsilon\alpha
i}\hat{c}_{\epsilon^{\prime}\alpha^{\prime}i^{\prime}},$ (34)

where $\delta T(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$ is the sum over
second order scattering amplitudes involving the creation of an intermediate
particle ($p$) with or hole ($h$) with energy $|\epsilon|\in[D-\delta D,D]$:

$\delta T(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})=\delta
T^{(p)}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})+\delta
T^{(h)}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime}).$ (35)

The Feynman diagrams associated to these scattering amplitudes are given in
Fig. S11. Evaluating them, we obtain

$\displaystyle m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime\prime},\epsilon^{\prime\prime}\alpha^{\prime\prime},i^{\prime\prime}\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime}$ (36)
$\displaystyle\delta T^{(p)}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$
$\displaystyle=\frac{1}{2}\rho
J_{\alpha\alpha^{\prime\prime}}S^{a}_{m_{0}m_{0}^{\prime\prime}}\sigma^{a}_{ii^{\prime\prime}}\frac{|\delta
D|}{-D}\frac{1}{2}\rho
J_{\alpha^{\prime\prime}\alpha^{\prime}}S^{a}_{m_{0}^{\prime\prime}m_{0}^{\prime}}\sigma^{b}_{i^{\prime\prime}i^{\prime}}$
$\displaystyle=-\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}(S^{a}S^{b})_{m_{0}m_{0}^{\prime}}(\sigma^{a}\sigma^{b})_{ii^{\prime}}\frac{|\delta
D|}{D}$

$m_{0}^{\prime}$$\epsilon\alpha
i$$m_{0}$$\epsilon^{\prime}\alpha^{\prime}i^{\prime}$$\epsilon^{\prime\prime}\alpha^{\prime\prime}i^{\prime\prime}$$m_{0}^{\prime\prime}$

$m_{0}$$\epsilon\alpha
i$$m_{0}^{\prime}$$\epsilon^{\prime}\alpha^{\prime}i^{\prime}$$\epsilon^{\prime}\alpha^{\prime}i^{\prime}$$m_{0}^{\prime\prime}$

Figure S11: Feynman diagrams for virtual processes with an intermediate
electron (left) and hole (left).

for the scattering through intermediate states with an excited electron, where
$\rho|\delta D|$ is the number of particle states in the $[D-|\delta D|,D]$
energy range. For the scattering through an intermediate state with the
incoming electron plus an electron-hole pair, we have

$\displaystyle m_{0},\epsilon\alpha i\leftarrow m_{0},\epsilon\alpha
i,\epsilon^{\prime\prime}\alpha^{\prime\prime}i^{\prime\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime}\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime}$ (37)
$\displaystyle\delta T^{(h)}(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$
$\displaystyle=-\frac{1}{2}\rho
J_{\alpha^{\prime\prime}\alpha^{\prime}}S^{a}_{m_{0}m_{0}^{\prime\prime}}\sigma^{a}_{i^{\prime\prime}i^{\prime}}\frac{|\delta
D|}{-D}\frac{1}{2}\rho
J_{\alpha\alpha^{\prime\prime}}S^{b}_{m_{0}^{\prime\prime}m_{0}^{\prime}}\sigma^{b}_{ii^{\prime\prime}}$
$\displaystyle=\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}(S^{a}S^{b})_{m_{0}m_{0}^{\prime}}(\sigma^{b}\sigma^{a})_{ii^{\prime}}\frac{|\delta
D|}{D}.$

Thus, the total amplitude is

$\displaystyle\delta T(m_{0},\epsilon\alpha i\leftarrow
m_{0}^{\prime},\epsilon^{\prime}\alpha^{\prime}i^{\prime})$
$\displaystyle=-\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}(S^{a}S^{b})_{m_{0}m_{0}^{\prime}}[\sigma^{a},\sigma^{b}]_{ii^{\prime}}\frac{|\delta
D|}{D}$ (38)
$\displaystyle=-\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}(S^{a}S^{b})_{m_{0}m_{0}^{\prime}}2i\epsilon_{abc}\sigma^{c}_{ii^{\prime}}\frac{|\delta
D|}{D}$
$\displaystyle=-\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}i\epsilon_{abc}[S^{a},S^{b}]_{m_{0}m_{0}^{\prime}}\sigma^{c}_{ii^{\prime}}\frac{|\delta
D|}{D}$
$\displaystyle=-\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}i\epsilon_{abc}i\epsilon_{abd}S^{d}_{m_{0}m_{0}^{\prime}}\sigma^{c}_{ii^{\prime}}\frac{|\delta
D|}{D}$
$\displaystyle=\frac{1}{4}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}2\delta_{cd}S^{d}_{m_{0}m_{0}^{\prime}}\sigma^{c}_{ii^{\prime}}\frac{|\delta
D|}{D}$
$\displaystyle=\frac{1}{2}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}S^{c}_{m_{0}m_{0}^{\prime}}\sigma^{c}_{ii^{\prime}}\frac{|\delta
D|}{D}.$

Introducing this in Eq. 34, we get

$\displaystyle\delta\hat{H}_{\text{K}}$
$\displaystyle=\sum_{m_{0},m_{0}^{\prime}}\int_{-D}^{D}d\epsilon\int_{-D}^{D}d\epsilon^{\prime}\frac{1}{2}\rho^{2}\left[J^{2}\right]_{\alpha\alpha^{\prime}}\frac{|\delta
D|}{D}S^{a}_{m_{0}m_{0}^{\prime}}\hat{c}^{\dagger}_{\epsilon\alpha
i}\sigma^{a}_{ii^{\prime}}\hat{c}_{\epsilon^{\prime}\alpha^{\prime}i^{\prime}}.$
(39)
$\displaystyle=\frac{1}{2}\rho\left[J^{2}\right]_{\alpha\alpha^{\prime}}\frac{|\delta
D|}{D}\hat{S}^{c}\hat{\psi}^{\dagger}_{\alpha
i}\sigma^{c}_{ii^{\prime}}\hat{\psi}_{\alpha^{\prime}i^{\prime}}.$

By comparison with Eq. 34, we find

$\delta
J_{\alpha\alpha^{\prime}}=\rho[J^{2}]_{\alpha\alpha^{\prime}}\frac{|\delta
D|}{D}.$ (40)

In the diagonal channel basis, this reduces to the simple expression

$\delta j_{a}=\rho j_{a}^{2}\frac{|\delta D|}{D},$ (41)

Therefore, as the energy scale is lowered, the magnitude $|j_{a}|$ of the
coupling will grow in antiferromagnetically coupled channels ($j_{a}>0$) until
reaching the characteristic perturbative divergence and slowly decrease in the
ferromagnetically coupled channels ($j_{a}<0$).

###  S4.3 Analysis of the coupling for $\bf 1$ and $\bf 2$

In order to examine the effective coupling of $\bf 1$ and $\bf 2$ to the
substrate in detail, we diagonalize the Kondo coupling. This allows us to
represent the coupling in an eigenstate basis where the molecule behaves as a
local moment coupled to several independent channels, meaning that conduction
electrons in one such channel cannot be scattered to another channel. We call
these new channels ”eigenchannels” and we label them with the letter $a$,
keeping $\alpha$ for the channels that couple to each MO, which we call ”MO
channels”. We shall determine (i) whether the coupling to the eigenchannels is
AFM or FM, (ii) the relative strength of these eigenchannels, and (iii) the
projection of these eigenchannels onto the MO channels.

We choose the possible energy values $\epsilon_{-}$ and $\epsilon_{+}$ by
approximating them as $\epsilon_{-}=\epsilon_{-}^{0}-\mu$ and
$\epsilon_{+}=\epsilon_{+}^{0}+\mu$, where $\epsilon_{-}^{0}$ and
$\epsilon_{+}^{0}$ are the energies of the excited multiplets when the
molecule is free, and $\mu$ is the chemical potential introduced to model the
coupling to the substrate. This leaves us with $\mu$ as a parameter, which we
allow to vary in the range where $\epsilon_{-},\epsilon_{+}>0$, ensuring that
the ground state does not change.

For $\bf 1$, using the data in Table I of the main paper, which corresponds to
the channels $\alpha=2,3,4$ relevant for this molecule. we find that the
coupling constant matrix is

$\mathbf{J}=V^{2}\begin{pmatrix}-\frac{0.40}{\epsilon_{-}}-\frac{0.29}{\epsilon_{+}}&-\frac{0.06}{\epsilon_{+}}&0\\\
-\frac{0.06}{\epsilon_{-}}&\frac{1.13}{\epsilon_{-}}+\frac{0.51}{\epsilon_{+}}&0\\\
0&0&-\frac{0.04}{\epsilon_{-}}-\frac{0.4}{\epsilon_{+}}\end{pmatrix}$ (42)

We index the matrix elements $\alpha=2,3,4$ so that it matches the MO
labeling. We immediately see that the $\alpha=4$ MO channel is an eigenchannel
with a FM coupling to the molecule. We discard it, as it is irrelevant at low
temperatures. In the upper block we find the eigenvalues $j_{2}$ and $j_{3}$,
which are shown Fig. S12. Using the free molecule energies
$\epsilon^{0}_{-}=-80$ meV and $\epsilon^{0}_{+}=490$ meV, we find that the
eigenvalue $j_{2}$ is negative in the whole $\mu$ window, so the coupling to
the $a=2$ eigenchannel is FM; $j_{3}$ is positive, so the coupling to $a=3$ is
AFM. In Fig. S13 we show the mixing of the MO channels in the eigenchannels in
terms of the quantities $|U_{\alpha a}|^{2}$, where $U$ are the diagonalizing
matrices that fulfil $\mathbf{j}=U^{\dagger}\mathbf{J}U$. As can be seen in
the graphs, the eigenchannel 2 (3) couples almost exclusively to the MO
channel 2 (3), hence the labeling of the eigenchannel. The coupling to the
secondary channel is of the order of $|U_{\alpha a}|\sim 10^{-3}$ for both
eigenchannels.

Figure S12: Eigenvalues $j^{\prime}_{2}=V^{-2}j_{2}$ and
$j^{\prime}_{3}=V^{-2}j_{3}$ for $\bf 1$.

(a)

(b)

Figure S13: For $\bf 1$, MO channel mixings $|U_{\alpha a}|^{2}$ for the
eigenchannels $a=2$ (a) and $a=3$ (b) and the primary (thick) and secondary
(dashed) MO channels $\alpha=2,3$.

For $\bf 2$, the coupling constant matrix calculated from the data on Table I
of the main paper is given by

$\mathbf{J}=V|^{2}\begin{pmatrix}\frac{1.24}{\epsilon_{-}}+\frac{1.24}{\epsilon_{+}}&\frac{0.1}{\epsilon_{-}}+\frac{0.83}{\epsilon_{+}}\\\
\frac{0.1}{\epsilon_{-}}+\frac{0.83}{\epsilon_{+}}&\frac{0.57}{\epsilon_{-}}+\frac{0.55}{\epsilon_{+}}\end{pmatrix}.$
(43)

Using the free molecule energies $\epsilon^{0}_{-}=-45$ meV and
$\epsilon_{+}^{0}=500$ meV, we find the eigenvalues $j_{3}$ and $j_{4}$ shown
in Fig. S14. Except for the near vicinity of the lower limit $\mu=-500$ meV,
which is well into the mixed valence regime, both eigenvalues are positive,
which means both eigenchannels develop an AFM coupling with the molecular
spin. The mixing of the MO channels in the eigenchannels is shown in Fig. S15.
The mixing is larger for $\bf 2$ than for $\bf 1$ as a result of the
relatively larger off-diagonal entries of $\mathbf{J}$, but it is still
possible to identify a primary MO channel for each eigenchannel in the whole
$\mu$ region: the eigenchannel 3 (4) has its largest component in the MO
channel 3 (4).

Figure S14: Eigenvalues $j^{\prime}_{3}=V^{-2}j_{3}$ and
$j_{4}^{\prime}=V^{-2}j_{4}$ for $\bf 2$.

(a)

(b)

Figure S15: For $\bf 2$, MO channel mixings $|U_{\alpha a}|^{2}$ for the
eigenchannels $a=3$ (a) and $a=4$ (b) and the primary (thick) and secondary
(dashed) MO channels $\alpha=3,4$.

###  S4.4 Relationship between the NTOs and the eigenchannels

Using the orbital indices $\alpha$, we can write the transition matrix
elements from which we compute the NTOs as

$\displaystyle T^{m_{0}}_{\alpha\alpha^{\prime}}$
$\displaystyle=\bra{M_{0},m_{0}-1}\hat{f}^{\dagger}_{\alpha\downarrow}\hat{f}_{\alpha^{\prime}\uparrow}\ket{M_{0}m_{0}}.$
(44)
$\displaystyle=\bra{M_{0}}|f^{\dagger}_{\alpha}f_{\alpha^{\prime}}|\ket{M_{0}}(1,-1;S_{0},m_{0}|S_{0},m_{0}-1),$

where we have applied the Wigner-Eckart theorem to decompose the transition
matrix element into a reduced part and a Clebsch-Gordan coefficient, taking
into account that $f^{\dagger}_{\alpha\downarrow}f_{\alpha^{\prime}\uparrow}$
is a tensor operator with total spin (irrep) $1$ and spin projection (partner
number) $-1$. Since the transition matrix elements in complete and reduced
form for a given $m_{0}$ differ only in a constant Clebsch-Gordan factor, they
have the same eigenvectors and therefore they give us the same NTOs.

We can also expand the transition matrix elements by introducing the identity
operator $\mathbf{I}=\sum_{M,m}\ket{M,m}\bra{M,m}$, which results in

$\displaystyle T^{m_{0}}_{\alpha\alpha^{\prime}}$
$\displaystyle=\sum_{M,m}\bra{M_{0},m_{0}-1}\hat{f}^{\dagger}_{\alpha\downarrow}\ket{Mm}\bra{Mm}\hat{f}_{\alpha^{\prime}\uparrow}\ket{M_{0}m_{0}}$
(45)
$\displaystyle=\sum_{M_{c}}\bra{M_{0},m_{0}-1}\hat{f}^{\dagger}_{\alpha\downarrow}\ket{M_{c},m_{0}-\tfrac{1}{2}}\bra{M_{c},m_{0}-\tfrac{1}{2}}\hat{f}_{\alpha^{\prime}\uparrow}\ket{M_{0}m_{0}}$
$\displaystyle=\sum_{M_{c}}\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}(\tfrac{1}{2},+\tfrac{1}{2};S_{c},m_{0}-\tfrac{1}{2}|S_{0},m_{0})(\tfrac{1}{2},-\tfrac{1}{2};S_{c},m_{0}-\tfrac{1}{2}|S_{0},m_{0}-1)$
$\displaystyle=S^{-}_{m_{0}-1,m_{0}}\sum_{M_{c}}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}}{2F(N_{0}-1,S_{c},S_{0})}$
$\displaystyle=\bra{S_{0}}|S^{-}|\ket{S_{0}}(1,-1;S_{0},m_{0}|S_{0},m_{0}-1)\sum_{M_{c}}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}}{2F(N_{0}-1,S_{c},S_{0})},$

where

$F(N_{c},S_{c},S_{0})=\begin{cases}(S_{0}+\frac{1}{2}),&\text{if
}N_{c}=N_{0}+1,S_{c}=S_{0}-\frac{1}{2},\\\ S_{0},&\text{if
}N_{c}=N_{0}-1,S_{c}=S_{0}-\frac{1}{2},\\\ -(S_{0}+\frac{1}{2}),&\text{if
}N_{c}=N_{0}+1,S_{c}=S_{0}+\frac{1}{2},\\\ -(S_{0}+1),&\text{if
}N_{c}=N_{0}-1,S_{c}=S_{0}+\frac{1}{2}\end{cases}$ (46)

contains the spin factor and the sign appearing in Eq. 33. The factors with
$N_{c}=N_{0}+1$ appear when we compute the matrix elements with the creation
and annihilation operators commuted,

$T^{m_{0}}_{\alpha\alpha^{\prime}}=-\bra{M_{0},m_{0}-1}f_{\alpha^{\prime}\uparrow}f^{\dagger}_{\alpha\downarrow}\ket{M_{0}m_{0}}.$
(47)

Putting together Eqs. 44 and 45, we arrive at the equivalent expressions

$\bra{M_{0}}|f^{\dagger}_{\alpha}f_{\alpha^{\prime}}|\ket{M_{0}}=\bra{S_{0}}|S^{-}|\ket{S_{0}}\sum_{M_{c}}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}}{2F(N_{0}+1,S_{c},S_{0})}$
(48)

and

$-\bra{M_{0}}|f_{\alpha}f^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}^{*}=\bra{S_{0}}|S^{-}|\ket{S_{0}}\sum_{M_{c}}\frac{\bra{M_{0}}|\hat{f}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|\hat{f}^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{2F(N_{0}-1,S_{c},S_{0})}$
(49)

for the reduced matrix elements computed by introducing intermediate excited
multiplets with positive and negative charge, respectively. The reduced matrix
element $\bra{S_{0}}|S^{-}|\ket{S_{0}}$ is constant, so we incorporate it in a
new transition matrix $\mathbf{\mathcal{T}}$ that contains only multiplet
quantum numbers and gives us the same NTOs as the original transition matrix
$\mathbf{T}^{m_{0}}$:

$\displaystyle\mathcal{T}_{\alpha\alpha^{\prime}}$
$\displaystyle:=\frac{2\bra{M_{0}}|f^{\dagger}_{\alpha}f_{\alpha^{\prime}}|\ket{M_{0}}}{\bra{S_{0}}|S^{-}|\ket{S_{0}}}=\sum_{M_{c}}^{N_{c}=N_{0}+1}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}}{F(N_{0}+1,S_{c},S_{0})}$
(50)
$\displaystyle=-\frac{2\bra{M_{0}}|f^{\dagger}_{\alpha}f_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{\bra{S_{0}}|S^{-}|\ket{S_{0}}}=\sum_{M_{c}}^{N_{c}=N_{0}-1}\frac{\bra{M_{0}}|\hat{f}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|\hat{f}^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{F(N_{0}+1,S_{c},S_{0})}.$

To establish a connection between the NTOs and the eigenchannels, we write the
coupling constant matrix as

$\mathcal{J}_{\alpha\alpha^{\prime}}=V^{-2}J_{\alpha\alpha^{\prime}}=\sum_{M_{c}}^{N_{c}=N_{0}-1}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})F(N_{0}-1,S_{c},S_{0})}+\sum_{M_{c}}^{N_{c}=N_{0}+1}\frac{\bra{M_{0}}|\hat{f}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}}{(\epsilon_{M_{c}}-\epsilon_{M_{0}})F(N_{0}+1,S_{c},S_{0})}.$
(51)

We observe two main differences between $\mathbf{\mathcal{T}}$ and
$\mathbf{\mathcal{J}}$: (i) $\mathbf{\mathcal{J}}$ contains energy factors
$(\epsilon_{M_{c}}-\epsilon_{M_{0}})^{-1}$ and (ii) $\mathbf{\mathcal{T}}$ can
be equivalently obtained from projecting onto intermediate states with either
$N_{c}=N_{0}-1$ or $N_{c}=N_{0}+1$, while $\mathbf{\mathcal{J}}$ includes both
types of terms. The similarity of the matrices is best observed in the case
where all the multiplets with $N_{c}=N_{0}\pm 1$ particles are degenerate and
have an energy $E=(\epsilon_{M_{c}}-\epsilon_{M_{0}})$ which yields

$\displaystyle\mathcal{J}_{\alpha\alpha^{\prime}}$
$\displaystyle=\frac{1}{E}\left[\sum_{M_{c}}^{N_{c}=N_{0}-1}\frac{\bra{M_{0}}|\hat{f}^{\dagger}_{\alpha}|\ket{M_{c}}^{*}\bra{M_{c}}|\hat{f}_{\alpha^{\prime}}|\ket{M_{0}}^{*}}{F(N_{0}-1,S_{c},S_{0})}+\sum_{M_{c}}^{N_{c}=N_{0}+1}\frac{\bra{M_{0}}|\hat{f}_{\alpha}|\ket{M_{c}}\bra{M_{c}}|\hat{f}^{\dagger}_{\alpha^{\prime}}|\ket{M_{0}}}{F(N_{0}+1,S_{c},S_{0})}\right]$
(52)
$\displaystyle=\frac{1}{E}\left[\mathcal{T}_{\alpha\alpha^{\prime}}^{*}+\mathcal{T}_{\alpha\alpha^{\prime}}^{*}\right]$
$\displaystyle=\frac{2}{E}\mathcal{T}_{aa}^{*}.$

This shows that in the limit where all the excited multiplets are degenerate,
both matrices contain equivalent information and therefore the same unitary
transformation gives us the NTOs and the eigenchannels when applied to the
orbital and conduction states, respectively. Furthermore, in Figs. S7 and S8
we show the numerical comparison between these orbitals and corresponding
dI/dV maps. In the case of 1, the agreement is good for the NTOs and the Kondo
orbitals. The main difference in the dI/dV maps is the brighter contrast in
the central part presented by the Kondo orbitals, which fits better to the
experimental evidence (see Fig. 1 in the main text). In the case of 2, we
found no apparent difference between both dI/dV maps.

##  S5 NRG calculations

###  S5.1 Method

For the NRG calculations we use a discretized Anderson Hamiltonian with the
conduction and hibridization terms given by

$\displaystyle\hat{H}_{\text{con}}=\sum_{n=1}^{\infty}\tau_{n,\alpha}\hat{c}^{\dagger}_{n,\alpha\sigma}\hat{c}_{n+1,\alpha\sigma}+\text{h.c.},$
(53)
$\displaystyle\hat{H}_{\text{hyb}}=\sum_{M_{i},M_{f}}\sum_{m_{i},m_{f}}\sqrt{2\Gamma/\pi}\bra{M_{f},m_{f}}\hat{f}^{\dagger}_{\alpha\sigma}\ket{M_{i},m_{i}}\ket{M_{f},m_{f}}\bra{M_{i},m_{i}}\hat{c}_{1,\alpha\sigma}+\text{h.c.},$

where $\Gamma=\pi\rho V^{2}$ is the constant hybridization; repeated indices
are to be summed over. The term $\hat{H}_{\text{con}}$ represents the
conduction channels as a semi-infinite tight-binding chains, which are
obtained using the interleaved discretization procedure [8]. Electrons in the
$n$-th site of this chain are annihilated (created) by the operators
$\hat{c}^{(\dagger)}_{n,\alpha\sigma}$. Neighboring sites are coupled by the
hopping constants $\tau_{n,\alpha}$. The hybridization term
$\hat{H}_{\text{hyb}}$ couples the molecule to the first site in the chain.
For all the calculations we use a half-bandwidth of $D=1\text{eV}$.

We use the PointGroupNRG code [9] to compute the zero-temperature spectral
functions and the thermodynamic curves (given below) following the procedure
presented in Ref. [10]. The spectral function, shown in Fig. 4 of the main
paper, is defined as

$A_{\alpha}(\epsilon)=\sum_{\sigma}\sum_{E}\left[|\bra{E}\hat{f}^{\dagger}_{\alpha\sigma}\ket{G}|^{2}\delta(\epsilon-(\epsilon_{E}-\epsilon_{G}))+|\bra{E}\hat{f}_{\alpha\sigma}\ket{G}|^{2}\delta(\omega+(\epsilon-\epsilon_{G}))\right],$
(54)

where $\ket{E}$ is an excited state, $\ket{G}$ is the ground state. The
calculation is performed by choosing for each NRG iteration $n$ a value
$\epsilon_{n}=2\omega_{n}$ for which the spectral function is computed, where
$\omega_{n}$ is the energy scale corresponding to the $n$-th step. The delta
peaks are broadened using a Gaussian kernel with

$\delta(\epsilon)\rightarrow\frac{1}{\eta_{n}\sqrt{\pi}}\exp\left[-\left(\epsilon/\eta_{n}\right)^{2}\right],$
(55)

where $\eta_{n}=\omega_{n}$ is chosen as the broadening factor. The results
are averaged over 4 values of the twisting parameter $z$ in order to improve
the smoothness of the curves. For the thermodynamic functions, smooth curves
are obtained by averaging over even and odd NRG steps. For both spectral and
thermodynamic calculations, we use a discretization parameter $\Lambda=3$ and
a cutoff of 1000 multiplets.

###  S5.2 Thermodynamic curves and Kondo temperatures

To provide further insight into the Kondo effect, we show the thermodynamic
curves for $\bf 1$ and $\bf 2$. These have been computed using the same
parameters $\pi\rho|V|^{2}$ and $\mu$ used in the calculation of the spectral
functions.

(a)

(b)

Figure S16: Calculated thermodynamic functions for molecules $\bf 1$ (a) and
$\bf 2$ (b): entropy (top panel) and magnetic susceptibility (bottom panel).
The parameters are $\Gamma=0.1$ eV for $\bf 1$, $\Gamma=0.15$ eV for $\bf 2$,
and $\mu=-0.2$ eV for both. The arrow in (a) indicates the Kondo temperature
$k_{B}T_{K}=4\times 10^{-4}$ eV. The arrow in (b) indicates the separation
between the two screening stages.

The calculated entropy and magnetic susceptibility curves for $\bf 1$ are
shown in Fig. 16(a). These curves show that the system develops a local moment
around $k_{B}T/D\sim 10^{-2}$, in the region where the entropy and magnetic
susceptibility are close to the local moment fixed point values
$S/k_{B}=\ln(2)$ and $k_{B}T\chi/(g\mu_{B})^{2}=1/4$. At lower temperatures
this local moment is screened, with the full screening corresponding to the
values $S=\chi=0$. The Kondo temperature $T_{K}$ obtained from the magnetic
susceptibility, calculated as the point where $d\chi/dT|_{T=T_{K}}=0$, is
$T_{K}\approx 5K$. This is close to the value $T_{K}\approx 8K$ calculated
from the full width at half-maximum $\Delta$ of the peak in the MO 3 spectral
function, $k_{B}T_{K}=\Delta$. Both values agree with the Kondo temperature
obtained by fitting to the experimental results, which lies in the
$T_{K}\in[1,10]K$ range.

Thermodynamic curves for $\bf 2$ are shown in Fig. 16(b). As the temperature
is decreased, the magnetic susceptibility curve starts approaching the $3/4$
value of the spin-$1$ local moment fixed point, which has an associated
entropy $S/k_{B}=\ln(3)$. The hybridization rapidly causes this tendency to
shift, screening the partially formed local moment. Two stages can be observed
in the screening, separated by the middle inflexion point in the magnetic
susceptibility curve, in agreement with the anistropy in the $j_{a}$
eigenvalues found in section S4.3. This is a characteristic feature of the
spin-$1$ local moment asisotropic Kondo effect[11].

It must be noted that, for both $\bf 1$ and $\bf 2$, the large hybridization
$\Gamma\sim\epsilon_{\pm}$ causes the system to evolve to the screened regime
without fully developing the local moment. If the hybridization were smaller,
we would observe a plateau close to the local moment fixed point values for a
large portion of the temperature range. In principle, this would mean that the
effective model is not a good approximation. Nevertheless, the results for the
spectral functions agree with the analysis based on the effective Kondo model
with respect to the nature of the magnetic coupling to the channels 2 and 3
(FM or AFM). Thus, the rationalization of the Kondo peaks based on the Kondo
model can be considered as qualitatively correct.

## References

  * [1] A. Cahlik, Kondo fit, https://github.com/cahlikales/Kondo_fit.
  * Frota [1992] H. Frota, Physical Review B 45, 1096 (1992).
  * Gruber  _et al._ [2018] M. Gruber, A. Weismann, and R. Berndt, Journal of Physics: Condensed Matter 30, 424001 (2018).
  * Kondo [1964] J. Kondo, Progress of theoretical physics 32, 37 (1964).
  * Martin [2003] R. L. Martin, The Journal of Chemical Physics 118, 4775 (2003).
  * Krejčí  _et al._ [2017] O. Krejčí, P. Hapala, M. Ondráček, and P. Jelínek, Physical Review B 95, 10.1103/physrevb.95.045407 (2017).
  * Anderson [1970] P. W. Anderson, J. Phys. C: Solid St. Phys 3 (1970).
  * Yoshida  _et al._ [1990] M. Yoshida, M. A. Whitaker, and L. N. Oliveira, Phys. Rev. B 41, 9403 (1990).
  * Calvo-Fernández  _et al._ [2023] A. Calvo-Fernández, M. Blanco-Rey, and A. Eiguren, The PointGroupNRG code for numerical renormalization group calculations with discrete point-group symmetries (2023), arXiv:2307.03658 [cond-mat.str-el] .
  * Bulla  _et al._ [2008] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
  * Jayaprakash  _et al._ [1981] C. Jayaprakash, H. R. Krishna-murthy, and J. W. Wilkins, Phys. Rev. Lett. 47, 737 (1981).