\section*{Supplementary Material}

\subsection{A. Renormalization group analysis}
We begin by reiterating that our RG procedure is carried out for the Hamiltonian defined in Eq.~\ref{eq:ham}. 
The effective actions 
$S$ and $S_{\chi}$ are introduced to describe the next-order fluctuations and calculate the correlations functions. 

The relevant correlation functions are:
\ba 
&&\langle \sum_\mu [\phi^\mu(\bm{r},\tau=0)]^2 \rangle _{imp}  \, ,
\nonumber \\
&&\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle _{imp} \, .
\label{eq:correlation}
\ea  
Here, $\langle \cdot \rangle _{imp}$ means we subtract the free-bulk contributions, and the bulk system has spacetime dimension $d=4-\epsilon$. Before performing calculations, we note that the exact formula of above correlations functions are 
\ba 
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}=\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }
\langle \sum_\mu S^\mu(\tau) S^\mu(\tau') \rangle_{imp} \nonumber \\
&& -\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} = K
\bigg\langle 
\bigg[ 
\delta_{\tau,\tau'}\delta_{k,k'}(\partial_{\tau'}  + \epsilon_{\bm{k}}) 
+ \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\sigma^\mu\bigg)
\bigg] 
\bigg\rangle_{imp} -G_0(\bm{r}_0,\tau),
\label{eq:corr_exact}
\ea  
where  
\ba 
h =-4A_\phi g^2\Gamma(-1+\epsilon)\sin(\frac{\epsilon \pi}{2}) \approx  [2\pi A_\phi +O(\epsilon) ]g^2.
\ea
$A_\phi$ is the normalization factor of boson propagator
and $G_0(\bm{r}_0,\tau)$ is the free fermion propagator. Above equations suggest the fermion and boson correlation functions can be represented as functional of $n$-point correlation functions of $S^\mu$ fields. Then, it's sufficient to evaluate the correlation functions of $S^\mu$. To do so, we first derive the effective action of $S^\mu$ by integrating out fermions and bosons:
\ba 
S_{eff}
=-iSS_B 
-K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg).
\ea 
Using $S_{eff}$ and Eq.~\ref{eq:corr_exact}, we're able to calculate fermion and boson correlation functions. The detailed calculations are shown in Supplementary Material (SM) B and the results are
\ba 
\langle \sum_\mu [ \phi^\mu(\bm{r},\tau=0)]^2\rangle_{imp}
&\propto &  \frac{1}{|\bm{r}|^2}
h\bigg[ 
1- 
\frac{1}{S\pi r}\int_0^J\frac{F(j)}{j}dj +\frac{1}{S\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
\bigg] \nonumber \\ 
\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} &\propto&
\frac{J^2\text{sgn}(\tau)}{(1+J^2)|\tau|}
\bigg[ 
1 
+ \frac{-1-J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj 
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&&
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj -\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg( 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)} 
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) 
\bigg] . \nonumber \\
 \label{eq:corr}
\ea 


