\section{\label{sec:DFT} Details of DFT and DFT+SO calculations}
Electronic properties of the Sn/Si(0001) system were simulated within generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional \cite{S_PBE} as implemented in the {\sc Quantum Espresso}~\cite{S_espresso} and Vienna ab-initio simulation package ({\sc VASP})~\cite{S_VASP,S_VASP_1}  with plane-wave basis set. The calculation parameters are as follows. The energy cutoff equals to 500 eV and the energy convergence criteria is 10$^{-8}$ eV. For the Brillouin zone integration a 20$\times$20$\times$1 Monkhorst-Pack mesh was used. A vacuum space of 16~\AA~between unit cell replicas in the vertical $z$ direction was introduced.  The optimized atomic positions and the band structure of the nonmagnetic ground state are in good agreement with previous works~\cite{S_Glass,S_Cho}. In particular, the corresponding Sn-Si and C-C distances are 2.05~\AA~and 0.18~\AA, respectively. The distance between nearest neighbor adatoms is 5.32~\AA, which is smaller than that for the Sn/Si(111) system, 6.61~\AA.

Fig.~\ref{fig:BANDS_SO} shows the calculated electronic spectrum of the Sn/Si(0001) system. The width of the energy band at the Fermi level is 0.29 eV, which is smaller than the corresponding value of 0.45~eV in Sn/Si(111). Taking into account spin-orbit coupling leads to the energy splitting, which is in agreement with Ref.~\onlinecite{S_Glass}. 

\begin{figure}[!h]
\includegraphics[width=0.7\textwidth]{FigS1.png}
\caption{ (Top) Band structures as obtained from DFT and DFT+SO calculations. The inset shows the Wannier function (top view) representing the band at the Fermi level. (Bottom) Energy splitting of the band in DFT+SO calculations due to spin-orbit coupling.}
\label{fig:BANDS_SO}
\end{figure} 

To construct the low-energy tight-binding model in the Wannier functions basis, we use the procedure of maximal localization \cite{S_Marzari,S_Marzari2} implemented in {\sc wannier90} package~\cite{S_wannier90}. The resulting Wannier function is well localized with the corresponding spread of 7.9~\AA $^2$, which is about twice as small as in Sn/Si(111), 16.8~\AA $^2$~\cite{S_Danis}. This difference can be explained by the fact that the Wannier function in Sn/SiC(0001) is less hybridized with the slab compared to Sn/Si(111). A significant delocalization of the electronic density (Fig.~\ref{fig:BANDS_SO}) suggests that magnetic properties cannot be reliably described only by using the atomic wave functions.



\section{\label{sec:Low-energy models }Construction of the low-energy model}
{\it Electronic model.}
In order to study the electronic and magnetic properties of the Sn/SiC(0001) system, we use an effective electronic model constructed in the Wannier functions basis taking into account spin-orbit coupling and electronic correlations: 

