\section{Conclusion}
There is an increasing demand for IoT systems and devices due to their potential in healthcare. IoMT technology has been aiding people with chronic illnesses in various ways, improving the quality of their lives. Making these systems smarter by introducing machine learning has great potential, with the focus shifting to obtaining faster real-time predictions over the years. Researchers have been experimenting with various distributed computing frameworks to achieve the highest possible prediction accuracy while being efficient in terms of resources at the same time.

This systematic literature review started with a broad scope of exploring how the healthcare industry had been combining machine learning and IoT. It then narrowed down to determining whether distributed machine learning could be applied to IoT healthcare data or not. Then narrowed further down to identifying and comparing popular distributed machine learning algorithms and frameworks being used in various state-of-the-art IoT healthcare systems.

Upon analyzing the primary studies, they were classified into three categories (scenarios) - monitoring health and vitals, diagnosing health conditions, and other innovative yet noteworthy systems being used in the healthcare industry today. Each of the primary studies touched upon distributed computing and distributed machine learning, with the latter being used by the majority. The studies discussing distributed computing were similar and utilized machine learning using the distributed frameworks proposed by the authors of the respective studies. Then the machine learning algorithms tested in the studies and their prediction accuracy were evaluated, giving us Random Forest as the algorithm that achieved the highest accuracy in a few cases. However, most studies proposed their own variant of a neural network tailored to their specific case studies and frameworks. This subject area is still fairly new with many contributions being made. It will take a few more years until we see state-of-the-art algorithms for IoMT scenarios.
\section{Discussion}
Collaborative (or distributed) machine learning algorithms were the focus of this review. The research goal was the identification of algorithms that performed the best in the IoT healthcare scenarios considered - live health monitoring systems and medical diagnostic systems. Completing the review gave us the answers we were looking for, but the results are not definitive i.e., Random Forest may not always achieve the highest accuracy when deployed in live monitoring systems. A few studies used other algorithms for their case studies and had different conclusions. Some even built their neural networks to achieve higher accuracy, the majority taking this approach for diagnostic systems.

Overall, no known algorithm performs consistently across all scenarios in healthcare as there are numerous variables involved. However, Random Forest could serve as a good starting point for distributed machine learning applications in IoMT. With the growing number of machine learning-driven IoT healthcare applications, more technical contributions to this subject area are expected with more conclusive results.
\section{Introduction}
Over the last few years, we have seen significant growth in the number of Internet of Things devices used in the healthcare industry. The widespread adoption of self-tracking smart wearable devices such as smartwatches has been a major contributor. Since such devices produce massive amounts of data, researchers and industry specialists have been exploring ways to make efficient use of them by building (artificial) intelligence into these devices via means of machine learning and making smart devices smarter. Given the small form factor of these devices, they are often limited in terms of hardware capabilities to carry out machine learning tasks on-device. Researchers have been exploring ways to tackle this problem through distributed machine learning techniques.

This systematic study reviews the trends of adoption of IoT systems in healthcare that utilize distributed machine learning and the state-of-the-art algorithms and frameworks used by those systems. We then identify the algorithms that perform best in terms of prediction accuracy in different healthcare scenarios in the subsequent sections. The study ends with a discussion of the findings and where this area of research could potentially be heading.

\section{Results}
In this section, the findings and results of the systematic literature review are reported. To make it easy to follow, this section has been split into three, each answering a research question.

\subsection{ML \& IoT in the Healthcare Industry Today}
After analyzing the primary studies for the ways in which ML and IoT are being combined in the healthcare industry today, a pattern was observed. The selected studies were then further classified based on categories created for each pattern as depicted in Figure \ref{fig:categories}.

\begin{figure}[h!]
    \centering
    \includegraphics[width=8.7cm]{charts/health_categories.PNG}
    \caption{Number of Studies in Each Category}
    \label{fig:categories}
\end{figure}