We are now in the position to carry out renormalization group study. We introduce the following relation between renormalized couplings $J,h$ and bare couplings $J_B,h_B$:
\begin{eqnarray}
&&h_B = Z_h\mu^{\epsilon}h \nonumber \\
&&J_B = Z_J\mu^{r}J 
\label{eq:bare_coupling_supp}
\end{eqnarray}
where $Z_J,Z_h$ are coupling renormalization factor, and $\mu$ is a renormalization energy scale. 
Without loss of generality, we let 
\ba 
Z_h = 1 + \frac{1}{S}x_h \nonumber \\
Z_J = 1 + \frac{1}{S}x_J \nonumber.
\ea 
We determine $x_J$ and $x_h$ by requiring all the poles as a function of $r,\epsilon$  in Eq.~\ref{eq:corr} are minimally removed. This leads to \ba 
x_h &=& \frac{1}{\pi r}\int_0^J\frac{F(j)}{j}dj 
-\frac{1}{\pi \epsilon }
\bigg[ 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg] +O(\frac{1}{S})\nonumber\\
x_J&=& \frac{1}{2\pi r}\int_0^J\frac{F(j)}{j}dj 
-\frac{1}{\pi r}\int_0^J\frac{jF(j)}{1+j^2} dj
+\frac{1}{r\pi}\int_0^J\frac{H(j)}{1+j^2}dj +\frac{-1+J^2}{2(1+J^2)\pi \epsilon}\bigg( 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
-\frac{iJ}{(1+J^2)\pi \epsilon} 
\bigg( 
-\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) +O(\frac{1}{S}) \nonumber \\ 
\ea 

Finally, we take $\mu$ derivative of Eq.~\ref{eq:bare_coupling_supp} at fixed $h_B$ and $J_B$ and obtain the following beta functions: 
\ba 
\beta_J &=&  \frac{J-J^3}{2S\pi (1+J^2)}F(J) + \frac{J^2}{S\pi(1+J^2)}H(J) +\bigg[ \frac{-(1-J^2)Jh+2iJ^2h}{2S(1+J^2)\pi c(J)(c(J)+h)}+\text{h.c.}
\bigg] +o(\frac{1}{S}) \nonumber  \, , \\
\beta_h& =& -\frac{1}{S}\tilde{\epsilon} h + \frac{h}{S\pi}F(J)
-\frac{1}{ S\pi }\frac{(c(J)+c(J)^*)h}{|c(J)|^2} 
-\frac{1}{S\pi }\frac{-2h^2-(c(J)+c(J)^*)h}{(c(J)+h)(c(J)^*+h)} +o(\frac{1}{S})  \, .
\ea 
where we also set $r=0$ at the last step.

Note that an alternative way to perform renormalization group study is to start from the action in Eq.~\ref{eq:action} and then expand in powers of $\chi$ fields~\cite{nahum2022fixed}. Combining the quadratic order term $S_{\chi}$ and fourth-order term $S_4$, we could find an interacting field theory of $\chi$: $S_{\chi}+S_4$. Then we are able to perform RG analysis based on this model. The resulting beta functions are expected to be consistent with current approach. 
{We reserve such an analysis for a future study.}


\subsection{B. Calculations of correlation functions}
Here we present the detailed calculation of Eq.~\ref{eq:corr_exact} using the effective action $S_{eff}
$.
We first introduce the spinor representation of $S^\mu$ with $z^\dag \frac{\sigma^\mu}{2} z$ and let $z=z_0+\delta z$ (see Eq.~3). This gives us
\ba 
S^\mu &= & S\delta_{\mu,z} + \delta S^\mu \nonumber \\
\delta S^x &=& \frac{i\sqrt{S}}{\sqrt{2}}(\chi^\dag -\chi) + O(\frac{1}{S^{1/2}})\nonumber \\
\delta S^y &=& \frac{\sqrt{S}}{\sqrt{2}}(\chi^\dag +\chi) + O(\frac{1}{S^{1/2}}) \nonumber  \\
\delta S^z &=& -|\chi|^2+O(\frac{1}{S}) 
. 
\label{eq:s_exp}
\ea 
We expand the action in Eq.~\ref{eq:action} to the quadratic order of $\chi$ fields:
\ba 
S_{eff}
&=& S_0 +S_\chi \nonumber \\
S_0&=&-\frac{Sg^2}{4}\int_{\tau,\tau'} \frac{\sum_\mu(z_0^\dag \sigma^\mu z_0)^2}{2}G_\phi(\tau-\tau')
-K\text{Tr}\log\bigg( 
\delta_{\tau,\tau'}\delta_{k,k'}(\partial_{\tau'}  + \epsilon_{\bm{k}}) 
+ \delta_{\tau,\tau'}J_K\sum_\mu (z_0^\dag \sigma^\mu z_0)\sigma^\mu\bigg) ,
\label{eq:simp_exp}
\ea 
with $S_0$ the saddle point contribution at order of $S$ and $S_\chi$ describing fluctuations at order $S^0$. We now derive the formula of $S_\chi$. It has contribution from three part: the Berry phase, bosonic coupling and Kondo
interaction.
The Berry phase term 
gives  
\ba 
\int i\Omega|\chi(\Omega)|^2 \frac{d\Omega}{2\pi} \, ,
\label{eq:berry_chi}
\ea 
{while}
 the bosonic coupling 
yields
\ba 
&&-\frac{g^2}{2}\int_{\tau,\tau'}G_\phi(\tau-\tau')[\chi^\dag(\tau)\chi(\tau') + \chi(\tau)\chi^\dag(\tau') -|\chi(\tau)|^2 -|\chi(\tau')|^2] \nonumber \\
&=&\int h|\Omega|^{1-\epsilon} |\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} \, ,
\label{eq:bose_chi}
\ea 
and, finally, the
 Kondo coupling gives (by expanding the fermion determinant)
\ba 
&&\frac{\kappa J_K^2}{4} \int_{\tau,\tau'}\bigg[ G_s(\tau-\tau')G_s(\tau'-\tau) -G_a(\tau-\tau')G_a(\tau'-\tau)\bigg]
\bigg[
\chi^\dag(\tau)\chi(\tau') + \chi^\dag(\tau')\chi(\tau) - |\chi(\tau)|^2 - |\chi(\tau')|^2
\bigg]\nonumber \\
&&+\frac{\kappa J_K^2}{2}\int_{\tau,\tau'}\bigg[ 
G_s(\tau-\tau')G_a(\tau'-\tau)
\bigg] 
\bigg[ 
-\chi^\dag(\tau)\chi(\tau') + \chi(\tau)\chi^\dag(\tau')
\bigg] \nonumber \\
&=&\int \frac{2\kappa J^2 |\Omega|^{1-2r} }{\pi(1+J^2 |\Omega|^{-2r})^2}|\chi(i\Omega)|^2 \frac{d\Omega}{2\pi}+\int \frac{i\Omega \kappa }{\pi }
\bigg[ 
\frac{J|\Omega|^{-r} -J^3 |\Omega|^{-3r}}{(1+J^2|\Omega|^{-2r})^2} - \arctan(J|\Omega|^{-r})
\bigg] |\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} ,
\label{eq:fermi_chi}
\ea 
where $G_s(\tau) = \sum_\sigma G_\sigma(\tau)$ and $G_a(\tau) = \sum_\sigma \sigma G_\sigma(\tau)$ and $G_\sigma(\tau)$ denotes the conduction electron Green's function at leading order $(S^1)$. The explicit formula of these functions are
\ba 
G_s(i\omega) &= \frac{-2i\rho_0\pi |\omega|^{-r}}{1+\rho_0^2\pi^2J_K^2 |\omega|^{-2r}} \text{sgn}(\omega) \nonumber \\
G_a(i\omega) &= \frac{-2\rho_0^2\pi^2J_K |\omega|^{-2r}}{1+\rho_0^2\pi^2J_K^2 |\omega|^{-2r}} .
\ea 
In addition, the coefficient before each $|\Omega|^{-r}$ should be a function of $r$. However, we only keep the leading order term ($r^0$) since we would set $r=0$ in the final step. Combining Eqs.~\ref{eq:berry_chi}, \ref{eq:bose_chi} and ~\ref{eq:fermi_chi}, we have the quadratic theory of the $\chi$ fields:
\ba 
&&S_{\chi} = \int  G^{-1}_\chi(i\Omega)|\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} \nonumber \\
&&G_\chi^{-1}(i\Omega) = i\Omega \bigg[ 1+p(J|\Omega|^{-r})\bigg] +|\Omega|\bigg[ q(J|\Omega|^{-r}) + h|\Omega|^{-\epsilon}\bigg].
\label{eq:Schisupp}
\ea 

For future reference, we introduce the series expansion of $G_\chi(i\Omega)$
\ba 
G_\chi(i\Omega) &=& \frac{1}{|\Omega|} 
\frac{1}{i(\text{sgn}(\Omega)+p(J|\Omega|^{-r})) + q(J|\Omega|^{-r}) 
} + \frac{1}{|\Omega|}\sum_{m=1}^\infty \frac{(-h)^m|\Omega|^{-m \epsilon} }
{[i(\text{sgn}(\Omega) +p(J)) + q(J)]^{m+1}}  
\ea 
Since we set $r=0$ in the last step, we drop the $r$ dependence in the summation by setting $r=0$. Including this contribution won't change the beta functions at the leading order. We can further expand the real and imaginary part of $G_\chi(i\Omega)$ as the following 
\ba 
&&G_\chi(i\Omega) + G_\chi(-i\Omega) = 
\sum_{n=0}^\infty a_n J^n |\Omega|^{-1-nr} + \sum_{m=1}^\infty 
(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) 
|\Omega|^{-1-m\epsilon} \nonumber \\
&&G_\chi(i\Omega)-G_\chi(i\Omega) = 
\sum_{n=0}^\infty ib_n J^n|\Omega|^{1-nr} +   \sum_{m=1}^\infty 
(-h)^m \bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) 
|\Omega|^{1-m\epsilon}
\ea 
where $a_n, b_n$ are the Taylor coefficients of $F(j)$ and $H(j)$ at $j=0$ and 
\ba 
F(j) &=& \frac{2q(j)}{(1+p(j))^2+q(j)^2}\nonumber \\ 
H(j) &= &\frac{-2(1+p(j))}{(1+p(j))^2 + q(j)^2}\nonumber \\
c(J) &= & i(1+p(J)) + q(J) \nonumber .
\ea 

Now, we combine Eqs.~\ref{eq:corr_exact}, \ref{eq:s_exp}, \ref{eq:simp_exp} and \ref{eq:Schisupp} and calculate the correlation functions of fermions and bosons. 
We first expand Eq.~\ref{eq:corr_exact} using Eq.~\ref{eq:s_exp}. For bosons, we have
\ba 
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}\nonumber \\ &=&\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }S^2 \nonumber \\
&&+ 
\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }\langle -S\delta S^z(\tau) -S\delta S^z(\tau') + \delta S^x(\tau) \delta S^x(\tau') + \delta S^y(\tau) \delta S^y(\tau') \rangle_{S_\chi}  + O(S^0) \nonumber \\
\ea 
where the last line denotes the contribution of the next-order fluctuations and are evaluated with respect to $S_{\chi}$. 
For the fermionic fields, it gives 
\ba 
&&-\langle  \sum_{a,\sigma} c_{\bm{r_0},a,\sigma}(\tau)
c_{\bm{r_0},a,\sigma}^\dag(0)\rangle \nonumber \\
&=&K\bigg[G_s(\tau) -G_0(\tau) \nonumber \\
&&
+\int_{\tau_1}\sum_\mu \frac{J_K\text{Tr}[G(\tau-\tau_1) \sigma^\mu G(\tau_1)]}{S}
 \langle \delta S^\mu(\tau_1)\rangle_{S_{\chi}}
+\int_{\tau_1,\tau_2}\sum_{\mu,\nu}\frac{J^2_K\text{Tr}[G(\tau-\tau_1)  \sigma^\mu G(\tau_1-\tau_2) \sigma^\nu G(\tau_2)]}{S^2} 
\langle
\delta S^\mu(\tau_1) \delta S^\nu(\tau_2)
\rangle_{S_{\chi}}
\bigg]\nonumber \\
&&+ O(S^0)
\ea 
where the third line denotes the contribution of the next-order fluctuations and are evaluated with respect to $S_{\chi}$ and $G(\tau) = \text{diag}[\frac{G_s(\tau)+G_a(\tau)}{2},\frac{G_s(\tau)-G_a(\tau)}{2}]$. 


To get the final results, we only need to calculate the correlation functions of $\delta S^\mu$ (or $\chi$) that appearing in the above equations. The boson correlation function becomes
\ba 
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}\nonumber \\
&=&\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }
S^2 \nonumber\\
&&\bigg[ 
1 + \frac{1}{S}\langle \chi^\dag(\tau)\chi(\tau') + \chi^\dag (\tau')\chi(\tau') - |\chi(\tau)|^2 -|\chi(\tau')|^2\rangle 
\bigg]+O(S^0)\nonumber \\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h 
+  
\frac{A_\phi Sh}{2\pi }
\int_0^\infty 
2(G_\chi(\Omega)+G_\chi(-\Omega))
\bigg( 
2\pi x^{-(1-\epsilon)}B[\frac{\epsilon-1}{2},x]^2|\Omega|^{2-2\epsilon} - \frac{\pi^2}{|\bm{r}|^2}
\bigg) 
\frac{d\Omega}{2\pi} \nonumber +O(S^0)\\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h 
+  
\frac{A_\phi Sh}{2\pi } 
\int_0^\infty 
2\bigg[ 
\sum_{n=0}^\infty  a_nJ^n|\Omega|^{1-nr}
+\sum_{m=1}^\infty (-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) 
|\Omega|^{-1-m\epsilon} 
\bigg] \nonumber \\
&&\frac{1}{|\bm{r}|^{2}}
\bigg(
2\pi x^{1-\epsilon}B[\frac{\epsilon-1}{2},x]^2-
\pi 2^{-\epsilon}\Gamma((1-\epsilon)/2)^2
\bigg) 
\frac{d\Omega}{2\pi} +O(S^0)\nonumber \\
&=&
\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h 
+  
\frac{A_\phi S h}{\pi |\bm{r}|^2} \bigg[ \frac{1}{r}
\sum_{n=0}^\infty \frac{-a_nJ^n\pi}{2n} 
+\frac{1}{\epsilon}\sum_{m=1}(-h)^m 
\bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) \frac{-\pi}{2n\epsilon} 
+O(r^0,\epsilon^0)\bigg] +O(S^0) \nonumber \\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}
h\bigg\{ 
1- \frac{1}{S}\bigg[\frac{1}{\pi r}\int_0^J\frac{F(j)}{j}dj 
+\frac{1}{\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
+O(\epsilon^0,r^0)\bigg]+o(\frac{1}{S})
\bigg\} 
\label{eq:chi_prop}
\ea 
where $x=|\Omega||\bm{r}|$ and $B(a,x)$ is the modified Bessel functions of the second kind. In the above expression, we only keep track of the poles as a function of $\epsilon$ and $r$. We also drop the $\frac{r}{\epsilon}$ poles which goes to zero as we set $r=0$. 

For the fermion correlation function,
we first perform Fourier transformation and separate the next-order contributions into four parts $I_1,I_2,I_3,I_4$:
\ba 
G_f(i\omega) = -\int_\tau\langle  \sum_{a,\sigma} c_{\bm{r_0},a,\sigma}(\tau)
c_{\bm{r_0},a,\sigma}^\dag(0)\rangle_{imp}
e^{i\omega\tau} = K\bigg[ \frac{2i\pi \rho_0 J^2}{1+J^2} \text{sgn}(\omega)
+I_1+I_2+I_3+I_4+o(\frac{1}{S})\bigg] 
\ea 
Each part is shown below:
\begin{eqnarray*}
I_1 &=& 
\frac{J_K^2}{4S}(G_s(i\omega)^2 +G_a(i\omega)^2)
\int G_\chi(i\Omega) (G_s(i\omega+i\Omega)+G_s(i\omega-i\Omega)-G_s(i\omega)-G_s(i\omega))\frac{d\Omega}{2\pi} \\
&=&\frac{-2i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^2S}
\int_0^\infty \bigg[ \sum_{n=0}^\infty a_nJ^n|\Omega|^{-1-nr} +\sum_{m=1}^{\infty}(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) 
|\Omega|^{-1-m\epsilon}\bigg]\\
&&\sum_{p=0}^\infty (-J^2)^p
\bigg[ \text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r} +
 \text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}
 -2
  \text{sgn}(\omega)|\omega|^{-2pr-r}