\begin{eqnarray}
\hat{\cal H}=\sum_{i j,\sigma\sigma'}t_{ij}^{\sigma\sigma'}\hat{a}_{i \sigma}^{+}\hat{a}_{j \sigma'} 
+ \frac{1}{2}\sum_{i j k l,\sigma\sigma'} K_{ijkl}\,\hat{a}_{i \sigma}^{+}\hat{a}_{j \sigma'}^{+}\hat{a}^{}_{l \sigma'}\hat{a}^{}_{k \sigma},
\label{eq:electron_ham}
\end{eqnarray}
where $i(j)$ and $\sigma (\sigma')$ are site and spin indices; $t_{ij}^{\sigma \sigma'}$ is the element of the hopping matrices with spin--orbit coupling;  $K_{ijkl}$  represents the element of the Coulomb interaction matrix. To define the model parameters, we used the state-of-the-art approaches based on the Wannier functions formalism, which are described below.

{\it Coulomb interactions.}
In the static limit ($\omega$ = 0), the bare value of the four-element Coulomb integral in the basis of Wannier functions $W_i(\boldsymbol{r})$  can be defined as:
\begin{equation}
K_{ijkl} = \int  \frac{ W^{*}_i(\boldsymbol{r}) W_k(\boldsymbol{r})  W^{*}_j(\boldsymbol{r}') W_l(\boldsymbol{r}') }{|\boldsymbol{r}-\boldsymbol{r}'|}  d\boldsymbol{r} d\boldsymbol{r}'.
\label{Coulomb_bare}
\end{equation}
Here, we consider the on-site and inter-site Coulomb  interactions, $V_{ii}= K_{iiii}$ and $V_{ij} = K_{ijij}$, respectively, as well as the nonlocal direct exchange interaction, $J^{F}_{ij} = K_{ijji}$. To capture the effects of screening, we adopt a numerical procedure previously applied to $sp$ systems \cite{S_Mazurenko,S_Rudenko}. Particularly, we use the constrained random phase approximation (RPA)~\cite{S_RPA}:
\begin{equation}
\kappa (\boldsymbol{q}) = [1 - K(\boldsymbol{q})P(\boldsymbol{q})]^{-1}K(\boldsymbol{q}),
\label{Ploarization}
\end{equation}
where $P(\boldsymbol{q})$ is the static single-particle polarizability function: 

\begin{equation}
P(\boldsymbol{q}) = \frac{1}{\Omega} \sum_{m \boldsymbol{k}}^{occ} \sum_{n \boldsymbol{k}'}^{unocc} \frac{| \langle  \Phi_{m \boldsymbol{k}} | e^{-i \boldsymbol{qr} } | \Phi_{n \boldsymbol{k}'} \rangle |^2 }    {\epsilon_{m \boldsymbol{k}} - \epsilon_{n \boldsymbol{k}'} +i \eta}.
\label{Polarization_RPA}
\end{equation}
In Eq.~\ref{Polarization_RPA}, $\boldsymbol{k}' = \boldsymbol{k} + \boldsymbol{q}$, summation runs over the Brillouin
zone including transitions from occupied $m$ to unoccupied $n$ states. Account of all possible transitions from occupied to empty states gives the fully screened value of the corresponding Coulomb integrals. At the same time, excluding ``self-screening'' transitions for the band at the Fermi level gives a partially screened value of $\kappa (\boldsymbol{q})$~\cite{S_silke2}. Here, $\Omega$ is the volume of the unit cell; $\epsilon_{m \boldsymbol{k}}$ and $\Phi_{m \boldsymbol{k}}$ stand for the eigenvalues and eigenvectors of the full DFT Hamiltonian, respectively; $\eta$ is a numerical smearing parameter chosen to be 10 meV. The bare $V$, partially screened $U$ and fully screened $W$ values of the on-site and inter-site Coulomb interactions were evaluated for Sn/SiC(0001), and are presented in Table~\ref{tab:Coulomb}. 

\begin{table}[!h]
\centering
\caption [Bset]{On-site (00) and inter-site (01) values (in eV) of the bare $V$, partially screened $U$ and fully screened $W$ Coulomb interactions, as obtained by using cRPA. See Fig.~\ref{fig:Model} for further details on the interaction paths. }
\begin{ruledtabular}
\begin {tabular}{c|cccc}
  
 path & $V$ & $U$ & $W$  \\
\hline
00 &  5.89 & 1.98 & 0.42 \\
01 &  2.45 & 0.80 & 0.34 \\
\end {tabular}
\end{ruledtabular}
\label{tab:Coulomb}
\end {table}

For comparison, we calculated the bare and screened values of the Coulomb interactions in Sn/Si(111). The obtained results are in good agreement with reported previously values~\cite{S_silke1}. For both systems, the bare parameters at long distances follow the classical Coulomb law $e^2/r$. For nearest and next-nearest neighbors, the deviation from the $e^2/r$ law can be explained by the fact that a particular shape of the constructed Wannier functions plays a dominant role when calculating their couplings. At the same time, at long distances the corresponding Wannier functions can be considered as point-like objects and their particular form is not important. Dependence of the Coulomb interactions on screening is almost the same for both the on-site and inter-site interactions (Fig.~\ref{fig:Coulomb}).  
\begin{figure}[!h]
\includegraphics[width=0.7\textwidth]{FigS2.png}
\caption{Calculated Coulomb interactions for Sn/SiC(0001) (top panel) and Sn/Si(111) (bottom panel) systems. The values denoted with tilde (green line) are taken from work \cite{S_silke1}.}
\label{fig:Coulomb}
\end{figure}

Generally, the direct exchange interaction parameter $J^{F}_{ij}$ is relatively small compared to the on-site and inter-site Coulomb integrals. As a consequence, the application of the cRPA scheme requires an extremely accurate integration over the Brillouin zone in order to calculate the polarizability function (Eq.~\ref{Polarization_RPA}), which is computationally challenging. We estimated the bare value of the direct exchange interaction,  $J^{F \, (bare)}_{01}$ = 5.20 meV, that defines the upper boundary for this coupling. To take into account the effects of screening, we used the renormalization factor of 3 obtained for the on-site and inter-site Coulomb interactions. Thus, one obtains $J^{F \, (screened)}_{01} = J^{F \, (bare)}_{01}/3$ = 1.73 meV. This value was used to define the total exchange interaction.   

{\it Hopping integrals.}
The calculated hopping parameters $t_{ij}^{\sigma \sigma'}$ in the Wannier functions basis are shown in Table~\ref{tab:hoppings}. Their diagonal parts are in excellent agreement with previously reported values obtained without spin-orbit coupling~\cite{S_Glass}.  As will be shown below, significant non-diagonal elements of the hopping matrix in spin space indicate that there are strong anisotropic magnetic excitations. 

\begin{table}[!h]
\centering
\caption [Bset]{Hopping integrals (in meV) between nearest and next-nearest neighbors as obtained from DFT+SO calculations. See Fig.~\ref{fig:Model} for further details on the interaction paths.}
\begin{ruledtabular}
\begin {tabular}{cccc}
   & $t_{01}$ & $t_{02}$ & \\
\hline
  & $\left(\begin{array}{cc} 27.94  &  7.39   \\ -7.39 & 27.94  \end{array} \right)$ &  $\left( \begin{array}{cc} -14.36 + 0.12 i &  -1.50 i \\  -1.50 i & -14.36 - 0.12 i \end{array}\right)$\\
\end {tabular}
\end{ruledtabular}
\label{tab:hoppings}
\end {table}

\section{Hartree-Fock solution}

The Hartree-Fock method provides the simplest approximation to the many-body problem \citep{S_Solovyev_HF}. Within this approach, the electronic model given by Eq.~\ref{eq:electron_ham} is solved as:

\begin{equation}
\left(\hat{t}_{\boldsymbol{k}}+\hat{\mathcal{V}}^{H}_{\boldsymbol{k}} + \hat{\mathcal{J}}^{H}_{\boldsymbol{k}} \right)|\varphi_{\boldsymbol{k}}\rangle=\varepsilon_{\boldsymbol{k}}|\varphi_{\boldsymbol{k}}\rangle,
\label{hf}
\end{equation}
\noindent where $\hat{t}_{\boldsymbol{k}}$ is the Fourier transform of the hopping parameters $\hat{t}_{ij}$, $\hat{\mathcal{V}}^{H}_{\boldsymbol{k}}$ and $\hat{\mathcal{J}}^{H}_{\boldsymbol{k}}$ are the Hartree-Fock potentials describing the Coulomb and direct exchange interactions, respectively, $\varepsilon_{\boldsymbol{k}}$ and $|\varphi_{\boldsymbol{k}}\rangle$ are the corresponding eigenvalues and eigenvectors in a given basis; a self-consistent solution of Eq.~(\ref{hf}) is achieved with respect to the density matrix:
\begin{equation}
\hat{n}=\sum\limits_{\boldsymbol{k}}|\varphi_{\boldsymbol{k}}\rangle \langle\varphi_{\boldsymbol{k}}|.
\end{equation}

\begin{figure}[!h]
\includegraphics[width=0.7\textwidth]{FigS3.png}
\caption{  Density of states (DOS) of Sn/SiC(0001) as obtained by using the Hartree-Fock approximation for the full model (Eq.~\ref{eq:electron_ham}). The ground-state magnetic order is shown in the inset.}
\label{fig:HF_DOS}
\end{figure} 

In our study, we simulated four magnetic structures: ferromagnetic (FM), 120$^\circ$-N\'eel (or 120$^\circ$-AFM), collinear row-wise (cRW or cAFM) and  120$^\circ$ row-wise (120$^\circ$-RW) structures for 3 $\times$ 3 unit cell. Due to strong spin-orbit couping, the FM and cRW configurations are found to be unstable. Having compared the energies of the noncollinear solutions, we conclude that the 120$^\circ$-RW state corresponds to the ground state of the Sn/SiC(0001) system (Fig.~\ref{fig:HF_DOS}). The calculated density of states (DOS) has the insulating gap of 2.0 eV, which is in excellent agreement with the results of photoemission experiment~\citep{S_Glass}. 

\section{Construction of the spin model}
To explore the excited magnetic states of the Sn/SiC(0001) system, we construct the spin model, given by Eq.~\ref{spinham} in the main text. The parameters of this model are estimated at the level of the superexchange theory~\cite{S_Anderson}. Particularly, the isotropic term consists of the antiferromagnetic kinetic part $J_{ij}^{kin}$ and the screened direct exchange interaction $J_{ij}^{F}$:
\begin{eqnarray}
J_{ij} =   \frac{1}{\widetilde U} {\rm Tr_{\sigma}} \{ \hat t_{ji} \hat t_{ij} \} -  J^{F}_{ij},
\label{exch}
\end{eqnarray}
where $\hat t_{ij}$ is the hopping integral with spin-orbit coupling, the effective local Coulomb interaction, $\widetilde U$ is estimated as $\widetilde U = U - V_{ij}$. The non-local Coulomb interaction plays a significant role for nearest neighbor interactions only. For the long-range interactions, its value is further reduced due to the exponential decay of the Wannier functions overlap and screening effects. The calculation gives $J_{01} = 1.42 - 1.73 = -0.31$ meV, whereas for the next-nearest interaction $J_{02}$ = 0.21 meV.

As for the anisotropic part of the magnetic model, the antisymmetric Dzyaloshinskii-Moriya (DM) and symmetric anisotropic exchange interactions are given by~\cite{S_Aharony1}:
\begin{equation}
\boldsymbol{D}_{ij} = \frac{i}{2\widetilde{U}}[{\rm Tr}(\hat {t}_{ij}){\rm Tr}(\hat{t}_{ji}\boldsymbol{\sigma})-{\rm Tr}(\hat{t}_{ji}){\rm Tr}(\hat{t}_{ij}\boldsymbol{\sigma})],
\label{eq:DM-vector}
\end{equation}
\begin{equation}
\overset{\leftrightarrow}{\Gamma}_{ij} = \frac{1}{2\widetilde{U}}[{\rm Tr}(\hat{t}_{ji}\boldsymbol{\sigma})\otimes {\rm Tr}(\hat{t}_{ij}\boldsymbol{\sigma})+{\rm Tr}(\hat{t}_{ij}\boldsymbol{\sigma})\otimes {\rm Tr}(\hat{t}_{ji}\boldsymbol{\sigma})],
\label{eq:Gamma-matrix}
\end{equation}
where $\boldsymbol{\sigma}$ are the Pauli matrices. The DMI vectors depicted in Fig.~2 paths of interactions are the following:  ${\bf D}_{01} =  (0, 0.70, 0)$, ${\bf D}_{02} = (0.043,  0, -0.003)$ meV. In accordance with the $C_{3v}$ symmetry, the nearest neighbor DMI vectors lie in $xy$-plane and are perpendicular to the corresponding bond, whereas the next-nearest neighbor DMI vectors alternate their $z$ component. 

The calculated symmetric anisotropic exchange interaction tensor has the following form (in meV):
\begin{eqnarray}
\Gamma_{01} = \left( \begin{array}{ccc} 0.00  &  0.00 & 0.00  \\ 0.00  & 0.19 & 0.00 \\ 0.00 & 0.00 & 0.00  \end{array}\right).
\end{eqnarray}
The tensors for the rest of the bonds can be obtained by 60$^{\circ}-$rotations around the z axis. 

\section{\label{sec:Magnetic characteristic}Analysis of the DFT magnetic moments and total energies }
In order to elucidate the magnetic structure and capture the strong delocalization effects, we performed the calculations by using an all-electron full-potential linearized augmented-plane wave method (FP-LAPW) implemented in the ELK code \cite{S_ELK}. Without spin-orbit coupling a significant interstitial contribution to the band at the Fermi energy (Fig.~\ref{fig:BANDS_SO}) was revealed. It is interesting to note that the magnetic moment has a similar composition. Spin-polarized calculations for the ferromagnetic state give 1 $\mu_B$, which contains only $\sim$50\% of the total moment at muffin-tin spheres, while the other half belongs to the interstitial area (Table.~\ref{tab:Magnetic_moments}). 

\begin{table}[!h]
\centering
\caption [Bset]{The total value and partial contributions to the magnetic moment obtained for the ferromagnetic ordering by using different first-principles methods (in $\mu_B$). Notation $MT$ and $I$ stands for sum of the moments in muffin-tin spheres and interstitial contribution, respectively.}
\label {basisset}
\begin{ruledtabular}
\begin{tabular}{lccc}
   & Quantum Espresso & VASP & ELK  \\
  \hline
 M$_{Sn}$  & 0.15 & 0.23 & 0.21   \\
 M$_{MT}$  & 0.35 & 0.52 & 0.51   \\
 M$_{Tot}$ & 1.00 & 1.00 & 0.96   \\
 M$_{I}$   &  -   &   -  & 0.45   \\

\end{tabular}
\end{ruledtabular}
\label{tab:Magnetic_moments}
\end {table}

\begin{figure}[!h]
\includegraphics[width=0.7\textwidth]{FigS4.png}
\caption{ Ferromagnetic (FM), collinear antiferromagnetic (cAFM) and ferrimagnetic (FI) configurations simulated in previous works \cite{S_Glass, S_Cho_Si_111}. Arrows represent the directions of local magnetic moments. Red circles denote the magnetic unit cell.}
\label{fig:cAFM}
\end{figure} 

In the previous calculations~\cite{S_Glass, S_Cho}, the authors considered ferromagnetic (FM), ferrimagnetic (FI) and two antiferromagnetic magnetic configurations: collinear (cAFM) and non-collinear 120$^\circ$ (120$^\circ$-AFM) ordering (Fig.~\ref{fig:cAFM}). It was reported that cAFM corresponds to the minimum of the total energy in Sn/SiC(0001). By using the calculated total energies of different magnetic configurations it is possible to estimate the individual magnetic interactions, between nearest and neighbors, $J_{01}$ and $J_{02}$, respectively. For that purpose, we consider the Heisenberg spin Hamiltonian $ \mathcal {H} = \sum_{ij} J_{ij} {\bf S}_{i} {\bf S}_{j}$. For convenience, a $3 \times 3$ unit cell is used. In this case, the energies for these magnetic configurations can be written in the following form:
\begin{equation}
 \begin{aligned}
    3E_{FM} = (18J_{01} + 18J_{02}) \times \frac{1}{4} \\        
    3E_{cAFM} = -(6J_{01} + 6J_{02}) \times \frac{1}{4} \\
    3E_{FI} = - (6J_{01} - 18J_{02}) \times \frac{1}{4} 
  \end{aligned}
\label{eq:Energies}
\end{equation}
Here, the energies $E_{FM}$, $E_{cAFM}$ and $E_{FI}$ are given per $\sqrt{3} \times \sqrt{3}$ unit cell. From these equations one can derive the expressions for nearest and next-nearest neighbors isotropic exchange interactions:

\begin{equation}
 \begin{aligned}
     J_{01} = \frac{1}{2} (E_{FM} - E_{FI}) \\
     J_{02} = \frac{1}{2} (E_{FI} - E_{cAFM}) \\
  \end{aligned}
\label{eq:J_calculation}
\end{equation}
The resulting interactions are given in Table~\ref{tab:Energies} in the main text.

\begin{table}[!h]
\centering
\caption [Bset]{Energies and corresponding values of isotropic exchange interactions (in meV) for Sn/SiC(0001) and Sn/Si(111) estimated by using the results of previous DFT calculations with different exchange-correlation (xc) functionals. Energies are given in meV per $\sqrt{3} \times \sqrt{3}$ unit cell.  Values in brackets are those evaluated with superexchange theory in this work and in Ref.~\onlinecite{S_Danis}. }
\label {basisset}
\begin{ruledtabular}
\begin{tabular}{ccccccc}
 &  xc-type &  $E_{FM}$ & $E_{cAFM}$ & $E_{FI}$ & $J_{01}$ & $J_{02}$ \\
  \hline
  {Sn/SiC(0001)} &  & & & & & \\
  \\
Glass {\it et al.}, Ref.~\onlinecite{S_Glass} & LDA & 3.03 & $-$2.93 & $-$2.03 &  2.53\,(1.42)  & 0.45\,(0.21) \\
   \hline
   \\
  {Sn/Si(111)} &  & & & & & \\
  \\
Cho {\it et al.}, Ref.~\onlinecite{S_Cho_Si_111} & hybrid & $-$47.9 & $-$65.1 & $-$64.2 &  8.15\,(7.69) & 0.45 \,(0.73) \\
 
\end{tabular}
\end{ruledtabular}
\label{tab:Energies}
\end {table}

\section{\label{sec:Monte Carlo simulation}Monte Carlo simulation}
Monte Carlo simulations were performed by using the heat-bath method combined with overrelaxation~\cite{S_Binder}. The corresponding model parameters are given up to next-nearest neighbors. In these calculations, supercells of various sizes from $N=96\times96$ to $150\times150$ spins with periodic boundary conditions were used, and a single run contained $(0.5-2.0)\cdot 10^{6}$ Monte Carlo steps. For initial relaxation the system is gradually cooled down from higher temperatures. 


\section{Scanning tunneling microscopy (STM) simulations }
Spin-polarized scanning tunneling microscopy technique \cite{S_Wiesendanger} is a powerful tool for exploring distinct spin textures on a surface, such as spin spirals or skyrmions. 
Theoretical background for such experiments was developed in Refs.~\cite{S_Fransson,S_Blugel}. According to these studies, the tunneling current can be written in the following form: 
\begin{eqnarray}
I(\boldsymbol{r},eV) \sim \int [f(\epsilon)-f(\epsilon-eV)]\times[n(\epsilon-eV)N(\boldsymbol{r},\epsilon)+  \boldsymbol{m}(\epsilon-eV)\boldsymbol{M}(\boldsymbol{r},\epsilon)] d \epsilon,
\label{eq:STM}
\end{eqnarray}
where $n(\epsilon)$ stands for the tip's DOS, $N(\boldsymbol{r},\epsilon)$ is the DOS of the surface. The magnetic part of the tunneling current contains the tip's and surface magnetization densities, $\boldsymbol{m}(\epsilon)$ and $\boldsymbol{M}(\boldsymbol{r},\epsilon)$, respectively. $V$ denotes the bias applied to the sample. To define the electronic and magnetization densities, we employed the calculated Wannier functions, $W(\boldsymbol{r})$. It allows us to take into account the delocalization of the electronic state and magnetic moment. Thus, one has $N(\boldsymbol{r},\epsilon) = |W(\boldsymbol{r})|^2 N_{loc}(\epsilon)$ and $\boldsymbol{M}(\boldsymbol{r},\epsilon) = |W(\boldsymbol{r})|^2 \boldsymbol{M}_{loc} N_{loc}(\epsilon)$, where  $N_{loc}(\epsilon)$ and $\boldsymbol{M}_{loc}$ is the local DOS and magnetic moment of adatom obtained from DFT calculations and Monte Carlo simulations, respectively.  Practically, we assumed that the tip's DOS $n(\epsilon)$  and magnetization $\boldsymbol{m}(\epsilon)$ do not depend on energy.   