\subsubsection{\textbf{Monitoring Health and Vitals}}
The most popular application combining ML and IoT in healthcare was monitoring the health condition and vitals of patients and quantified selves. IoT technology, although fairly new, has increasingly been gaining popularity due to its potential in a world with a rising number of people with chronic illnesses. This along with the aging population has been identified to be placing significant stress on the healthcare system due to rising resource demands from them. The standard model for IoT healthcare systems proposed by \cite{Seven} suggested a rehabilitation system for people suffering from knee injuries which would enable the calculation of the position and angle of knees, and assist the patient with knee-related motor movements, all while monitoring the knee position and angle. The system could also be adapted for various other monitoring purposes including but not limited to the management of hypertension and monitoring critical health metrics such as pulse, blood pressure, body temperature, etc \cite{Seven}\cite{Eight}. Similarly, several parameters such as physiological parameters, speech, posture, skin, motion, and environmental conditions were monitored in \cite{ThirtySix} using various machine learning techniques. \par

According to the study performed by \cite{Sixteen}, poor medication adherence was identified as another leading factor significantly stressing the healthcare system. Medication adherence was measured on how closely patients followed their prescribed medication regimens including the time and dosage. To combat this by capitalizing on IoT and ML, \cite{Sixteen} developed a distributed cloud-based medication intake monitoring system paired with an Android smartwatch application that collected sensor data. The holistic multi-layer architecture presented by \cite{Thirteen} took things a step further by enabling real-time insights, helping improve the decision-making abilities of patients as well as healthcare providers. It not only employs machine learning on the cloud but also utilizes an artificial neural network for performing ML on-chip, all while ensuring less computation intensity to compensate for the limitations of low-end battery-powered devices. \par

In the typical cloud, fog, edge computing IoT stack, fully distributed fog-based systems were observed to offer unbound and limit-free operations but at the cost of accuracy due to computation limitations at the edge layer. To combine the best of both fog and cloud, a novel hierarchical fog-assisted computing architecture (HiCH) was presented by \cite{TwentyEight}, aimed at hierarchical partitioning of ML algorithms across all layers. A linear Support Vector Machine (SVM) classifier was chosen given its low computation cost and minimal performance and accuracy loss, to detect arrhythmia. The SVM classifier returned good results, but they certainly were not at the level and quality offered by deep learning algorithms. Since IoT-based health systems offer the ability to monitor patients beyond general clinical settings, results of the highest quality are desired as even a minor flaw in the results could prove fatal for patients. The authors of \cite{ThirtyTwo} presented a system that capitalized on the HiCH architecture by \cite{TwentyEight} and deployed a Convolutional Neural Network (CNN) classifier to enable deep learning. The proposed system was used to classify ECG signals to better monitor and detect arrhythmia. The federated learning framework proposed by \cite{Twenty} follows a similar approach but without HiCH. It partitions and allocates deep neural networks to IoT devices and a centralized server while significantly minimizing computation load on both. However, this approach resulted in a small accuracy loss when deployed on a real-world arrhythmia detection system. Comparing the three systems proposed above, the one proposed by \cite{TwentyEight} was observed to offer the best mix of both performances as well as accuracy. \par

Health monitoring solutions were not only limited to patients and quantified selves. Studies \cite{Six} and \cite{Eighteen} took things to the battlefield to track the location and monitor the health of soldiers who got lost, injured on the battlefield, or are deployed in adverse environmental conditions. This was achieved by the use of Wireless Body Area Networks (WBAN) comprising various sensors. The systems not only monitored the health metrics of soldiers, but also tracked their position, helped detect nearby bombs, and predicted potential war zones, helping army control units minimize operation times and plan their actions quicker \cite{Six}. \par

The real-time medical emergency response system proposed by \cite{Fifty} also used a WBAN comprising of a plethora of sensors, data generated by which was analyzed, and necessary actions were taken accordingly. The proposed system provides direct connectivity to the end-user to avail various kinds of medical help such as hospitals and ambulances, and also police stations. The connected entities are automatically contacted in the event of an emergency such as a car accident.

