
\section{Continuum Model and Symmetries}






\section{The Kondo lattice Hamiltonian}

In this section we detail the connection between the tight-binding model and the Schrieffer-Wolff transformation to obtain the Kondo lattice Hamiltonian given in the main text.

\subsection{The tight-binding Hamiltonian}
Assuming $W$-layer weakly correlated the tight-binding Hamiltonian that describes the low-energy properties of the system reads $H=H_0+H_W+H_t$: 
\begin{equation}
\label{H_tb}
\begin{split}
    H_0=&-\frac{\Delta}{2} (N_{Mo}-N_W)+U_{}\sum_{\mathbf{r}\in Mo}n_{\mathbf{r}\uparrow}n_{\mathbf{r}\downarrow}\\
    H_W=&- t_{W}\sum_{\langle \mathbf{r},\mathbf{r}'\rangle\in W} e^{-i\nu_{\mathbf{r},\mathbf{r}'}2\pi s_\sigma/3} c^\dagger_{\mathbf{r}\sigma} c_{\mathbf{r}'\sigma},\\
    H_t=&-t_{Mo} \sum_{\langle \mathbf{r},\mathbf{r}'\rangle\in Mo} c^\dagger_{\mathbf{r}\sigma} c_{\mathbf{r}'\sigma}-t_h\sum_{\langle\mathbf{r},\mathbf{r}'\rangle}c^\dagger_{\mathbf{r}}c_{\mathbf{r}'}.
\end{split}
\end{equation}
We notice that differently from $H_0$ and $H_W$ the contribution $H_t$ changes the local configuration in the $Mo$ layer.    

\subsection{The Schrieffer-Wolff transformation}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{hybr.pdf}
    \caption{Evolution of $\Gamma(\epsilon_F)/\Delta$ as a function of the filling $x$ of the conduction band.}
    \label{fig:hybrd}
\end{figure}
We assume that $U$ is the largest energy scale, $U\gg t_{Mo}$, and $\Gamma(\epsilon_F)/\Delta<1$ where $\Gamma(\epsilon_F)$ is the hybridization function:
 \begin{equation}
     \Gamma(\epsilon_F)=\pi t^2_h\rho(\epsilon_F)\langle V^*_\mathbf{k} V_\mathbf{k}\rangle_{FS}=\pi t^2_h\left(\frac{3\sqrt{3}a^2_M}{8\pi^2}\right)\oint_{FS} dk_t\frac{V^*_\mathbf{k} V_\mathbf{k}}{|\nabla_\mathbf{k} \epsilon_{\mathbf{k}\sigma}|}.
 \end{equation}
Fig.~\ref{fig:hybrd} show the evolution of $\Gamma(\epsilon_F)$ as a function of the filling $x$ for the value of $t^2_h/\Delta$ corresponding to $J_K/t_W=1$ value used for the numerical calculations shown in the manuscript. Within these assumptions valence fluctuations in $Mo$ layer are suppressed. We now observe that the tunneling $H_t$ can be decomposed as:
\begin{equation}
    H_t=\sum_{ q=-1}^1\sum_{ d=-1}^1 T_{ q, d},
\end{equation}
where $T_{ q, d}$ gathers all tunneling events that change the number of electrons unbalance between $Mo$ and $W$ by $ q$ and the double occupancies in $Mo$ layer by $d$. We list the $T_{q,d}$ operators below: \begin{equation}
    \begin{split}
        &T_{+1,+1}=-t_h\sum_{\mathbf{r}\in Mo}\sum_\sigma\sum_{j=1}^3 n_{\mathbf{r}\bar\sigma}f^\dagger_{\mathbf{r}\sigma}c_{\mathbf{r}+\bm\delta_j\sigma},\\
        &T_{+1,0}=-t_h\sum_{\mathbf{r}\in Mo}\sum_\sigma\sum_{j=1}^3 h_{\mathbf{r}\bar\sigma}f^\dagger_{\mathbf{r}\sigma}c_{\mathbf{r}+\bm\delta_j\sigma},\\
        &T_{0,+1}=-t_{Mo}\sum_{\mathbf{r}\in Mo}\sum_\sigma\sum_{j=1}^3\left[ n_{\mathbf{r}\bar\sigma}f^\dagger_{\mathbf{r}\sigma}f_{\mathbf{r}+\bm\gamma_j\sigma}h_{\mathbf{r}+\bm\gamma_j\bar\sigma}+ h_{\mathbf{r}+\bm\gamma_j\bar\sigma}f^\dagger_{\mathbf{r}+\bm\gamma_j\sigma}f_{\mathbf{r}\sigma}n_{\mathbf{r}\bar\sigma} \right],\\
        &T_{0,0}=-t_{Mo}\sum_{\mathbf{r}\in Mo}\sum_\sigma\sum_{j=1}^3\left[ n_{\mathbf{r}\bar\sigma}f^\dagger_{\mathbf{r}\sigma}f_{\mathbf{r}+\bm\gamma_j\sigma}n_{\mathbf{r}+\bm\gamma_j\bar\sigma}+ h_{\mathbf{r}\bar\sigma}f^\dagger_{\mathbf{r}\sigma}f_{\mathbf{r}+\bm\gamma_j\sigma}h_{\mathbf{r}+\bm\gamma_j\bar\sigma} + h.c. \right],
    \end{split}
\end{equation}
where we have introduced the operator $f_{\mathbf{r}\sigma}$ for electrons in $Mo$-layer, $T_{-q,-d}=T^\dagger_{q,d}$, $n_{\mathbf{r}\sigma}=f^\dagger_{\mathbf{r}\sigma} f_{\mathbf{r}\sigma}$ and $h_{\mathbf{r}\sigma}=1-n_{\mathbf{r}\sigma}$. Before moving on, we notice that the intralayer hopping $H_{Mo}$ does not change the charge imbalance $N_{Mo}-N_{W}$ ($q=0$), i.e. $H_{Mo}$ commutes with $N_{Mo}-N_W$ ($[N_{Mo}-N_{W},H_{Mo}]=0$). Moreover, $[N_{Mo}-N_W,T_{ q, d}]=2q T_{ q, d}$ and
\begin{equation}
    \sum_{\mathbf{r}\in Mo}[n_{\mathbf{r}\uparrow}n_{\mathbf{r}\downarrow},T_{ q, d}]= d T_{q, d}.
\end{equation}
We seek  an unitary transformation $U=\exp(-iS)$ which eliminates hops between states with different numbers of doubly occupied sites and interlayer charge imbalance~\cite{PhysRev.149.491,PhysRevB.37.9753}:
\begin{equation}
    \bar H\equiv e^{iS} H e^{-iS}=H+[iS,H]+\frac{[iS,[iS,H]]}{2!}+\cdots,
\end{equation}
where we applied the Baker-Campbell-Haussdorf formula. We now notice that to the lowest order in the expansion we have: 
\begin{equation}
\label{SW_1st}
    [H_0,iS]=\sum_{q,d}^{(q,d)\neq0}T_{q,d},
\end{equation}
where in previous sum we exclude the term with $q=0$ and $d=0$ in shorthand notation $(q,d)\neq0$. We readily realize that the solution of Eq.~\eqn{SW_1st} reads: 
\begin{equation}
    S=-i\sum_{q,d}^{(q,d)\neq0}\frac{T_{q,d}}{d U-q\Delta}.
\end{equation}  
From the latter expression we find:
\begin{equation}
\begin{split}    
   
   
    \bar H=& H_0+T_{0,0}+H_W+\frac{1}{2}\sum_{q,d}^{(q,d)\neq0}\sum_{q',d'}^{(q',d')\neq0}\frac{[T_{q,d},T_{q',d'}]}{-q\Delta + dU}+\sum_{q,d}^{(q,d)\neq0}\frac{[T_{q,d},H_W+T_{0,0}]}{-q\Delta + dU}+O([H_W+T_{0,0}]S^2).
\end{split}
\end{equation}
Projecting the model in the low-energy subspace with one-electron in $Mo$ layer, i.e. $n_{\mathbf{r}}=1$ for $\mathbf{r}\in Mo$, we find that the latter term vanishes. 
Furthermore, the projection constraints $q'=-q$ and $d'=-d$:  
\begin{equation}
\label{sm:Hsd_0}
    \bar H=- t_{W}\sum_{\langle \mathbf{r},\mathbf{r}'\rangle\in W} e^{-i\nu_{\mathbf{r},\mathbf{r}'}2\pi s_\sigma/3} c^\dagger_{\mathbf{r}\sigma} c_{\mathbf{r}'\sigma}-\mu N_W+\frac{1}{2}\sum_{q,d}^{(q,d)\neq0}\frac{[T_{q,d},T_{-q,-d}]}{ dU-q\Delta},
\end{equation}
where the constant energy term $\Delta N_{W}/2$ has been absorbed in the chemical potential shift and $T_{0,0}$ vanishes in the low-energy subspace. 
By performing straightforward calculations we find that the commutator in Eq.~\eqn{sm:Hsd_0} gives: 
\begin{equation}
\begin{split}
    \frac{1}{2}\sum_{q,d}^{(q,d)\neq0}\frac{[T_{q,d},T_{-q,-d}]}{ dU-q\Delta}=&t^2_h\left(\frac{1}{\Delta}+\frac{1}{U-\Delta}\right)\sum_{\mathbf{r}\in Mo}\sum_{jl=1}^3\bm S_\mathbf{r}\cdot c^\dagger_{\mathbf{r}+\bm\delta_j \sigma} \,\bm \sigma_{\sigma\sigma'}\, c_{\mathbf{r}+\bm\delta_l\sigma'}+\frac{4t^2_{Mo}}{U}\sum_{\langle\mathbf{r},\mathbf{r}'\rangle\in A}\bm S_\mathbf{r} \cdot\bm S_{\mathbf{r}'}\\
    &+t^2_h\left(\frac{1}{\Delta}-\frac{1}{U-\Delta}\right)\sum_\mathbf{k} |V_\mathbf{k}|^2 c^\dagger_\mathbf{k} c_\mathbf{k}. 
\end{split}
\end{equation}
Expanding close to the bottom of the conduction band [$\kappa(\kappa')$ for $\uparrow(\downarrow)$] $|V_{\kappa/\kappa'+\mathbf{k}}|^2\simeq 9(a_Mk)^2/4$ we realize that the latter terms can be simply accounted by a redefinition of the bare mass $m_W\to m_W/[1+(t_h/\Delta-t_h/(U-\Delta))t_h/(2t_W)]$.
Finally, we find the effective spin-fermion Hamiltonian:
\begin{equation}
\label{sm:Hsd}
    \bar H = \sum_{\mathbf{k}}  \xi_{\mathbf{k}\sigma}c^\dagger_{\mathbf{k}\sigma} c_{\mathbf{k}\sigma}+ J_H\sum_{\langle\mathbf{r},\mathbf{r}'\rangle\in A}\bm S_\mathbf{r} \cdot\bm S_{\mathbf{r}'}+\frac{1}{2N}\sum_{\mathbf{r}\in \text{Mo}}\sum_{\mathbf{k},\mathbf{p}} e^{-i(\mathbf{k}-\mathbf{p})\cdot\mathbf{r}}
    J_{\mathbf{k},\mathbf{p}}
    \bm S_\mathbf{r}\cdot c^\dagger_{\mathbf{k} \sigma} \,\bm \sigma_{\sigma\sigma'}\, c_{\mathbf{p}\sigma'},
\end{equation}
where $J_H=4t^2_{Mo}/U$, $J_{\mathbf{k},\mathbf{p}}=J_K V^*_\mathbf{k} V_\mathbf{p}$, $J_K=2t^2_h[1/\Delta+1/(U-\Delta)]$, $\xi_{\mathbf{k}\sigma}=\epsilon_{\mathbf{k}\sigma}-\epsilon_F$, $\epsilon_{\mathbf{k}\sigma}$ is the electron dispersion and $\epsilon_F$ fixes the number of electron in the conduction band. The dispersion relation reads $\epsilon_{\mathbf{k}\sigma}=-2t_W\sum^3_{j=1}\cos(\mathbf{k}\cdot\bm a_i+2\pi s_\sigma/3)$, $\bm a_1=\sqrt{3}a_{M}(1,0)$, $\bm a_{2,3}=\sqrt{3}a_{M}\left(-1/2,\pm\sqrt{3}/2\right)$ are the lattice vectors with $a_{M}=5$nm the moir\'e cell lattice constant. The form factor is $V_\mathbf{k}=\sum_{j=1}^3 e^{i\mathbf{k}\cdot\bm\delta_j}$. Despite $U$ is large and $J_H$ small the particular form of the exchange interaction $J_{\mathbf{k},\mathbf{p}}=J_K V^*_\mathbf{k} V_\mathbf{p}$ gives rise to a non-trivial competition between a low-density magnetic phase and a paramagnetic heavy Fermi liquid.

\subsubsection{On the higher-order corrections}

Higher order corrections in the expansion lower the symmetry of the spin interaction in the $Mo$-layer. This can be simply realized noticing that the phase factor in the Hamiltonian $H_W$ lower the spin-symmetry of the model to $U(1)$ around $z$. Including higher-order hopping processes mediated by $t_h$ gives an XXZ spin model with Dzyaloshinskii–Moriya interactions:   
\begin{equation}
    H_S=J_H\sum_{\langle\mathbf{r},\mathbf{r}'\rangle\in \text{Mo}}\left(S^z_\mathbf{r} S^z_{\mathbf{r}'}+\gamma S^+_\mathbf{r} S^-_{\mathbf{r}'}+h.c.\right)+D\sum_{\langle\mathbf{r},\mathbf{r}'\rangle\in \text{Mo}}\left(\bm S_\mathbf{r}\times \bm S_{\mathbf{r}'}\right)_z,
\end{equation}
we refer to Ref.~\cite{Devakul_magic_2021} for additional details. For simplicity we consider the isotropic limit $\gamma=1/2$ and $D=0$ in our calculations.

\section{The mean-field approach}

In this section we detail the mean-field of Abrikosov fermions discussed in the main text. We perform the decomposition of the spin-$1/2$ into spinons $\bm S_\mathbf{r}=\chi^\dagger_{\mathbf{r}\alpha}\bm\sigma_{\alpha\beta}\chi_{\mathbf{r}\beta}/2$ and performing the mean-field decomposition in the magnetic and excitonic channels we obtain the mean-field Hamiltonian: 
\begin{equation}
\label{mean_field_Ham}
\begin{split}    
    \bar H_{\rm mf}=&\lambda\sum_{\mathbf{k}} \chi^\dagger_\mathbf{k} \chi_\mathbf{k}+\frac{h^z}{2}\sum_{\mathbf{k}} \chi^\dagger_\mathbf{k} \sigma^z \chi_\mathbf{k} +\frac{h^\parallel}{2}\sum_{\mathbf{k}}\left[ \chi^\dagger_{\mathbf{k}+\mathbf{Q}} \sigma^- \chi_{\mathbf{k}}+h.c. \right]+\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}\sigma}c^\dagger_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}\\
    &+\frac{ M^z}{2}\sum_\mathbf{k} J_{\mathbf{k},\mathbf{k}} c^\dagger_\mathbf{k}\sigma^z c_\mathbf{k}+\frac{M^\parallel}{2}\sum_\mathbf{k}\left[ J_{\mathbf{k}+\mathbf{Q},\mathbf{k}} c^\dagger_{\mathbf{k}+\mathbf{Q}}\sigma^- c_\mathbf{k}+h.c.\right]\\
    &+J_K\sum_{\mathbf{k}}\left[V^*_\mathbf{k}\Phi^* c^\dagger_{\mathbf{k}} \chi_{\mathbf{k}}+h.c.\right].
\end{split}
\end{equation}
The mean-field free-energy reads: 
\begin{equation}
\begin{split}
\label{meafield_free}
    F_{\rm mf}=&-k_BT\sum_{\mathbf{k}}\sum_\lambda\log\left[1+e^{-E_{\mathbf{k}\lambda}/k_BT}\right]+2J_K N |\Phi |^2+\frac{3J_H N}{2} {M^\parallel}^2\\
    &-3J_H N {M^z}^2-J_K M^z N m^z-J_K N M^\parallel m^\parallel-\lambda N+N \mu x,
\end{split}
\end{equation}
where $E_{\mathbf{k}\lambda}$ is the mean-field dispersion relation. The  mean-field solution is obtained by minimizing $F_{\rm mf}$ with respect to the variational parameters $M^z, M^\parallel, \Phi, m^z, m^\parallel$. The Lagrange multiplier $\lambda$ imposes the Gutzwiller constraint $\sum_\mathbf{k}\langle \chi^\dagger_\mathbf{k}\chi_\mathbf{k}\rangle/N=1$, while the chemical potential $\mu$ fixes the number of particle in conduction band $\sum_\mathbf{k} \langle c^\dagger_{\mathbf{k}}c_{\mathbf{k}}\rangle/N=x$. Taking the saddle point of Eq.~\eqn{meafield_free} with respect to the variational parameters gives a set of self-consistency equations that are solved by find-root algorithm. 

\subsection{Energetics of the HFL and AFM states}
\label{energetics}

In this section we detail the weak coupling expansion to determines the characteristic energy scales of the AFM and HFL phases that are the RKKY energy and the Kondo temperature, respectively. 
We will also introduce the effective model describing the quasiparticle excitations in the two different regimes.   

\subsubsection{AFM}
In the magnetic regime the local moments $\bm S_\mathbf{r}$ form a 120$^\circ$ AFM order with $\langle\bm S_\mathbf{r}\rangle=\left(M^\parallel\cos\mathbf{Q}\cdot\mathbf{r},M^\parallel \sin\mathbf{Q}\cdot\mathbf{r},M^z\right)$, $\mathbf{Q}=\kappa-\kappa'$,  correspondingly the conduction electron Hamiltonian reads: 
\begin{equation}
\begin{split}
    \bar H^c_{\rm mf} =&\sum_\mathbf{k} \xi_{\mathbf{k} \sigma}c^\dagger_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}+\frac{J_K M^z}{2}\sum_\mathbf{k} |V_{\mathbf{k} }|^2 c^\dagger_\mathbf{k}\sigma^z c_\mathbf{k}+\frac{J_K M^\parallel}{2}\sum_{\mathbf{k}}\left[V^*_{\mathbf{k} }V_{\mathbf{k}+\mathbf{Q} }c^\dagger_\mathbf{k}\sigma^+ c_{\mathbf{k}+\mathbf{Q}}+h.c.\right].
