
\section{Heuristic Derivation of the Generalized Widom Formula}
In $d$ dimensions we consider a specific real-space partition in which the boundary between the two subsystems is a plane whose normal direction is $\hat{n}_d$. This partition preserves the translational symmetries in $d-1$ dimensions perpendicular to $\hat{n}_d$, and one can perform partial Fourier transformation for all the physical degrees of freedom along these $d-1$ axes, since the momenta $k_{1,2,\cdots, d-1}$ are good quantum numbers. We thus view the momentum space as consisting of arrays of parallel 1D chains with spacings $\delta k_{1,2,...,d-1} = 2\pi/L_\bot$, where $L_\bot$ is the linear size of these transverse directions.

Using well-established results for free fermions and coupled harmonic lattice systems with critical surfaces \cite{Wolf, GioevKlich, dingprx12, Lai_EE4EBL}, we assume that each 1D chain in momentum space intersecting the critical Fermi surfaces (critical points) contributes a 2nd Renyi entanglement entropy (REE) $(c_{eff}/4) \ln \mathcal{L}_\|$ [or a von Neumann EE ($\nu$EE) $(c_{eff}/3) \ln \mathcal{L}_\|$] to the total leading REE~\cite{Calabrese2004, Calabrese2009}, where $\mathcal{L}_\|$ is the linear size of the (smaller) subsystem along $\hat{n}_d$ and $c_{eff}$ represents the effective number of free gapless modes for each 1D chain. Note that we only consider the universal part of the leading terms in REE \cite{Calabrese2004, Calabrese2009}, \textit{i.e.}, for each 1D chain the leading REE explicitly should be $\frac{c_{eff}}{4} \ln \mathcal{L}_\parallel + c_2$ \cite{Calabrese2009}, where $c_2$ is a non-universal constant that we ignore. For free fermions or Fermi liquids (FL) \cite{GioevKlich, Zhang11, dingprx12, Lai_EEreconstruct}, $c_{eff} = c_{F}=1$; for coupled harmonic lattice models realizing the lattice version of the Exciton Bose liquid phase (EBL)~\cite{Paramekanti02}, $c^{EBL}_{eff} =2$~\cite{Lai_EE4EBL}. For the Gutzwiller-projected wave function in 2D, the $c_{eff}$ is not known. Nevertheless, the leading REE can be obtained by counting the total number of chains (in momentum space) intersecting the critical surface, which corresponds to the critical surface cross-sectional area divided by the $(d-1)$ dimensional area spacing between the chains, {\it i.e.,} $(2\pi/L_\bot)^{d-1}$. Explicitly, the leading universal part of REE is
\begin{eqnarray}
S_{dD} &\dot{=}& \frac{c_{eff}}{4} \ln \mathcal{L}_\|  \frac{1}{2}  \frac{\int_{\partial \Gamma} \left| d\hat{S}_\Gamma \cdot \hat{n}_d\right|}{(2\pi/L_\bot)^{d-1}} \label{Eq:Lai}\\
\nonumber & = & \frac{c_{eff}}{8}\ln \mathcal{L}_\| \left( \frac{L_\bot}{2\pi}\right)^{d-1} \frac{\int_{\partial \Gamma}\int_{\partial A} \left | d\hat{S}_\Gamma \cdot d\vec{S}_A\right|}{2L_\bot^{d-1}} \\
& = & \frac{c_{eff}}{16} \frac{\ln \mathcal{L}_\|}{(2\pi)^{d-1}} \int_{\partial A} \int_{\partial \Gamma} \left| d\vec{S}_A \cdot d\hat{S}_\Gamma \right|,\label{supp:eq_WGK}
\end{eqnarray}
where $\dot{=}$ represents the leading contribution. $\mathcal{L}_\|$ is the linear size of the (smaller) subsystem along $\hat{n}_d$. The factor $1/2$ in the first line is due to the over counting of the cross-section. In the second line, we rewrite $\hat{n}_d$ as a real-space partition surface integral (with $d\vec{S}_A$ being the corresponding oriented area element whose direction is along the local normal direction) divided by the partition surface area in $d-1$ dimensions, $2L_\bot^{d-1}$. $\int_{\partial \Gamma}$ represents the surface integral along the critical surface in momentum space (with $d\hat{S}_\Gamma$ being the corresponding oriented area element). While we arrived at Eq.~\eqref{supp:eq_WGK} by considering a special partition, it is actually the correct formula for the free-fermion state for arbitrary cuts \cite{GioevKlich} if we set $c_{eff}=1$ and replace $\mathcal{L}_\|$ by the generic linear size of the smaller subsystem. For the Gutzwiller-projected wave function, if we consider an $L_y$-legged chain system with infinite length along the $x$-axis, the total number of gapless modes is $2L_y -1$.  In real 2D systems, taking periodic boundary conditions along the $y$ direction, we expect $c_{eff}$ for each line in momentum space to be $c_{eff} = 2 - 1/L_y\big{|}_{L_y\rightarrow \infty} \rightarrow 2$.