\subsubsection{\textbf{Diagnosing Health Conditions}}
The second most common application observed was the health diagnosis of patients as well as predicting health conditions based on the results of the diagnosis or otherwise. Not to intertwine with health monitoring which is mainly aimed at providing real-time insights, health diagnosis takes into account a plethora of data points to analyze and comprehensively assess the health of an individual and predict potential diseases. Early detection of diseases would lead to quicker actions being taken to treat them. The cloud-centric IoT-based m-healthcare disease diagnosing framework proposed by \cite{Nine} does exactly this as well as identifying the level of the severity, enabling healthcare providers to provide the correct treatment. A smart student healthcare prototype was designed to test this framework, which employed a Decision Tree (DT), k-Nearest Neighbor (k-NN), Naïve Bayes (NB), and SVM for disease classification along with the corresponding levels of severity. If the results of a diagnosis indicated an emergency, the nearest hospital was immediately notified to handle the emergency \cite{Nine}. In the journal paper by \cite{Eleven}, a similar but scalable cloud-based teleophthalmology architecture was proposed for predicting age-related macular degeneration (AMD) disease. The retinal images of patients were diagnosed by a modified ResNet CNN (AMD-ResNet) in the cloud to identify the level of AMD severity with high precision. Precision is critical in this scenario as the medication dosage would have to be prescribed accordingly. \par

The ML algorithms employed in the works mentioned above provided good results and predictions, but as discussed in Section 3.1.1, deep learning networks were identified to offer the highest quality and the most accurate predictions. The HealthFog framework by \cite{ThirtyFive} combines the power of complex deep learning networks with cloud-based edge computing devices to detect heart-related diseases with high accuracy and very low latency. The framework was built for efficient heart-related patient data management and diagnosing heart health while also identifying disease severity. \par

The rise in the number of wearable pieces of technology has made it easy for people to get access to their health data of all sorts. The large amount of sensor data generated by these wearables could be capitalized effectively by feeding it to ML models. However, the problem of data islands and data security and privacy come into play when doing this. FedHealth, a federated transfer learning framework was designed by \cite{FiftyOne} to overcome these challenges. By aggregating the sensor data via means of federated learning and building custom models by transfer learning, it provided quality results without hampering security and privacy. Given how general this framework is, it could easily be built upon and applied to various healthcare applications and scenarios. \par

Another approach to processing healthcare data is via means of big data technologies. These can handle extensive amounts of data efficiently. The works by \cite{FortyOne} and \cite{FortyFour} make use of various big data technologies such as Spark and Kafka to apply distributed machine learning to healthcare data to perform diagnoses and drive predictions. While these technologies ace when it comes to processing large amounts of data, they are often resource-hungry and need powerful computing hardware. These will be discussed in-depth in the coming sections.

\subsubsection{\textbf{Other Innovative Systems}}
Finally, a handful of new and unique applications and systems were also identified. The majority of the solutions discussed in the previous sections were built on the cloud, fog, edge IoT computing stack and employed machine learning on it. The architectures and frameworks proposed were tailored towards handling one kind of a problem - diagnosing a particular type of disease, for instance. While these contributions were breakthroughs in their specific niches, a more complex system comprising various IoT components was not something these were aimed at or thought of. This is where \cite{TwentyNine} studied and developed a health-related data sharing system between IoT devices and a central system, using the distributed ledger technology (DLT). The system benefitted from DLT by making it secure and highly scalable, two traits that are highly desired in a data-sharing system. This architecture could enable more comprehensive IoT systems comprising of various sensors to be designed for all kinds of diagnoses, and eventually, a single system that could diagnose all possible health conditions. \par

While the concept is a single system for everything sounds both interesting and complex, \cite{TwentyFive} built on this idea and presented a holistic architecture for an IoT eHealth ecosystem based on the cloud, fog, edge computing stack. Such an ecosystem would comprise various connected IoT-based medical devices monitoring a patient at the hospital for instance. This architecture would also enable a seamless transition from clinic-centric treatment, while the patient is at the clinic, to patient-centric, when the patient is discharged but continues the treatment at home or elsewhere, all while connected to each other at all times. The condition of the patient could then be monitored by healthcare providers from the clinic without having to visit the patient in person. Having an intelligent assisted living environment for home-based healthcare as presented by \cite{ThirtySeven} would also enable caregivers to tend to specific needs at any time. Given the criticality of the IoT-based medical devices in such a scenario, it is important that these devices work correctly at all times and provide correct insights. To ensure that these devices work correctly, they would have to undergo periodic preventive maintenance. The authors of \cite{ThirtyThree} investigated this and proposed autonomous integrity monitoring architecture that would provide complete visibility into the medical devices as well as predict potential failures before they actually happen. The architecture utilizes big data technologies to autonomously diagnose the components for wear and tear. \\