\end{split}
\end{equation}
In the limit of low-doping the Fermi surface is a small electron pocket around $\kappa(\uparrow)$ and $\kappa'(\downarrow)$ that are folded into the origin of the magnetic Brillouin zone $\gamma_m$, depicted in Fig.~\ref{fig:fig_continuum}(b). Expanding close to quadratic order around $\gamma_m$ and keeping only the two lowest energy bands we find the continuum model: 
\begin{equation}
\label{conduction_electrons}
    \bar{\mathcal H}^c_{\rm mf}=\sum_{\mathbf{k}}c^\dagger_{\mathbf{k}}\left[\left(\frac{\hbar^2 k^2}{2m_W}-\epsilon_F\right)\sigma_0+\frac{9J_Ka^2_M }{8}\sum_{a}d_a(\mathbf{k})\sigma_a \right]c_\mathbf{k},
\end{equation}
where $d_z(\mathbf{k})=M^z k^2$, $d_x(\mathbf{k})=-M^\parallel(k^2_x-k^2_y)$ and $d_y(\mathbf{k})=2M^\parallel k_xk_y$.
We readily realize that due to the SOC term the spin is no longer a good quantum number, the eigenstates $\ket{u_{\mathbf{k}\lambda}}$ are labeled by $\lambda=\pm$ and the corresponding eigenvalues are $\epsilon_{\mathbf{k}\lambda}=(\hbar^2k^2/2m_W-\epsilon_F)+\lambda 9J_Ka^2_Mk^2|M|/8$ with $|M|=\sqrt{(M^z)^2+(M^\parallel)^2}$.
We observe that the theory is O$(3)$ invariant under rotation of the magnetization $M^a\to R_{ab}M^b$. The resulting stabilization energy does not depend on the orientation of the local moments. 
The Fermi momentum obtained by setting $\epsilon_{\mathbf{k}\lambda}=0$ reads $k^{\pm}_F=\sqrt{2m_\lambda\epsilon_F/\hbar^2}$ where $m_\lambda=m_W/[1+\lambda (J_K/8t_W)]$ and the Fermi energy is $\epsilon_F\simeq x/\left(\sum_\lambda \rho_\lambda\right)$ with $\rho_\lambda=\rho_0 m_\lambda/m_W$. We notice that in the Kondo regime $J_K$ is smaller than the bandwidth of conduction electrons $9t_W$ so that the mass is always positive. The kinetic energy variation with respect to the normal state reads:
\begin{equation}
    \delta \epsilon^c_{\rm mf}=\epsilon^c_{\rm mf}-\epsilon^c_{J_K=0}=-\rho_0\bar J^2_K|M|^2/2,
