\section*{Appendices}
In these appendices we provide supporting documents for some of the statements in the paper. Appendix \ref{refA} shows that the problem of $N$ multi-channel Kondo impurities immersed in 2D and 3D can be reduced to a set of one-dimensional problems, including the superconducting proximity pairing term. Appendix \ref{refB} summarizes our convention for the Klein factors that are used in bosonization.
\subsection{Dimensional reduction}\label{refA}
For completeness, here we show how a problem of $M$ impurities in two dimension can be reduced to a sum of 1D metallic systems terminated at the spin impurities, appropriate for field theory analysis. This is a generalization of $M=2$ case from Ref.\,\cite{Affleck95}. We start from a two-dimensional system described by $H=H_0+H_K+H_\Delta$ defined in Eqs.\,\ref{eq1} and \ref{eq2}, and generalize it to the case that there are $M$ Kondo impurities located at positions $\vec d_n$ coupled to the conduction band.

\subsubsection{Kondo coupling}
We consider $M$ impurities located on a 2D plane at positions $\vec d_n$ shown in Fig.\,3. We have
\begin{equation}
H_K=\int{\frac{d^2kd^2k'}{(2\pi)^4}}\psi^{\dagger}_{\vec k}\frac{\vec\sigma}{2}\psi^{\vphantom{\dagger}}_{k'}\cdot\sum_{n=1}^MJ_n\vec S_n e^{-i(\vec k-\vec k').\vec d_n}.
\end{equation}
We start by defining the fermions
\begin{equation}
\psi_{nE}=\int{\frac{d^2k}{(2\pi)^2}}\delta(\epsilon_k-E)e^{i\vec k\cdot\vec d_n}\psi^{\vphantom{\dagger}}_{\vec k},
\end{equation}
in terms of which the Kondo Hamiltonian is
\begin{equation}
H_K=\int{dEdE'}\sum_{n=1}^MJ_n\psi^{\dagger}_{nE}\frac{\vec\sigma}{2}\psi^{\vphantom{\dagger}}_{nE'}\cdot\vec S_n.
\end{equation}
The problem is that $\psi_{nE}$ do not obey standard anticommutation relations. Rather
\begin{equation}
\acom{\psi^{\vphantom{\dagger}}_{nE}}{\psi^{\dagger}_{mE'}}=g_{nm}(E)\delta(E-E'),
\end{equation}
where
\begin{equation}
g_{nm}(E)=\int{\frac{d^2k}{(2\pi)^2}}\delta(\epsilon_k-E)e^{i\vec k\cdot(\vec d_n-\vec d_m)}.
\end{equation}
Representing the vector $\vec d_{nm}\equiv \vec d_n-\vec d_m$ in polar coordinates by $\vec d_{nm}=d_{nm}(\cos\phi_{nm},\sin\phi_{nm})$ we can write
\begin{eqnarray}
g_{nm}(E)&=&\frac{k_E}{2\pi\partial_k\epsilon_E}\int_0^{2\pi}\frac{d\phi}{2\pi}e^{ik_Ed_{nm}\cos(\phi-\phi_{nm})}\nonumber\\
&=&\frac{k_E}{2\pi\partial_k\epsilon_E}\sum_pe^{-ip\phi_{nm}}J_p(k_Ed_{nm})i^p\int_0^{2\pi}\frac{d\phi}{2\pi}e^{ip\phi}\nonumber\\
&=&g_EJ_0(k_Ed_{nm}),
\end{eqnarray}
where
\begin{equation}
g_E=\frac{k_E}{2\pi\partial_k\epsilon_E}.
\end{equation}
Remarkably, $g_{nm}$ depends only on mutual distances of the impurities, measured by the corresponding wavelength $k_Ed_{nm}$. The matrix $g_{nm}$ is real and symmetric and it has 
real eigenvalues $\lambda_n$ and can be diagonalized by orthogonal eigenvectors $\vec u_n$ where $g\vec u_n=\lambda_n\vec u_n$. So we can write $g_{nm}=\sum_p u_{np}\lambda_{p}u^*_{mp}$ and the orthogonality is $\sum_pu_{np}u^*_{mp}=\delta_{nm}$. Defining new operators with
\begin{equation}
\psi_{nE}=\sum_i u_{ni}(E)\sqrt{\lambda_i(E)}\tilde{\psi}_{iE},
\end{equation}
and
\begin{equation}
\tilde{\psi}_{iE}\equiv\frac{1}{\sqrt{\lambda_i(E)}}\sum_n{u^*_{ni}}(E)\psi^{\vphantom{\dagger}}_{nE}
\end{equation}
we find that they are orthonormal
\begin{equation} \acom{\tilde{\psi}_{iE}}{\tilde{\psi}_{jE'}}=\delta_{ij}\delta({E-E'}),
\end{equation}
and the Kondo Hamiltonian becomes
\begin{equation}
H_K=\int{dEdE'}\sum_{nij}J_{nij}(E,E')\vec S_n\cdot \tilde{\psi}^{\dagger}_{iE}\frac{\vec\sigma}{2}\tilde{\psi}^{\vphantom{\dagger}}_{jE'}, 
\end{equation}
where
\begin{equation}
J_{nij}(E,E')=J_n\sqrt{\lambda_i(E)\lambda_j(E')}u_{ni}(E)u^*_{nj}(E').
\end{equation}
Although the $J_{nij}$ couplings are complex in general, the relation $J^*_{nij}(E,E')=J_{nji}(E',E)$ ensures the hermiticity of the Hamiltonian. Whether or not a spin $n$ couples the channels $i$ and $j$ depend on the product of wavefunctions of $i$ and $j$ at the site $n$, which can be tuned by moving potential scatterings induced by scanning tips (Fig.\,3). Next, we do a Taylor expansion of $J_{nij}(E,E')$ function around Fermi energy and keep only the leading relevant term. 
$H_K$ can be written in the matrix form $[\bb J_n]_{ij}=J_{nij}$,
\begin{equation}
H_K=\sum_n  \Psi^{\dagger}(0)[\bb J_n\vec S_n\cdot\frac{\vec \sigma}{2}]\Psi^{\vphantom{\dagger}}(0).
\end{equation}
where $\Psi$ is a vector in $m=1...M$ index and the spin is implicit. As we see the impurities talk to all the channels and scatter electron between all the channels. A result of the Taylor expansoin of $J_{nij}(E,E')$ is that RKKY interactions will be induced between the spins
\begin{equation}
H_K\to H_K+\sum_{ij}J_H^{ij}\vec S_i\cdot\vec S_j.
\end{equation}
\subsubsection{Mapping to 1D}
The mapping to (unfolded) left-movers follows Ref.\,\cite{Affleck95}
\begin{equation}
\psi_L(x,t)=\frac{1}{\sqrt v}\int_{-D}^D dE e^{-iE(t+x/v)}\psi_E,
\end{equation}
for $x\in(-\infty,+\infty)$ 
with commutation relations
\begin{equation}
\{\psi^{\dagger}_L(x),\psi^{\vphantom{\dagger}}_L(y)\}=2\pi\delta(x-y).
\end{equation}
Alternatively, we can work in the folded space of left and right-movers defined as
\begin{equation}
\psi_R(x,t)=\psi_L(-x,t), \qquad (x>0).
\end{equation}
\subsubsection{Kinetic term}
As we saw above, the Kondo interaction involves $M$ 1D conduction bands pulled out of the 2D conduction band. The natural guess for their kinetic energy  is
\begin{eqnarray}
H_0(\vec d)&\equiv&\int{dE E\psi^{\dagger}_{nE}\psi^{\vphantom{\dagger}}_{nE}}\\
&=&\int\frac{d^2kd^2k'}{(2\pi)^4}\epsilon_k\delta(\epsilon_k-\epsilon_k')e^{i(\vec k'-\vec k)\cdot\vec  d}\psi^{\dagger}_{\vec k}\psi^{\vphantom{\dagger}}_{\vec k'}.\nonumber
\end{eqnarray}
Summing over all the positions we have
\begin{eqnarray}
\int{d^2r}H_0(\vec r)
&=&\delta(0)\int\frac{d^2kd^2k'}{(2\pi)^2}\epsilon_k\psi^{\dagger}_{\vec k}\psi^{\vphantom{\dagger}}_{\vec k},
\end{eqnarray}
which is the total Hamiltonian up to a $\delta(0)$, i.e. the total volume. Transformation from $\psi_{nE}$ to $\tilde{\psi}_{iE}$ is a unitary transformation at each energy and it doesn't change the kinetic part. Therefore, we have
\begin{equation}
H_0\to\sum_{n=1}^M\int{dE E\tilde{\psi}^{\dagger}_{nE}\tilde{\psi}^{\vphantom{\dagger}}_{nE}}
\end{equation}
\subsubsection{Pairing term}
Without loss of generality, we consider a channel-diagonal s-wave singlet pairing
\begin{equation}
H_\Delta=\int{\frac{d^2k}{(2\pi)^2}}(\Delta\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{-k\downarrow}+h.c.).
\end{equation}
Let us look at 
\begin{eqnarray}
H_\Delta^+(\vec d_n)&\equiv&\Delta\int{dE}\psi^{\dagger}_{nE\uparrow}\psi^{\dagger}_{n,E\downarrow}\nonumber\\
&&\hspace{-1.75cm}=\Delta\int{dE}\int{\frac{d^2kd^2k'}{(2\pi)^4}}\delta(\epsilon_k-E)\delta(\epsilon_{k'}-E)\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{q\downarrow}e^{-i(\vec k+\vec k')\cdot \vec d_n}\nonumber
\end{eqnarray}
Summing over all the points gives a $(2\pi)^2\delta(\vec k+\vec k')$ term that  gives
\begin{equation}
\int{d^2r}H^+_\Delta(\vec r)=\delta(0)\Delta\int{\frac{d^2k}{(2\pi)^2}}\psi^{\dagger}_{k\uparrow}\psi^{\dagger}_{-k\downarrow},
\end{equation}
which is our initial Hamiltonian up to the system volume. Therefore,
\begin{eqnarray}
H_{\Delta}&\to &\Delta\sum_{n=1}^M\int{dE}(\psi^{\dagger}_{nE\uparrow}\psi^{\dagger}_{nE\downarrow}+h.c.)\nonumber\\
&=&\Delta\sum_{ij}\int{dE}(\tilde{\psi}^{\dagger}_{iE\uparrow}\tilde{\psi}^{\dagger}_{jE\downarrow}+h.c.).
\end{eqnarray}
\subsection{Klein factors}\label{refB}
Here, we briefly review the treatment of Klein-factors following Ref.\,\cite{Zarand00}. Denoting the total number of fermions $\mu$ with
\begin{equation}
N_\mu=\intoinf{dx}c^{\dagger}_{\mu}(x)c^{\vphantom{\dagger}}_{\mu}(x), \quad \mu\in\{1\uparrow,1\downarrow,2\uparrow,2\downarrow\},
\end{equation} 
the Klein factors have the property (no summation):
\begin{eqnarray}
&&F^{\vphantom{\dagger}}_\mu F^{\dagger}_\mu=F^{\dagger}_\mu F^{\vphantom{\dagger}}_\mu=1\\
&&\{F^{\vphantom{\dagger}}_\mu,F^{\dagger}_{\nu}\}=2\delta_{\mu\nu}\\
&&\{F^{\vphantom{\dagger}}_\mu,F^{\vphantom{\dagger}}_\nu\}=0, \quad \forall \mu\neq \nu\\
&&\com{F_\mu}{N_\nu}=\delta_{\mu\nu}F_\nu\\
&&\com{F_\mu}{\phi_\nu(x)}=\com{F_\mu}{\theta_\nu(x)}=0.
\end{eqnarray}
In the paper we transform into collective boson basis:
\begin{eqnarray}
\mat{{\cal N}_c \\ {\cal N}_s \\ {\cal N}_f \\ {\cal N}_{sf}}=\frac{1}{2}\matc{rrrr}
{1 & 1 & 1 & 1 \\
1 & -1 & 1 & -1 \\ 
1 & 1 & -1 & -1 \\ 
1 & -1 & -1 & 1}
\mat{N_{1\uparrow} \\ N_{1\downarrow} \\N_{2\uparrow} \\ N_{2\downarrow}}.
\end{eqnarray}
In accord with this transformation, we define new Klein factors ${\cal F}_\mu$ for $\mu=c,s,f,sf$ that satisfy 
\begin{eqnarray}
&&{\cal F}^{\vphantom{\dagger}}_\mu {\cal F}^{\dagger}_\mu={\cal F}^{\dagger}_\mu {\cal F}^{\vphantom{\dagger}}_\mu=1\\
&&\{{\cal F}^{\vphantom{\dagger}}_\mu,{\cal F}^{\dagger}_{\nu}\}=2\delta_{\mu\nu}\\
&&\{{\cal F}^{\vphantom{\dagger}}_\mu,{\cal F}^{\vphantom{\dagger}}_\nu\}=0, \quad \forall \mu\neq \nu\\
&&\com{{\cal F}_\mu}{{\cal N}_\nu}=\delta_{\mu\nu}{\cal F}_\nu\\
&&\com{{\cal F}_\mu}{\phi_\nu(x)}=\com{{\cal F}_\mu}{\theta_\nu(x)}=0.
\end{eqnarray}
It can be seen that these any off-diagonal relation combination, like $F^{\dagger}_\mu F^{\vphantom{\dagger}}_\nu$ change the population of the collective modes ${\cal N}_\zeta$ by integers and thus can be represented by a combination of new Klein factors\,\cite{Zarand00}. For example,
\begin{equation}
F^{\dagger}_{1\uparrow}F^{\vphantom{\dagger}}_{1\downarrow}={\cal F}^{\dagger}_{sf}{\cal F}^{\dagger}_s, \qquad F^{\dagger}_{2\uparrow}F^{\vphantom{\dagger}}_{2\downarrow}={\cal F}_{sf}{\cal F}^{\dagger}_s.
\end{equation} 
In the paper, we repeatedly use these and similar relations.