The number of applications and how ML and IoT in healthcare are being combined is growing every day. Looking at all the applications discussed in this section, we get a good overview and an answer to RQ1.

\subsection{Collaborative Machine Learning \& IoMT}
While we explored the different ways in which ML and IoT in healthcare are being combined today in section 3.1, we observed that all the works mentioned some form of collaborative computing. Both sound very similar but have a few differences. To be precise, 67\% of the primary studies discussed or made use of distributed machine learning, and the remaining 33\% proposed solutions involving distributed data processing and employing machine learning at a later stage. With the majority of the primary studies directly addressing distributed machine learning, it is evident that distributed machine learning can be applied to IoT healthcare data, answering RQ2.

\begin{figure}[H]
    \centering
    \includegraphics[width=8.7cm]{charts/dml_dc.PNG}
    \caption{Distributed ML vs Distributed Processing}
    \label{fig:dml_dc}
\end{figure}

\subsection{Distributed Machine Learning Algorithms}
Having explored the various applications of ML in healthcare IoT in section 3.1 and confirmed the possibility of applying distributed machine learning techniques to the same in section 3.2, this section evaluates the distributed machine learning algorithms and distributed computing frameworks that have been proposed in the primary studies. The best performing algorithms, in terms of prediction accuracy, are then identified for each healthcare scenario - live health monitoring systems, and medical diagnostic systems.

\subsubsection{\textbf{Live Health Monitoring Systems}}
The growing number of people suffering from chronic illnesses has led to the development and innovation of several live health monitoring systems. IoT-based live health monitoring systems are essential for such people to keep their health in check, and this could be a life or death scenario. Considering how critical such systems and their applications are, they need to produce accurate predictions in real-time. For benchmarking the standard model for IoT healthcare systems proposed by \cite{Seven}, a few ML algorithms were trained in the cloud on blood oxygen level data to monitor people's health and identify underlying conditions. The predictions produced were reasonably accurate compared to the actual results of patients undergoing surgery in a clinical setting.

The evaluation of ML algorithms, namely K-NN, Support Vector Machine, Decision Trees, Random Forest, and MLP, on publicly available datasets for various diseases by \cite{Eight} offered insights into their performances. Random Forest outperformed the rest in most scenarios, achieving an accuracy of 96.42\% for breast cancer, 97.26\% in dermatology, and 81.16\% for diabetes. The accuracy could be further improved by working with larger training datasets. To test the performance of distributed processing framework for the smartwatch app developed by \cite{Sixteen}, the authors tested four algorithms - Random Forest, Gradient-boosted Tree, Logistic Regression, and Support Vector Machine and compared their F1 scores. Gradient-boosted tree performed the best with an F1 score of 0.983 but took the most time to train, with Random Forest coming second at 0.977 and taking the least time to train.

The authors of \cite{Thirteen} discussed the applicability of several ML algorithms on different kinds of healthcare datasets and concluded that none of them are tailored for IoT applications. This led to the development of a custom shallow Artificial Neural Network for ML on-chip for low-end devices. The ANN achieved an accuracy of 95\%, which was comparable to the performance of a Convolutional Neural Network trained on the same ECG data on an edge node.

Study \cite{ThirtyTwo} put the HiCH architecture, built by \cite{TwentyEight}, to test by training a CNN classifier on ECG data from MIT's Arrhythmia dataset first to set a baseline. Then it was retrained on incoming ECG signals from the patient to be monitored to tailor the accuracy to the patient. With this method, the authors managed to achieve an accuracy of 96\%. While tailoring algorithms to the end-users yields higher accuracy, this approach is not feasible for time-critical applications. Similarly, the authors of \cite{Six} built upon the three-layer distributed IoT architecture proposed by \cite{Eighteen} and came up with a system to track the location and vitals of soldiers on the battlefield in real-time. K-Means Clustering was the algorithm of choice which classified different activities. The model produced results well enough to align with real-world scenarios, but accuracy numbers were not reported. \\

Looking at the overall findings from all the studies on live health monitoring systems, the majority of the studies concluded that Random Forest performed the best in terms of overall accuracy.