\bigg]\frac{d\Omega}{2\pi} +O(r^0)\\
&=&\frac{-i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^2\pi S}\bigg[ \sum_{n=0}^\infty a_nJ^n\sum_{p=0}^\infty (-J^2)^P\frac{-2}{nr} 
+
\sum_{m=1}^\infty (-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) \frac{-2}{m\epsilon} 
\bigg] \text{sgn}(\omega)  \\
&=&\frac{2i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^3\pi S}
\bigg[ 
\frac{1}{r} \int_0^J \frac{F(j)}{j}dj
-\frac{1}{\epsilon}\bigg( 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
+O(r^0,\epsilon^0)\bigg) 
\bigg] \text{sgn}(\omega) .
\end{eqnarray*}

\begin{eqnarray*}
I_2 &=& 
-\frac{J_K^2}{2S}G_s(i\omega)G_a(i\omega)
\int G_\chi(i\Omega)
(G_a(i\omega+i\Omega)+G_a(i\omega-i\Omega)-G_a(i\omega)-G_a(i\omega))\frac{d\Omega}{2\pi} \\
&=&-\frac{2iJ^3(-2J\rho_0\pi)}{(1+J^2)^2S}\text{sgn}(\omega)
\int_0^\infty 
\bigg[ 
\sum_{n=0}^\infty a_nJ^n|\Omega|^{-1-nr} +\sum_{m=1}^{\infty}(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) 
|\Omega|^{-1-m\epsilon}
\bigg] \\
&&
\sum_{p=0}^\infty
(-J^2)^p \bigg[ 
|\omega+\Omega|^{-2pr-2r}
+|\omega-\Omega|^{-2pr-2r}
-2|\omega|^{-2pr-2r}
\bigg] 
\frac{d\Omega}{2\pi}\\
&=&-\frac{2iJ^3(-2J\rho_0\pi)}{2\pi (1+J^2)^2S}\text{sgn}(\omega)
\bigg[ 
\sum_{n=0}^\infty a_n
\sum_{p=0}^\infty (-J^2)^p
J^n(\frac{2}{2pr+2r+nr} - \frac{2}{nr})
+O(r^0,\epsilon^0)
\bigg] 
 \\