\end{equation}
where $\bar J_K$ is the average over the FS of the Kondo exchange.
We conclude that the total energy per site in the magnetic regime is given by: 
\begin{equation}
    \epsilon_{\rm mf}=-3J_H |M|^2/2-\rho_0\bar J^2_K|M|^2/2.
\end{equation}
The energy gain from the coupling between the conduction band and the local moments goes quadratically in the electron density $x$.  

We conclude observing that in the large $J_K$ regime corresponding to a very small charge transfer energy one of the two band in Eq.~\eqn{conduction_electrons} develops a negative curvature around $k=0$. Including higher-order corrections in momentum regularize the theory in the UV allowing for a finite doping in conduction band. It is known~\cite{Sun_2009} that in this regime the model is unstable towards ferromagnetism which opens a gap at the quadratic band touching. The analysis requires including interaction effects in the conduction electrons and goes beyond the regime of validity of the Kondo lattice model we investigate here.   

\subsubsection{HFL}

In the paramagnetic regime electrons are described by the mean-field Hamiltonian reads: 
\begin{equation}
\label{hfl_Hamiltonian}
    \bar H_{\rm mf}=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}\sigma}c^\dagger_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}+\lambda\sum_\mathbf{k} \chi^\dagger_\mathbf{k}\chi_\mathbf{k}+J_K\sum_{\mathbf{k}\sigma}\left[V^*_\mathbf{k}\Phi^*c^\dagger_\mathbf{k}\chi_\mathbf{k}+h.c.\right].
