#  Supplementary Material
Many-body quantum interference route to the two-channel Kondo effect:
Inverse design for molecular junctions

Sudeshna Sen  Department of Physics, IIT(ISM) Dhanbad, Dhanbad-826004,
Jharkhand, India  Andrew K. Mitchell  School of Physics, University College
Dublin, Belfield, Dublin 4, Ireland  Centre for Quantum Engineering, Science,
and Technology, University College Dublin, Ireland

##  S-I Microscopic models of
molecular junctions

We describe molecular junctions (and coupled quantum dot devices constituting
artificial molecular junctions) in terms of generalized quantum impurity
models of the type,

$\hat{H}=\hat{H}_{\rm mol}+\hat{H}_{\rm leads}+\hat{H}_{\rm hyb}+\hat{H}_{\rm
gate}$ (S-1)

where the part of the Hamiltonian $\hat{H}_{\rm mol}$ describing the isolated
molecule is given by an extended Hubbard model,

$\hat{H}_{\rm
mol}=\sum_{\sigma=\uparrow,\downarrow}\>\sum_{m,n}t_{mn}^{\phantom{\dagger}}d_{m\sigma}^{\dagger}d_{n\sigma}^{\phantom{\dagger}}+\tfrac{1}{2}\sum_{m,n}U_{mn}\hat{n}_{m}\hat{n}_{n}$
(S-2)

as given in Eq. 1 of the main text. Here $d_{m\sigma}^{(\dagger)}$ annihilates
(creates) an electron on molecule orbital $m$ with spin $\sigma$ and
$\hat{n}_{m}=\sum_{\sigma}\hat{n}_{m\sigma}$ is the total number operator for
orbital $m$, where
$\hat{n}_{m\sigma}=d_{m\sigma}^{\dagger}d_{m\sigma}^{\phantom{\dagger}}$. The
single-particle tunneling matrix elements between orbitals $m\neq n$ are given
by $t_{mn}$ ($=t_{nm}^{*}$) whereas the onsite potentials are
$\epsilon_{m}=t_{mm}+\tfrac{1}{2}U_{mm}$. The term $U_{mn}$ ($=U_{nm}$)
embodies electronic interactions on the molecule. The standard local Hubbard
interaction term $U_{m}\hat{n}_{m\uparrow}\hat{n}_{m\downarrow}$ between up
and down spin electrons on a given site $m$ is already included in the
diagonal part of the second term of Eq. S-2, with $U_{m}\equiv U_{mm}$. The
off-diagonal parts $U_{mn}$ for $m\neq n$ correspond to capacitive (Coulomb)
repulsion between different sites. In principle one can also add spin-orbit
coupling terms, anisotropies, and exchange interactions. The formalism
presented below can be straightforwardly adapted to treat any $\hat{H}_{\rm
mol}$.

The source ($\alpha=s$) and drain ($\alpha=d$) leads are taken to be thermal
reservoirs of non-interacting conduction electrons, $\hat{H}_{\rm
leads}=\sum_{\alpha}\hat{H}_{\rm leads}^{\alpha}$ with

$\hat{H}_{\rm
leads}^{\alpha}=\sum_{k,\sigma}\epsilon_{k}^{\phantom{\dagger}}c_{\alpha
k\sigma}^{\dagger}c_{\alpha k\sigma}^{\phantom{\dagger}}$ (S-3)

where for simplicity we take the dispersion $\epsilon_{k}$ to be independent
of spin $\sigma$ and equivalent for both leads $\alpha$ (although again this
can be straightforwardly generalized without affecting the following).

The molecule frontier orbital $d_{r_{\alpha}\sigma}$ couples to the conduction
electrons of lead $\alpha$ according to the hybridization term

$\hat{H}_{\rm hyb}=\sum_{\alpha,k,\sigma}V_{\alpha
k}^{\phantom{\dagger}}(d_{r_{\alpha}\sigma}^{\dagger}c_{\alpha
k\sigma}^{\phantom{\dagger}}+{\rm H.c.})\;.$ (S-4)

We now define the local lead orbitals at the molecule position as
$c_{\alpha\sigma}=\tfrac{1}{V_{\alpha}}\sum_{k}V_{\alpha k}c_{\alpha\sigma k}$
with $V_{\alpha}^{2}=\sum_{k}V_{\alpha k}^{2}$ such that $\hat{H}_{\rm
hyb}=\sum_{\alpha,\sigma}V_{\alpha}(d_{r_{\alpha}\sigma}^{\dagger}c_{\alpha\sigma}^{\phantom{\dagger}}+{\rm
H.c.})$.

The free lead density of states at the junction position (taken to the same
for both leads and independent of spin) is $\rho(\omega)=-\tfrac{1}{\pi}{\rm
Im}~{}G_{\rm lead}^{0}(\omega)$ where $G_{\rm lead}^{0}(\omega)=\langle\langle
c_{\alpha\sigma}^{\phantom{\dagger}};c_{\alpha\sigma}^{\dagger}\rangle\rangle$
is the free Green’s function for the local lead orbital $c_{\alpha\sigma}$ to
which the molecule couples. For simplicity we use
$\rho(\omega)=\rho_{0}\theta(D-|\omega|)$ such that the lead density of states
is constant $\rho_{0}=1/2D$ in a band of half-width $D$. Throughout this work
we set $D=1$.

The gate voltage $V_{g}$ allows to tune the number of electrons on the
molecule via $\hat{H}_{\rm gate}=V_{g}\sum_{m}\hat{n}_{m}$.

###  S-I.1 Interacting nanowire (Hubbard chains)

In the main text we discuss a particular limiting case of Eq. S-2
corresponding to an $M$-site 1d Hubbard chain,