\begin{figure*}
   \centering
   \includegraphics[width=\linewidth]{ee_fs_uni_lattice.pdf}
   \caption{(a) Illustration of the different subsystem geometries on the triangular lattice that we consider for the REE calculations. We set $t\equiv 1$, and consider $t'/t = 1.0$ for the isotropic model. The colored regions represent the subsystems that we considered. The blue and red regions represent subsystems with size $L_A \times \alpha L_A$, where $\alpha =1/2, 1$. The green region represents the subsystem preserving translation along the $\bm{a}_2$ direction. (b) The REE, $S_2$, for different subsystem geometries for the isotropic case. We choose two setups to extract REE. For red squares, green diamonds and blue circles, we fix the whole system size and choose the subsystem size to be $L_A \times \alpha L_A$ with $\alpha = 1/2$ (blue circles), $1$ (red squares and green diamonds). For purple triangles, we fix the ratio of the linear size of the subsystem and that of the whole system to be $1/2$ ($L_A / L = 1/2$) and vary the value of $L$ to extract the REE. Based on Eq.~\eqref{supp:eq_WGK}, we plot $S_2/((1+\alpha)L_A)$ or $\gamma S_2/L_A$ vs $\ln (L_A)$. The slopes of the lines give the prefactor of the leading REE. }
\label{Fig:ee_fs_uni}
\end{figure*}

\section{REE by Variational Monte Carlo}
The variational wave function used here for a critical spin liquid with a spinon Fermi surface (SFS) is defined as
\begin{eqnarray}
|\psi\rangle=\mathcal{P}_{G}|\psi_{0}\rangle,
\end{eqnarray}
where $\mathcal{P}_{G}=\prod_{i}(1-n_{i\uparrow}n_{i\downarrow})$ is the Gutzwiller projector,
which enforces no double occupation on each site.
$|\psi_{0}\rangle$ is the ground state of the following mean-field Hamiltonian on the triangular lattice:
\begin{equation}\label{eq:SL}
{\cal H}_{\rm MF} = \sum_{\langle i,j\rangle,\sigma} t_{ij}c^{\dag}_{i,\sigma}c_{j,\sigma} + h.c.
\end{equation}
In a system with subsystems $A$ and $B$, the REE of order $n$ in $A$ is defined as
\begin{eqnarray}
S_{n}=\frac{1}{1-n}log[Tr\rho_{A}^{n}],
\end{eqnarray}
with the reduced density matrix associated with $A$, $\rho_{A}=Tr_{B}|\psi\rangle\langle\psi|$. We will focus on $n=2$ REE, $S_2$. To compute $S_2$, as in Ref.~\cite{Zhang11} we introduce an identical copy of the original system:
We divide the original system into two subsystems $a$ and $b$ and, likewise, the replica into $a'$ and $b'$. 
The operator $Swap$ is introduced as $Swap|a,b\rangle|a',b'\rangle=|a,b'\rangle|a',b\rangle$, and the REE is
\begin{eqnarray}
e^{-S_{2}}=\frac{\langle\Psi|Swap|\Psi\rangle}{\langle\Psi|\Psi\rangle},
\end{eqnarray}
with $|\Psi\rangle$ the wave function of the product between the original system and its replica. The operator
\begin{eqnarray}
\langle Swap\rangle=\underset{a,b;a',b'}{\sum}P(a,b,a',b')\frac{\psi(a',b)\psi(a,b')}{\psi(a,b)\psi(a',b')}
\end{eqnarray}
is calculated according to the weight $P(a,b,a',b')=\frac{|\psi(a,b)|^{2}|\psi(a',b')|^{2}}{\underset{a,b}{\sum}|\psi(a,b)|^{2}\underset{a',b'}{\sum}|\psi(a',b')|^{2}}$ by variational Monte Carlo (VMC). In this method, we generate the Markov chain for each configuration of both the original system and its replica, and use the Metropolis algorithm \cite{metropolis} to update the configurations according to the probability distribution $P(a,b,a',b')$. 

\begin{figure}
    \subfigure[von Neumann entanglement entropy]{\label{Fig:FFnuE_1to1}\includegraphics[width=\linewidth]{FFnuE_121.pdf}}
    \subfigure[Renyi entanglement entropy]{\label{Fig:FFRE_1to1}\includegraphics[width=\linewidth]{FFRE_121.pdf}}
    \caption{The Fermi surface size of free fermions is fixed to be equal to the one obtained from the spin structure factor. Here we choose the subsystem with $L_A \times L_A$ lattice sites along the $\bm{a}_1$ and $\bm{a}_2$ directions and calculate the von Neumann entanglement entropy ($\nu$EE) (a) and the REE (b).  (a) We find that the fitting line is $y\simeq 0.278 x + 0.355$. (b) We observe a stronger oscillating behavior in the Renyi entropy, which makes it harder to obtain a conclusive fitting line. Based on the current data, we get the fitting line to be $y\simeq 0.194 x + 0.331$.}
\label{Fig:FFEE_121}
\end{figure}

\begin{figure}
    \subfigure[von Neumann entanglement entropy]{\label{Fig:FFnuE_bipartite}\includegraphics[width=\linewidth]{FFnuE_bipartite.pdf}}
    \subfigure[Renyi entanglement entropy]{\label{Fig:FFRE_bipartite}\includegraphics[width=\linewidth]{FFRE_bipartite.pdf}}
    \caption{The Fermi surface size of free fermions is fixed to be equal to the one determined from the spin structure factor result. Here we choose the subsystem to preserve the translation symmetry along the $\bm{a}_2$ direction to obtain the $\nu$EE (a) and REE (b). In both cases, we find stronger oscillating behaviors compared with those in Fig.~\ref{Fig:FFEE_121}, but we still can obtain a consistent result for the $\nu$EE. (a) We find the fitting line to be $y = 0.282 x + 0.424$, whose slope is comparable to that in Fig.~S\ref{Fig:FFnuE_1to1} obtained in a different subsystem geometry. (b) The REE shows a very strong oscillating behavior, which makes it very difficult to find a fitting line. A linear fitting line based on the current data is $y = 0.260 x + 0.199$.}
\label{Fig:FFEE_bipartite}
\end{figure}

\section{REE of the Isotropic Case by Variational Monte Carlo}
The results for REE in the isotropic case $(t=t')$, with $t/t'$ representing the fermion hopping amplitudes in the Hamiltonian~Eq.~(\ref{eq:SL}) (prior to Gutzwiller projection), are shown in Fig.~\ref{Fig:ee_fs_uni}. Panel (a) illustrates three different subsystems that we use for calculating the REE. The red and blue regions represent, in general, the subsystem with $L_A \times \alpha L_A$ sites, while the green region represents the subsystem with $L_A \times L$ sites with periodic boundary condition along the $\bm{a}_2$ direction, where $L$ represents the full system length. In panel (b), the slopes of the lines give the value of $c_{eff} A_{sf}$. The red squares are consistent with the data in Ref.~\cite{Zhang11} on a $18 \times 18$ triangular lattice with periodic boundary condition along both the $\bm{a}_1$ and $\bm{a}_2$ directions. The green diamonds represent our numerical data on a system with $20 \times 20$ sites along the $\bm{a}_1$ and $\bm{a}_2$ directions. We remark that the consistency between these two results shows the convergence of REE from $L=18$ to $L=20$ clusters, and then the $L=18$ cluster is sufficient to illustrate our VMC results. The open purple triangles represent the REE associated with the subsystem, which is half of the whole system in our setup, with periodic boundary conditions along $\bm{a}_2$. Within the error bars, we can see that the slope of the latter is \textit{equal} to the former obtained using a different subsystem geometry. The equality can be explained using the formula, Eq.~\eqref{supp:eq_WGK}, which we elaborate below.

In the first setup with $\alpha=1$, the REE is associated with the subsystem with four boundaries with equal surface area $L_A$ on a 2D triangular lattice with $L\times L$ sites. For each surface in real space, the integral along the real-space surface contributes $2 L_A A_{sf}$ to REE leading to $S^I_2\doteq c_{eff} A_{sf} /(4\pi)L_A \ln (L_A)$ based on Eq.~\eqref{supp:eq_WGK}. On the other hand, if the subsystem only preserves translational symmetry along the $\bm{a}_2$ direction, in real space there are only two surfaces with surface area $L$. For each boundary the integral gives $2 L A_{sf}$ that results in REE as $S^{II}_2 \doteq c_{eff} A_{sf} /(8\pi)L \ln (L_A) = (2\gamma)^{-1} c_{eff} A_{sf} /(4\pi) L_A \ln (L_A)$, where we rewrite $L = L_A/\gamma$. Setting $\gamma = 1/2$ gives $S^I_2 \doteq S^{II}_2$. The consistency between the leading terms of $S^{I}_2$ and $S^{II}_2$ as shown in Fig.~\ref{Fig:ee_fs_uni}(b) suggests the applicability of the Widom formula for the gapless spin liquid with a SFS described by the Gutzwiller-projected wave function. 

To further illustrate the applicability of the Widom formula to the gapless spin liquid states, we calculate the REE associated with a subsystem with different lengths along the $\bm{a}_1$ and $\bm{a}_2$ directions, which is chosen to be $L_A \times L_A/2$ in our work. If the Widom formula is applicable, the slopes of $S_2/((1+\alpha) L_A)$ vs $\ln(L_A)$ should be the same for different values of $\alpha$. The blue open circles in Fig.~\ref{Fig:ee_fs_uni}(b) are the data obtained in $L_A \times 1/2 L_A$ subsystems, and we can see that the slope is comparable to the previous results using different subsystem setups, which again suggest the validity of the Widom formula in this system.

\section{$\nu$EE and REE for Spinless Free Electrons in the Isotropic Case}
For free-fermion systems, we use the correlation function method \cite{Peschel2003} to obtain the $\nu$EE and REE in Figs.~\ref{Fig:FFEE_121} and \ref{Fig:FFEE_bipartite}. We use two kinds of subsystem setups illustrated in Fig.~\ref{Fig:ee_fs_uni}(a):  Setup 1: $L_A\times L_A$ subsystem [red region in Fig.~\ref{Fig:ee_fs_uni}(a)]; Setup 2: $L_A \times L$ subsystem with a periodic boundary along the $\bm{a}_2$ direction [green region in Fig.~\ref{Fig:ee_fs_uni}(a)]. The results are illustrated in Figs.~\ref{Fig:FFEE_121}-\ref{Fig:FFEE_bipartite}. In Fig.~\ref{Fig:FFEE_121}, we plot $S^F_\nu/((1+\alpha) L_A)$ $[S^F_2/((1+\alpha) L_A]$ vs $\ln(L_A)$, with $\alpha = 1$. In Fig.~\ref{Fig:FFEE_bipartite}, we plot $\gamma S^F_\nu/L_A$ $(\gamma S^F_2/L_A)$ vs $\ln(L_A)$, with $\gamma$ being the ratio between the subsystem length and that of the whole system, $\gamma = L_A/L = 1/2$. The numerical calculations have been performed up to the triangular lattice with $L=72$. In both cases, we observe stronger oscillating behaviors for the $S_2$ data in the free fermion systems, which makes it difficult to obtain the fitting lines. Focusing on the von Neumann entropy data, $S_\nu^F$, we find that the two different setups give comparable results. If we average the two slopes obtained in these two setups, we get the slope to be $\sim 0.28$. Utilizing the theoretical understandings that $S_\nu^F \dot{=} 4/3 S_2^F$ (the universal part of the entanglement entropy), we find that the theoretical slope for the REE data to be $0.21$, which is comparable to the average of the slopes in the two REE data in Figs.~\ref{Fig:FFEE_121}(b) and \ref{Fig:FFEE_bipartite}(b). The slope we obtained in this spinless free-fermion system is half of the value obtained in the Gutzwiller-projected wave function system, which again suggests that $c_{eff}\simeq 2$.

\begin{figure*}[t]
   \centering
   \includegraphics[width=\linewidth]{sq_ani07.pdf}
   \caption{(a) Side view of the spin structure factor within BZ for the anisotropic case $t'/t = 0.7$, where the black hexagon represents the BZ. There is a sharp and singular peak at ${\bf q}=0$ corresponding to the uniform real-space power-law decaying behavior. The much weaker singular lines near the boundary of the BZ correspond to the oscillating real-space behavior caused by the presence of the SFS. (b) Top view of the spin structure factor. The weak singular lines on the surface are due to the presence of the emergent SFS and the theoretical locations of these lines are at $\bm{k}_{FR} - \bm{k}_{FL}$, \textit{i.e.}, for an observation direction we can define a momentum $\bm{k}_{FR/L}$ associated with a right-moving (left-moving) patch, represented by the blue lines. The theoretical blue lines match with remarkable accuracy  with the singular lines obtained by numerical calculations on a triangular lattice consisting of $30\times 30$ sites.
 }
   \label{Fig:sq_ani07}
 \end{figure*}

\section{Long Wavelength Analysis of the Spin Structure Factor}
To establish the precise formula for the leading entanglement entropy of the Gutzwiller-projected wave function describing the gapless spin liquid with a SFS, we can numerically determine $c_{eff}$ by comparing the result of the spin structure factor $\mathcal{D}_{\bf q}$ and the REE results, where
\begin{eqnarray}\label{Eq:Sq}
\mathcal{D}_{\bf q} \equiv \sum_j \chi^s_j e^{-i {\bf q} \cdot {\bf r}_j},
\end{eqnarray}
and $\chi^s_j $ is the real-space spin correlation functions defined as
\begin{eqnarray} \label{Eq:spincorr}
\chi^s_j \equiv \sum_{\mu =x, y, z} \langle S^\mu_j S^\mu_{i=0}\rangle.
\end{eqnarray}
The emergent SFS in momentum space can be viewed as consisting of patches of critical surfaces which at the low-energy descriptions are independent to each other \cite{SSLee2009}. It is expected that these independent patches contribute \textit{equally} to REE with the same $c_{eff}$. To determine $c_{eff}$, we focus on the isotropic critical surface case, \textit{i.e.}, a circular critical surface. For an arbitrary observing direction $\hat{n}$, we can determine two momenta corresponding to the right-moving patch, $\bm{k}^n_{FR}$, and the left-moving patch, $\bm{k}^n_{FL}$, respectively. For an isotropic convex critical surface with an inversion center, the vector ${\bf A}_s^n \equiv \bm{k}^n_{FR} - \bm{k}^n_{FL}$ must pass through the center of the critical surface and, therefore, the cross-section of the isotropic convex critical surface can be determined by the length $A_s^n = |{\bf A}_s^n|= |\bm{k}^n_{FR} -\bm{k}^n_{FL}|$. 

The wave vectors $\bm{k}^n_{FR} - \bm{k}^n_{FL}$ can be extracted by examining the spin structure factor. The power-law correlations in real space correspond to singularities in momentum space that can be revealed in the spin structure factor. At mean-field level, the power-law behavior can be explicitly determined. In general, the low-energy description dictates that the spin structure factor should show singularities at ${\bf q} = {\bf 0},~\bm{k}^n_{FR } - \bm{k}^n_{FL}$. Identifying $\bm{k}^n_{FR} - \bm{k}^n_{FL}$ along different directions can determine the cross-sections of the SFS at different directions. The exact result for the isotropic case is illustrated in Fig.~2 in the main text, which shows wiggling lines on the surface of the spin structure factor corresponding to the weak singularities of $\bm{k}^n_{FR} - \bm{k}^n_{FL}$. To numerically fit  the exact $\bm{k}^n_{FR} - \bm{k}^n_{FL}$, we assume that the Gutzwiller projection does not dramatically change the locations of the singularities, while only the exponents of the power-law behaviors of the singularities are modified. We then adopt a mean-field fermionic state with a Fermi surface at $1/2$-filling. Extracting $\bm{k}^n_{m,FR} - \bm{k}^n_{m,FL}$ of the mean-field ansatz, we find that the $\bm{k}^n_{m,FR} - \bm{k}^n_{m, FL}$ can fit the exact $\bm{k}^n_{FR} - \bm{k}^n_{FL}$, which indeed suggests that the Gutzwiller projection does not change dramatically the geometric information of the SFS in this system. 

At mean-field level, the Hamiltonian for an isotropic fermionic tight-binding model without enlargement of unit cells is
\begin{eqnarray}
H_{mf}  =  \sum_\alpha \sum_{\langle {\bf r}, {\bf r}' \rangle} f^{\alpha \dagger}({\bf r}) A({\bf r} - {\bf r}') f^\alpha ({\bf r}'),
\end{eqnarray}
where $\langle...\rangle$ represents nearest-neighbors and $\alpha = 1,2$ represents the spin flavor of fermions. $f^{\alpha \dagger}({\bf r}) (f^\alpha({\bf r}))$ represents the fermion creation (annihilation) operator with flavor $\alpha$ at location ${\bf r}$. We assume translational invariance so that $A_{{\bf r} {\bf r}'} = A({\bf r} - {\bf r}')$, where the matrix $A({\bf r} - {\bf r}')$ represents the hopping matrix between sites at ${\bf r}$ and ${\bf r}'$.  The mean-field Hamiltonian can be diagonalized by the complex fermions in Fourier space as
\begin{eqnarray}
f^\alpha({\bf r}) = \sqrt{\frac{1}{N_{s}}} \sum_{{\bf k}\in BZ} e^{i {\bf k}\cdot {\bf r}}\mathfrak{f}^{\alpha}({\bf k}),
\end{eqnarray}
where $N_{s}$ is the number of sites, and the complex fermion field $f$ satisfies the usual anti-commutation relations, $\{ \mathfrak{f}^{\alpha \dagger} ({\bf k}),\mathfrak{f}^{\alpha'}({\bf k}') \} = \delta^{\alpha \alpha'} \delta_{{\bf k} {\bf k}'}$. The component of the spin-$1/2$ spin operator is
\begin{eqnarray}
S^\mu({\bf r}) = \sum_{\alpha, \beta=1,2} f^{\alpha\dagger}({\bf r}) \left( \frac{\sigma^\mu_{\alpha \beta}}{2} \right) f^{\beta}({\bf r}),
\end{eqnarray}
which can be expressed in Fourier space as
\begin{eqnarray}
&& S^{\mu}({\bf r})\simeq \sum_{{\bf k},{\bf k'}\in {\bf B.Z.}}\sum_{\alpha,\beta}  \frac{\sigma^\mu_{\alpha \beta}}{2N_s} \mathfrak{f}^{\alpha \dagger}({\bf k}) \mathfrak{f}^\beta ({\bf k'}) e^{-i({\bf k}-{\bf k'})\cdot {\bf r}}.
\end{eqnarray}
In order to determine the long-distance behavior at separation ${\bm r}$, we focus on patches near the Fermi surface where the group velocity is parallel or antiparallel to the observation direction $\hat{\bm n} = {\bm r}/|{\bm r}|$, because at large separation $|{\bm r}| \gg k_F^{-1}$, the main contributions to the correlations arise precisely from such patches.  Specifically, we introduce Right(R) and Left(L) Fermi patch fields and the corresponding energies
\begin{eqnarray}\label{fermi patches:momentum space}
&& f^{\alpha,(\hat{\bm n})}_P(\delta {\bm k}) = \mathfrak{f}^{\alpha}({\bm k}^{(\hat{\bm n})}_{FP} + \delta {\bm k}) ~,\\
&& \epsilon^{(\hat{\bm n})}_P (\delta {\bm k}) = |{\bm v}^{(\hat{\bm n})}_{FP}| \left( P \delta k_\parallel + \frac{\mathfrak{C}^{(\hat{\bm n})}_P}{2} \delta k_\perp^2 \right) ~, \label{fermi patches:energy}
\end{eqnarray}
where the superscript $(\hat{\bm n})$ refers to the observation direction and $P = R/L = +/-$; ${\bm v}^{(\hat{\bm n})}_{FP}$ is the corresponding group velocity (parallel to $\hat{\bm n}$ for the Right patch and anti-parallel for the Left patch); $\mathfrak{C}_{P = R/L}$ is the curvature of the Fermi surface at the Right/Left patch; $\delta k_\parallel$ and $\delta k_\perp$ are respectively components of $\delta {\bm k}$ parallel and perpendicular to $\hat{\bm n}$.  It is convenient to define fields in real space
\begin{eqnarray}\label{fermi patches:real space}
f^{\alpha,(\hat{\bm n})}_P({\bm r}) \sim \sum_{\delta {\bm k} \in {\rm Fermi~Patch}} f^{\alpha,(\hat{\bm n})}_P(\delta {\bm k}) e^{i\delta{\bm k} \cdot {\bm r}} ~,
\end{eqnarray}
which vary slowly on the scale of the lattice spacing [from now on, we will drop the superscript $(\hat{\bm n})$].  In this long-wavelength analysis, the relevant terms in the spin operator are
\begin{eqnarray}\label{Eq:spincorr_longwl}
 S^{\mu}({\bf r})\sim \sum_{P,P'}\sum_{\alpha, \beta} \sigma^{\mu}_{\alpha \beta} f^{\alpha\dagger}_{P}({\bf r}) f^{\beta}_{ P'}({\bf r}) e^{-i({\bf k}_{FP}-{\bf k}_{FP'})\cdot {\bf r}}.
\end{eqnarray}
The above long-wavelength expression for the $S^{\mu}$ operator implies that the corresponding correlation function defined in Eq.~\eqref{Eq:spincorr} contains contributions with ${\bf q}={\bf 0}$ and  ${\bf q}_{-}\equiv\bm{k}_{FR}- \bm{k}_{FL}$. More explicitly, for a patch specified by $\epsilon_P (\delta {\bm k})$ in Eqs.~\eqref{fermi patches:momentum space}-\eqref{fermi patches:energy}, we can derive the Green's function for the continuum complex fermion fields as
\begin{eqnarray}\label{long-wavelength:green function}
\langle f^{\alpha\dagger}_{R/L}({\bm 0}) f^{\alpha}_{R/L}({\bm r}) \rangle = \frac{\exp[\mp i \frac{3\pi}{4}]}{2^{3/2} \pi^{3/2} \mathfrak{C}_{R/L}^{1/2} |{\bm r}|^{3/2}} ~.
\end{eqnarray}
Using this and Eq.~\eqref{Eq:spincorr_longwl}, we can obtain the spin correlation
\begin{eqnarray}
\chi_s({\bf r}) &\sim& -\frac{1}{\mathfrak{C}_{R} |{\bf r}|^3} - \frac{1}{\mathfrak{C}_{L}|{\bf r}|^3}~\label{uniform:spin}\\
&+& \frac{2 \sin [ ( \bm{k}_{FR}-\bm{k}_{FL})\cdot {\bf r}] }{\mathfrak{C}_{R}^{1/2}\mathfrak{C}_{L}^{1/2} |{\bf r}|^3}. ~\label{krminuskl:spin}
\end{eqnarray}
Focusing on the structure factors $\mathcal{D}_{\bf q}$ defined in Eq.~\eqref{Eq:Sq}, we expect that there should be a cone-shaped singularity at ${\bf q}=0$, based on Eq.~\eqref{uniform:spin}:
\begin{eqnarray}
\mathcal{D}_{{\bf q}\sim 0} \sim |{\bf q}|,
\end{eqnarray}
which can be seen straightforwardly by performing Fourier transform exactly or by scaling analysis with ${\bf q} \sim {\bf r}^{-1}$. Furthermore, the spin structure factor should also reveal the singular surface at ${\bf Q}_{-} $, as expected from Eq.~\eqref{krminuskl:spin}.  At the long wavelength analysis at the mean-field level, we note that the singularities are expected to be one-sided,
\begin{eqnarray}
&& \mathcal{D}_{{\bf Q}_{-} + \delta {\bf q}}\sim |\delta q_{||}|^{3/2} \Theta (-\delta q_{||}).\end{eqnarray}

Fitting the exact $\bm{q} = \bm{k}_{FR} - \bm{k}_{FL}$ in the spin structure factor data illustrated in Fig.~2(b) in the main text, we extract the cross sections of the emergent SFS in the isotropic case to be the $5.24 \pm 0.05$ (where we set the lattice constant $a$ to be 1 and $\hbar = 1$), which allows us to obtain $c_{eff}\simeq2.01 \pm 0.02$.

\section{Spin Structure Factor of the Anisotropic Case}
The Fig.~\ref{Fig:sq_ani07} shows the 3D spin structure factor obtained on a triangular lattice with $30 \times 30$ sites for the anisotropic case $t'=0.7t$. The spin structure factor is defined as $D_{\bf q} \equiv \sum_j \chi^s_j e^{-i {\bf q} \cdot {\bf r}_j}$ with the real-space spin correlation function $\chi^s_j \equiv \sum_{\mu =x, y, z} \langle S^\mu_j S^\mu_{i \equiv0}\rangle$. It is known that for an observation direction $\hat{n}$, $D_{\bm{q}}$ should show singular peaks at $\bm{q} = \bm{0},~\bm{k}^{\hat{n}}_{FR} - \bm{k}^{\hat{n}}_{FL}$ associated with forward and backward scattering process. The Fig.~\ref{Fig:sq_ani07}(a) gives a side view of the $D_{\bm{q}}$ in the hexagonal Brillouin zone (B.Z.) where we observe a sharp singular point at $\bm{q} = 0$ and the wiggle lines on the surface correspond to the singular lines located at $\bm{k}^n_{FR} - \bm{k}^n_{FL}$. Unlike the isotropic case, we can see that the structure factor breaks the $C_6$ rotation. The Fig.~\ref{Fig:sq_ani07}(b) shows the top view of $D_{\bm{q}}$. To determine the location of the $\bm{k}^n_{FR} - \bm{k}^n_{FL}$, we adopt a mean-field fermionic state with a Fermi surface at $1/2$-filling. Extracting $\bm{k}^n_{m,FR} - \bm{k}^n_{m,FL}$ of the mean-field ansatz, where the subscript $m$ means mean-field, we find that the $\bm{k}^n_{m,FR} - \bm{k}^n_{m, FL}$ can fit the exact $\bm{k}^n_{FR} - \bm{k}^n_{FL}$ quite well for an arbitrary observation direction $\hat{n}$, which gives the locations of the weak singular lines on the surface of the 3D spin structure factor. In Fig.~\ref{Fig:sq_ani07}(b) the blue lines are obtained by examining the mean-field SFS spin liquid state, which overlap the weak singular lines on the surface of the 3D spin structure factor obtained by exact numerical calculations.