&=& \frac{2iJ^4\rho_0\pi}{\pi (1+J^2)^2S}\text{sgn}(\omega)
 \bigg[ 
 \frac{1}{rJ^2}\int_0^J \frac{2jF(j)}{1+j^2}dj 
 -\frac{2}{r(1+J^2)}\int_0^J\frac{F(j)}{j}dj
 +O(r^0,\epsilon^0)\bigg].
\end{eqnarray*}

\begin{eqnarray*}
I_3 &=&\frac{J_K^2}{4S}
\bigg[ 
G_a(i\omega)^2 +G_s(i\omega)^2 
\bigg] \int G_\chi(i\Omega) [
G_a(i\Omega+i\omega) - G_a(i\omega-i\Omega)]
\frac{d\Omega}{2\pi} \\
&=& \frac{J^2(-1+J^2)}{2\pi(1+J^2)^2S} 
\int_0^\infty \bigg[ 
\sum_{n=0}ib_nJ^n|\Omega|^{1-nr} +\sum_{m=1}(-h)^m
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)|\Omega|^{-1-m\epsilon} \bigg]
\\
&&\bigg[ -2J\rho_0\pi \sum_{p=0}^\infty (-J^2)^p 
\bigg( 
|\Omega+\omega|^{-2r-2pr} -
|-\Omega+\omega|^{-2r-2pr}
\bigg) 
\bigg] 
d\Omega \\
&=&\frac{1}{S}[0+O(r^0,\epsilon^0)].
\end{eqnarray*} 