\end{equation}
We easily realize that the Green's function of the problem reads:
\begin{equation}
    \bm G^{-1}_\sigma(\mathbf{k},i\epsilon)=\left(\begin{matrix} i\epsilon-\lambda & -J_K\Phi^* V^*_\mathbf{k} \\  -J_K\Phi V_\mathbf{k} & i\epsilon-\xi_{\mathbf{k}\sigma} \end{matrix}\right),
\end{equation}
so that the saddle-point equation for the bosonic amplitude $\Phi$ can be written as: 
\begin{equation}
\label{sce_phi}
    \Phi=-\frac{T}{2N}\sum_{\mathbf{k}\sigma}\sum_{i\epsilon}V^*_\mathbf{k} \mathcal G_{\chi,\sigma}(\mathbf{k},i\epsilon)J_K V_\mathbf{k} \Phi G_{c,\sigma}(\mathbf{k},i\epsilon),
\end{equation}
where $\mathcal G_{\chi,\sigma}(\mathbf{k},i\epsilon)=(i\epsilon-\lambda)^{-1}$ is the bare Green's function and $G_{c,\sigma}(\mathbf{k},i\epsilon)=[i\epsilon-\xi_{\mathbf{k}\sigma}-\Sigma_{c,\sigma}(\mathbf{k},i\epsilon)]^{-1}$ with 
\begin{equation}
\label{self_energy_c}
    \Sigma_{c,\sigma}(\mathbf{k},i\epsilon)=J^2_K|\Phi|^2|V_\mathbf{k}|^2/(i\epsilon-\lambda).