$\hskip-5.69046pt\hat{H}_{\rm
mol}=t\sum_{m=1}^{M-1}\sum_{\sigma}\left(d_{m\sigma}^{\dagger}d_{m+1\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)+\sum_{m=1}^{M}\left(\epsilon\hat{n}_{m}+U\hat{n}_{m\uparrow}\hat{n}_{m\downarrow}\right)$
(S-5)

which describes an interacting nanowire (for example the $\pi$-system of a
conjugated hydrocarbon polymer [1]). We take constant nearest-neighbour
tunneling $t$, constant onsite potentials $\epsilon$ and constant local
Hubbard repulsion $U$. The source and drain leads are connected to the ends of
the nanowire at sites $r_{s}=1$ and $r_{d}=M$, such that

$H_{\rm
hyb}=\sum_{\sigma}\left(V_{s}^{\phantom{\dagger}}c_{s\sigma}^{\dagger}d_{1\sigma}^{\phantom{\dagger}}+V_{d}^{\phantom{\dagger}}c_{d\sigma}^{\dagger}d_{M\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\;.$ (S-6)

The low-energy effective model for this system with odd chain length $M$ is
discussed below in Sec. S-II.2.

###  S-I.2 Symmetric 5-site cluster

Introducing additional tunnel-couplings between different sites of the above
1d Hubbard chain clearly produces greater complexity and, in particular, more
possibilities for QI effects.

However, if we want exact $ph$ symmetry then the molecule connectivity must be
compatible with a bipartite graph – which in the present context simply means
that there should not be any odd-membered loops. The simplest system with an
odd total number of sites with $sd$-symmetry to contain multiple paths through
the molecule, while also realizing exact particle hole symmetry at $V_{g}=0$,
is the $M=5$ chain with next-next-nearest-neighbour tunneling connecting sites
$1\leftrightarrow 4$ and $2\leftrightarrow 5$,

$\begin{split}H_{\rm
mol}=&\sum_{m=1}^{5}\left(-\tfrac{1}{2}U\hat{n}_{m}+U\hat{n}_{m\uparrow}\hat{n}_{m\downarrow}\right)\\\
&+t\sum_{m=1}^{4}\sum_{\sigma}\left(d_{m\sigma}^{\dagger}d_{m+1\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\\\
&+t^{\prime}\sum_{\sigma}\left(d_{1\sigma}^{\dagger}d_{4\sigma}^{\phantom{\dagger}}+d_{2\sigma}^{\dagger}d_{5\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\end{split}$ (S-7)

which is equivalent to Eq. 3 of the main text. We attach leads to sites $1$
and $5$ via Eq. S-6. We discuss the low-energy effective model for this system
below in Sec. S-II.3

###  S-I.3 Low-symmetry 4-site cluster

Finally we consider a general low-symmetry structure with just $M=4$
interacting sites. We simplify to local Hubbard interactions:

$\hat{H}_{\rm
mol}=\sum_{m,n}\sum_{\sigma}t_{mn}^{\phantom{\dagger}}d_{m\sigma}^{\dagger}d_{n\sigma}^{\phantom{\dagger}}+\sum_{m}U_{m}\hat{n}_{m\uparrow}\hat{n}_{m\downarrow}$
(S-8)

Leads are attached to sites $1$ and $4$ via Eq. S-6. We discuss gradient
descent optimization of this system to satisfy the QI-2CK conditions in Sec.
S-VII.1.

The parameters used in Fig. 5 of the main text are $t_{11}\simeq 0.697$,
$t_{22}\simeq 0.642$, $t_{33}\simeq 0.480$, $t_{44}\simeq 0.253$,
$V_{g}=0.42$, $t_{12}\simeq 0.102$, $t_{13}\simeq 0.053$, $t_{14}\simeq 0.1$,
$t_{23}\simeq 0.109$, $t_{24}\simeq 0.046$, $t_{34}\simeq 0.169$.

##  S-II Generalized numerical Schrieffer-Wolff transformation

The Schrieffer-Wolff transformation (SWT) was first applied to the mapping
between the Anderson impurity model (AIM) and the Kondo model [2]. For the AIM
it involves projecting onto the impurity spin states by eliminating virtual
excitations and folding in their contribution to effective model parameters.
The SWT is in essence a single-step renormalisation group (RG) procedure.
Typically the mapping is done perturbatively to second order in the impurity-
bath hybridization, which is then equivalent to performing second-order
Brillouin-Wigner perturbation theory (BWPT) [3]. For the AIM-to-Kondo mapping
this can be done straightforwardly and analytic expressions for the Kondo
model parameters are easily obtained.

However, in the present context of complex molecular junctions, an exact
analytic treatment of the SWT is in general not possible, and we must resort
to a numerical evaluation of the effective model parameters [4, 5]. Such a
generalized numerical SWT is an ideal tool for efficiently describing many-
body quantum interference (QI) effects. It can also be integrated with machine
learning (ML) algorithms for inverse design, by identifying candidate
molecular junctions that have specific desired effective model properties, as
shown below.

Before proceeding we make a comment on the validity of the SWT. The RG
character of quantum impurity systems means that a second-order SWT typically
yields the correct operator structure of the low-energy effective model, even
if the precise value of the derived parameters is not exact [6]. However, even
infinite-order SWT [7] neglects retardation effects from the conduction
electrons resulting in bandwidth renormalization [8], and more sophisticated
non-perturbative techniques such as the numerical renormalization group [9]
(NRG) are ultimately required. On the other hand, in this work we demonstrate
that QI effects are very well captured by SWT (even quantitatively) when
compared with NRG.

###  S-II.1 Formalism

In this work we are interested in deriving generalized two-channel Kondo (2CK)
models of the type Eq. 2 from the main paper. Such a description arises in
two-lead molecular junctions as the low-energy effective model when the
isolated molecule ground state manifold comprises a unique spin-doublet [4,
5]. The structure of the effective model is constrained by $U(1)$ charge
conservation and $SU(2)$ spin symmetry. The spin-doublet ground state of the
molecule implies that the ground state hosts an odd number $N$ of electrons,
and that tunneling to and from the leads involves $N\pm 1$ even-electron
sectors with integer spin.

With $\hat{H}^{\prime}_{\rm mol}=\hat{H}_{\rm mol}+\hat{H}_{\rm gate}$ taken
to include gate voltage effects on the isolated molecule, we define a set of
exact eigenstates $\\{|\psi^{n;S^{z}}_{j}\rangle\\}$ satisfying the
Schrödinger equation $\hat{H}^{\prime}_{\rm
mol}|\psi^{n;S^{z}}_{j}\rangle=E_{j}^{n;S^{z}}|\psi^{n;S^{z}}_{j}\rangle$,
where $n$ is the number of electrons on the isolated molecule, $S^{z}$ is the
total molecule spin projection, and $j$ is an index distinguishing molecule
states with the same $n$ and $S^{z}$. The ground spin doublet states are
specified by $|\psi^{N;\sigma}_{0}\rangle$ with $\sigma=\uparrow,\downarrow$,
and the corresponding ground state energy is $E_{\rm gs}\equiv
E_{0}^{N;\uparrow}=E_{0}^{N;\downarrow}$.

Eigenstates of $\hat{H}^{\prime}_{\rm mol}$ can be expressed in the
(computational) many-body product-state basis
$\\{|\phi^{n;S^{z}}_{k}\rangle\\}$, which are labelled by the same quantum
numbers, viz:

$|\psi^{n;S^{z}}_{j}\rangle=\sum_{k}U^{n;S^{z}}_{k,j}|\phi^{n;S^{z}}_{k}\rangle\;.$
(S-9)

For notational convenience we define product basis vectors
$\vec{\phi}^{n,S^{z}}=\left(|\phi^{n;S^{z}}_{0}\rangle,~{}|\phi^{n;S^{z}}_{1}\rangle,~{}|\phi^{n;S^{z}}_{2}\rangle,~{}...\right)$
and eigenbasis vectors
$\vec{\psi}^{n,S^{z}}=\left(|\psi^{n;S^{z}}_{0}\rangle,~{}|\psi^{n;S^{z}}_{1}\rangle,~{}|\psi^{n;S^{z}}_{2}\rangle,~{}...\right)$.
These are related by Eq. S-9 as
$\vec{\psi}^{n,S^{z}}=\vec{\phi}^{n,S^{z}}\bm{U}^{n;S^{z}}$ in terms of the
transformation matrices $\bm{U}^{n;S^{z}}$ with elements $U^{n;S^{z}}_{k,j}$.
The Hamiltonian matrix for the isolated but interacting molecule in the
computational product basis of a given quantum number subspace can then be
written as
$\bm{H}_{\phi}^{n;S^{z}}=\left(\vec{\phi}^{n,S^{z}}\right)^{\dagger}\hat{H}^{\prime}_{\rm
mol}~{}\vec{\phi}^{n,S^{z}}$, whereas in the eigenbasis we have
$\bm{H}_{\psi}^{n;S^{z}}=\left(\vec{\psi}^{n,S^{z}}\right)^{\dagger}\hat{H}^{\prime}_{\rm
mol}~{}\vec{\psi}^{n,S^{z}}=\left(\bm{U}^{n;S^{z}}\right)^{\dagger}\bm{H}_{\phi}^{n;S^{z}}\bm{U}^{n;S^{z}}\equiv{\rm
diag}\left(E_{0}^{n,S^{z}},E_{1}^{n,S^{z}},E_{2}^{n,S^{z}},...\right)$. These
expressions hold separately in each subspace.

Finally, we will find it useful to define the vector of expansion coefficients
for the $N$-electron ground state with spin $\sigma$ as
$\vec{U}^{N;\sigma}_{\rm
gs}=\left(U^{N;\sigma}_{0,0},~{}U^{N;\sigma}_{1,0},~{}U^{N;\sigma}_{2,0},~{}...\right)^{T}$
such that
$|\psi^{N;\sigma}_{0}\rangle=\vec{\phi}^{N;\sigma}\cdot\vec{U}^{N;\sigma}_{\rm
gs}$. We may also denote $\vec{\phi}_{\rm gs}\equiv\vec{\phi}^{N;\uparrow}$
and $|\psi_{\rm gs}\rangle\equiv|\psi^{N;\uparrow}_{0}\rangle$ and
$\vec{U}_{\rm gs}\equiv\vec{U}^{N;\uparrow}_{\rm gs}$ for notational
compactness where applicable.

With the above preliminaries established, we now use BWPT to project the full
Hamiltonian $\hat{H}=\hat{H}_{\rm leads}+\hat{H}_{\rm
mol}^{\prime}+\hat{H}_{\rm hyb}$ onto the low-energy manifold containing the
spin-doublet molecule ground states, by perturbatively eliminating excitations
to higher-lying molecule states. We define the projector
$\hat{P}=\sum_{\sigma=\uparrow,\downarrow}|\psi_{0}^{N,\sigma}\rangle\langle\psi_{0}^{N,\sigma}|$
for the ground states to be kept and $\hat{Q}=\hat{1}-\hat{P}$ for the excited
states to be discarded in our effective model. To second order in
$\hat{H}_{\rm hyb}$ we obtain an effective Hamiltonian,

$\hat{H}_{\rm eff}=\hat{P}\hat{H}\hat{P}+\hat{P}\hat{H}_{\rm
hyb}\hat{Q}(E_{\rm gs}-\hat{H}^{\prime}_{\rm mol})^{-1}\hat{Q}\hat{H}_{\rm
hyb}\hat{P}\;.$ (S-10)

With the specific form of the hybridization Hamiltonian,

$\hat{H}_{\rm hyb}=\sum_{\begin{subarray}{c}s=\uparrow,\downarrow\\\
\alpha=s,d\end{subarray}}V_{\alpha}^{\phantom{\dagger}}\left(c_{\alpha
s}^{\dagger}d_{r_{\alpha}s}^{\phantom{\dagger}}+d_{r_{\alpha}s}^{\dagger}c_{\alpha
s}^{\phantom{\dagger}}\right)\;,$ (S-11)

the effective model takes the form,

$\hat{H}_{\rm eff}=\hat{H}_{\rm
leads}+\sum\left(V_{\alpha}V_{\beta}~{}A{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\ \alpha,\beta\end{subarray}}\right)c_{\beta
s^{\prime}}^{\dagger}c_{\alpha
s}^{\phantom{\dagger}}|\psi_{0}^{N;\sigma^{\prime}}\rangle\langle\psi_{0}^{N;\sigma}|\;,$
(S-12)

where a sum over all indices of the tensor $A$ is implied. Charge
conservation, $SU(2)$ spin symmetry and Hermiticity impose constraints on the
the tensor $A$. Together with the definitions,

$\displaystyle\hat{\mathbb{S}}^{+}$
$\displaystyle=|\psi_{0}^{N;\uparrow}\rangle\langle\psi_{0}^{N;\downarrow}|\;$
(S-13) $\displaystyle\hat{\mathbb{S}}^{-}$
$\displaystyle=|\psi_{0}^{N;\downarrow}\rangle\langle\psi_{0}^{N;\uparrow}|\;$
(S-14) $\displaystyle\hat{\mathbb{S}}^{z\;}$
$\displaystyle=\tfrac{1}{2}\left[|\psi_{0}^{N;\uparrow}\rangle\langle\psi_{0}^{N;\uparrow}|-|\psi_{0}^{N;\downarrow}\rangle\langle\psi_{0}^{N;\downarrow}|\right]\;$
(S-15)

for the retained molecule states and

$\hat{\vec{s}}_{\alpha\beta}=\tfrac{1}{2}\sum_{ss^{\prime}}c_{\beta
s^{\prime}}^{\dagger}\vec{\bm{\sigma}}_{s^{\prime}s}~{}c_{\alpha
s}^{\phantom{\dagger}}\;,$ (S-16)

for the conduction electrons, we obtain the effective generalized 2CK model
[4, 5], Eq. 2 of the main paper:

$\hat{H}_{\rm eff}=\hat{H}_{\rm
leads}+\sum_{\alpha\beta}\left[J_{\alpha\beta\;}\hat{\vec{\mathbb{S}}}\cdot\hat{\vec{s}}_{\alpha\beta}+W_{\alpha\beta\;}^{\phantom{\dagger}}\sum_{\sigma}c_{\beta\sigma}^{\dagger}c_{\alpha\sigma}^{\phantom{\dagger}}\right]\;.$
(S-17)

In the following we consider the rescaled couplings $j_{\alpha\beta}$ and
$w_{\alpha\beta}$ defined by
$J_{\alpha\beta}=V_{\alpha}V_{\beta}j_{\alpha\beta}$ and
$W_{\alpha\beta}=V_{\alpha}V_{\beta}w_{\alpha\beta}$, which follow from the
factorized form of the couplings up to second order in Eq. S-12. It
immediately follows that:

$\displaystyle j_{\alpha\beta}$
$\displaystyle=2\left(A{\begin{subarray}{c}\uparrow\uparrow\\\
\uparrow\uparrow\\\
\alpha,\beta\end{subarray}}-A{\begin{subarray}{c}\uparrow\uparrow\\\
\downarrow\downarrow\\\ \alpha,\beta\end{subarray}}\right)\;$ (S-18)
$\displaystyle w_{\alpha\beta}$
$\displaystyle=\tfrac{1}{2}\left(A{\begin{subarray}{c}\uparrow\uparrow\\\
\uparrow\uparrow\\\
\alpha,\beta\end{subarray}}+A{\begin{subarray}{c}\uparrow\uparrow\\\
\downarrow\downarrow\\\ \alpha,\beta\end{subarray}}\right)\;$ (S-19)

in terms of the amplitudes appearing in Eq. S-12. There are other equivalent
forms which follow from the symmetry transformations but we shall use these
ones.

The final step is of course to calculate the amplitudes $A$. This is done by
inserting Eq. S-11 into Eq. S-10 and comparing with Eq. S-12. The amplitudes
decompose into particle and hole tunneling processes

$A{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}=p{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\
\alpha,\beta\end{subarray}}-h{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\ \alpha,\beta\end{subarray}}$ (S-20)

where

$p{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}=\sum_{j}\frac{\langle\psi^{N;\sigma^{\prime}}_{0}|d_{r_{\beta}s^{\prime}}^{\phantom{\dagger}}|\psi^{N+1;\sigma+s}_{j}\rangle\langle\psi^{N+1;\sigma+s}_{j}|d_{r_{\alpha}s}^{\dagger}|\psi^{N;\sigma}_{0}\rangle}{E_{\rm
gs}-E^{N+1;\sigma+s}_{j}}$ (S-21)

for the particle process and

$h{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}=\sum_{j}\frac{\langle\psi^{N;\sigma^{\prime}}_{0}|d_{r_{\beta}s^{\prime}}^{\dagger}|\psi^{N-1;\sigma-s}_{j}\rangle\langle\psi^{N-1;\sigma-s}_{j}|d_{r_{\alpha}s}^{\phantom{\dagger}}|\psi^{N;\sigma}_{0}\rangle}{E_{\rm
gs}-E^{N-1;\sigma-s}_{j}}$ (S-22)

for the hole process. Note that spin conservation requires
$\sigma+s=\sigma^{\prime}+s^{\prime}$ for particle tunneling and
$\sigma-s=\sigma^{\prime}-s^{\prime}$ for hole tunneling.

Numerical calculations can be made dramatically more efficient by exploiting
the following observation: we may express particle and hole amplitudes in
matrix form, viz:

$\displaystyle p{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}$
$\displaystyle=\left[\left(\bm{X}^{N;\sigma^{\prime}}_{\beta
s^{\prime}+}\right)^{\dagger}\left[E_{\rm
gs}\bm{I}-\bm{H}_{\psi}^{N+1;\sigma+s}\right]^{-1}\left(\bm{X}^{N;\sigma}_{\alpha
s+}\right)\right]_{00}$ (S-23) $\displaystyle
h{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}$
$\displaystyle=\left[\left(\bm{X}^{N;\sigma^{\prime}}_{\beta
s^{\prime}-}\right)^{\dagger}\left[E_{\rm
gs}\bm{I}-\bm{H}_{\psi}^{N-1;\sigma-s}\right]^{-1}\left(\bm{X}^{N;\sigma}_{\alpha
s-}\right)\right]_{00}$ (S-24)

where $\bm{X}^{N;\sigma}_{\alpha s\pm}=\left(\vec{\psi}^{N\pm 1,\sigma\pm
s}\right)^{\dagger}d_{r_{\alpha}s}^{(\dagger)}~{}\vec{\psi}^{N;\sigma}$ is
evaluated in the eigenbasis. Switching now to the product basis, we write
$\bm{X}_{\alpha s\pm}^{N;\sigma}=\left(\bm{U}^{N\pm 1;\sigma\pm
s}\right)^{\dagger}\bm{\mathbb{X}}^{N;\sigma}_{\alpha
s\pm}\;\bm{U}^{N;\sigma}$, where $\bm{\mathbb{X}}^{N;\sigma}_{\alpha
s\pm}=\left(\vec{\phi}^{N\pm 1,\sigma\pm
s}\right)^{\dagger}d_{r_{\alpha}s}^{(\dagger)}~{}\vec{\phi}^{N;\sigma}$.
Inserting into Eqs. S-23 and S-24 we finally obtain,

$\displaystyle p{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}$
$\displaystyle=\left(\vec{U}^{N;\sigma^{\prime}}_{\rm
gs}\right)^{\dagger}\bm{\mathbb{M}}{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\ \alpha,\beta\end{subarray}}(+)~{}\vec{U}^{N;\sigma}_{\rm gs}$
(S-25) $\displaystyle h{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\ \alpha,\beta\end{subarray}}$
$\displaystyle=\left(\vec{U}^{N;\sigma^{\prime}}_{\rm
gs}\right)^{\dagger}\bm{\mathbb{M}}{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\ \alpha,\beta\end{subarray}}(-)~{}\vec{U}^{N;\sigma}_{\rm gs}$
(S-26)

with

$\displaystyle\bm{\mathbb{M}}{\begin{subarray}{c}\sigma\sigma^{\prime}\\\
s,s^{\prime}\\\
\alpha,\beta\end{subarray}}(\pm)=\left(\bm{\mathbb{X}}^{N;\sigma^{\prime}}_{\beta
s^{\prime}\pm}\right)^{\dagger}\left[E_{\rm gs}\bm{I}-\bm{H}_{\phi}^{N\pm
1;\sigma\pm s}\right]^{-1}\left(\bm{\mathbb{X}}^{N;\sigma}_{\alpha
s\pm}\right)$ (S-27)

and where $\vec{U}^{N;\sigma}_{\rm gs}$ is the vector of expansion
coefficients for the ground state wavefunction only, as above.

The beauty of Eqs. S-25, S-26, S-27 is that only the ground eigenstate
$|\psi^{N;\sigma}_{0}\rangle$ and the ground state energy $E_{\rm gs}$ of the
isolated molecule are required: one does not need to fully diagonalize the
entire molecule. The matrix $\mathbb{M}$ is simply evaluated in the original
product basis. This is highly efficient since one can use eigensolvers
targeting just the ground state for the calculation, meaning that large
systems can be dealt with cheaply. This in turn permits high throughput
calculations which can be integrated within ML routines for inverse design.

Although the above formulation is entirely general, we note from Eqs. S-18,
S-19 that in the present context, the 2CK mapping requires only a few elements
of the full tensor $A$. Specifically, we may set
$\sigma=\sigma^{\prime}=\uparrow$, which we assume from now on. Furthermore,
spin conservation then implies that $s=s^{\prime}$. Therefore we abbreviate
our full $A$, $p$ and $h$ tensors to,

$A{\begin{subarray}{c}\uparrow\uparrow\\\ s,s\\\
\alpha,\beta\end{subarray}}\to
A^{s\alpha\beta}=p^{s\alpha\beta}-h^{s\alpha\beta}\;.$ (S-28)

Likewise we abbreviate,

$\bm{\mathbb{M}}{\begin{subarray}{c}\sigma\sigma^{\prime}\\\ s,s^{\prime}\\\
\alpha,\beta\end{subarray}}(\pm)\to\bm{\mathbb{M}}^{s\alpha\beta}_{\pm}\;,$
(S-29)

as used in the main paper. Thus we may write,

$\displaystyle p^{s\alpha\beta}$
$\displaystyle=\left(\vec{U}^{N;\uparrow}_{\rm
gs}\right)^{\dagger}\bm{\mathbb{M}}^{s\alpha\beta}_{+}~{}\left(\vec{U}^{N;\uparrow}_{\rm
gs}\right)$ (S-30) $\displaystyle h^{s\alpha\beta}$
$\displaystyle=\left(\vec{U}^{N;\uparrow}_{\rm
gs}\right)^{\dagger}\bm{\mathbb{M}}^{s\alpha\beta}_{-}~{}\left(\vec{U}^{N;\uparrow}_{\rm
gs}\right)$ (S-31)
$\displaystyle\bm{\mathbb{M}}^{s\alpha\beta}_{\pm}=\left(\bm{\mathbb{X}}^{N;\uparrow}_{\beta
s\pm}\right)^{\dagger}\left[E_{\rm gs}\bm{I}-\bm{H}_{\phi}^{N\pm
1;\tfrac{1}{2}\pm s}\right]^{-1}\left(\bm{\mathbb{X}}^{N;\uparrow}_{\alpha
s\pm}\right)$ (S-32) $\displaystyle\bm{\mathbb{X}}^{N;\uparrow}_{\alpha
s\pm}=\left(\vec{\phi}^{N\pm 1,\tfrac{1}{2}\pm
s}\right)^{\dagger}d_{r_{\alpha}s}^{(\dagger)}~{}\left(\vec{\phi}^{N;\uparrow}\right)$
(S-33)

###  S-II.2 SWT results for odd Hubbard chains

As a simple application, we consider the SWT and resulting effective 2CK model
parameters as a function of gate voltage $V_{g}$ for chains of $M$ interacting
sites with on-site Coulomb repulsion $U$, connected by nearest-neighbour
hoppings $t$. The Hamiltonian is given by Eq. S-5, and we take $\epsilon=-U/2$
here so that the system possesses exact particle-hole symmetry when $V_{g}=0$.
Leads are connected to either end of the chain, Eq. S-6. For odd $M$, the
molecule supports a unique spin-doublet ground state for gate voltages around
$V_{g}=0$ (the Coulomb blockade transitions to $M\pm 1$ electron states depend
on $U$, and $t$). When the molecule hosts a spin-doublet ground state, the
effective 2CK model, Eq. S-17, describes the low-energy behavior.

Figure S-1: Effective model parameters $j_{\alpha\beta}$ and $w_{\alpha\beta}$
as a function of gate voltage $V_{g}$ from SWT calculations on chains of
length $M=1,3,5,7$ sites, with local interaction $U$ and nearest-neighbour
hopping $t$. In this example we have taken $U/t=1$ and set $t=0.5$. Source and
drain leads are connected to sites $1$ and $M$. Kondo exchange interactions
$j_{ss}$ and $j_{dd}$ (blue dashed lines) are equal and antiferromagnetic
throughout the entire gate voltage range shown where the many-body molecule
ground state is a spin-doublet. Potential scatterings $w_{ss}$, $w_{dd}$, and
$w_{sd}$ vanish at the particle-hole symmetric point $V_{g}=0$ in all cases.
The sign of $j_{sd}$ is seen to alternate, being antiferromagnetic for
$M=4m+1$ and ferromagnetic for $M=4m+3$.

Fig. S-1 shows the effective 2CK model parameters $j_{\alpha\beta}$ and
$w_{\alpha\beta}$ as a function of gate voltage $V_{g}$ for odd $M=1,3,5,7$ as
derived using the numerical SWT with $U/t=1$. We note that the prescription
developed above allows the calculations to be performed cheaply even for the
$M=7$ chain using exact diagonalization, with the entire gate-voltage sweep
obtained in a few minutes on a standard desktop computer.

The $M=1$ case is equivalent to a two-lead version of the AIM, and is a
special limit since both source and drain lead couple to the same molecule
orbital ($r_{s}=r_{d}=1$). This is the only ‘molecule’ to satisfy
‘proportionate coupling’, which allows us to reduce the system exactly to a
single-channel model. This is achieved by canonical transformation of the lead
operators from source-drain basis to even-odd basis,

$\displaystyle c_{e\sigma}$
$\displaystyle=\tfrac{1}{\tilde{V}}\left(V_{s}c_{s\sigma}+V_{d}c_{d\sigma}\right)\;,$
$\displaystyle c_{o\sigma}$
$\displaystyle=\tfrac{1}{\tilde{V}}\left(V_{d}c_{s\sigma}-V_{s}c_{d\sigma}\right)\;,$
(S-34)

where $\tilde{V}^{2}=V_{s}^{2}+V_{d}^{2}$, such that the odd lead is
decoupled. This exact decoupling on the level of the bare Hamiltonian is not
possible for $M>1$. For $M=1$ however, the model reduces to the 1-channel AIM
and the usual SWT results then pertain on mapping to an effective 1-channel
Kondo (1CK) model involving only the even lead [3]:

$\displaystyle\hat{H}_{1CK}$ $\displaystyle=\hat{H}_{\rm
lead}^{e}+J_{ee}\hat{\vec{\mathbb{S}}}\cdot\hat{\vec{s}}_{ee}+W_{ee}\sum_{\sigma}c_{e\sigma}^{\dagger}c_{e\sigma}^{\phantom{\dagger}}$
(S-35) $\displaystyle J_{ee}$
$\displaystyle=2\tilde{V}^{2}\left[\frac{1}{U+\epsilon+V_{g}}-\frac{1}{\epsilon+V_{g}}\right]\;,$
(S-36) $\displaystyle W_{ee}$
$\displaystyle=-\tilde{V}^{2}\left[\frac{1}{U+\epsilon+V_{g}}+\frac{1}{\epsilon+V_{g}}\right]\;.$
(S-37)

On inverting the transformation Eq. S-34 and returning to the physical source-
drain lead basis, the effective 1-channel Kondo model Eq. S-35 becomes
precisely our original generalized 2CK model, Eq. S-17. However, the
underlying single-channel description for $M=1$ implies additional relations
between the coupling parameters that are not present in general [5]. For $M=1$
we have $j_{ss}=j_{dd}=j_{sd}=j_{ds}$ and likewise
$w_{ss}=w_{dd}=w_{sd}=w_{ds}$ (where, as before, we have divided out the
hybridization dependence $J_{\alpha\beta}=V_{\alpha}V_{\beta}j_{\alpha\beta}$
and $W_{\alpha\beta}=V_{\alpha}V_{\beta}w_{\alpha\beta}$). The important
conclusion is that the source-drain mixing terms $j_{sd}$ and $w_{sd}$ are
equal to the local terms when there is an underlying 1CK description. Although
the single-channel limit is certainly not generic for multi-orbital
nanostructures, $j_{sd}$ and $w_{sd}$ are typically finite: these are the
terms that mediate a source-drain current through the molecular junction. By
contrast, for the 2CK critical point we require $j_{sd}=w_{sd}=0$, which
corresponds to a conductance node.

Figure S-2: Effective model parameters $j_{\alpha\beta}$ and $w_{\alpha\beta}$
as a function of gate voltage $V_{g}$ from SWT calculations for the $M=5$ site
molecular junction. We take local on-site interaction $U$ and nearest-
neighbour hopping $t$, with $U/t=1$ and $t=0.5$ as usual. Here, the next-next-
nearest-neighbour hopping $t^{\prime}$ (connecting sites $1\leftrightarrow 4$
and sites $2\leftrightarrow 5$ of the chain) is now also included. The
inclusion of $t^{\prime}$ does not break exact particle-hole symmetry at
$V_{g}=0$, but does continuously interpolate between $M=3$ and $M=5$ limiting
cases of the pure chains shown in Fig. S-1 by providing a “shortcut” route
through the molecule. This implies the existence of a node in $j_{sd}$ at
$V_{g}=0$, which we confirm at a critical value of
$t^{\prime}=t^{\prime}_{c}\simeq 0.4283$, as shown in the middle panel.

Returning now to Fig. S-1, we see that for $M=1$ the 1CK condition on the
couplings is indeed satisfied. However, more interestingly, we find that
$j_{ss}=j_{dd}=j_{sd}=j_{ds}$ and $w_{ss}=w_{dd}=w_{sd}=w_{ds}$ very
accurately holds (although not exactly) for $M=5$, suggesting that (perhaps
surprisingly) the odd lead combination is almost decoupled in the $M=5$ case
at energies or temperatures $\ll U$ where the SWT effective model is expected
to hold. This holds also for $M=9$ (not shown) and we conjecture is the case
for all $M=4m+1$ (note however that the range of gate over which the doublet
ground state occurs narrows on increasing $M$). Furthermore, for $M=3$ and
$M=7$ (and again we suspect $M=4m+3$ in general), we have
$j_{ss}=j_{dd}\simeq-j_{sd}=-j_{ds}$ and $w_{ss}=w_{dd}\simeq-w_{sd}=-w_{ds}$.
This also implies an effective 1CK description at low energies or
temperatures, but this time involving the odd lead combination (and with the
even lead decoupling).

###  S-II.3 SWT results for 5-site cluster

The effective 1CK behavior of the odd-length Hubbard chains shown above is
particular to the chain geometry around half-filling. Introducing couplings
between different sites clearly produces greater complexity on the molecule.
Here we consider the 5-site molecular cluster, Eq. S-7, which has $sd$
symmetry, and also an exact $ph$ symmetry at $V_{g}=0$.

We expect interesting crossover behavior as a function of the ratio
$t^{\prime}/t$. For small $t^{\prime}/t$, the $t^{\prime}$ connections are a
small perturbation and we expect physics similar to the pure 5-site chain. On
the other hand, for large $t^{\prime}/t$, we may regard $t$ as the
perturbation . Source and drain leads then have a dominant connecting path
that is of length 3; the $t^{\prime}$ bonds offer a kind of ‘shortcut’ through
the $M=5$ chain. Tuning $t^{\prime}/t$ therefore allows us to continuously
deform between the $M=5$ and $M=3$ limits. Since the sign of $j_{sd}$ changes
on going from $M=5$ to $M=3$ we expect a critical value
$t^{\prime}=t^{\prime}_{c}$ where $j_{sd}=0$ just vanishes at $V_{g}=0$. As
shown in Fig. S-2, this is precisely what we find. Here we plot the SWT model
parameters for the model Eq. S-7, with source and drain leads connected to
sites 1 and 5, and for $U/t=1$ and $t=0.5$ as a function of gate voltage
$V_{g}$. The ratio $t^{\prime}/t$ is increased from top to middle to bottom
panels. In the top panel, we see behavior similar to that of the pure $M=5$
chain shown in Fig. S-1, but with the effective single-channel condition now
violated ($j_{ss}=j_{dd}\neq j_{sd}$ and $w_{ss}=w_{dd}\neq w_{sd}$). In the
bottom panel we see behavior instead more similar to the pure $M=3$ chain,
with negative $j_{sd}$. The middle panel shows the ‘sweet spot’ condition
where $j_{sd}$ precisely vanishes at $V_{g}=0$ due to many-body QI effects of
the competing pathways through the 5-site molecular moiety.

##  S-III QI Classes

Ref. [4] made a detailed study of many-body QI effects in molecular junctions,
focusing on the transport problem within perturbation theory. Alternant
hydrocarbon molecules were studied, whose structures can be accommodated on a
bipartite graph with no odd loops, and to which the Coulson-Rushbrooke pairing
theorem applies [10, 11]. Exchange cotunneling and potential scattering terms
equivalent to our $J_{sd}$ and $W_{sd}$ were identified for systems with a net
spin-doublet ground state.

A key insight from Ref. [4] is that QI effects in such systems can be
classified according to the number of nodes found by varying the gate voltage
all the way across a given charge state. In our notation, for a particular
molecular junction they define,

$\mathcal{Q}_{W}=\frac{\langle\psi_{0}^{N;\tfrac{1}{2}}|d_{r_{s}\uparrow}^{\dagger}|\psi_{0}^{N-1;0}\rangle\langle\psi_{0}^{N-1;0}|d_{r_{d}\uparrow}|\psi_{0}^{N;\tfrac{1}{2}}\rangle}{\langle\psi_{0}^{N;\tfrac{1}{2}}|d_{r_{d}\uparrow}|\psi_{0}^{N+1;1}\rangle\langle\psi_{0}^{N+1;1}|d_{r_{s}\uparrow}^{\dagger}|\psi_{0}^{N;\tfrac{1}{2}}\rangle}\;,$
(S-38)

and

$\mathcal{Q}_{J}=\frac{\langle\psi_{0}^{N;-\tfrac{1}{2}}|d_{r_{s}\downarrow}^{\dagger}|\psi_{0}^{N-1;0}\rangle\langle\psi_{0}^{N-1;0}|d_{r_{d}\uparrow}|\psi_{0}^{N;\tfrac{1}{2}}\rangle}{\langle\psi_{0}^{N;-\tfrac{1}{2}}|d_{r_{d}\uparrow}|\psi_{0}^{N+1;0}\rangle\langle\psi_{0}^{N+1;0}|d_{r_{s}\downarrow}^{\dagger}|\psi_{0}^{N;\tfrac{1}{2}}\rangle}\;,$
(S-39)

where a sum over possible degeneracies of the $N\pm 1$ ground states is
implied.

Ref. [4] showed that for bipartite molecules $\mathcal{Q}>0$ describes the odd
QI class with an odd number of QI nodes found upon scanning $V_{g}$, whereas
$\mathcal{Q}<0$ describes the even QI class with an even number of QI nodes.
Indeed, $\mathcal{Q}_{W}>0$ if the source and drain leads connect to the same
sublattice, whereas $\mathcal{Q}_{W}<0$ if the source and drain leads connect
to different sublattices. Furthermore, bipartite molecules with an odd number
of sites (which therefore have a spin-doublet ground state at their $ph$
symmetric point) satisfy the rule $\mathcal{Q}_{W}=-\mathcal{Q}_{J}$, meaning
that $W$ and $J$ always necessarily belong to different QI classes. In a given
molecule charge sector, if $W$ shows an odd number of nodes on scanning
$V_{g}$ then $J$ must show an even number. This is borne out by our SWT
results in Fig. S-1, where the odd Hubbard chains (which have leads coupled to
the same sublattice) have a single node in $W$ but no nodes in $J$. Ref. [4]
concludes that a true conductance node is therefore ‘generally not possible’
in such systems since $J_{sd}$ and $W_{sd}$ do not generally have a node at
the same $V_{g}$.

However, it is also possible that two QI nodes in the even QI class can occur
at the same gate voltage (degenerate case). This is precisely what we see in
Fig. S-2. $W$ is in the odd QI class with a node in $W_{sd}$ at $V_{g}=0$. On
the other hand $J$ is in the even QI class and has no nodes for
$t^{\prime}<t^{\prime}_{c}$ but two nodes for $t^{\prime}>t_{c}$. The two
nodes coalesce and are exactly degenerate at $V_{g}=0$ when precisely at
$t^{\prime}=t^{\prime}_{c}$. This is our 2CK critical point. The QI
classification therefore provides a topological perspective on the 2CK phase
transition in terms of coalescing and annihilating QI nodes in the 2CK model
parameters.

Finally, we note that this analysis holds only for bipartite molecules.
Therefore more general structures with odd loops (such as the one considered
in Fig. 5 of the main text) need not be so constrained.

##  S-IV Numerical Renormalization Group

As noted above, the SWT is a perturbative technique and more sophisticated
methods are required to perform the exact mapping [6]. In particular, the SWT
does not take into account renormalization effects from the conduction
electrons in the leads [8]. It is also confined to the perturbative regime of
large $U/V$. And of course, the solution of the effective model is still
required to understand the low-energy physics of the system and compute
physical observables of interest relevant to experiments, such as the series
conductance.

The numerical renormalization group [9] (NRG) is considered the gold-standard
method of choice for solving generalized quantum impurity problems, such as
those describing the low-energy physics of molecular junctions [5].
Furthermore, experimentally-relevant transport signatures can be accurately
calculated using NRG down to low temperatures [12], including renormalization
effects from strong electron correlations, such as Coulomb blockade, Kondo
effect, and many-body quantum interference phenomena. NRG is a numerically-
exact, non-perturbative technique for solving quantum impurity type problems.
The NRG method involves the following steps [9].

(i) The source and drain lead density of states $\rho(\omega)$ is divided up
into intervals on a logarithmic grid according to the discretization points
$\pm D\Lambda^{-n}$, where $D$ is the bare conduction electron bandwidth,
$\Lambda>1$ is the NRG discretization parameter, and $n=0,1,2,3,...$. The
continuous electronic density in each interval is replaced by a single pole at
the average position with the same total weight, yielding $\rho^{\rm
disc}(\omega)$. In this work, as is conventional, we assume metallic leads
with a constant density of states $\rho_{0}=1/2D$ within a flat band of half-
width $D$.
(ii) The conduction electron part of the Hamiltonian $H_{\rm leads}$ is then
mapped onto semi-infinite tight-binding chains for each lead $\alpha=s,d$:

$\displaystyle H_{\rm leads}\to H_{\rm leads}^{\rm
disc}=\sum_{\alpha,\sigma}\sum_{n=0}^{\infty}$ $\displaystyle
t_{n}^{\phantom{\dagger}}\left(f_{\alpha n,\sigma}^{\dagger}f_{\alpha
n+1,\sigma}^{\phantom{\dagger}}+{\rm H.c.}\right)\;,$

where the Wilson chain coefficients $\\{t_{n}\\}$ are determined such that the
density of states at the end of each chain reproduces exactly the discretized
lead density of states, that is $-\tfrac{1}{\pi}{\rm Im}~{}\langle\langle
f_{\alpha 0,\sigma}^{\phantom{\dagger}};f_{\alpha
0,\sigma}^{\dagger}\rangle\rangle=\rho^{\rm disc}(\omega)$. Due to the
logarithmic discretization, the Wilson chain parameters decay roughly
exponentially down the chain, $t_{n}\sim\Lambda^{-n/2}$. However the detailed
form of the $t_{n}$ encode the specific lead density of states [9].
(iii) The hybridizing molecule orbitals are coupled to site $n=0$ of the two
Wilson chains. We define a sequence of Hamiltonians $H_{N}$ comprising the
molecule and the first $N+1$ Wilson chain sites,

$\begin{split}H_{N}=&H_{\rm mol}+H_{\rm
hyb}+\sum_{\alpha,\sigma}\sum_{n=0}^{N-1}t_{n}^{\phantom{\dagger}}\left(f_{\alpha
n,\sigma}^{\dagger}f_{\alpha n+1,\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\;,\end{split}$ (S-40)

We now define the recursion relation,

$H_{N+1}=H_{N}+\sum_{\alpha,\sigma}t_{N}^{\phantom{\dagger}}\left(f_{\alpha
N,\sigma}^{\dagger}f_{\alpha N+1,\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\;,$ (S-41)

such that the full (discretized) model is obtained as $H^{\rm
disc}=\lim_{N\to\infty}H_{N}$.
(iv) Starting from the molecule, the chain is built up successively by adding
Wilson chain sites using the recursion, Eq. S-41. At each step $N$, the
intermediate Hamiltonian $H_{N}$ is diagonalized, and only the lowest $N_{s}$
states are retained to construct the Hamiltonian $H_{N+1}$ at the next step,
with the higher energy states being discarded. With each iteration we
therefore focus on progressively lower energy scales. Furthermore, the
iterative diagonalization and truncation procedure can be viewed as an RG
transformation [9], $H_{N+1}=\mathcal{R}[H_{N}]$.
(v) The sequence of $H_{N}$ can be viewed as coarse-grained versions of the
full model, which faithfully capture the physics at progressively lower and
lower temperatures. Useful information is therefore extracted from each step,
and physical observables can be obtained at essentially any temperature or
energy scale [9].

In this work we use a discretization parameter $\Lambda=3$ and retain
$N_{s}=5000$ states at each NRG step. Total charge and spin projection abelian
quantum numbers are exploited in the block diagonalization procedure.

In the main text we present results for the linear response series
differential conductance $G_{c}=dI_{sd}/dV_{b}$ through the molecular junction
due to a small source-drain bias voltage $V_{b}$. Since anything more complex
than the $M=1$ single-impurity limit does not satisfy the proportionate
coupling condition, transport coefficients cannot be simply related to
molecule Green’s functions. However, the Kubo formula can still be used [13].
Here we use the ‘improved Kubo’ formulation given in Ref. [12], which produces
much more accurate results in the context of NRG calculations for the linear
response electrical conductance.

##  S-V 2CK critical point

By tuning the ratio $t^{\prime}/t$ in our 5-site molecular moiety Eq. S-7, the
SWT results indicate that we can simultaneously eliminate both $j_{sd}$ and
$w_{sd}$ terms of the effective 2CK model. The $sd$ symmetry of the model
means that the local Kondo couplings $j_{ss}=j_{dd}$ are equal. Since they are
also antiferromagnetic, we expect to be able to realize the QI-2CK effect in
this prototype molecular junction. In Fig. S-3 we demonstrate that this is
indeed the case, using full NRG results. Furthermore, the QI-2CK effect is
found to be remarkably robust, holding beyond the perturbative regime of
applicability of the SWT.

Figure S-3: NRG results for the molecule contribution to the entropy $S_{\rm
mol}$ as a function of temperature $T$ for the 5-site system tuned to the
critical point at $t^{\prime}=t^{\prime}_{c}$, for different interaction
strengths $U/t$ in panels (a), (b), and (c). As usual we set $t=0.5$ and the
conduction electron bandwidth is $D=1$. In each case, the molecule-lead
hybridization $V$ is varied, showing that the many-body QI-driven 2CK critical
point with finite residual ($T=0$) entropy $S_{\rm mol}=\tfrac{1}{2}\ln(2)$
can be realized in all cases. The 2CK Kondo temperature $T_{\rm K}$ is seen to
depend sensitively on $V$ (see Fig. 2 of the main paper); but we note that a
remarkably large $T_{\rm K}\sim 10^{-2}U$ can be achieved in all cases, where
no incipient local moment regime is observed and RG flow proceeds directly to
the 2CK fixed point. For the black, red, blue, green and magenta lines we have
$V=0.03,0.04,0.05,0.06,0.08$ for $U/t=0.1$ in panel (a);
$V=0.1,0.125,0.15,0.2,0.3$ for $U/t=1$ in panel (b); and
$V=0.3,0.4,0.5,0.6,0.75$ for $U/t=10$ in panel (c). The magenta lines show the
maximum possible $T_{\rm K}$ achievable for each ratio of $U/t$.

In Fig. S-3 we plot the molecule contribution to the total entropy $S_{\rm
mol}(T)$ as a function of temperature $T$ for systems tuned precisely to the
QI-2CK critical point ($V_{g}=0$ and $t^{\prime}=t^{\prime}_{c}$). Panels
(a,b,c) show different interaction strengths $U/t$ (with $t=0.5$ and $D=1$
fixed), over a very wide range from perturbative ($U/t=0.1$) to
intermediate/non-perturbative strength ($U/t=1$), to strong coupling
($U/t=10$). We have also tested other ratios of $U/t$ and find equivalent
results. The different lines plotted in each panel are for different molecule-
lead hybridization strengths $V$. Small $V$ and large $U$ corresponds to the
regime of applicability of the SWT calculation, but NRG allows to interrogate
the non-perturbative regimes as well – as shown in the figure. In particular,
the magenta line in each panel is for a comparatively large $V$ which yields
the largest 2CK Kondo temperature $T_{\rm K}$, while the black lines are for
smaller $V$ where the corresponding $T_{\rm K}$ is then smaller.

Interestingly, we find that for all $U/t$ considered, we can access the QI-2CK
effect for any $V$. The hallmark of the 2CK effect is the residual ($T=0$)
molecule entropy of $\tfrac{1}{2}\ln(2)$ [14], which we here see in all cases.
For smaller $V$ the Kondo scale $T_{\rm K}$ is small and we have a window over
which the molecule entropy is $\ln(2)$. This corresponds to the local moment
regime where the molecule hosts an essentially free spin-doublet state (this
is the ground molecule doublet state that we project onto using the SWT).
However, the most remarkable observation from our NRG results is that large
$T_{\rm K}$ can be realized by tuning $V$. The largest $T_{\rm K}$ possible
for each $U/t$ is in fact of order $10^{-2}U$, which is a high-energy scale
when compared with the usual exponentially-small Kondo temperatures [3] found
in regular single molecule junctions or semiconductor quantum dot devices.
This is promising for experimental realization at accessible temperature
scales. Finally, we note that as $U\to 0$, the maximum $T_{\rm K}$ also
vanishes, meaning that as the non-interacting limit is approached, the 2CK
effect also disappears. For typical interaction strengths in molecules or
quantum dots, however, we expect the QI-2CK effect to be accessible at
unusually high temperatures.

Figure S-4: To confirm that the many-body QI-driven critical point in the
5-site molecular junction indeed realizes the classic 2CK critical point of
the pure 2CK model, we show the ‘impurity’ entropy $S_{\rm imp}(T)$ vs
$T/T_{\rm K}$ (top panel) and the $T=0$ scattering t-matrix $t(\omega)$ vs
$\omega/T_{\rm K}$ (bottom panel) rescaled in terms of the Kondo temperature
$T_{\rm K}$ in the universal regime where $T_{\rm K}$ itself is $\ll U,D,t$.
Full NRG results for the molecular junction (black solid lines) are compared
with those for the pure 2CK model (dashed red lines), and agree perfectly.
Shown for $U/t=10$, $t=0.5$, $D=1$, $V=0.3$, $V_{g}=0$,
$t^{\prime}=t^{\prime}_{c}$.

In Fig. S-4 we confirm the SWT prediction that at the critical point of our
5-site molecule we realize the standard 2CK effect. For this, we select a
molecule realization in the universal regime with a small $T_{\rm K}$. When
scaled in terms of $T_{\rm K}$, physical observables should collapse to a
universal curve characteristic of the pure 2CK model. The top panel of Fig.
S-4 shows NRG results for the entropy, rescaled in terms of $T/T_{\rm K}$
(black solid line), compared with standard results for the pure 2CK model. In
the lower panel we show the energy-resolved conduction electron scattering
t-matrix $t(\omega)$ at $T=0$ vs $\omega/T_{\rm K}$ for the molecule, compared
with universal results of the pure 2CK model. Both agree perfectly in the
universal regime considered. On the other hand, we note that the standard 2CK
model presupposes the existence of an ‘impurity’ spin-$\tfrac{1}{2}$ degree of
freedom. Our results show that for larger $V$, the molecule hosts a 2CK
critical regime that is not contained within the pure 2CK model since there is
no incipient local moment regime.

##  S-VI FL crossover

Figure S-5: NRG results for the Fermi liquid crossover in the 5-site molecular
junction near the QI-2CK critical point, with $t^{\prime}$ successively
approaching the critical point at $t^{\prime}_{c}$. Main panel shows the
entropy crossover, in which the critical $\tfrac{1}{2}\ln(2)$ entropy is
quenched below $T^{*}$. Inset shows the dependence of $T^{*}$ itself on the
perturbation $|t^{\prime}_{c}-t^{\prime}|$. Dashed line fit shows the robust
power-law behavior $T^{*}\sim|t^{\prime}_{c}-t^{\prime}|^{2}$. We take
$U/t=10$, $V=0.75$, $t=0.5$, $D=1$ and $V_{g}=0$.  Figure S-6: NRG results for
the Fermi liquid crossover in the 5-site molecular junction, with
$t^{\prime}=t^{\prime}_{c}$ but varying the gate voltage $V_{g}$ to
destabilize the QI-driven 2CK critical point. Main panel shows the
temperature-dependence of the entropy for each gate voltage setting. Inset
shows the resulting crossover scale $T^{*}$ as a function of $V_{g}$,
demonstrating with the black dashed line the power-law behavior
$T^{*}\sim|V_{g}|^{4}$. We take $U/t=10$, $V=0.75$, $t=0.5$, $D=1$.

Having established that the 2CK critical point is accessible in molecular
junctions by exploiting the many-body QI effect, we now wish to study in more
detail the vicinity of the critical point and the relevant perturbations that
drive the system along the NFL-FL crossover towards the stable fully-Kondo-
screened ground state. We consider here two perturbations that destabilize the
2CK critical point in our 5-site molecular junction: deviations in the
couplings away from $t^{\prime}=t^{\prime}_{c}$ (see Fig. S-5) and deviations
in the gate voltage away from particle hole symmetry at $V_{g}=0$ (see Fig.
S-6).

When the perturbations are sufficiently small, the RG flow for the system
proceeds first to the 2CK critical point (with characteristic
$\tfrac{1}{2}\ln(2)$ entropy) on the scale of $T_{\rm K}$, and then it flows
away from the critical fixed point and towards the FL ground state (with zero
entropy) on the emergent low-energy scale $T^{*}$ [14, 15, 16]. For good scale
separation $T^{*}\ll T_{\rm K}$ we expect a simple power-law dependence of the
$T^{*}$ scale on the strength of the destabilizing perturbation. The power-
laws reflect the scaling dimension of the leading RG relevant operators at the
fixed point induced by the perturbation, and are a characteristic fingerprint
of the critical fixed point. In the present context of the generalized 2CK
model, we expect [14] $T^{*}\sim aJ_{sd}^{2}+bW_{sd}^{2}$. But $J_{sd}$ and
$W_{sd}$ themselves depend on $t^{\prime}$ and $V_{g}$. In Fig. S-5 we show
that finite $|t^{\prime}-t^{\prime}_{c}|$ does indeed generate a finite FL
scale $T^{*}$ and corresponding entropy crossover to a fully Kondo screened
state. The inset shows specifically that

$T^{*}\sim|t^{\prime}-t^{\prime}_{c}|^{2}\qquad:~{}V_{g}=0$ (S-42)

Likewise in In Fig. S-6 we see that the FL crossover is also generated by
finite $V_{g}$ even at $t^{\prime}=t^{\prime}_{c}$ with,

$T^{*}\sim|V_{g}|^{4}\qquad~{}~{}~{}~{}~{}:~{}t^{\prime}=t^{\prime}_{c}$
(S-43)

The difference is because although $J_{sd}$ depends quadratically on both
$t^{\prime}$ and $V_{g}$, the critical value of $t^{\prime}$ is at a finite
value $t^{\prime}_{c}$ while the critical value of $V_{g}$ is at $V_{g}=0$.

##  S-VII ML optimization and
gradient descent

Although we were able to identify conditions for the QI-2CK effect in
molecular junctions in the preceding sections using a combination of symmetry
arguments and an observation about the SWT for odd chains, a more systematic
approach is required for inverse design in general. Indeed, even with the
topology of the model Eq. S-7 set, some fine-tuning of the ratio
$t^{\prime}/t$ is still required. On the level of the SWT, one can fix all
parameters except $t^{\prime}$ and then use a simple root-finding algorithm to
determine the value of $t^{\prime}$ that yields $j_{sd}=0$. The critical value
of $t^{\prime}$ so determined will however only be an approximation to the
true value of $t^{\prime}_{c}$ found in full NRG calculations, since the SWT
is of course only a perturbative approximation (albeit a rather good one). For
the QI-2CK problem, one can again fix all parameters except $t^{\prime}$ and
run the NRG multiple times, varying $t^{\prime}$, and search for the critical
value $t^{\prime}_{c}$ that minimizes the FL scale $T^{*}$. With only a single
parameter to optimize, such a search is very fast and efficient, requiring
only a handful of ‘function evaluations’ to converge exponentially quickly
towards the critical point (here ‘function evaluation’ means either doing the
SWT mapping numerically to find $j_{sd}$ or performing an NRG run to determine
the $T^{*}$ scale). Furthermore, the SWT result provides a good starting point
guess for the NRG calculation.

However, inverse design of molecular junctions is in general much more
complicated than the simplest constrained single-parameter optimization
described above. We illustrate this in the main paper with the example of
unconstrained optimization of the 4-site molecular cluster, demonstrating that
machine learning (ML) methodologies can be applied to the numerical SWT
calculation to realize the QI-2CK condition. We envision that this approach
will be useful beyond the specific example of the 2CK problem, and will be
applicable to the inverse design of molecular junctions for functionalities
such as switching, rectification, thermoelectric properties etc, where the
desired low-energy effective model parameters are known, but the optimal
parent microscopic systems to realize them have yet to be identified.

###  S-VII.1 Inverse design using SWT

For molecular junctions, we imagine fixing the lead and hybridization
Hamiltonians, and searching through the large configuration space of molecular
structures that yield the desired effective model parameters from the SWT. For
demonstration here we do not yet consider realistic molecules whose structure
and chemistry are determined from ab initio methods or from experimental
analysis, but rather simply take the parameters of our generalized
$\hat{H}_{\rm mol}$ molecular model to be freely tunable. Consider an abstract
model for the isolated molecule of the form,

$\hat{H}_{\rm mol}=\sum_{i}\theta_{i}\hat{h}_{i}$ (S-44)

with coupling constants $\theta_{i}$ (which we assemble into a vector
$\vec{\theta}$) and operators $\hat{h}_{i}$ that are defined on the molecule
Hilbert space. Provided that the molecule ground state over the range of gate
voltages of interest is a unique spin-double state, we can apply the numerical
SWT prescription described above in Sec. S-II. The desired effective model
parameters are obtained by evaluating Eqs. S-18-S-20 and S-25-S-27.

The desired functionality of the molecular junction places conditions on the
effective model parameters $j_{\alpha\beta}$ and $w_{\alpha\beta}$. We encode
these into a loss function
$\mathcal{L}(\vec{\theta})=f(j_{ss},j_{dd},j_{sd},w_{ss},w_{dd},w_{sd})$ where
$f$ is some objective function to be minimized that defines the optimization
landscape, whereas $j_{\alpha\beta}$, $w_{\alpha\beta}$ depend implicitly on
$\vec{\theta}$ through the SWT calculation. Despite the typically high
complexity of $\hat{H}_{\rm mol}$, the ML optimization of the molecule can be
done efficiently by gradient descent (GD). This involves calculating the
derivative of the loss function $\vec{\nabla}_{\theta}\mathcal{L}$ with
respect to the optimization parameters $\vec{\theta}$ and taking a finite step
in the direction of the gradient vector from an initial guess
$\vec{\theta}_{0}$ to a new guess $\vec{\theta}_{1}$ such that
$\mathcal{L}(\vec{\theta}_{1})<\mathcal{L}(\vec{\theta}_{0})$. This is
repeated iteratively until convergence within some tolerance to the minimum in
$\mathcal{L}$, whence we have obtained the optimal set of coupling constants
$\vec{\theta}$. We note however that care must be taken since the optimization
landscape may contain many local minima or indeed degenerate global minima.

To make this more concrete, consider the generalized 4-site molecule model,
Eq. S-8, with $U_{m}\equiv U=1$. We will treat the set of single-particle
coupling constants $t_{mn}$ and the overall gate voltage $V_{g}$ as free
tuning parameters that comprise the vector $\vec{\theta}$ to be optimized.
Focusing on the QI-2CK effect, we demand a spin-doublet ground state giving an
effective model for which $j_{sd}=0$ and $w_{sd}=0$, as well as
$j_{ss}=j_{dd}>0$. A simple starting point for the loss function is,

$\mathcal{L}(\vec{\theta})=a(j_{sd})^{2}+b(w_{sd})^{2}+c(j_{ss}-j_{dd})^{2}+d\sum_{\alpha}(|j_{\alpha\alpha}|-j_{\alpha\alpha})$
(S-45)

where $a,b,c,d$ are optimization hyperparameters that affect the optimization
landscape and hence the learning curve. We also add a term that punishes gate
voltages close to the Coulomb blockade steps, since the SWT breaks down in
these regimes and the effective parameters $j_{\alpha\beta}$ and
$w_{\alpha\beta}$ spuriously diverge. In practice, in a given iteration of the
optimization procedure we fix a set of $t_{mn}$ and scan over the entire gate
voltage range to identify regions with a spin-doublet ground state, and find
the minimum of $\mathcal{L}$ for that gate sweep. Designing a good loss
function is the art behind any complex optimization procedure.

In order to perform the GD, we would like the analytic gradient of the loss
function. A numerical evaluation is costly when the number of parameters to be
optimized is large, and finite-difference methods are inaccurate in practice.
Although the analytic gradient with respect to the effective model parameters
$j_{\alpha\beta}$ and $w_{\alpha\beta}$ is trivial, the desired gradient with
respect to the bare model parameters $\vec{\theta}$ is much more challenging.
The key step is to work out the derivative of the particle and hole tunneling
amplitudes in Eqs. S-25, S-26 with respect to a given tuning parameter
$\theta_{i}$.

An important implication of our SWT formulation is that the only non-trivial
information needed to compute the gradient of the loss function analytically
is the derivative of the ground eigenstate $\partial_{\theta_{i}}\vec{U}_{\rm
gs}$ and the derivative of the ground state energy
$\partial_{\theta_{i}}E_{\rm gs}$. Suppressing indices to simplify the
notation in the following, from Eq. S-25,

$\begin{split}\partial_{\theta_{i}}p=&\left(\partial_{\theta_{i}}\vec{U}_{\rm
gs}\right)^{\dagger}\mathbb{M}_{+}\vec{U}_{\rm gs}+\vec{U}_{\rm
gs}^{\dagger}\mathbb{M}_{+}\left(\partial_{\theta_{i}}\vec{U}_{\rm
gs}\right)\\\ &+\vec{U}_{\rm
gs}^{\dagger}\Big{(}\partial_{\theta_{i}}\mathbb{M}_{+}\Big{)}\vec{U}_{\rm
gs}\end{split}$ (S-46)

and similarly for $\partial_{\theta_{i}}h$. Then from Eq. S-27 it follows,

$\begin{split}&\partial_{\theta_{i}}\mathbb{M}_{\pm}=\\\
&\mathbb{X}_{\pm}^{\dagger}\left[E_{\rm
gs}\bm{I}-\bm{H}_{\phi}\right]^{-1}\left[\bm{h}^{i}_{\phi}-\partial_{\theta_{i}}E_{\rm
gs}\bm{I}\right]\left[E_{\rm
gs}\bm{I}-\bm{H}_{\phi}\right]^{-1}\mathbb{X}_{\pm}\end{split}$ (S-47)

where $\bm{h}^{i}_{\phi}=\vec{\phi}^{\dagger}\hat{h}_{i}~{}\vec{\phi}$ such
that $\bm{H}_{\phi}=\sum_{i}\theta_{i}\bm{h}^{i}_{\phi}$. Note that
$\mathbb{X}_{\pm}$, $\bm{H}_{\phi}$ and $\bm{h}^{i}_{\phi}$ are all evaluated
in the product basis and are therefore independent of the optimization
procedure. Only $\partial_{\theta_{i}}\vec{U}_{\rm gs}$ and
$\partial_{\theta_{i}}E_{\rm gs}$ must be recomputed at each new optimization
step. For the latter, we may simply use the Hellmann–Feynman theorem,

$\partial_{\theta_{i}}E_{\rm gs}=\langle\psi_{\rm gs}|\hat{h}_{i}|\psi_{\rm
gs}\rangle\;,$ (S-48)

whereas the former follows from a standard result in non-degenerate
perturbation theory [17],

$\partial_{\theta_{i}}\vec{U}_{\rm gs}=\vec{\phi}^{\>\dagger}_{\rm
gs}\sum_{j\neq{\rm gs}}\frac{\langle\psi_{j}|\hat{h}_{i}|\psi_{\rm
gs}\rangle}{E_{\rm gs}-E_{j}}|\psi_{j}\rangle\;,$ (S-49)

where the sum runs over states in the ground quantum number subspace
($N,\sigma$) of the isolated molecule only.

These expressions allow us to compute accurately and cheaply the analytic
gradient of the effective model parameters from the SWT, and hence the
analytic gradient of any optimization loss function that depends on them. In
turn, we may perform GD optimization of our effective model efficiently.

In the main text we demonstrate this with the example of a 4-site molecular
cluster given by Eq. S-8 using the loss function Eq. S-45. Starting from a
random guess for the parameters $\vec{\theta}_{0}$, we rapidly converge to a
solution satisfying the QI-2CK conditions, as shown in Fig. 5 of the main
paper. The parameters used in Fig. 5 resulting from our optimization procedure
are $t_{11}\simeq 0.697$, $t_{22}\simeq 0.642$, $t_{33}\simeq 0.480$,
$t_{44}\simeq 0.253$, $V_{g}=0.42$, $t_{12}\simeq 0.102$, $t_{13}\simeq
0.053$, $t_{14}\simeq 0.1$, $t_{23}\simeq 0.109$, $t_{24}\simeq 0.046$,
$t_{34}\simeq 0.169$.

###  S-VII.2 Inverse design using NRG

A more sophisticated treatment is provided using NRG which is a non-
perturbative method [9], and does not rely on the faithful derivation of a
low-energy effective model [6], but can work directly with the bare model (at
least for modest molecular complexity [5]). Inverse design of molecular
junctions using NRG affords the possibility of engineering loss functions for
the optimization that work directly with the physical observables of interest.

For example, in the present case of the QI-2CK effect, the quantum critical
point is associated with a vanishing $T^{*}$ scale. This scale can be read off
from the temperature dependence of the entropy, as shown in the main paper and
in Figs. S-5, S-6. Defining $T^{*}$ through the relation $S_{\rm
mol}(T=T^{*})=\ln(2)/4$ and a loss function $\mathcal{L}(\vec{\theta})=T^{*}$
allows us to use ML methodologies to optimize the model parameters
$\vec{\theta}$ to minimize $T^{*}$ by calculation of the molecule entropy.
This is the approach we adopt in this work to find the critical value of
$t^{\prime}$ in the 5-site system. However, alternative approaches can be
used. Fig. 4 of the main paper shows that the QI-2CK critical point is also
associated with a conductance node at $T=0$. Therefore an alternative loss
function might be $\mathcal{L}(\vec{\theta})=G_{c}(T=0)$.

A numerical (finite-difference) approximation to the gradient of the chosen
loss function can be computed by running several standard NRG runs. NRG can
then be integrated within a standard GD routine for inverse design model
optimization. For constrained models with only a few optimization parameters,
this might be sufficient.

However, for more complex problems (especially those for which each NRG run is
computationally expensive and optimization efficiency is thus a priority), it
is clearly advantageous to have the analytic gradient of the loss function
from the NRG calculation. This is more challenging than for the SWT
calculation because the NRG algorithm itself is much more complex. Also the
loss function might comprise physical observables rather than simply effective
model parameters. Fortunately this problem has been overcome recently by using
differentiable programming techniques [18] with the NRG method. This
‘automatic differentiable NRG’ method allows exact gradients of NRG
observables to be calculated quickly and efficiently [19]. This brings inverse
design of molecular junctions using an NRG solver within reach.

##  S-VIII Non-interacting limit

To better understand the QI-2CK critical point of the 5-site molecular
junction studied above, here we consider the non-interacting $U\to 0$ limit of
the model,

$\begin{split}H_{\rm
mol}=&t\sum_{\sigma}\left(d_{1\sigma}^{\dagger}d_{2\sigma}^{\phantom{\dagger}}+d_{2\sigma}^{\dagger}d_{3\sigma}^{\phantom{\dagger}}+d_{3\sigma}^{\dagger}d_{4\sigma}^{\phantom{\dagger}}+d_{4\sigma}^{\dagger}d_{5\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\\\
&+t^{\prime}\sum_{\sigma}\left(d_{1\sigma}^{\dagger}d_{4\sigma}^{\phantom{\dagger}}+d_{2\sigma}^{\dagger}d_{5\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\\\ H_{\rm
hyb}=&V\sum_{\sigma}\left(c_{s\sigma}^{\dagger}d_{1\sigma}^{\phantom{\dagger}}+c_{d\sigma}^{\dagger}d_{5\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\end{split}$ (S-50)

which possesses $sd$-symmetry and particle-hole symmetry. Since $U=0$ the
$\sigma=\uparrow$ and $\downarrow$ sectors are completely independent and can
be treated separately. The two sectors are identical due to $SU(2)$ spin
symmetry. Here we focus on the $\sigma=\uparrow$ sector for concreteness.

Figure S-7: Through-molecule Green’s function $G_{15}(\omega)$ for the 5-site
molecular junction in the non-interacting ($U=0$) limit, plotted for $t=1$ and
$\Gamma_{s}=\Gamma_{d}=0.1$ for different $t^{\prime}$. From Eq. S-56 the
zero-temperature series conductance $G_{c}(0)=2e^{2}/h$ is seen to take the
maximum value for a single electron transistor for any $t^{\prime}\neq t$, but
precisely vanishes due to a single-particle QI node at $t^{\prime}=t$.

The first question is: what is the analogue of the QI-2CK critical point when
$U=0$? The answer is provided by the SWT calculations shown in Fig. 2 of the
main paper. For small but finite $U/t$ we see that the critical value of
$t^{\prime}$ tends to the value of $t$. This is supported by full NRG
calculations. Thus, we deduce that $t^{\prime}/t\to 1$ as $U/t\to 0$. Below,
we therefore examine Eq. S-50 in the vicinity of the point $t^{\prime}/t=1$.
Indeed, we find that $t^{\prime}/t=1$ is a special point.

We start by performing a rotation to an even/odd orbital basis, viz:

$\displaystyle d_{1\begin{subarray}{c}e\\\ o\end{subarray}\sigma}$
$\displaystyle=\tfrac{1}{\sqrt{2}}[d_{1\sigma}\pm\ d_{5\sigma}]\;$ (S-51)
$\displaystyle d_{2\begin{subarray}{c}e\\\ o\end{subarray}\sigma}$
$\displaystyle=\tfrac{1}{\sqrt{2}}[d_{2\sigma}\pm\ d_{4\sigma}]\;$ (S-52)
$\displaystyle c_{\begin{subarray}{c}e\\\ o\end{subarray}\sigma}$
$\displaystyle=\tfrac{1}{\sqrt{2}}[c_{s\sigma}\pm\ c_{d\sigma}]\;.$ (S-53)

In this basis, the model exactly decouples into even and odd sectors,

$\begin{split}H_{\rm
mol}=&\sum_{\sigma}\left[(t+t^{\prime})d_{1e\sigma}^{\dagger}d_{2e\sigma}^{\phantom{\dagger}}+\sqrt{2}td_{2e\sigma}^{\dagger}d_{3\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right]\\\
&~{}~{}~{}+(t-t^{\prime})\sum_{\sigma}\left[d_{1o\sigma}^{\dagger}d_{2o\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right]\\\ H_{\rm
hyb}=&V\sum_{\sigma}\left(c_{e\sigma}^{\dagger}d_{1e\sigma}^{\phantom{\dagger}}+c_{o\sigma}^{\dagger}d_{1o\sigma}^{\phantom{\dagger}}+{\rm
H.c.}\right)\end{split}$ (S-54)

Importantly, we see that the orbital $d_{2o\sigma}$ exactly decouples from the
rest of the system when $t^{\prime}=t$, and sits at zero energy. This directly
implies that the $T=0$ contribution from the molecule to the total entropy of
the system is $S_{\rm mol}(0)=\ln(4)$ at this special point, since all other
molecule degrees of freedom are screened by the leads.

By contrast, we know from our NRG results at small but finite $U/t$ where the
critical point arises at $t^{\prime}\simeq t$ that the $T=0$ residual molecule
entropy is $\tfrac{1}{2}\ln(2)$. The critical point is therefore a genuine
non-Fermi liquid (NFL) in the sense that it is not adiabatically connected to
the non-interacting limit. However, the energy scale $T_{\rm K}$ corresponding
to the formation of the NFL state also vanishes as $U/t\to 0$ as argued above.
We speculate that it is the decoupled single-particle state $d_{2o\sigma}$ of
the non-interacting molecule that seeds the NFL state when interactions are
turned on.

Finally, we consider dynamical properties of the $U=0$ molecule. The molecule
Green’s functions of the lead-coupled system can be simply determined,

$\bm{G}_{\rm
mol}^{U=0}(\omega)=\begin{bmatrix}\omega^{+}+i\Gamma&~{}-t~{}&0&~{}-t^{\prime}~{}&0\\\
-t&\omega^{+}&~{}-t~{}&0&-t^{\prime}\\\ 0&-t&\omega^{+}&-t&0\\\
-t^{\prime}&0&-t&\omega^{+}&-t\\\
0&-t^{\prime}&0&-t&\omega^{+}+i\Gamma\end{bmatrix}^{-1}$ (S-55)

where $\left[\bm{G}_{\rm mol}\right]_{ij}=\langle\langle
d_{i\sigma}^{\phantom{\dagger}};d_{j\sigma}^{\dagger}\rangle\rangle$,
$\omega^{+}=\omega+i0^{+}$ and $\Gamma=\pi V^{2}\rho_{0}$, with $\rho_{0}$ the
lead density of states (we have assumed the wide-flat-band limit for
simplicity).

Of particular interest is the through-molecule Green’s function that connects
the source and drain leads, $G_{15}(\omega)$, since this controls series
transport in the non-interacting limit via the Landauer formula [20, 12],

$G_{c}(T)=\left(\frac{2e^{2}}{h}\right)\int d\omega\frac{-\partial
f(\omega)}{\partial\omega}\times|2\Gamma G_{15}(\omega)|^{2}\;,$ (S-56)

where $G_{c}(T)$ is the series electrical conductance $dI_{sd}/dV_{b}$
resulting when a current $I_{sd}$ flows due to a bias voltage $V_{b}$ between
source and drain leads in linear response, and
$f(\omega)=1/[1+\exp(\omega/T)]$ is the equilibrium fermi function. The $T=0$
conductance then follows as $G_{c}(0)=(2e^{2}/h)\times|2\Gamma G_{15}(0)|^{2}$
and depends on the fermi level behavior of the through-molecule Green’s
function $G_{15}(0)$. Analysis of Eq. S-55 shows that,

$2\Gamma G_{15}^{U=0}(\omega\to 0)=\begin{cases}-i&~{}~{}:~{}t^{\prime}\neq
t\\\ 0&~{}~{}:~{}t^{\prime}=t\end{cases}$ (S-57)

The full energy-dependence is shown in Fig. S-7, and establishes that for
$t^{\prime}\neq t$ a narrow peak pinned at the fermi energy to $|2\Gamma
G_{15}(0)|^{2}=1$ collapses at $t^{\prime}=t$. Thus the $T=0$ conductance
$G_{c}(0)=2e^{2}/h$ for any $t^{\prime}\neq t$ but has an exact conductance
node $G_{c}(0)=0$ for $t^{\prime}=t$. This is a result of the single-particle
QI in the system at $t^{\prime}=t$, giving rise to destructively interfering
paths through the molecule. We note that in the interacting case $U>0$ the
$T=0$ conductance also saturates to its maximum $G_{c}(0)=2e^{2}/h$ away from
the QI-2CK critical point $t^{\prime}\neq t^{\prime}_{c}$, but series
conductance vanishes $G_{c}(0)=0$ at the critical point
$t^{\prime}=t^{\prime}_{c}$. However, in this case the conductance node is a
true many-body QI effect.

## References

  * Lafferentz et al. [2009] L. Lafferentz, F. Ample, H. Yu, S. Hecht, C. Joachim, and L. Grill, Conductance of a single conjugated polymer as a continuous function of its length, Science 323, 1193 (2009).
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Relation between the anderson and kondo hamiltonians, Physical Review 149, 491 (1966).
  * Hewson [1993] A. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, Cambridge, England, 1993).
  * Pedersen et al. [2014] K. G. Pedersen, M. Strange, M. Leijnse, P. Hedegård, G. C. Solomon, and J. Paaske, Quantum interference in off-resonant transport through single molecules, Physical Review B 90, 125413 (2014).
  * Mitchell et al. [2017] A. K. Mitchell, K. G. Pedersen, P. Hedegård, and J. Paaske, Kondo blockade due to quantum interference in single-molecule junctions, Nature communications 8, 1 (2017).
  * Rigo and Mitchell [2020] J. B. Rigo and A. K. Mitchell, Machine learning effective models for quantum systems, Physical Review B 101, 241105 (2020).
  * Chan and Gulácsi [2004] R. Chan and M. Gulácsi, The exact schrieffer–wolff transformation, Philosophical Magazine 84, 1265 (2004).
  * Haldane [1978] F. Haldane, Theory of the atomic limit of the anderson model. i. perturbation expansions re-examined, Journal of Physics C: Solid State Physics 11, 5015 (1978).
  * Bulla et al. [2008] R. Bulla, T. A. Costi, and T. Pruschke, Numerical renormalization group method for quantum impurity systems, Reviews of Modern Physics 80, 395 (2008).
  * Coulson and Rushbrooke [1940] C. Coulson and G. Rushbrooke, Note on the method of molecular orbitals, in Mathematical Proceedings of the Cambridge Philosophical Society, Vol. 36 (Cambridge University Press, 1940) pp. 193–200.
  * Tsuji et al. [2018] Y. Tsuji, E. Estrada, R. Movassagh, and R. Hoffmann, Quantum interference, graphs, walks, and polynomials, Chemical reviews 118, 4887 (2018).
  * Minarelli et al. [2022] E. L. Minarelli, J. B. Rigo, and A. K. Mitchell, Linear response quantum transport through interacting multi-orbital nanostructures, arXiv preprint arXiv:2209.01208 (2022).
  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, Journal of the physical society of Japan 12, 570 (1957).
  * Affleck and Ludwig [1993] I. Affleck and A. W. Ludwig, Exact conformal-field-theory results on the multichannel kondo effect: Single-fermion green’s function, self-energy, and resistivity, Physical Review B 48, 7297 (1993).
  * Sela et al. [2011] E. Sela, A. K. Mitchell, and L. Fritz, Exact crossover green function in the two-channel and two-impurity kondo models, Phys. Rev. Lett. 106, 147202 (2011).
  * Mitchell and Sela [2012] A. K. Mitchell and E. Sela, Universal low-temperature crossover in two-channel kondo models, Phys. Rev. B 85, 235127 (2012).
  * Feynman [1939] R. P. Feynman, Forces in molecules, Physical review 56, 340 (1939).
  * Bartholomew-Biggs et al. [2000] M. Bartholomew-Biggs, S. Brown, B. Christianson, and L. Dixon, Automatic differentiation of algorithms, Journal of Computational and Applied Mathematics 124, 171 (2000).
  * Rigo and Mitchell [2022] J. B. Rigo and A. K. Mitchell, Automatic differentiable numerical renormalization group, Physical Review Research 4, 013227 (2022).
  * Landauer [1957] R. Landauer, Spatial variation of currents and fields due to localized scatterers in metallic conduction, IBM Journal of research and development 1, 223 (1957).