\begin{eqnarray*} 
I_4 &=& 
\frac{-J_K^2}{2}G_a(i\omega)G_s(i\omega)
\int G_\chi(i\Omega) \bigg(
G_s(i\Omega+i\omega)-G_s(i\omega-i\Omega)
\bigg) \frac{d\Omega}{2\pi} \\
&=&\frac{-2J^3\rho_0\pi}{S(1+J^2)^2}\int_0^\infty 
\bigg[ 
\sum_{n=0}ib_nJ^n|\Omega|^{1-nr} +\sum_{m=1}(-h)^m
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)|\Omega|^{-1-m\epsilon} \bigg]\\
&&
\sum_{p=0}(-J^2)^p
\bigg[ 
\text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r}
-\text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}
\bigg] 
\frac{d\Omega}{\pi} \\
&=&\frac{-2J^3\rho_0\pi}{S(1+J^2)^2} 
\bigg[ 
\sum_{n=0}ib_nJ^n\sum_{p=0} (-J^2)^p \frac{2}{\pi(nr+2pr+r)} +\sum_{m=1}(-h)^m\frac{1}{1+J^2}
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)\frac{2}{\pi m\epsilon} \bigg]\text{sgn}(\omega)  \\
&=&\frac{-2J^3\rho_0\pi}{(1+J^2)^2} 
\bigg[ \frac{2}{ r\pi J}\int_0^J\frac{iH(j)}{(1+j^2)}dj
+\frac{2}{\pi \epsilon(1+J^2)}\bigg(
-\frac{\log(1+h/c(J))}{c(J)}
+\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) +O(r^0,\epsilon^0)
\bigg]\text{sgn}(\omega)
.
\\
\end{eqnarray*}
In the above derivation, we only keep terms that diverge as $r\rightarrow 0$ and $\epsilon \rightarrow 0$, with $r$ going to zero first.