\end{equation}
Discarding the $\Phi=0$ solution the equation reduces to
\begin{equation}
    \frac{1}{J_K}=\frac{1}{2N}\sum_{\mathbf{k}\sigma}\frac{f(E_{\mathbf{k}-\sigma})-f(E_{\mathbf{k}+\sigma})}{\sqrt{(\xi_{\mathbf{k}\sigma}-\lambda)^2+4J^2_K|\Phi|^2|V_\mathbf{k}|^2}}=\frac{1}{2N}\sum_{\mathbf{k}\sigma}\frac{\theta(-E_{\mathbf{k}-\sigma})}{\sqrt{(\xi_{\mathbf{k}\sigma}-\lambda)^2+4J^2_K|\Phi|^2|V_\mathbf{k}|^2}},
\end{equation}
where $E_{\mathbf{k}\pm}=(\xi_{\mathbf{k}\sigma}+\lambda)/2\pm\sqrt{(\xi_{\mathbf{k}\sigma}-\lambda)^2+4J^2_K|\Phi|^2|V_\mathbf{k}|^2}/2$ and the RHS is obtained taking the we took the zero temperature limit and considering the case $x<1$ where only the lower band is filled. In addition we also have the self-consistent equation for $\lambda$: 
\begin{equation}
\label{sce_lambda}
\frac{T}{N}\sum_{\mathbf{k}\sigma}\sum_{i\epsilon}\frac{1}{i\epsilon-\lambda-\Sigma_{\chi,\sigma}(\mathbf{k},i\epsilon)}=1\implies\frac{1}{2N}\sum_{\mathbf{k}\sigma}\theta(-E_{\mathbf{k}-\sigma})\left(1+\frac{\xi_{\mathbf{k}\sigma}-\lambda}{\sqrt{(\xi_{\mathbf{k}\sigma}-\lambda)^2+4J^2_K|\Phi|^2|V_\mathbf{k}|^2}}\right)=1.
\end{equation}
The onset of the HFL instability can be determined by determining the instability condition of the normal state to hybridization. In this case the solution of Eq.~\eqn{sce_lambda} is $\lambda=0$, i.e. the local moments are pinned at the Fermi level, and Eq.~\eqn{sce_phi} becomes: 
\begin{equation}
\label{Kondo_critical}
    \frac{1}{J_K}+\Pi^0_{\chi c}(\mathbf{q}\to0,i\Omega=0)=0,