\subsubsection{\textbf{Medical Diagnostic Systems}}
Diagnosing health and vitals is as important as monitoring them. The authors of the study \cite{Nine} designed an architecture for smart healthcare focused on diagnosing heart-related diseases. Four different classification algorithms - Naive Bayes, k-Nearest Neighbor, Decision Tree, and Support Vector Machine,  were tested. Decision Tree came out at the top-yielding an accuracy of 92.8\%, and also yielded the highest values for specificity, sensitivity, and F-measure.

The cutting-edge cloud-based teleophthalmology framework built by \cite{Eleven} was designed to utilize AMD ResNet, which is a modified ResNet Convolutional Neural Network and a long–short term memory-based neural network in one. Upon testing, the model achieved an accuracy of 79.4\%, which is significantly higher than the accuracy of an expert ophthalmologist. The HealthFog framework by \cite{ThirtyFive}, built for diagnosing heart health and disease severity, used an ensemble deep neural network for analysis and binary classification. The neural network achieved an accuracy of 89\%, courtesy of how ensemble learning works.

The FedHealth framework built by \cite{FiftyOne} took the personalized healthcare model approach similar to the one taken by \cite{TwentyEight} and brought it to wearable devices, utilizing a CNN and optimizing with Stochastic Gradient Descent. The authors tested the framework's performance through fine-tuning and a transfer with maximum mean discrepancy - the first being a tuning approach and the second transfer learning. Transfer learning outperformed by 4\% in terms of average accuracy, with the authors claiming FedHealth produced higher accuracy for wearable activity recognition when compared to traditional methods \cite{FiftyOne}. \\

Looking at the overall findings from all the studies on medical diagnostic systems, the majority of the studies proposed their modified versions of Neural Networks which, as a result, helped attain the best accuracy numbers for each. No existing state-of-the-art machine learning algorithm stood out.

\section{Study Design}\label{sec:design}
This study was carried out by setting a research goal and determining research questions to help achieve it. The query string was then run on Google Scholar to obtain preliminary results on the studies available on the topic. To narrow down the scope of the search, inclusion and exclusion criteria were determined based on the research goal and questions, which helped remove studies that were irrelevant and select primary studies. The primary studies were then deeply analyzed and relevant data from them was extracted to help answer the research questions framed.

\label{sec:design}

\subsection{Research Goal}
This research was designed to create a path for researchers and practitioners to build upon and potentially introduce new technical contributions based on the results. We identify and compare popular distributed machine learning algorithms and frameworks used in state-of-the-art IoT healthcare systems based on their performance in terms of prediction accuracy. Since there is no one-solution-fits-all framework or algorithm available, we evaluate the performance of each in different scenarios and identify those that produce the best prediction results for each healthcare scenario.

\subsection{Research Questions}
To clarify the direction of this literature review, research questions were determined along with their rationales. The research questions drove the whole study and had a special influence on
\begin{enumerate}
    \item the search and selection of primary studies,
    \item the data extraction process,
    \item data synthesis and analysis.
\end{enumerate}

\label{sec:questions}
\begin{enumerate}
    \item[\textbf{{RQ}1:}] \textit{How has the healthcare industry been combining Machine Learning and the Internet of Things in recent years?}
    \item[{ }] The number of connected IoT devices generating health-related sensor data has grown exponentially in recent years. To make effective and efficient use of the massive amounts of data produced by these devices, various forms of AI - machine learning techniques, in particular, are being used to make smart personalized predictions. To get better insights, we review some of the state-of-the-art solutions being used in real-world applications and their evolution trends in the healthcare industry in the period 2016-2020.
\end{enumerate}

\begin{enumerate}
    \item[\textbf{{RQ}2:}] \textit{Can Distributed Machine Learning be applied to IoT Healthcare data?}
    \item[{ }] The most popular machine learning algorithms are only capable of handling small to medium-sized datasets. As the amount of data generated by IoT healthcare devices is bound to grow in volume over time, processing this data will become more and more computationally challenging. To continue driving quality predictions, all of the data must be leveraged. To mine data efficiently, the machine learning workload will have to be distributed across multiple processing agents.
\end{enumerate}