Summing over all the contributions, we have 
\ba 
&&G_f(i\omega) \nonumber \\
&=& 2iK\pi \rho_0\text{sgn}(\omega) \frac{J^2}{1+J^2} 
\bigg[ 
1
+ \frac{-1+J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj 
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&& -\frac{2J^2}{(1+J^2)^2\pi Sr}
\int_0^J\frac{F(j)}{j}dj
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj \nonumber \\ 
&&-\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg( 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) 
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) 
\bigg] 
\nonumber \\
\label{eq:g_prop}
\ea 

In summary, we have 
\ba 
\langle \sum_\mu [ \phi^\mu(\bm{r},\tau=0)]^2\rangle_{imp}
&\propto &  \frac{1}{|\bm{r}|^2}
h\bigg[ 
1- 
\frac{1}{S\pi r}\int_0^J\frac{F(j)}{j}dj +\frac{1}{S\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
\bigg] \nonumber \\ 
\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} &\propto&
\frac{J^2\text{sgn}(\tau)}{(1+J^2)|\tau|}
\bigg[ 
1 
+ \frac{-1-J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj 
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&&
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj -\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg( 
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)} 
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) 
\bigg] . \nonumber 
\ea 


\subsection{C. Beta function analysis} 
We numerically evaluate the beta functions and the resulting RG flows are shown in Fig.~\ref{fig:rg_num}.
Besides numerical evaluations, we analyze the beta equations in certain extreme cases. 