\end{equation}
where $\Pi^0_{\chi c}(\mathbf{q}\to0,\tau)=-\langle T_\tau(\hat V_{\mathbf{q}=0}(\tau)\hat V^\dagger_{\mathbf{q}=0})\rangle/N$, $\hat V_{\mathbf{q}=0}=\sum_\mathbf{k} V^*_\mathbf{k} c^\dagger_\mathbf{k}\chi_\mathbf{k} /\sqrt{2}$ and $\Pi^0_{\chi c}(\mathbf{q}\to0,i\Omega)=\int^\beta_0 e^{i\Omega\tau}\Pi^0_{\chi c}(\mathbf{q}\to0,\tau)$. Expanding close to the bottom of the band Eq.~\eqn{Kondo_critical} becomes: 
\begin{equation}
    \frac{1}{J_K}-\frac{\rho_0}{2}\int^{E_\Lambda}_{\epsilon_F}\frac{|V_\epsilon|^2}{\epsilon-\epsilon_F}-\frac{\rho_0}{2}\int^{\epsilon_F}_{0}\frac{|V_\epsilon|^2}{\epsilon_F-\epsilon}=0,
\end{equation}
where $|V_\epsilon|=9a^2_M k^2_\epsilon/4$ with $k_\epsilon=\sqrt{2m_W\epsilon/\hbar^2}$.The integral is characterized by a log singularity at $\epsilon_F$. We introduce the IR cutoff $T_K$, $\epsilon\to\epsilon\pm T_K$, that regularize the divergence. Finally, by performing simple calculations we find: 
\begin{equation}
    T_K\simeq \epsilon_Fe^{-1/(\rho_0\bar J_K)}.
\end{equation}
Since the average over the Fermi surface $\bar J_K$ goes linearly with the doping $x$ in conduction band we find that the Kondo temperature is exponentially suppressed in the limit $x\to0^+$.






\section{Transport properties}
 
 In this section we detail the evaluation of the transverse and longitudinal conductivities in the various phases of the phase diagram. 
 
\subsection{Charge transport in the HFL}
 
 
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{Chern_flux.pdf}
    \caption{  Spin Hall effect (top panel) and anomalous Hall effect (bottom panel) in the HFL as a function of the Zeeman field.}
    \label{fig:sHE_aHE}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{circle.pdf}
    \caption{Ratio between the numerical evaluation of Eq.~\eqn{longitudinal} and $\sigma_{xx}=e^2\tau(1-x)/\bar m^*$ obtained assuming a circular Fermi surface with average velocity.}
    \label{fig:circular}
\end{figure}

In the paramagnetic regime the local moments fractionalize giving rise to a finite density of holes in MoTe$_2$ layer. In this regime the field $\Phi_{\mathbf{r}}$ is equivalent to the holon operator carrying physical charge $-1$~\cite{Chowdhury_2018}, i.e. opposite to the electron charge. As a result both conduction electrons $c$ and spinon $\chi$ contribute to the charge current. Within the semiclassical Boltzmann equation approach and in the relaxation time approximation the transport properties are simply obtained as integrals over the quasiparticle Fermi surface~\cite{PhysRevB.43.193,Haldane_2004}: 
\begin{equation}
\label{longitudinal}
    \sigma_{xx}=\frac{3\sqrt{3} a^2_Me^2\tau}{8\pi^2}\sum_\sigma\oint_{FS} dk_t \frac{\bar v^x_{F,\sigma}\bar v^x_{F,\sigma}}{\hbar |\bar{\bm v}_{F,\sigma}|},
\end{equation}
and
\begin{equation}
\label{transverse}
    \sigma_{xy}\equiv \sigma^{\rm Ohm}_{xy}+\sigma^{\rm AH}_{xy}=\frac{3\sqrt{3} a^2_M}{8\pi^2}\frac{e^3\tau^2 B}{\hbar}\sum_\sigma\oint_{FS} dk_t \frac{\bar v^y_{F,\sigma}\left[\bar {v}^y_{F,\sigma}\partial_{k_x}-\bar {v}^x_{F,\sigma}\partial_{k_y}\right]\bar v^x_{F,\sigma}}{\hbar |\bar{\bm v}_{F,\sigma}|}+\frac{e^2}{2\pi h}\sum_\sigma\oint_{FS} dk_t\, \bm t\cdot\bm A_{\sigma}(\mathbf{k}),