\begin{enumerate}
    \item[\textbf{{RQ}3:}] \textit{Which Distributed Machine Learning algorithms are the most suitable for processing IoT data in each healthcare scenario?}
    \item[{ }] Since no machine learning algorithm performs consistently and fits all scenarios in healthcare, we compare algorithms and distributed learning frameworks most used in state-of-the-art IoT healthcare systems based on their performance and the quality of the predictions they produce. The identification of such works and cutting-edge techniques would create a path for future technical contributions that build upon and improve existing solutions.
\end{enumerate}

\subsection{Initial Search}
This stage aimed to perform an automatic search on an electronic database and retrieve potential publications touching upon the topic of this study. Google Scholar was chosen for this as it searches across several academic and scientific databases. Google Scholar has gained credibility over the years as it been observed to have coverage advantages over controlled databases \cite{GoogleScholar}, making it a good choice for this time-boxed literature review. Next, keywords and query strings were identified and extracted from the research questions. The query string helped narrow down the scope and area of the search. The following query string was used for finding literature for this study. \\

\noindent \texttt{(Machine Learning) \textbf{AND} (Internet of Things \textbf{OR} IoT \textbf{OR} IoMT) \textbf{AND} (Health \textbf{OR} Healthcare \textbf{OR} Medicine \textbf{OR} Medical) \textbf{AND} (Distributed \textbf{OR} Federated) \textbf{AND} (Algorithm)} \\

After running the aforementioned query string and obtaining 115 results, publications that were outside the scope of this literature study were filtered out. Results such as articles, books, workshop proceedings, and the likes were removed and only 61 research papers were retained. From the remaining results, duplicates were then identified and removed, leaving us with 51 studies to apply the selection criteria to, discussed in the next section.

\begin{figure}[h!]
    \centering
    \includegraphics[width=8.7cm]{charts/years_all.PNG}
    \caption{Distribution of Selected Studies Over Years}
    \label{fig:years_all}
\end{figure}

\subsection{Application of Selection Criteria}
As mentioned in \cite{SysReviewProc}, the selection criteria are defined when defining the execution protocol of the systematic review. The 11 quality assessment criteria methodology by \cite{ApplyingSystematicReviews} was used to define the inclusion \textit{(I)} and exclusion \textit{(E)} criteria to enable us to narrow the scope further down and only select studies that were relevant.

\begin{enumerate}
\item[\textbf{{I}1:}] Studies that address machine learning in the field of IoT.
\item[\textbf{{I}2:}] Studies that address IoT-based healthcare systems.
\item[\textbf{{I}3:}] Studies that propose a new framework, system, architecture, tool, algorithm, or methodology involving collaborative machine learning in healthcare.
\item[\textbf{{I}4:}] Studies published in the domain of artificial intelligence and healthcare.
\item[\textbf{{I}5:}] Studies developed by either academics or industry.
\item[\textbf{{I}6:}] Studies written in English.
\item[\textbf{{E}1:}] Studies that do not consider the use of machine learning with IoT device data.
\item[\textbf{{E}2:}] Studies mainly targeting the hardware architectural side of IoT in healthcare or distributed machine learning.
\item[\textbf{{E}3:}] Studies not available as full-text.
\end{enumerate}

After filtering out selected papers using the defined inclusion and exclusion criteria, we were left with 28 publications. However, in order to carry out a valid systematic literature review, according to \cite{ApplyingSystematicReviews}, 15-20 publications should be selected as primary studies. Looking into and assessing the abstract, introduction, results, and conclusion of each, 21 publications were selected as primary studies for this review.

\begin{figure}[h!]
    \centering
    \includegraphics[width=8.7cm]{charts/years.PNG}
    \caption{Distribution of Primary Studies Over Years}
    \label{fig:years}
\end{figure}

\subsection{Data Extraction and Synthesis}
In this stage, the data from the primary studies were prepared for in-depth analysis to help answer the research questions. The type of publication was noted, and the studies were summarized based on the IoMT scenarios covered, distributed computing frameworks proposed, and the distributed machine learning algorithms used.

\begin{figure}[h!]
    \centering
    \includegraphics[width=8.7cm]{charts/study_type.PNG}
    \caption{Distribution of Primary Study Types}
    \label{fig:study_type}
\end{figure}

Based on the results of the data extraction process, each research question was answered extensively, limitations were identified, and conclusions were derived.


\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}