At $h=0$, the model becomes a standard Kondo impurity model. The beta function of $J$ is 
\ba 
\beta_J\bigg|_{h=0} &=& \frac{J-J^3}{2S\pi (1+J^2)}F(J) + \frac{J^2}{S\pi(1+J^2)}H(J)  \nonumber \\
&=&-\frac{2J^2}{S\pi} +\frac{2\kappa J^3}{S\pi^2}  +O(J^4/S). 
\ea
Numerically, we find $J$ is always relevant at $\kappa=2$ and the system flows to a strong coupling fixed point $K$ at $J=\infty$. The first two terms in the second line are also consistent with the result obtained from perturbation theory as shown in Ref.~\cite{affleck1995conformal} (note that we've scaled the Kondo coupling).


At $J=0$, the model becomes a Bose-Kondo model with the following beta function of $h$:
\ba 
\beta_h\bigg|_{J=0} = -\frac{1}{S}\tilde{\epsilon} h + \frac{2h^2}{S\pi (1+h^2)}.
\ea 
This is consistent with previous results obtained in Refs.~\cite{cuomo2022spin,nahum2022fixed}. It gives two stable fixed points at $h = (1-\sqrt{1-\tilde{\epsilon}^2})/ {\tilde{\epsilon}},\infty$ (L and L$'$) and one critical point at $h=(1+\sqrt{1-\tilde{\epsilon}^2})/ {\tilde{\epsilon}}$ (LC). After introducing the Kondo coupling, we find $J$ is irrelevant at both L and LC, but marginal at L$'$. 

At $h=\infty$, we have a line of fixed points (L$'$) where $J$ is always marginal. The beta function at the large $h$ is:
\ba 
\beta_J = \frac{J-J^3}{\pi S (1+J^2)}\frac{1}{h} + O(\frac{1}{h^2S}).
\ea 
Clearly, $\beta_J \rightarrow 0 $ as $h\rightarrow \infty$.

At $J=\infty$, besides the K fixed point, there is an unstable fixed point KD$'$ at $(h_c=\frac{2}{\tilde{\epsilon}\pi},J=\infty)$. To see the existence of such a critical point, we give the beta function of $h$ at large $J$:
\ba 
\beta_h= -\frac{\tilde{\epsilon} h}{S}+\frac{2}{S\pi} + O(\frac{1}{SJ}). 
\ea 
It suggests that $h$ flows to $0$ ($\infty$) when $h<(>) h_c$ at $J\rightarrow \infty$. 

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{rg.pdf}
    \caption{Renormalization group flow at $\tilde{\epsilon}<\tilde{\epsilon}_1$ (a), $\tilde{\epsilon}_1<\tilde{\epsilon}<\tilde{\epsilon}_2$ (b) and $\tilde{\epsilon}_2<\tilde{\epsilon}$ (c) obtained from numerical evaluation. For purpose of illustration, we normalize all the vectors to have the same length. }
    \label{fig:rg_num} 
\end{figure*}


\subsection{D. Topological contribution from the Kondo coupling} 
In this section, we'll show that the Kondo coupling 
contributes an additional topological term to the impurity spin at sufficiently large values of $J_K$. For the purpose of discussion, we consider pure Kondo model, in other words, the BFK model at $g=0$ or $h=0$. We then integrate out the fermionic fields, which gives the following action
\ba 
S_{K}
=-iSS_B 
-K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg).
\ea 
We now evaluate the fermion determinant contribution at sufficiently large $J_K$ via  a Trotter decomposition
\ba 
\exp\bigg ( K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg) \bigg) = \bigg[\lim_{N\rightarrow \infty }\text{Tr}[\prod^N_i \exp[-H(\tau_i) d\tau]]\bigg]^K
\ea 
where
\ba 
H(\tau_i) = J_K\sum_\mu c_{\bm{r}_0}\frac{\sigma^\mu}{2} c_{\bm{r}_0}S^\mu(\tau_i) + \sum_{\bm{k} \sigma } \epsilon_{\bm{k}}c_{\bm{k},\sigma}^\dag c_{\bm{k},\sigma}.
\ea 
Since the $K$ orbitals are
degenerate, it's sufficient to consider one orbital and multiply the contribution by $K$. When $J_K$ is sufficiently large, the hopping between site $\bm{r}_0$ and the other sites is a high order effect ($\frac{t^2}{J_K}$). Here, we focus on the leading order contribution, and it's sufficient to take $H(\tau_i) = J_K/2\sum_\mu c_{\bm{r}_0}\sigma^\mu c_{\bm{r}_0}S^\mu(\tau_i)$. The eigenstates and eigenvalues of $H(\tau_i)$ are 
\ba 
&&|v_{+}(\tau_i)\rangle = \cos(\frac{\theta(\tau_i)}{2})c_{\bm{r_0},\up}^\dag 
|0\rangle +\sin(\frac{\theta(\tau_i)}{2})e^{i\phi(\tau_i)}c_{\bm{r_0},\dn}^\dag |0\rangle \nonumber \\
&&|v_{-}(\tau_i)\rangle = -\sin(\frac{\theta(\tau_i)}{2})c_{\bm{r_0},\up}^\dag 
|0\rangle +\cos(\frac{\theta(\tau_i)}{2})e^{i\phi(\tau_i)}c_{\bm{r_0},\dn}^\dag |0\rangle \nonumber \\
&&|v_0(\tau_i)\rangle = |0\rangle \nonumber  \\
&&|v_d(\tau_i)\rangle = |\up\dn\rangle \nonumber \\
&&E_{\pm}(\tau_i) = \pm \frac{J_KS}{2},\quad E_0(\tau_i) = 0 , \quad E_d(\tau_i) =0
\ea 
where we let $\bm{S}(\tau_i) =S (\sin(\theta(\tau_i))\cos(\phi(\tau_i)), \sin(\theta(\tau_i))\sin(\phi(\tau_i)), \cos(\theta(\tau_i)))$. We have the following decomposition via the eigenvalues and eigenstates
\ba 
e^{-H(\tau_i)d\tau} &=& e^{-E_+(\tau_i)d\tau}|v_{+}(\tau_i)\rangle\langle v_{+}(\tau_i)|+ e^{-E_-(\tau_i)d\tau}|v_{-}(\tau_i)\rangle\langle v_{-}(\tau_i)|
\nonumber \\
&&
+e^{-E_0(\tau_i)d\tau} |v_0(\tau_i)\rangle \langle v_0(\tau_i)| 
+e^{-E_d(\tau_i)d\tau} |v_d(\tau_i)\rangle \langle v_d(\tau_i)|
.
\ea 
Combining the above identity with a Trotter decomposition of the fermion determinant, we obtain the following leading order contribution
\ba 
&&\exp(K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg)) \nonumber \\
&\propto& \exp\bigg(K\int_\tau \langle \partial_\tau v_-(\tau)|v_-(\tau)\rangle d\tau +   O(\frac{1}{J_K}) \bigg)\nonumber \\
&\propto& \exp\bigg(i\frac{K}{2}\int_\tau (1-\cos(\theta))\partial_\tau \phi d\tau +   O(\frac{1}{J_K}) \bigg)\nonumber \\
&=& \exp\bigg( -i\frac{K}{2}S_B+
O(\frac{1}{J_K})\bigg).
\ea 
The effective action  then
becomes 
\ba 
S_{K}
= -i(S-\frac{K}{2})S_B + O(\frac{1}{J_K})
\ea 
which indicates the non-trivial topological contribution from the coupling to the gapless fermions. 
Qualitatively,
 the spin has been Kondo-screened, with the effective spin becoming $S-\frac{K}{2}$. 
 {For the}
 perfect screening case
 {that we consider}, $S-\frac{K}{2}=0$, and we obtain a Kondo singlet.


\end{document}