\end{equation}
where the anomalous Hall contribution $\sigma^{\rm AH}_{xy}$ arises from the circulation of the Berry connection $\bm A_{\sigma}(\mathbf{k})=i\braket{u_{\mathbf{k}\sigma}}{\partial_{k_a} u_{\mathbf{k} \sigma}}$, with $\ket{u_{\mathbf{k} \sigma}}$ occupied eigenstate of the Hamiltonian in Eq.~\eqn{hfl_Hamiltonian}, along the FS. 
In Eqs.~\eqn{longitudinal} and~\eqn{transverse} $k_t$ is the component of $\mathbf{k}$ along the tangent $\bm t$ to the FS curve. 
We notice that the FS of the heavy quasiparticle is obtained by the set of $\mathbf{k}$ points solution of the equation:
\begin{equation}
    {\text{heavy FS}}: \xi_{\mathbf{k}\sigma}-J^2_K|\Phi|^2 |V_\mathbf{k}|^2/\lambda=0,
\end{equation}
where the second term comes from the conduction electron self-energy $\Sigma_{c,\sigma}(\mathbf{k},0)$ introduced in Eq.~\eqn{self_energy_c} computed at $i\omega=0$. The Fermi velocity $\bm v_{F,\sigma}$ is obtained replacing $\mathbf{k}\to \mathbf{k}+e\mathbf{A}(t)/\hbar$ in $\bar\xi_{\mathbf{k}\sigma}=Z_{\mathbf{k}\sigma}[\xi_{\mathbf{k}\sigma}-J^2_K|\Phi|^2 |V_\mathbf{k}|^2/\lambda]$ and expanding to linear order around the FS we find $\bar \xi_{\mathbf{k}_F+\frac{e}{\hbar}\mathbf{A}(t)\sigma}\simeq e \bar{\bm v}_{F}\cdot\mathbf{A}(t)$ where: 
\begin{equation}
    \bar{\bm v}_{F,\sigma} = \frac{Z_{\mathbf{k}_F,\sigma}}{\hbar}\nabla_\mathbf{k}[\xi_{\mathbf{k} \sigma}-J^2_K|\Phi|^2 |V_\mathbf{k}|^2/\lambda]\Big|_{\mathbf{k}_F},
\end{equation}
and $Z_{\mathbf{k},\sigma}=[1-\partial_z\Sigma_{c,\sigma}(\mathbf{k},z)|_{z=0}]^{-1}$ is the quasiparticle weight. 
The evaluation of $\sigma_{xx}$ and $\sigma^{\rm Ohm}_{xy}$ is considerably simplified observing that the Fermi surface consists of a hole-pocket around $\bm\kappa'$ for spin $\uparrow$ and $\bm\kappa$ for spin $\downarrow$, respectively. Assuming a circular hole-like Fermi surface with average mass $\bar m^*$ the longitudinal contribution becomes 
\begin{equation}
\label{longitudinal_sigma}
    \sigma_{xx}=2\frac{e^2\tau}{\bar m^*}\left(\frac{\sqrt{3}a^2_M}{8\pi^2}\right)k^2_F\int^{2\pi}_0d\phi \cos^2\phi=\frac{e^2\tau}{\bar m^*}\left(2\pi k^2_F\frac{\sqrt{3}a^2_M}{8\pi^2}\right)=\frac{e^2\tau(1-x)}{\bar m^*}.
\end{equation}
As a sanity check we show in Fig.~\ref{fig:circular} the ratio between Eq.~\eqn{longitudinal_sigma} and the numerical evaluation of Eq.~\eqn{longitudinal} for various concentrations of electrons in the $W$-layer. By following the same line of reasoning we obtain $\sigma^{\rm Ohm}_{xy}=-e^3\tau^2 B(1-x)/\bar m^*{}^2$. We now look at the anomalous contribution $\sigma^{\rm AH}_{xy}$ which can be conveniently written as $\sigma^{\rm AH}_{xy}=e^2( W_\uparrow+ W_\downarrow)/h$ where $W_\sigma=\int d^2\mathbf{k}  f(\epsilon_{\mathbf{k}\sigma}) \Omega_{\sigma}(\mathbf{k})/(2\pi)$ and $\Omega_\sigma(\mathbf{k})=\partial_{k_x}A^y_\sigma(\mathbf{k})-\partial_{k_y}A^x_\sigma(\mathbf{k})$. Due to the opposite winding of spin $\uparrow$ and $\downarrow$ we find that in the absence of a magnetic field $W_\uparrow=- W_\downarrow$. We observe that the difference $ W_\uparrow- W_\downarrow$ gives a finite spin Hall (SH) conductivity~\cite{Kane_2005} even in the absence of the external field. Fig.~\ref{fig:sHE_aHE} shows $ W_\uparrow\pm  W_\downarrow$ as a function of the magnetic field at doping $x=0.46$ and $J_K/t_W=1$. The small value of $W_\uparrow+W_\downarrow$ at finite $B$ follows from the Berry curvature distribution which is peaked around the bare FS of the conduction electrons. As a result we find a small AH contribution and a large SH one. 








\subsubsection{The magnetic regime}
\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{fig3.pdf}
    \caption{(a) Dispersion of the electrons close to the bottom of the conduction band. Dashed lines are obtained with the continuum Hamiltonian~\eqn{} while the solid ones are the lowest energy bands of the mean-field Hamiltonian~\eqn{}. (b) Berry curvature around the origin of the magnetic Brillouin zone. Solid lines show the two Fermi surfaces. The calculations has been performed at doping $x=0.08$, $J_K/t_W=1$, $J_H/t_W=0.025$ and magnetic field $h_W/t_W=0.08$.
    }
    \label{fig:fig_continuum}
\end{figure}
We now turn our attention in the magnetic regime where only the $c$-electrons contribute to charge transport. The low-energy Hamiltonian is given in Eq.~\eqn{conduction_electrons} and describes electrons with dispersion relation $\epsilon_{\mathbf{k}\lambda}=\hbar^2 k^2/2m_\lambda$ shown in Fig.~\ref{fig:fig_continuum}(a) and group velocity $\bm v_{\mathbf{k}\lambda}=\hbar \mathbf{k} /m_\lambda$ with $m_\lambda=m_W/[1-\lambda J_K|M|/(4t_W)]$ and $\lambda=\pm$. Under the assumption of a single transport time, i.e. momentum- and band-independent, we apply Eqs.~\eqn{longitudinal} and \eqn{transverse} to find $\sigma_{xx}=e^2\tau x/\left(\sum_\lambda m_\lambda/2\right)$ and $\sigma^{\rm Ohm}_{xy}=e^3\tau^2 B x/\left(m_W\sum_\lambda m_\lambda/2\right)$. We now conclude our analysis considering the AHE contribution coming from the the Berry phase winding introduced by the $d_x(\mathbf{k})$ and $d_y(\mathbf{k})$ terms in Eq.~\eqn{conduction_electrons}. In the absence of an external magnetic field $M^z=0$ and the eigenstates are simply $\ket{u_{\mathbf{k}\pm}}=(1,\pm e^{-2i \varphi_\mathbf{k}})\sqrt{2}$ with $2\pi$-Berry phase around the origin of the magnetic Brillouin zone. A small Zeeman term opens a gap in the band structure and gives rise to a finite Berry curvature $\Omega_{\pm}(\mathbf{k})$:
\begin{equation}
\label{Berry_curvature}
    \Omega_{\pm}(\mathbf{k})=\pm\frac{64 h^2k^2\eta^2}{\left[16k^4\eta^2+\left(2h+\frac{m_{\uparrow}-m_{\downarrow}}{m_{\uparrow}m_{\downarrow}}k^2\right)^2\right]^{3/2}},
\end{equation}
where $\eta=9J_K M^\parallel/8$ and $m_{\uparrow/\downarrow}=m_W/[1\pm J_K M^z/(4t_W)]$.
We notice that the Berry curvature is an even function of $k$, vanishes quadratically at $k=0$ and takes its maximum value at finite momentum $k$. The maximum is located at $k^4_{\rm max}=h^2/(8\eta^2)$ for $M^z=0$. The momentum space distribution in the magnetic moir\'e~ Brillouin zone is given in Fig.~\ref{fig:fig_continuum}(b). We observe that from Eq.~\eqn{Berry_curvature} the anomalous Hall conductivity is obtained as: 
\begin{equation}
    \sigma^{\rm AH}_{xy}=\frac{e^2}{2\pi h}\sum_\lambda \int d^2\mathbf{k} \Omega_{\lambda}(\mathbf{k}) f(\epsilon_{\mathbf{k}\lambda}-\mu),
\end{equation}
which in the zero temperature limit becomes: 
\begin{equation}
    \sigma^{\rm AH}_{xy}=\frac{e^2}{h}\int^{k^+_F}_{k^-_F} dk\,k \,\Omega_{+}(k)=\frac{e^2}{h}\frac{h[k^2(m_\downarrow-m_\uparrow)-2hm_\uparrow m_\downarrow]}{\sqrt{16k^4\eta^2(m_\uparrow m_\downarrow)^2+[k^2(m_\uparrow-m_\downarrow)+2hm_\uparrow m_\downarrow]^2}}\Bigg|^{k^+_F}_{k^-_F}.
\end{equation}
Fig.~\ref{fig:ahe_m} shows the evolution of the anomalous Hall conductivity in the low-density regime. 
Finally, Fig.~\ref{fig:ahe_2D} gives bird's-eye view over the anomalous Hall effect in the plane of the $W$-layer concentration and the Zeeman field.  

\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{ahe_magnetic.pdf}
    \caption{Top panel: evolution of the in-plane and out-of-plane magnetization as a function of the applied field. Bottom panel: anomalous Hall conductivity. The calculations has been performed at doping $x=0.08$, $J_K/t_W=1$, $J_H/t_W=0.025$. }
    \label{fig:ahe_m}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{AHE_2d.pdf}
    \caption{Intrinsic contribution to AHE conductivity in the plane of $W$-layer concentration $x$ and $h_W$ Zeeman field.}
    \label{fig:ahe_2D}
\end{figure}








\bibliographystyle{apsrev}
