\section{Supplementary Materials}

\subsection{Momentum-space tight-binding characterization}

The Bloch states of the hexagonal lattice $V(\mathbf{r})=-V_{0}|%
\sum_{j=1}^{3}\mathbf{\epsilon }_{j}\exp [i\mathbf{k}_{L,j}\cdot (\mathbf{r}-%
\mathbf{r_{0}})]|^{2}$ can be written as $\phi _{l\mathbf{k}}(\mathbf{r}%
)=e^{i\mathbf{k}\cdot \mathbf{r}}u_{l\mathbf{k}}(\mathbf{r})$. The periodic
part can be expanded as $u_{l\mathbf{k}}(\mathbf{r})=\sum_{\mathbf{p}}c_{l%
\mathbf{k}}^{\mathbf{p}}e^{i\mathbf{p}\cdot \mathbf{r}}$, where $\mathbf{p}%
=p_{1}\mathbf{e}_{1}+p_{2}\mathbf{e}_{2}$ with $\mathbf{e}_{i}$ the
primitive reciprocal lattice vectors and $p_{i}$ integers. By substituting
expansion of $\phi _{l\mathbf{k}}(\mathbf{r})$ into the Schr\"{o}dinger
equation $[-\frac{\mathbf{\nabla }^{2}}{2m}+V(\mathbf{r})]\phi _{l\mathbf{k}%
}=\mathcal{E}_{l\mathbf{k}}\phi _{l\mathbf{k}}$, the Bloch states (i.e., the
coefficients $c_{l\mathbf{k}}^{\mathbf{p}}$) and bands can be obtained. The
bare Bloch bands $\mathcal{E}_{sl\mathbf{k}_{s}}$ and states $\phi _{sl%
\mathbf{k}_{s}}$ of the two pseudospin states can be obtained similarly, as
shown in Fig.~\ref{fig:S1}. Here we keep the expansion coefficients up to $%
p_{i}=\pm 6$ in the calculation. We see that the two lowest bands (i.e., $%
\mathcal{E}_{sl\mathbf{k}_{s}}$ with $l=1,2$), which have a gap from higher
bands, form two Dirac cones at valleys $K_{s}$ and $K_{s}^{\prime }$~\cite%
{science.aad5812S}. For typical lattice depths, the two Dirac bands are
asymmetric with respect to Dirac-point energy $E=0$ due to the long-range
tunnelings beyond nearest neighbors that break the sublattice symmetry. Such
effect is more significant for shallower lattices.

\begin{figure}[b]
\includegraphics[width=0.85\linewidth]{FigS1.pdf}
\caption{Band structures of a hexagonal lattice with (a) $V_{0}=6E_{R}$, (b)
$V_{0}=4E_{R}$ and (c) $V_{0}=2E_{R}$. We have shifted the Dirac-point
energy to $E=0$. For a shallower lattice, the tunnelings (both short and
long range tunnelings) are stronger, therefore, the Dirac velocity is larger
and the nonlinearity of the dispersion around the Dirac points is stronger.
In addition, the two Dirac bands become more asymmetric with respect to the
Dirac-point energy $E=0$ due to the enhanced long-range tunnelings that
break the sublattice symmetry.}
\label{fig:S1}
\end{figure}

We are interested in the low-energy physics near the Dirac points, therefore
we only retain the two-lowest bare Dirac bands (i.e., $\mathcal{E}_{sl%
\mathbf{k}_{s}}$ with $l=1,2$) and drop all higher bands with $l>2$ safely.
That is, we keep only one Wannier orbital at each hexagonal lattice site,
and these Wannier orbitals form the full tight-binding basis set. We find
that it is more convenient to work in the Bloch basis $\{\phi _{sl\mathbf{k}%
_{s}}\}$, which is equivalent to the Wannier basis up to a Fourier
transformation. For commensurate twists $\cos (\theta )=\frac{n^{2}+m^{2}+4mn%
}{2(n^{2}+m^{2}+mn)}$, one can fold the bare Brillouin zone of each spin
states to the Moir\'{e} Brillouin zone. The inter-spin coupling coefficient
can be obtained as $\langle \phi _{\uparrow l\mathbf{k}_{\uparrow }}|\Omega
|\phi _{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}\rangle =\delta _{%
\mathbf{q},\mathbf{q}^{\prime }}J_{\mathbf{g}_{\uparrow }\mathbf{g}%
_{\downarrow }}^{ll^{\prime }}(\mathbf{q})$, where $\mathbf{k}_{\uparrow }=%
\mathbf{q}+\mathbf{g}_{\uparrow }$ and $\mathbf{k}_{\downarrow }=\mathbf{%
q^{\prime }}+\mathbf{g}_{\downarrow }$, with $\mathbf{g}_{s}$ the Moir\'{e}
reciprocal lattice vectors and $\mathbf{q}$ the Moir\'{e} Bloch momentum
that is a good quantum number. We obtain the Hamiltonian Eq.~(1) in the main
text, which incorporates all range real-space tunnelings with high accuracy.
Another advantage of this momentum-space approach is that we only need to
keep $\mathbf{g}_{s}$ around the Dirac cones to correctly characterize the
low-energy physics at small twist angles, leading to rather rapid
convergence of the basis set.

Now we show how to evaluate the inter-spin coupling coefficients
\begin{eqnarray}
\langle \phi _{\uparrow l\mathbf{k}_{\uparrow }}|\Omega |\phi _{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle &=&\Omega \int d^{2}\mathbf{r}%
e^{-i\mathbf{k}_{\uparrow }\cdot \mathbf{r}}u_{\uparrow l\mathbf{k}%
_{\uparrow }}^{\ast }(\mathbf{r})e^{i\mathbf{k}_{\downarrow }\cdot \mathbf{r}%
}u_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}(\mathbf{r})  \nonumber \\
&=&\Omega \sum_{\mathbf{M}}e^{i(\mathbf{q^{\prime }}-\mathbf{q})\cdot
\mathbf{M}}\langle u_{\uparrow l\mathbf{k}_{\uparrow }}|e^{i(\mathbf{k}%
_{\downarrow }-\mathbf{k}_{\uparrow })\cdot (\mathbf{r-M})}|u_{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle _{\mathbf{M}}  \nonumber \\
&=&\Omega \delta _{\mathbf{q},\mathbf{q}^{\prime }}\langle u_{\uparrow l%
\mathbf{k}_{\uparrow }}|e^{i(\mathbf{g}_{\downarrow }-\mathbf{g}_{\uparrow
})\cdot \mathbf{r}}|u_{\downarrow l^{\prime }\mathbf{k}_{\downarrow
}}\rangle .
\end{eqnarray}%
Here $\mathbf{M}$ denotes the Moir\'{e} lattice vectors, and the term $%
\langle \cdots \rangle _{\mathbf{M}}$ in the second line (with integral over
the $\mathbf{M}$-th Moir\'{e} unit cell) is independent of $\mathbf{M}$. The
coefficients $J_{\mathbf{g}_{\uparrow }\mathbf{g}_{\downarrow }}^{ll^{\prime
}}(\mathbf{q})$ is
\begin{eqnarray}
J_{\mathbf{g}_{\uparrow }\mathbf{g}_{\downarrow }}^{ll^{\prime }}(\mathbf{q}%
) &=&\Omega \langle u_{\uparrow l\mathbf{k}_{\uparrow }}|e^{i(\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{r}}|u_{\downarrow
l^{\prime }\mathbf{k}_{\downarrow }}\rangle  \nonumber \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow }}c_{\uparrow l%
\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast }c_{\downarrow l^{\prime
}\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow }}\int d^{2}\mathbf{r}%
e^{i(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow
}-\mathbf{g}_{\uparrow })\cdot \mathbf{r}}  \nonumber \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow },\mathbf{M}%
}c_{\uparrow l\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast
}c_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow
}}\frac{e^{i(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{L}_{1}}-1}{(\mathbf{p}%
_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow }-\mathbf{g}%
_{\uparrow })\cdot \mathbf{L}_{1}}\frac{1-e^{i(\mathbf{p}_{\downarrow }-%
\mathbf{p}_{\uparrow }+\mathbf{g}_{\downarrow }-\mathbf{g}_{\uparrow })\cdot
\mathbf{L}_{2}}}{(\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow }+\mathbf{g}%
_{\downarrow }-\mathbf{g}_{\uparrow })\cdot \mathbf{L}_{2}}  \label{eq:J} \\
&=&\Omega \sum_{\mathbf{p}_{\uparrow },\mathbf{p}_{\downarrow },\mathbf{M}%
}c_{\uparrow l\mathbf{k}_{\uparrow }}^{\mathbf{p}_{\uparrow }\ast
}c_{\downarrow l^{\prime }\mathbf{k}_{\downarrow }}^{\mathbf{p}_{\downarrow
}}\delta _{\mathbf{p}_{\downarrow }-\mathbf{p}_{\uparrow },\mathbf{g}%
_{\uparrow }-\mathbf{g}_{\downarrow }}.  \nonumber
\end{eqnarray}%
To obtain the last step, we have used $\mathbf{p}_{s}\cdot \mathbf{L}%
_{1}=2\pi (mp_{s1}+np_{s2})$ (similarly for $\mathbf{L}_{2}$), since $%
\mathbf{p}_{s}=p_{s1}\mathbf{e}_{s\mathbf{k}_{1}}+p_{s2}\mathbf{e}_{s\mathbf{%
k}_{2}}$ and $\mathbf{L}_{1}=m\mathbf{e}_{s1}+n\mathbf{e}_{s2}$, with $%
\mathbf{e}_{s1},\mathbf{e}_{s2}$ ($\mathbf{e}_{s\mathbf{k}_{1}},\mathbf{e}_{s%
\mathbf{k}_{2}}$) the primitive (reciprocal) lattice vectors of trap $V_{s}$.

\begin{figure}[b]
\includegraphics[width=1.0\linewidth]{FigS2.pdf}
\caption{(a) and (b) Typical density distributions of the Moire flat band
Bloch states $\protect\psi _{j\mathbf{q}}(\mathbf{r},\downarrow )$ in one
Moir\'{e} unit cell for $\mathbf{q}=K$ and $\mathbf{q}=\Gamma $,
respectively. Atoms are mainly distributed around the $AA$ region for $%
\mathbf{q}=K$ and become more uniformly distributed for $\mathbf{q}=\Gamma $%
. Other parameters are $\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$ and $%
\Omega =0.116$. (c) Moir\'{e} bands at incommensurate twist angle $\protect%
\theta =4^{\circ }$ based on the generalized one-valley continuum model,
with $V_{0}=6E_{R}$ and $\Omega =0.0915$. The red and blue lines are bands
from $K_{s}$ and $K_{s}^{\prime }$ valleys, respectively. We have retained a
shell (containing 37 MBZs up to the third Moir\'{e} reciprocal lattice
vector shell) around the valley to construct the basis set. (d) Moir\'{e}
bands at commensurate twist $\protect\theta =5.086^{\circ }$ based on the
generalized one-valley continuum model, with $V_{0}=6E_{R}$ and $\Omega
=0.116$. The green thin dashed lines are Moir\'{e} bands obtained using the
full tight-binding basis.}
\label{fig:S2}
\end{figure}

With these coupling coefficients, we can diagonalize the single particle
Hamiltonian $H_{0}=H_{\uparrow }+H_{\downarrow }+H_{\uparrow \downarrow }$.
In the Moir\'{e} Bloch eigenbasis $\{\psi _{j\mathbf{q}}\}$, it reads
\begin{equation}
H_{0}=\sum_{j,\mathbf{q}}E_{j\mathbf{q}}\beta _{j\mathbf{q}}^{\dag }\beta _{j%
\mathbf{q}}.
\end{equation}%
The typical distributions of $\psi _{j\mathbf{q}}(\mathbf{r},s)$ are shown
in Figs.~\ref{fig:S2}a and \ref{fig:S2}b. Atoms are mainly distributed
around the $AA$ region for all $\mathbf{q}$ except a small area near $%
\mathbf{q}=\Gamma $, where atoms become more uniformly distributed. As a
result, the interaction and thereby the pairing is weak at $\Gamma $ point.

The above results and the Hamiltonian Eq.~(1) in the main text apply for any
commensurate twist angles. For incommensurate twist angles, there are no
well-defined Moir\'{e} patterns and Moir\'{e} bands if the twist angle is
too large. However, if the twist angle is small enough, Moir\'{e} patterns
can form even for incommensurate twist. In this case, the low-energy physics
is mainly determined by the states around the Dirac valleys, and the
inter-valley coupling is also negligible. Therefore, we can adopt a similar
approach as the continuum model by only keeping $\mathbf{g}_{s}$ around one
valley in the Hamiltonian, and the inter-spin couplings can be obtained
using Eq.~\ref{eq:J}. In Fig.~\ref{fig:S2}c, we plot the Moir\'{e} bands at
a small incommensurate twist angle using the generalized one-valley
continuum model mentioned above. We also plot the Moir\'{e} bands at a small
commensurate twist angle using the one-valley continuum model approach and
compare it with the results based on the full tight-binding basis (see Fig.~%
\ref{fig:S2}d). Their agreement confirms the validation of the one-valley
continuum model at small twist angles.

As we discussed in the main text, the energy scales of the magic physics is
small and at the order of $10^{-3}E_{R}$. Therefore, in the Wannier basis,
long-range (intra- and inter-spin) tunnelings beyond nearest neighbors that
are $\gtrsim 10^{-3}E_{R}$ should be taken into account to obtain correct
magic flat bands~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S,PNAS.108.12233S,PhysRevB.86.155449S}%
. In particular, the inter-spin site separations take various values and are
nearly continuously distributed for small twists, therefore the on-site, NN,
next-NN, etc., inter-spin tunnelings should be considered even when the
tunneling cutoff distance is small. The inter-spin tunnelings are also
highly anisotropic, which depend on both the relative orientation and
distance between the corresponding sites due to the three-fold rotational
symmetry of the Wannier orbitals. A small deviation in the tunnelling
coefficients may result in significant change in the band structures near
the `magic angle' due to the narrow bandwidths and approximate degeneracy of
the flat bands. The simplified real-space tight-binding model in~\cite%
{PhysRevA.100.053604S}, which assumes a simply isotropic Gaussian Wannier
function and includes only the nearest-neighbor intra-spin tunneling and
on-site inter-spin couplings, has significant deviation in determining the
bands at small twist angles.

To elaborate more on the effects of long-range real-space tunnelings, we
first solve for the maximally localized Wannier functions $w_{s,j}(\mathbf{r}%
)$ (with $j$ the site index) of the two spin-dependent lattices, then
calculate the band structures using the real-space tight-binding Hamiltonian
\[
H_{s}=\sum_{j,j^{\prime }}t_{jj^{\prime }}\alpha _{s,j}^{\dag }\alpha
_{s,j^{\prime }}
\]%
and
\[
H_{\uparrow \downarrow }=\sum_{j,j^{\prime }}J_{jj^{\prime }}\alpha
_{\uparrow ,j}^{\dag }\alpha _{\downarrow ,j^{\prime }}+h.c.,
\]%
where $\alpha _{s,j}^{\dag }$ is the creation operator of the Wannier
orbitals, and the tunneling coefficients are $t_{jj^{\prime }}=\int d^{2}%
\mathbf{r}w_{s,j}^{\ast }(\mathbf{r})[-\nabla _{\mathbf{r}}^{2}+V_{s}(%
\mathbf{r})]w_{s,j^{\prime }}(\mathbf{r})$ and $J_{jj^{\prime }}=\int d^{2}%
\mathbf{r}w_{\uparrow ,j}^{\ast }(\mathbf{r})\Omega w_{\downarrow ,j^{\prime
}}(\mathbf{r})$. For the lattice depth considered in this paper, we have $%
(t_{0};t_{1};t_{2};t_{3};t_{4};t_{5})=(-4.3458;-0.1967;0.0166;-0.0047;-0.0025;0.0009)E_{R}
$ for $V_{0}=6E_{R}$ and $%
(t_{0};t_{1};t_{2};t_{3};t_{4};t_{5})=(-2.2498;-0.2351;0.0245;-0.0089;-0.0045;0.002)E_{R}
$ for $V_{0}=4E_{R}$, with $t_{jj^{\prime }}=t_{\tilde{\imath}}$ if site $j$
is the $\tilde{\imath}$-th neighbor of site $j^{\prime }$ (where $t_{0}$ is
the on-site energy). Therefore, the cutoff distance for tunneling should be
up to the 5-th nearest neighbors of the bare hexagonal lattice. The
inter-spin coupling $J_{jj^{\prime }}$ is nearly continuously distributed
which decays exponentially with site separation (the on-site tunneling rate
is close to $\Omega $). In Figs.~\ref{figS:real_space_TB}a and ~\ref%
{figS:real_space_TB}b, we show the band structures for different cutoff
distance. We see that the flat band energy is largely shifted (by $\sim
0.25t_{1}$) for a cutoff only to the 1st nearest neighbors of the bare
hexagonal lattice. Nevertheless, magic flat bands are still observed with a
larger width (which can be reduced further by tuning $\Omega $), indicating
the robustness of the magic-angle physics.

\begin{figure}[b]
\includegraphics[width=0.95\linewidth]{FigS_real_space_TB_r.pdf}
\caption{(a) and (b) Moir\'{e} bands based on real-space tight-binding
Hamiltonian for $V_0=4E_R$ and $V_0=6E_R$, respectively. For cutoff distance
$I_c$, we keep the (intra- and inter-spin) tunnelings with site separation
not greater than the $I_c$-th neighbor separation of the bare hexagonal
lattice (e.g., we keep $t_{\tilde{\imath}}$ with $\tilde{\imath}\leq I_c$). The blue,
green and red solid lines correspond to the cutoff distance up to the $1$st,
$2$nd and $5$-th neighbors of the bare hexagonal lattice. The thin dashed
lines are the results of the momentum-space tight-binding Hamiltonian. Other
parameters are $\protect\theta =5.086^{\circ }$ and $\Omega =\Omega_f$. (c)
Schematic illustration of the pairing (indicated by dashed lines) in the
Brillouin zones, between $\uparrow$ states (green dots at $+\mathbf{k}$)
around valley $K^{\prime }_\uparrow$ and $\downarrow$ states (red dots at $-%
\mathbf{k}$) around valley $K_\downarrow$. The two pairing states are
located at the same side of the corresponding valleys (as indicated by the
black arrows), and thus the relative phase between $A$ and $B$ sublattice
sites (which is related to the chirality of the valley) leads to the
staggered pairing order.}
\label{figS:real_space_TB}
\end{figure}

\subsection{Superfluid orders}

The mean-filed interaction reads
\begin{eqnarray}
H_{\text{int}} &=&U_{0}\int d^{2}\mathbf{r}\hat{\Psi}_{\uparrow }^{\dag }(%
\mathbf{r})\hat{\Psi}_{\downarrow }^{\dag }(\mathbf{r})\hat{\Psi}%
_{\downarrow }(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})  \nonumber \\
&=&\int d^{2}\mathbf{r}[\hat{\Psi}_{\uparrow }^{\dag }(\mathbf{r})\hat{\Psi}%
_{\downarrow }^{\dag }(\mathbf{r})\Delta (\mathbf{r})+\hat{\Psi}_{\downarrow
}(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})\Delta ^{\ast }(\mathbf{r})-%
\frac{|\Delta (\mathbf{r})|^{2}}{U_{0}}]
\end{eqnarray}%
with local pairing amplitude $\Delta (\mathbf{r})=U_{0}\langle \hat{\Psi}%
_{\downarrow }(\mathbf{r})\hat{\Psi}_{\uparrow }(\mathbf{r})\rangle =\sum_{%
\mathbf{g}}\Delta _{\mathbf{g}}e^{i\mathbf{g}\cdot \mathbf{r}}$. We expand
the field operator in the Moir\'{e} Bloch basis $\hat{\Psi}_{s}(\mathbf{r}%
)=\sum_{j,\mathbf{q}}\beta _{j\mathbf{q}}\psi _{j\mathbf{q}}(\mathbf{r},s)$
and obtain the gap equation
\begin{equation}
\Delta _{\mathbf{g}}=U_{0}\sum_{j^{\prime },j,\mathbf{q}}\chi _{j^{\prime
}j}^{\mathbf{q}}(\mathbf{g})C_{\mathbf{q}}^{j^{\prime }j}[\Delta (\mathbf{r}%
)],  \label{eq:GapS}
\end{equation}%
where $\chi _{j^{\prime }j}^{\mathbf{q}}(\mathbf{g})=\frac{1}{L}\int d^{2}%
\mathbf{r}e^{-i\mathbf{g}\cdot \mathbf{r}}\psi _{j^{\prime }-\mathbf{q}}(%
\mathbf{r},\downarrow )\psi _{j\mathbf{q}}(\mathbf{r},\uparrow )$. The
correlation $C_{\mathbf{q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }-%
\mathbf{q}}\beta _{j\mathbf{q}}\rangle $ can be obtained by solving the
Bogoliubov-de Gennes (BdG) Hamiltonian%
\begin{equation}
H_{\text{BdG}}=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}}-\mu)\beta _{j\mathbf{q}}^{\dag
}\beta _{j\mathbf{q}}+\sum_{j,j^{\prime },\mathbf{q}}[\bar{\Delta}%
_{jj^{\prime }}(\mathbf{q})\beta _{j\mathbf{q}}^{\dag }\beta _{j^{\prime }-%
\mathbf{q}}^{\dag }+h.c.]
\end{equation}%
with $\bar{\Delta}_{jj^{\prime }}(\mathbf{q})=L\sum_{\mathbf{g}}\Delta _{%
\mathbf{g}}\chi _{j^{\prime }j}^{\mathbf{q}\ast }(\mathbf{g})$. We solve the
gap equation self-consistently by retaining only the four flat bands which
have much larger DOS than nearby bands.

As we discussed in the main text, at finite temperature, the relevant
physics is the Berezinskii-Kosterlitz-Thouless (BKT) transition since there
is no long-range superfluid order in two dimensions due to phase
fluctuations, and free vortex-antivortex pairs are formed spontaneously
below a critical temperature $T_{\text{BKT}}$~\cite{PhysRevLett.109.105302S,PhysRevLett.114.110401S}. The critical temperature $%
T_{c}$ by mean-field theory is often overestimated. Nevertheless, the
mean-field critical temperature $T_{c}$ provides an upper bound of $T_{\text{%
BKT}}$. Furthermore, the true BKT transition temperature could be obtained
from the mean-field superfluid weight $\rho _{s}$ which can be calculated
based on current response~\cite{PhysRevLett.112.086401S,julku2018superfluidS}%
. On the Moir\'{e} superlattice scale, our system corresponds to an $s$-wave
superfluid, we introduce a slowly varying superfluid phase $\Delta (\mathbf{r%
})\rightarrow \Delta (\mathbf{r})e^{i\Theta (\mathbf{r},\tau)}=\sum_{\mathbf{g}%
}\Delta _{\mathbf{g}}e^{i\mathbf{g}\cdot \mathbf{r}}e^{i\Theta (\mathbf{r}%
,\tau)}$, and the phase $\Theta (\mathbf{r},\tau)=\mathbf{Q}\cdot \mathbf{r}$
(with $\mathbf{Q}\rightarrow 0$) leads to a supercurrent $\mathbf{j}%
_{s}=\rho _{s}\mathbf{Q}$. The supercurrent can also be calculated using the
variation of the free energy $\mathbf{j}_{s}=\frac{1}{L}\langle \frac{\delta
H_{\text{BdG}}}{\delta \mathbf{Q}}\rangle $, the superfluid weight is
obtained by $\rho _{s}=\mathbf{j}_{s}/\mathbf{Q}$, and the BKT transition
temperature is given by $T_{\text{BKT}}=\frac{\pi }{2}\rho _{s}(T_{\text{BKT}%
})$~\cite{PhysRevLett.112.086401S,julku2018superfluidS}. Here the BdG
Hamiltonian reads
\begin{equation}
H_{\text{BdG}}(\mu ,\Delta ,\mathbf{Q})=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}%
}-\mu)\beta _{j\mathbf{q}}^{\dag }\beta _{j\mathbf{q}}+\sum_{j,j^{\prime },%
\mathbf{q}}[\bar{\Delta}_{jj^{\prime }}(\mathbf{q,Q})\beta _{j\mathbf{q+Q/2}%
}^{\dag }\beta _{j^{\prime }\mathbf{-q+Q/2}}^{\dag }+h.c.]
\end{equation}%
where $\bar{\Delta}_{jj^{\prime }}(\mathbf{q,Q})=L\sum_{\mathbf{g}}\Delta _{%
\mathbf{g}}\chi _{j^{\prime }j}^{\mathbf{q,Q}\ast }(\mathbf{g})$ and $\chi
_{j^{\prime }j}^{\mathbf{q,Q}}(\mathbf{g})=\frac{1}{L}\int d^{2}\mathbf{r}%
e^{-i\mathbf{g}\cdot \mathbf{r}}\psi _{j^{\prime }\mathbf{-q+Q/2}}(\mathbf{r}%
,\downarrow )\psi _{j\mathbf{q+Q/2}}(\mathbf{r},\uparrow )$. The gap
equation is $\Delta _{\mathbf{g}}=U_{0}\sum_{j^{\prime },j,\mathbf{q}}\chi
_{j^{\prime }j}^{\mathbf{q,Q}}(\mathbf{g})C_{\mathbf{q,Q}}^{j^{\prime }j}$
with $C_{\mathbf{q,Q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }\mathbf{%
-q+Q/2}}\beta _{j\mathbf{q+Q/2}}\rangle $.

Typically, $\rho _{s}(T)$ is
nearly a constant at low temperature $T$, and therefore we use the following
approximation $T_{\text{BKT}}\simeq \frac{\pi }{2}\rho _{s}(T=0)$. We
numerically calculate $T_{\text{BKT}}$ (the results are shown in Fig.~3a in
the main text) and find that $T_{\text{BKT}}\simeq 0.4T_{c}$. In particular,
we first set $\Theta (\mathbf{r},\tau)=0$ and solve for the order parameter $%
\Delta _{0}(\mathbf{r})$ with fixed chemical potential $\mu $, the
corresponding zero-temperature free energy is $\mathcal{F}(\mathbf{Q}%
=0)\equiv \frac{1}{L}\langle H_{\text{BdG}}(\mu ,\Delta _{0},\mathbf{Q}%
=0)\rangle $. Then we introduce the superfluid phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$ by setting the order parameter as $\Delta _{0}(%
\mathbf{r})e^{i\mathbf{Q}\cdot \mathbf{r}}$, and we calculate the free
energy $\mathcal{F}(\mathbf{Q})\equiv \frac{1}{L}\langle H_{\text{BdG}}(\mu
,\Delta _{0},\mathbf{Q})\rangle $. Since $\mathbf{j}_{s}=\rho _{s}\mathbf{Q}=%
\frac{\delta F(\mathbf{Q})}{\delta \mathbf{Q}}$ for $\mathbf{Q}\rightarrow 0$
at $T=0$, we have $F(\mathbf{Q})-F(0)=\frac{1}{2}\rho _{s}\mathbf{Q}^{2}=%
\frac{1}{2}\rho _{s}(\nabla \Theta )^{2}$ and the superfluid weight is
calculated as $\rho _{s}(T=0)=2\times \frac{F(\mathbf{Q})-F(0)}{\mathbf{Q}%
^{2}}$. The BKT critical temperature is obtained as $T_{\text{BKT}}\simeq
\pi \times \frac{F(\mathbf{Q})-F(0)}{\mathbf{Q}^{2}}$.

It can be verified that, the results based on current response are
consistent with those obtained through the effective action based on quantum
field theory. The partition function can be written as $Z=\text{Tr}%
(e^{-H/T})=\int \mathcal{D}[\bar{\Psi},\Psi ]e^{-S_{\text{eff}}[\bar{\Psi}%
,\Psi ]}$. The effective action is
\begin{equation}
S_{\text{eff}}[\bar{\Psi},\Psi ]=\int_{0}^{1/T}d\tau \left[ \int d^{2}%
\mathbf{r}\sum_{s}\bar{\Psi}_{s}\partial _{\tau }\Psi _{s}+\mathcal{H}(\bar{%
\Psi},\Psi )\right] ,
\end{equation}%
with $\mathcal{H}(\bar{\Psi},\Psi )=\int d^{2}\mathbf{r}\left\{ \sum_{s}[%
\bar{\Psi}_{s}(-\nabla ^{2}+V_{s}-\mu)\Psi _{s}+\bar{\Psi}_{s}\Omega \Psi _{\bar{%
s}}]+U_{0}\bar{\Psi}_{\uparrow }\bar{\Psi}_{\downarrow }\Psi _{\downarrow
}\Psi _{\uparrow }\right\} $ obtained by replacing field operator $\hat{\Psi}%
_{s}^{\dag }(\mathbf{r},t)$ and $\hat{\Psi}_{s}(\mathbf{r},t)$ with Grassman
field number $\bar{\Psi}_{s}$ and $\Psi _{s}$ in the Hamiltonian, and $\tau $
is the imaginary time. We can integrate out the quartic interaction term by
Hubbard-Stratonovich transformation, where the order parameter $\Delta $ is
defined. By further integrating out fermion fields, the partition function
reads $Z=\int \mathcal{D}[\Delta ^{\ast },\Delta ]e^{-S_{\text{eff}}[\Delta
^{\ast },\Delta ]}$ with
\begin{equation}
S_{\text{eff}}[\Delta ^{\ast },\Delta ]=\int_{0}^{1/T}d\tau \int d^{2}%
\mathbf{r}\frac{-|\Delta |^{2}}{U_{0}}-\frac{1}{2}\text{Tr ln}\left( i\omega
_{m}+M_{\text{BdG}}\right) +\frac{1}{T}\sum(\mathcal{E}_{sl\mathbf{k}_{s}}-\mu ).
\end{equation}%
Where $M_{\text{BdG}}$ is the BdG matrix in the Nambu basis, and $\omega _{m}
$ are Matsubara frequencies. For a static phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$, the effective action induced by the phase
fluctuation is
\begin{equation}
S_{\text{fluc}}=\frac{1}{2}\text{Tr ln}[i\omega _{m}+M_{\text{BdG}%
}(0)]-\frac{1}{2}\text{Tr ln}[i\omega _{m}+M_{\text{BdG}}(\mathbf{Q})].
\end{equation}
We can
further sum over Matsubara frequencies and obtain
\begin{equation}
S_{\text{fluc}}=\frac{1}{2%
}\sum_{j,\mathbf{q}}\text{ln cosh}\frac{E_{\text{BdG}}(j,\mathbf{q,0})}{2T}-\frac{%
1}{2}\sum_{j,\mathbf{q}}\text{ln cosh}\frac{E_{\text{BdG}}(j,\mathbf{q,Q})}{2T},
\end{equation}
with $E_{\text{BdG}}(j,\mathbf{q,Q})$ the BdG eigenfrequencies, which can be calculated numerically.
It is hard to proceed further to obtain an analytical
expression of $S_{\text{fluc}} (\mathbf{Q})$ due to the complicated setup.
Instead, we
assume  $S_{\text{fluc}}$ takes the form
\begin{equation}
S_{\text{fluc}%
}=\int_{0}^{1/T}d\tau \int d^{2}\mathbf{r}\frac{1}{2}\rho _{s}(\nabla\Theta)^{2}=\int_{0}^{1/T}d\tau \int d^{2}\mathbf{r}\frac{1}{2}\rho _{s}\mathbf{Q}^{2}
\end{equation}
for $\mathbf{Q}\rightarrow 0$ and a static phase $\Theta (\mathbf{r},\tau)=%
\mathbf{Q}\cdot \mathbf{r}$~\cite{PhysRevLett.109.105302S,PhysRevLett.114.110401S},
leading to $\rho _{s}=\frac{2T}{
\mathbf{Q}^{2}L}S_{\text{fluc}}$.
In the low temperature limit $T\rightarrow 0
$, we have
\begin{equation}\lim_{T\rightarrow 0}\frac{T}{L}S_{\text{fluc}}=\frac{1}{2L}\sum_{j,\mathbf{q}%
,E<0}[E_{\text{BdG}}(j,\mathbf{q,Q})-E_{\text{BdG}}(j,\mathbf{q,0})]=F(%
\mathbf{Q})-F(0).\end{equation}
Therefore, the results based on current response are
consistent with those obtained through the effective action. We have
numerically verified that, $\rho _{s}(T)=\frac{2T}{L\mathbf{Q}^{2}}S_{\text{fluc%
}}$ is nearly a constant at low temperature with $\rho _{s}(T_{\text{BKT}})\simeq\rho _{s}(T=0)$,
and $T_{\text{BKT}}=%
\frac{\pi }{2}\rho _{s}(T_{\text{BKT}})\simeq \frac{\pi }{2}\rho _{s}(T=0)$
is a very good approximation.

The attractive $s$-wave interaction pairs atoms from opposite valleys of the
bare hexagonal lattice and the intra-sublattice pairing is dominant. Note
that the pairing is between Moir\'{e} momentum $\pm \mathbf{q}$, which is
mainly determined by the bare Bloch states $\phi _{\uparrow l\mathbf{k}}$
and $\phi _{\downarrow l-\mathbf{k}}$ nearest to the valleys (which
contribute most to the flat bands). Due to the relative twist, $\pm \mathbf{k%
}$ are located at the same side of the corresponding valleys (as illustrated
in Fig.~\ref{figS:real_space_TB}c). Therefore, we have $\phi _{\uparrow l%
\mathbf{k}}\propto \lbrack 1,e^{i\gamma _{\uparrow \mathbf{k}}}]^{T}$ and $%
\phi _{\downarrow l-\mathbf{k}}\propto \lbrack 1,e^{i\gamma _{\downarrow -%
\mathbf{k}}}]^{T}$ in the $A$ and $B$ sublattice basis, with the relative
phases $\gamma _{\uparrow \mathbf{k}}\simeq -\gamma _{\downarrow -\mathbf{k}%
}+\pi $ (which are related to the chirality of the valleys, i.e., the Berry
phase on loops surrounding the valley). The relative phase $\gamma _{s%
\mathbf{k}}$ between $A$ and $B$ sublattice sites are responsible for the
Larkin-Ovchinnikov phase with staggered pairing order $\Delta (\mathbf{r}%
)\propto \langle \phi _{\uparrow l\mathbf{k}}\phi _{\downarrow l-\mathbf{k}%
}\rangle \propto \lbrack 1,-1]^{T}$.

\begin{figure}[tb]
\includegraphics[width=1.0\linewidth]{FigS3_r.pdf}
\caption{(a)-(c) The correlations $C_{\mathbf{q}}^{j^{\prime }j}$ for
different $j^{\prime },j$ with $U_{0}=-0.5$. (d) The superfluid band
structures (two middle particle/hole bands are not shown). (e) Superfluid
order $\Delta (\mathbf{r}=0)$ as a function of $U_{0}$ obtained by retaining
12 bands (circles) and 4 bands (plus signs) with $\protect\mu =0.00485$. (f)
Superfluid order $\Delta (\mathbf{r})$ at $U_{0}=-0.5$. Common parameters: $%
\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$ and $\Omega =0.116$.}
\label{fig:S3}
\end{figure}

\begin{figure}[b]
\includegraphics[width=1.0\linewidth]{FigS4.pdf}
\caption{(a) and (b) Moir\'{e} bands for $AB$ stacking with twist axis at
one coinciding site (where $A$ site of $V_{\uparrow }$ coincides with $B$
site of $V_{\downarrow }$). All other parameters in (a) and (b) are the same
as Figs.~2a and 2b in the main text, respectively. The bands are almost
identical with that in Figs.~2a and 2b in the main text, except that the two
valence (conduction) flat bands become degenerate at $K$ and $K^{\prime }$
points and a tiny gap develops between the valence and conduction flat
bands. (c) The normalized pairing order $\Delta (\mathbf{r})$ (maximum is
normalized to 1) starting from $AA$ stacking with twist axis at one hexagon
center. Other parameters are $\protect\theta =5.086^{\circ }$, $V_{0}=6E_{R}$%
, $\Omega =\Omega _{f}=0.116$ and $U_{0}=-0.5$.}
\label{fig:S4}
\end{figure}

The correlations $C_{\mathbf{q}}^{j^{\prime }j}=\langle \beta _{j^{\prime }-%
\mathbf{q}}\beta _{j\mathbf{q}}\rangle $ for different $j,j^{\prime }$ are
shown in Figs.~\ref{fig:S3}a-c. The pairing is mainly between bands from
different valleys, and the conduction bands from different valleys become
degenerate along the high symmetric $\Gamma $-$K$ lines with avoided
crossing (a tiny gap) due to inter-valley couplings. As we discussed in the
main text, $C_{\mathbf{q}}^{11}$ changes from characterizing $K_{s}$-$K_{%
\bar{s}}^{\prime }$ to $K_{s}^{\prime }$-$K_{\bar{s}}$ correlations;
therefore, $C_{\mathbf{q}}^{11}$ changes the sign across the $\Gamma $-$K$
lines where $C_{\mathbf{q}}^{12}$ is mainly distributed.
This means that though $C_{\mathbf{q}}^{j^{\prime }j}$ show $f$-wave
structures, their combined effects lead to the nearly uniform superfluid gap
(see Fig.~\ref{fig:S3}d) and the pairing is $s$-wave. The small superfluid
gap (weak pairing) at $\Gamma$ is due to the nearly uniform Moir\'{e}
wavefunction there (see Fig.~\ref{fig:S2}b). Moreover, we also calculate the
results by including 8 nearby Moir\'{e} bands (12 bands in total with 4 flat
bands, 4 higher and 4 lower bands) for comparison, and find that these
nearby bands have very little effects on the phase diagram and superfluid
order (as shown in Figs.~\ref{fig:S3}e and \ref{fig:S3}f).

In the above mean-field approach, we have assumed a real-space pairing
order. Alternatively, we can assume the momentum-space pairing order. In
particular, we first write the interaction Hamiltonian as
\[
H_{\text{int}}=\sum U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{Q;q,q^{\prime }}%
}\beta _{j_{1}\mathbf{Q+q}}^{\dag }\beta _{j_{2}\mathbf{Q-q}}^{\dag }\beta
_{j_{3}\mathbf{Q-q^{\prime }}}\beta _{j_{4}\mathbf{Q+q^{\prime }}}
\]%
with $U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{Q;q,q^{\prime }}}=U_{0}\int d^{2}%
\mathbf{r}\psi _{j_{1}\mathbf{Q+q}}^{\ast }(\mathbf{r},\uparrow )\psi _{j_{2}%
\mathbf{Q-q}}^{\ast }(\mathbf{r},\downarrow )\psi _{j_{3}\mathbf{Q-q^{\prime
}}}(\mathbf{r},\downarrow )\psi _{j_{4}\mathbf{Q+q^{\prime }}}(\mathbf{r}%
,\uparrow )$. Here $\mathbf{Q}$, $\mathbf{q}$ and $\mathbf{q^{\prime }}$ are
superlattice momenta in the Moir\'{e} BZ. We restrict the interaction to the
$\mathbf{Q}=0$ BCS channel, and assume the momentum-space order $\bar{\Delta}%
_{j_{1}j_{2}}(\mathbf{q})=\sum_{j_{3},j_{4},\mathbf{q}^{\prime
}}U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{0;q,q^{\prime }}}C_{\mathbf{q}^{\prime
}}^{j_{3}j_{4}}$. The correlation $C_{\mathbf{q}^{\prime }}^{j_{3}j_{4}}$
can be obtained by solving the BdG Hamiltonian%
\begin{equation}
H_{\text{BdG}}=\sum_{j,\mathbf{q}}(E_{j\mathbf{q}}-\mu)\beta _{j\mathbf{q}}^{\dag
}\beta _{j\mathbf{q}}+\sum_{j_{1},j_{2},\mathbf{q}}[\bar{\Delta}%
_{j_{1}j_{2}}(\mathbf{q})\beta _{j_{1}\mathbf{q}}^{\dag }\beta _{j_{2}-%
\mathbf{q}}^{\dag }+h.c.],
\end{equation}%
which allows us to obtain the superfluid order $\bar{\Delta}_{j_{1}j_{2}}(%
\mathbf{q})$ self-consistently. Using this approach, we calculate superfluid
order by keeping only the four flat bands and find that the orders $\bar{%
\Delta}_{j_{1}j_{2}}(\mathbf{q})$ are the same (up to tiny numeric errors)
as those obtained by assuming a real-space pairing order [i.e., $%
\sum_{j_{3},j_{4},\mathbf{q}^{\prime }}U_{j_{1}j_{2};j_{3}j_{4}}^{\mathbf{%
0;q,q^{\prime }}}C_{\mathbf{q}^{\prime }}^{j_{3}j_{4}}\simeq L\sum_{\mathbf{g%
}}\Delta _{\mathbf{g}}\chi _{j_{1}j_{2}}^{\mathbf{q}\ast }(\mathbf{g})$].
The two approaches lead to the same superfluid phase, correlation $C_{%
\mathbf{q}}^{j_{3}j_{4}}$, as well as $\Delta (\mathbf{r})$ (which is
determined by $C_{\mathbf{q}}^{j_{3}j_{4}}$).

\subsection{Effects of different stackings and twist axes}

We have focused on the twists starting from $AA$ stacking with the twist
axis at one sublattice site. Like the magic behaviors in TBG~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S,PNAS.108.12233S,PhysRevB.86.155449S}%
, here the twist axis or stacking position do not affect the appearance of
magic flat bands at small twist angles, as shown in Figs.~\ref{fig:S4}a and %
\ref{fig:S4}b. When $\theta $ is small, the Moir\'{e} bands for different
stackings and twist axes are almost identical. For $AB$ stacking with twist
axis at one coinciding site (where $A$ site of $V_{\uparrow }$ coincides
with $B$ site of $V_{\downarrow }$), the two valence (conduction) flat bands
become degenerate at $K$ and $K^{\prime }$ points and a tiny gap develops
between the valence and conduction flat bands. We find that the Moir\'{e}
bands for $AB$ ($AA$) stacking with twist axis at one hexagon center is the
same as that for $AA$ ($AB$) stacking with twist axis at one coinciding
site. For different stackings and twist axes, the superfluid orders and the
phase diagrams are similar, where the pairing order is staggered and
distributed mainly around the $AA$ region (see Fig.~\ref{fig:S4}c).

\subsection{Difference with twisted bilayer graphene}

In this section, we reemphasize some of the differences between our system
and the TBG system.

(1) The two twisted lattice potentials are state dependent, and one
potential does not affect atoms trapped by the other. This is different from
the electrons in TBG, where electrons in one layer can feel the potential of
atoms in the other layer.

(2) Our system is physically a single-layer system and we twist the lattice
for atomic (pseudo-)spin states (i.e., atomic internal energy levels). The $%
z $-direction is tightly confined by an additional state-independent
potential using the so-called magic-wavelength lasers. Therefore, the two
spin states have identical Wannier orbital along the $z$ direction. The
inter-spin tunnelings, realized by additional lasers, are different from the
inter-layer tunnelings between $p_{z}$ orbitals~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S} in TBG where a large
inter-layer distance exist. The existence of magic behaviors in our system
is not a straightforward derivative of TBG.

(3) The optical lattice potential here takes a cosine form which is much
simpler comparing to the atomic potential in graphene. Therefore the bare
bands and inter-spin couplings can be obtained accurately by directly
solving for the Bloch states $\phi _{sl\mathbf{k}_{s}}$ in our system. While
for the TBG, real-space tight-binding approximation based on Slater-Koster
parameters is usually adopted~\cite%
{PhysRevLett.99.256802S,PhysRevB.81.165105S}.

(4) Long-range tunnelings are more significant in our system because the
optical lattices considered here are relatively shallow. A shallow lattice
not only improves the atomic lifetime (through reducing the atomic decay
rate), but also increases the bare Dirac velocity $v_{D}$ (a larger $v_{D}$
leads to larger gaps and better isolation of the flat bands).

(5) The interactions are dominated by the $s$-wave scattering between atoms
in different spin states that are coupled and relatively twisted. In TBG,
the electronic interactions are more complex and include both Coulomb
repulsive interaction and/or phonon-mediated attractive interactions, which
mainly involve electrons in the same layer with no relative twist~\cite%
{PhysRevB.98.220504S,PhysRevLett.122.257002S,PhysRevLett.121.257001S,
PhysRevX.8.041041S,PhysRevLett.121.217001S,PhysRevLett.122.026801S}. The
unique interaction in our system can lead to interesting Larkin-Ovchinnikov
superfluid orders that do not exist in TBG.

(6) Finally, the advantage of cold atom system is that the parameters (e.g.,
inter-spin tunnelings, lattice depth, lattice constant, interactions, etc.)
are highly tunable. This not only leads to magic behaviors in a wide range
of parameter space, but also opens various possibilities for exploring novel
twistronics in cold atom systems.


\section{Introduction}
Although different researches are fighting against DDoS attacks, even today, the deployment or development of effective methods and mechanisms of these researches cannot resist DDoS attacks. It is a hard problem to solve because it is multifaceted. There are multifarious ways of creating attacks that are extremely effective. Section~\ref{sec:classification} and ~\ref{sec:ddos} give an overview of the classification and types of the different DDoS attacks. Section~\ref{sec:defense} gives an overview of the commonly used defense mechanisms. I am interested in working on types of attacks that so far we have not been able to handle effectively. These are mainly low-rate attacks and their related works are mentioned in Section~\ref{sec:noteffective}.

DDoS attacks possess immense threat to the current Internet community, which comprises over 3.5 billion users \cite{IU}. Figure \ref{fig:ddos-growth} shows the year-wise growth statistics of DDoS attacks, collated from Arbor Networks: Worldwide Infrastructure Security Report \cite{Arbor}, Calyptix Security Blog \cite{Caplyptix}, DDoS attack on Dyn\cite{Dyn} in 2016 and the largest volumetric DDoS attack, ever recorded, on Github \cite{Github}.
The 2016 Dyn cyberattack\cite{Dyn}, that took place on October 21, 2016 included multiple DDoS attacks targeting systems operated by Domain Name System (DNS) provider Dyn. It caused major Internet platforms and services to be unavailable to a large number of users spread across Europe and North America. And the magnitude of these attacks peaked to 1.2 Tbps.

The Github DDoS attack\cite{Github} generated a flood of internet traffic that peaked at 1.35Tbps, making it the largest on record. Interestingly, no botnets were involved in this attack. Instead, the hackers went with an amplification attack. They spoofed GitHub's IP address, and sent queries to several memcached servers that are typically used to speed up database-driven sites. The memcached systems then sent amplified responses from those requests to GitHub.

A computer or networked device under the control of an intruder is known as a bot. In general, in a DDoS attack, an attacker begins by exploiting a vulnerability in a computer system, making it the DDoS master bot. The attack master system identifies other vulnerable systems and gains control over them by either infecting the systems with malware or through bypassing the authentication controls. These systems are called DDoS slave bots. The attacker creates a command-and-control server to command the network of bots, referred to as a botnet.\cite{journals/ejwcn/LiuXGDZ09} Then, the attacker uses the traffic generated by the compromised devices to inundate the target so as to get its services down. Figure \ref{fig:ddos-arch} represents a typical DDoS architecture \cite{journals/cn/DouligerisM04}.

Also, the attackers aim that the DDoS attack packets do not reveal any obvious characteristic to segregate malicious and legitimate traffic. With minimal effort, the attackers aim to create the maximum-possible catastrophic impact.

\begin{figure}
\includegraphics[width=\columnwidth,scale=1]{DDoS-Stats}
\caption{Growth in DDoS Volume over the years \label{fig:ddos-growth}}
\end{figure}
The reasons or motivations behind DDoS attacks can be specific. But, in general, the major causes for DDoS attacks are as follows:

\begin{enumerate}
\item Competition: DDoS attacks can be intended to cripple company operations, damage reputation and devastate sales, which may directly benefit competitors.
\item Financial Gains: DDoS attacks can be used to achieve financial gains. Attackers can be hired, paid well or can even ask for ransom.
\item Politics: DDoS attacks have the potential to digitally silence opposition  parties. They can be used by political parties and terrorists.
\item Revenge: Current employees, ex-employees, angry customers or anyone with a dispute may have a motive for a DDoS attack.
\end{enumerate}

\begin{figure}
\includegraphics[width=\columnwidth,scale=1]{DDoS-Arch}
\caption{DDoS Attack Architecture \label{fig:ddos-arch}}
\end{figure}

\section{Classification of DDoS Attacks}
\label{sec:classification}

There are a wide variety of possible DDoS attacks. And it is challenging to represent all of them with a single classification. Different literatures have classified DDoS attacks in different ways \cite{journals/ccr/MirkovicR04, conf/ISCApdcs/SpechtL04, journals/comsur/ZargarJT13,journals/nca/GuptaB17}. Arbor networks have classified DDoS attacks as volumetric attacks and application layer attacks. Volumetric attacks send a huge amount of traffic, or request packets, to a targeted network so as to overwhelm its bandwidth capabilities. Typically request sizes are in the 100's of Gbps. However, recent attacks have scaled to over 1Tbps. Application-layer attacks generally require a lot less packets and bandwidth to get a website down. They are silent and small, when compared to network-layer attacks, but extremely disruptive.

Reference \cite{journals/ccr/MirkovicR04} presents the taxonomy of DDoS attacks and the different classifications include degree of automation, exploited weaknesses to deny service, source address validity, attack rate dynamics, possibility of characterization, persistence of agent set, victim type and impact on victim. Reference \cite{conf/ISCApdcs/SpechtL04} portrays a taxonomy of a few major DDoS attacks, shown in Figure \ref{fig:ddos-classification}. It includes two major classes of DDoS attacks: bandwidth depletion and resource depletion attacks. Bandwidth depletion attacks flood the victim
network with undesired traffic, preventing legitimate traffic from reaching the victim. And resource depletion attacks tie up the resources of a victim system making the victim unable to process legitimate requests.

\begin{figure}
\includegraphics[width=\columnwidth,scale=1]{DDoS-Classification}
\caption{Classification of DDoS Attacks \label{fig:ddos-classification}}
\end{figure}

\section{Different DDoS Attacks}
\label{sec:ddos}

The different DDoS attacks are enumerated below:
\subsection{Volumetric DDoS Attacks}

The volumetric DDoS attacks generate huge volumes of attack request or response traffic, usually measured in Gigabits per second (Gbps). The concept of a volumetric attack is simple, which is to send as much traffic as possible to a site to overwhelm its bandwidth. Some of the well known volumetric DDoS attacks are as follows:

\begin{enumerate}

\item DNS Amplification Attack: A DNS amplification attack is a DDoS attack in which the attacker exploits vulnerabilities in domain name system (DNS) servers to turn initially small queries into much larger payloads, which are used to bring down the victim’s servers. It is a reflection attack that manipulates publicly-accessible domain name systems, making them flood a target with large quantities of UDP packets. Using various amplification techniques, perpetrators can “inflate” the size of these UDP packets, making the attack so potent as to bring down even the most robust Internet infrastructure.

\item CharGEN Amplification Attack: A CharGEN amplification attack is usually carried out by sending small packets carrying a spoofed IP address of the target to the internet enabled devices running CharGEN. These spoofed requests to such devices are then used to send UDP floods as responses from these devices to the target. Internet-enabled printers and copiers have this protocol enabled by default and can be used to execute a CharGEN attack. This can be used to flood a target with UDP packets on port number 19. The server eventually exhausts its resources and goes down.

\item SNMP Amplification Attack: Simple Network Management Protocol(SNMP) is used for network management. SNMP amplification attack is carried out by sending small packets carrying spoofed IP address of the target to the internet enabled devices running SNMP. These spoofed requests to such devices are then used to send UDP floods as responses from these devices to the target. However, amplification effect in SNMP can be greater when compared with CHARGEN and DNS attacks. When the target tries to make sense of this flood of requests, it will end up exhausting its resources and will go down.

\item NTP Amplification Attack: Network Time Protocol(NTP) is used to synchronize the clocks of computers to some time reference. The NTP protocol is another publicly accessible network protocol. NTP amplification is essentially a reflection attack. Reflection attacks involve eliciting a response from a server to a spoofed IP address. The attacker sends a packet with a forged IP address, which belongs to the victim and the server replies to this address. The response is considerably larger than the request, amplifying the amount of traffic directed at the target server and ultimately leading to a degradation of service for legitimate requests.

\item SSDP Amplification Attack: The Simple Service Discovery Protocol(SSDP) is based on the Internet Protocol Suite for advertisement and discovery of network services and presence information. SSDP allows universal plug and play devices to send and receive information using UDP on port number 1900.  SSDP is attractive to DDoS attackers because of its open state that allows spoofing and amplification. The attacker conducts a scan looking for plug-and-play devices that can be utilized as amplification factors. As the attacker discovers networked devices, they create a list of all the devices that respond. The attacker creates a UDP packet with the spoofed IP address of the targeted victim. The attacker then uses a botnet to send a spoofed discovery packet to each plug-and-play device with a request for as much data as possible by setting certain flags, specifically ssdp:rootdevice or ssdp:all. As a result, each device will send a reply to the targeted victim with huge response.
The target then receives a large volume of traffic from all the devices and becomes overwhelmed and results in denial-of-service to legitimate traffic.

\item UDP Flood Attack: User Datagram Protocol(UDP) is a connectionless protocol that uses datagrams embedded in IP packets for communication without creating a session between two
devices. In a UDP Flood attack, a lot of UDP packets are sent to either random or specific ports on the victim system. In order to determine the requested application, the victim system processes the incoming data. In case of absence of the requested application on the requested port, the victim system sends a "Destination unreachable" message to the sender. Also, to conceal the identity, the attacker often spoofs the source IP address of the attacking packets. UDP flood attacks often deplete the bandwidth of network around the victim's system. Thereby, the systems close to the victim are also impacted due to the UDP flooding attack.

\end{enumerate}
\subsection{Protocol Attacks}

An internet protocol is basically a discrete set of rules to exchange information on the internet. These rules can be exploited by a bad actor. Some of the well known protocol attacks are as follows:

\begin{enumerate}
\item TCP SYN Flood Attack: The Transmission Control Protocol(TCP) SYN Flood attack is a DDoS attack that exploits the design of the three way TCP communication process between a client, host, and server. First, a client initiates a new session by sending a SYN packet to the server. The server responds with a SYNACK packet and then the client sends an ACK packet to the server before a connection is established. In a SYN flood attack, the attacker sends repeated SYN packets to every port on the targeted server, often using fake IP addresses. The server, unaware of the attack, receives multiple, apparently legitimate looking requests to establish communication. It responds to each of these with a SYN-ACK packet from each open port. Before the connection can time out, another SYN packet arrives. This leaves an increasingly large number of half-open connections. Ultimately, the server’s connection overflow tables get flooded and the service to legitimate clients are denied. Finally, the server may even malfunction or crash.

\item TCP SYN-ACK Flood Attack: The second step of the three-way TCP communication process is exploited by this DDoS attack. In a TCP SYN-ACK DDoS attack, attackers send spoofed SYN packets to a large numbers of servers and proxies on the Internet that generate large numbers of SYN-ACK packets in response to incoming SYN requests from the spoofed attackers. This SYN-ACK flood is not directed back to the botnet, but instead, is directed back to victim’s network and often exhausts the victim’s firewalls by forcing state-table lookups for every incoming SYN-ACK packet.  This denial of service attack can render stateful devices inoperable and can also consume excessive amounts of resources on routers, servers and IPS/IDS devices.

\item PUSH+ACK Flood Attack: When connecting with a server, the client can ask for confirmation that the information was received by setting the ACK flag, or it can force the server to process the information in the packet by setting the PUSH flag. Both requests require the server to do more work than with other types of requests. By flooding a server with spurious PUSH and ACK requests, an attacker can prevent the server from responding to valid traffic.

\item RST/FIN Flood Attack: After a successful three or four-way TCP-SYN session, RST or FIN packets are exchanged by servers to close the TCP-SYN session between a host and a client machine. In an RST or FIN Flood attack, a target server receives a large number of spoofed RST or FIN packets that do not belong to any session on the target server. This attack tries to exhaust a server’s resources as the server tries to process these invalid requests. The result is a server unavailable to process legitimate requests due to resources exhaustion.

\item IGMP Flood Attack: Internet Group Management Protocol(IGMP) is yet another connectionless protocol, used by IP hosts to report or leave their multicast group memberships for adjacent routers. An IGMP flood attack involves a large number of IGMP message reports being sent to a network or router, significantly slowing down and eventually preventing legitimate traffic from being transmitted across the target network.

\item ICMP Flood Attack/Ping Flood Attack: Internet Control Message Protocol (ICMP) is another connectionless protocol used for IP operations, diagnostics, and errors. Just as with a UDP flood, an ICMP flood does not rely on any specific vulnerability to achieve denial-of-service. An ICMP flood attack can involve an ICMP echo request. Once enough ICMP traffic is sent to a target server, it gets overwhelmed with requests, resulting in a denial-of-service.

\item Multi-Vector Attack: Attackers can also combine more than one DDoS attack mechanisms to fire DDoS attacks keeping the engineers dealing with the DDoS attack confused. These attacks are the hard to discern and are capable of taking down some of the best-protected servers and networks. For example, HTTP Flood Attack can be combined with
Recursive HTTP GET Flood Attack and fired together.

\item DNS Flood Attack:  In a Domain Name System(DNS) flood attack, the attacker targets one or more DNS servers belonging to a given zone, attempting to hamper resolution of resource records of that zone and its sub-zones. DNS servers help requesters find the servers they seek. A DNS zone is a distinct portion of the domain name space in the Domain Name System. The administrative responsibility is delegated to a single server cluster for each zone. In a DNS flood attack the offender tries to overbear DNS servers with apparently valid traffic, overwhelming server resources and impeding the servers' ability to direct legitimate requests to zone resources. The attacker sends multiple DNS requests to the victim’s DNS server directly or via a botnet. The DNS servers become flooded with requests,cannot process all of its incoming requests, eventually crashes.

\item Smurf Attack: The Smurf attack is a DDoS attack in which a large number of ICMP packets with the intended victim's spoofed source IP are broadcast to a computer network using an IP broadcast address. Most devices on a network will, by default, respond to this by sending a reply to the source IP address. If the number of machines on the network that receive and respond to these packets is very large, the victim's computer will be flooded with traffic. This can slow down the victim's computer to the point where it becomes impossible to work on.

\item Fraggle Attack: A Fraggle Attack is a DDoS attack that involves sending a large amount of spoofed UDP traffic to a router’s broadcast address within a network. It is very similar to a Smurf Attack.

\item SIP Flood Attack: The Session Initiation Protocol(SIP) flood attack is a DDoS attack that exploits the application layer protocol SIP, used in VOIP call setup. An attack can be made using different types of SIP request messages like SIP REQUEST, SIP INVITE or the SIP call control messages like SIP INFO, SIP NOTIFY, SIP RE-INVITE. The goal of this attack is to flood the proxy server or the SIP registration server and to consume all of its resources. Here, any army from the botnet sends thousands of messages to the SIP registrar server which is responsible to accept REGISTRAR requests as well as to keep records of the addresses and parameters of the user agents, unless it receives any server error message. As a result, the server overwhelms the legitimate users experience service outage and cannot reach the server.

\item IP Null Attack: Packets contain IPv4 headers which carry information about which transport protocol is being used. When attackers set the value of this field to zero, these packets can bypass security measures designed to scan TCP, IP, and ICMP. When the target server tries to put process these packets, it will eventually exhaust its resources and reboot.

\item HTTP Flood Attack: The real IP of the bots is used to avoid suspicion. The number of bots used to execute the attack is same as the source IP range for this attack. Since the IP addresses of the bots are not spoofed, there is no reason for defense mechanisms to flag these valid HTTP requests. A bot can be used to send a large number of GET, POST or other HTTP requests to execute an attack. Many bots can be collated in an HTTP DDoS attack to completely cripple the target server.

\item Single Session HTTP Flood: An attacker can exploit a loophole in HTTP 1.1 to send several requests from a single HTTP session. This allows attackers to send a large number of requests from a handful of sessions. In other words, attackers can bypass the limitations imposed by DDoS defense mechanisms on the number of sessions allowed. Single Session HTTP flood attacks target a server’s resources resulting in complete system shutdown or poor performance of the server.

\item Single Request HTTP Flood Attack: Single request HTTP Flood attack is based on the idea of sending several HTTP requests in a single HTTP session by masking these requests within one HTTP packet. When defense mechanisms evolved to block multiple incoming packets for a single session, attacks like Single Packet HTTP Flood were designed with workarounds to defeat these defenses. This evolution of an HTTP flood exploits another loophole in the HTTP technology. Several HTTP requests can be made by a single HTTP session by masking these requests within one HTTP packet. Keeping packet rates within the allowed limits, this technique allows an effective attack, exhausting a server’s resources.

\item Recursive HTTP GET Flood Attack: For an attack to be highly successful, it must remain undetected for as long as possible. This is yet another kind of HTTP flood attack wherein the attacker requests multiple website pages, analyses replies, and then recursively requests each object of the website. This method can be used together with any kind of HTTP to make the attack as unnoticeable as possible. It is hard to detect as the recursive requests of the website objects resemble legitimate requests.

\item Random Recursive GET Flood Attack: This attack is a purpose built variation of Recursive GET attack. It is designed for forums, blogs and other websites that have pages in a sequence. Like Recursive GET it also appears to be going through pages. Since page names are in a sequence, to keep up appearance as a legitimate user, it uses random numbers from a valid page range to send a new GET request each time. Random Recursive GET also aims to deflate its target’s performance with a large number of GET requests leading to denial-of-service for actual legitimate users.

\item IP Packet Option Field Attack: This attack targets the optional fields of an IP packet and randomizes its value. For example, if in an IP packet, all qualities of service bits are set to 1, it causes the victim to apply some additional time to analyze the packet. Thus, it can inundate the processing ability of the victim if a flood of packets arrives with such deformation.

\item Ping of Death Attack: In a ping of death attack, an attacker intentionally forms a data packet that exceeds the maximum packet size which causes the victims to freeze or crash. This attack can be initiated by the attacker without requiring a botnet. It is an IP fragmentation attack that exploits the inherent size limitation that a packet can be transmitted in.

\item Teardrop Attack: In this attack, the attacker  manipulates the offset value which in turn generates errors in fragmentation and reassembly of packets. The attacker sends fragmented packets with overlapping offset numbers. Thus, during the time of the packet re-assembly, invalid packets are created and crash the targeted server.

\item Session Attack: To bypass defenses, instead of using spoofed IP addresses, this attack uses the real IP address of the bots being used to carry out an attack. This attack is executed by creating a TCP-SYN session between a bot and the target server. This session is then stretched out until it times out by delaying the ACK packets. Session attacks try to exhaust a server's resources through these empty sessions. That, in turn, results in a complete system shutdown or unacceptable system performance.

\item UDP Fragmentation Flood Attack: It is another one of those cleverly masked DDoS attacks that are not easily detected. The activity generated by this attack resembles valid traffic and all of it is kept within limits. This version of the UDP Flood attack sends larger yet fragmented packets to exhaust more bandwidth by sending fewer fragmented UDP packets.When a target server tries to put these unrelated and forged fragmented UDP packets together, it will fail to do so. Eventually, all available resources are exhausted and the server goes down.

\item VoIP Flood Attack: This version of application specific UDP flood attack targets VoIP servers. An attacker sends a large number of spoofed VoIP request packets from a very large set of source IP. When a VoIP server is flooded with spoofed requests, it exhausts all available resources while trying to serve the valid and invalid requests.

This reboots the server or takes a toll on the server’s performance and exhausts the available bandwidth. VoIP floods can contain fixed or random source IP. Fixed source IP address attack is not easy to detect as it masks itself and looks no different from legitimate traffic.

\item ICMP Fragmentation Flood Attack: This version of ICMP Flood attack sends larger packets to exhaust more bandwidth by sending fewer fragmented ICMP packets. When the target server tries to put these forged fragmented ICMP packets with no correlation together, it will fail to do so. The server eventually exhausts its resources and reboots.
\end{enumerate}
\subsection{Low Rate DDoS Attacks}
\label{sec:low}
Detection of low rate DDoS attacks is extremely hard as they
possess the ability of concealing their traffic, because the traffic is
very much like the normal, legitimate traffic. They are able to elude the current
anomaly-based detection schemes. Some of the well known low rate DDoS attacks are as follows:

\begin{enumerate}
\item Flash-Crowd DDoS Attacks: A Flash Crowd Attack(FCA) is a DDoS attack that consumes the resources of a targeted service with legitimate looking requests generated by numerous bots. It is extremely hard to detect as bots may request legitimate content. An attacker may further employ several bots, each sending requests at a low rate.

\item TCP Shrew Attack: These are DDoS attacks which generate periodic, short bursts of high-volume traffic and create congestion. This forces legitimate TCP connections to drastically reduce their sending rate. Shrew attacks exploit, the deficiencies in the retransmission
time-out (RTO) mechanism of TCP flows. They throttle legitimate
TCP flows by periodically sending burst pulses with high peak
rate in a low frequency. As such, the TCP flows see congestion
on the attacked link every time it recovers from RTO. Indeed,
such a shrew attack may reduce the throughput of TCP applications
down to almost zero.

\item Hash Collisions DDoS Attacks: Most application servers create hash tables to index POST session parameters and are sometimes required to manage hash collisions when similar hash values are returned. Collision resolutions are
resource intensive, as they require an additional amount of CPU to
process the requests. In Hash Collision DDoS attacks, the
attacker sends a specially crafted POST message with a multitude
of parameters. These parameters are built in a way that causes hash
collisions on the server side. This slows down the response processing
dramatically. Hash Collisions DDoS attacks are very effective and can
be launched from a single attacker computer, slowly exhausting the
application server’s resources.

\item Costly Requests Attack: These are DDoS attacks wherein the attacker requests for content that requires costly database queries, which deplete bandwidth and connection pool between a front-end server and a database. The requests that involve long processing times, create a load on database processing and consume ample amount of server resources are selected by the attackers to launch the attack.

\item Slowloris Attack: A Slowloris DDoS attack works by opening multiple connections to the targeted web server. And it keeps them open as long as possible. This is achieved by continuously sending partial HTTP requests, none of which are ever completed. The attacked servers open more and connections open, waiting for each of the attack requests to be completed. As the requests are not completed, the targeted server’s maximum concurrent connection pool gets filled, and additional legitimate connection attempts are denied.

\end{enumerate}
\subsection{Vulnerability-Based DDoS Attacks}

Servers or protocols may have are unknown or unpatched  vulnerabilities. Attackers exploit these vulnerabilities to launch attacks. Examples of vulnerability-based DDoS attacks are:
\begin{enumerate}
\item Zero Day DDoS Attack: These are DDoS attacks that exploit new vulnerabilities. These ZERO Day DDoS vulnerabilities do not have patches or effective defensive mechanisms at the time of the attack.

\item LAND Attack: A LAND(Local Area Network Denial) attack is a DDoS attack that consists of sending a special poison spoofed packets to a computer, causing it to lock up.  Whenever a system receives this type of packets, it replies back to itself which in turn creates an infinite loop. As a result of this, the system crashes eventually.

\end{enumerate}

\section{DDoS Mitigation}
\label{sec:defense}

DDoS attacks pose immense threat to the resources of the victim and to the network bandwidth and infrastructure. Although for each attack, or a group of attacks, described in Section~\ref{sec:ddos}, there are specific defense mechanisms available in literature, some of the commonly used ways to mitigate DDoS attacks are as follows:

\begin{enumerate}

\item Hop-count filtering:
A commonly used defense is hop-count filtering \cite{jin2003hop} that filters packets that contain spoofed IP addresses. In this approach, the authors have used the concept that it is not possible to alter the number of hops of an IP packet when it travels from a source to a destination. Therefore, the authors have used this count to determine the validity of a packet. In this method, time to live value is used to count the number of hops. This hop counts are stored in a mapping table against each source address. Upon receiving a packet, the number of hops required for this packet is calculated and matched against the mapping table. A packet is detected as spoofed packet if a mismatch is found in this comparison. If a flow of packets is identified as a flow of spoofed packets, the filter discards those packets as a prevention of an attack. The deployment of such a technique is easier as it requires implementation in the victim's system. But, it has a limitation in the process of hop count. As this method counts the number of hops based on time to live, the number of false positive is larger in this method. This is because the initial time to live value is usually different for different operating systems. Additionally, the attacker can forge valid hop counts in their packets which allows the packets to pass the filter. Finally, for a flood of malicious packets, the system cannot perform the calculation and comparison. Hence it becomes the victim of the DDoS attack.

\item Source Address Validity Enforcement protocol:
Source Address Validity Enforcement (SAVE) protocol \cite{li2002save} enforces the routers to send messages containing updated source information to each destination routers connected to a source. Then, each router updates its forwarding table with current information and uses it to filter the packets based on the methods of RPF. Thus, it overcomes the problem encountered in the RPF for the asymmetric and dynamic nature of the Internet. However, implementation of such protocol requires a change in the existing routing protocol which is a time-consuming process. Also, a partial deployment of the protocol does not guarantee full success in filtering spoofed IP addresses.

\item Rate Limiting:
Rate limiting is one of common ways to defend the DDoS attacks. Generally, if the results of a detection mechanism are found to be partially successful, that is, if it produces large false negatives or cannot identify a precise distinction between the legitimate and malicious traffic, it is reasonable to apply rate limiting rather than filtering.

\item Ingress/Egress filtering:
Ingress/Egress filtering is a commonly used technique to prevent traffic with spoofed IP addresses to enter into a protected network. Basically, ingress filtering filters the malicious traffic destined to a local network and egress filtering discards the malicious traffic leaving a local network. Ingress filtering defined in RFC 226768 allows the traffic to enter the network which matches with a predefined range of domain prefixes of the network. If an attacker uses spoofed IP address that do not match with the prefixes, it is discarded at the routers. Thus, this filtering technique ensures prevention from DDoS attacks where spoofed IP addresses are used. However, it is not a useful mechanism in the cases where the valid IP addresses of the botnets are used as a source IP addresses during the attack. Also, the success of these filtering depends on the knowledge of the range of expected IP addresses for a port which is not always achievable for the complicated topologies used in different networks. Moreover, if an attacker uses some spoofed IP addresses that fall into the valid address range, the filters in the routers cannot detect the malicious traffic in such scenarios. Also, for these filters, itunneling is required for mobile IP users so that they are not filtered by the routers using ingress/egress filtering. Moreover, as it does not ensure incentives to the Internet Service Providers, it is partially deployed in the networks. Additionally, many Internet Service Providers do not enforce this approach.

\item Load balancing:
The key function of a load balancer is to spread workloads across multiple servers to prevent overloading servers, optimize productivity, and maximize uptime. Load balancers also add resiliency by re-routing live traffic from one server to another if a server is under a DDoS attack. Therefore, load balancers help to eliminate single points of failure and reduce the attack surface. Load Balancing is an approach which tries to balance the loads of different systems so that no one system gets overloaded. The result of the load balancing helps to gain the optimal productivity and the maximum uptime. In order to ensure the maximum load balancing, a bandwidth increase is required on all critical connections. A good number of replicated servers and data centers are required to guarantee elimination of single point of failure.

\item Martian address filtering and source address validation:
Martian address filtering is defined in RFC 1812. It works for filtering spoofed IP addresses that are generated from a limited set of addresses. This filtering ensures that a router must not forward any packet which has an invalid source or destination IP address. These invalid IP addresses range from reserved or special IP addresses as well as the unallocated range of IP addresses. Additionally, it ensures that any packet with the destination IP address 255.255.255.255/32 must be discarded at the router.

Source address validation is also specified in RFC 1812 and implemented in the routers. In this filtering mechanism, the router compares the source address of a packet with the logical interface of the packet where it is received. If this interface does not match with the interface where the packet needs to be destined to reach the specified source address, the router discards the packet. Thus, it can filter packets with spoofed source IP addresses. However, the rate of the false positive may become high for the asymmetric routes of the Internet. This asymmetric nature does not guarantee the match of the interfaces upon receiving or return of a packet from a specific source address. Thus, this filtering may discard a large number of legitimate traffics. However, the essential issue here is that not all routers in the Internet implement these approaches.

\item Capability-based response:
This mechanism can help prevent flooding related DDoS attacks. The receiver cannot control how much traffic it would receive. There exists flow control and congestion control mechanisms, but a misbehaving sender does not care about it. Therefore, the techniques involved in capability-based DDoS response work to discern a solution to control such misbehaving sender. Stateless inter flow filter (SIFF)134 and traffic validation architecture (TVA)135 are two example methods of this type.

\item Route-based Packet Filtering:
Route-based Packet Filtering(RPF) filters packets with spoofed source IP addresses. This filtering technique enhances the scope of the ingress filtering by providing service to the core routers. It does provide filtering based on the route information of a packet in each link of a core router. It depends on the principle that each link of the core router accepts traffic from only a limited number of source addresses. Thus, an IP packet which has a different source address than this set of addresses is discarded by the core routers as it appears to be spoofed to the router. In order to implement this technique, it requires information of the Border Gateway Protocol routing topology. According to the simulation of Park and Lee,\cite{park2001effectiveness} a significant success of the technique will be achieved if 18\% of the autonomous systems implement this filtering technique. However, this number is found impractical in current Internet scenario. Additionally, in order to include the source addresses in the BGP message, it is required to modify the scheme used in BGP messaging. This increases the message size and processing time of the BGP messages. Moreover, if the routers do not maintain updated information, this technique can discard legitimate packets for unwanted root change. Also, as RPF filters packets are based on BGP messaging information, an attacker can deceive root information and filtering rules by stealing BGP sessions. Finally, the attacker can carefully choose IP addresses that do not resemble the spoofed IPs. This can make the method ineffective in protecting DDoS from spoofed IP. Reference \cite{duan2008controlling} extends the idea of Park and Lee\cite{park2001effectiveness} by designing a packet filtering mechanism which considers update messages of local BGP to filter out spoofed IP addresses. This method is easy to deploy on the current architecture which relies on the BGP routing protocol. It also reduces the rate of false positive and simplifies IP traceback process.

\item History-based filtering:
History-based filtering is a packet marking–based filtering technique where the history of the normal traffic is used to filter out the malicious traffic\cite{peng2003protection}. Here, the destination of an attack maintains a database of IP addresses. This database comprises IP addresses which are commonly found in the destination. Therefore, when a bandwidth attack is targeted to this system, the system only allows those IP addresses that appear in their database and discards all other IP packets. However, this filtering technique cannot successfully detect and discard the malicious flow if an attacker can simulate its attack traffic as a normal traffic.

\item Path identifier:
The path identifier method \cite{yaar2003pi} eliminates packets based on a path identifier that identifies the path of the attacker. Identifying the path is a deterministic approach where each packet is stamped with an identifier based on the path it has traveled. The packets that travel the same path contain the same identifier. Thus, if the victim can identify a packet traveling from the attacker, it can filter all the subsequent packets sent by the attacker. This approach works well when half of the routers get involved to mark the packets. Because it works with a small sized identification field, there remains the possibility that different paths can show the same path information. Therefore, it can give false-positives and false-negatives. Reference \cite{yaar2006stackpi} proposes an improved version of the path identifier approach called StackPi that improves path identifier's performance in terms of incremental deployment. Additionally the improved filtering mechanism is capable of identifying malicious flows based on just a single packet. As per the paper, this method can provide a reasonable amount of DDoS protection only if 20\% of the routers implement this marking scheme. However, this method does not consider the detection of spoofed IP address packets rather it marks and filters a malicious packet based on the deterministic packet marking mechanism. Moreover, a host or router looking for the path information through this method may need to have some supportive software as well as expense of processing making this method difficult to deploy, as mentioned in reference \cite{beverly2009understanding}.

\item PacketScore:
PacketScore \cite{kim2006packetscore} is a proactive filtering technique that uses Bayes theorem to compute conditional legitimate probability(CLP). The conditional legitimate probability is used to determine the likelihood of a legitimate packet based on the baseline profile value and the attribute value of a packet. The packet filtering works based on this CLP value and a dynamic threshold. As the filtering takes into account the statistical analysis, this method works well for new attack signatures as well as non spoofed attack traffics. However, this method requires a large amount of storage to deal with the increasing number of attack attributes. It introduces a DDoS defense scheme that supports automated online attack characterizations and accurate attack packet discarding based on statistical processing. The key idea is to prioritize a packet based on a score which estimates its legitimacy given the attribute values it carries. Once the score of a packet is computed, this scheme performs score-based selective packet discarding where the dropping threshold is dynamically adjusted based on the score distribution of recent incoming packets and the current level of system overload.

\item Secure overlay:
The idea behind secure overlay \cite{adkins2003towards} is to build up an overlay network on top of the IP network.
It uses three principles: (i) enabling end-hosts to communicate without revealing their IP address, (ii) giving end-hosts control to defend against Denial-of-Service (DoS) attacks at the overlay level, and (iii) making sure that the added functionality does not introduce vulnerabilities not present in the Internet. This overlay network is the entry point for the outside network to establish a communication to the protected network. It is assumed that the isolation can be achieved if a protected network hides its IP addresses or uses a distributed firewall. This firewall ensures that only trusted traffic from the nodes of the overlay network can get entry to the protected network. Although this mechanism ensures prevention from the DDoS attack, it is applicable to a private network.

\item Honeypots:
A honeypot \cite{weiler2002honeypots} attracts attackers to attack them. Therefore, the actual system remains protected. Also, honeypots  can be used to extract important information like records of attack activity, tools, and software used for the attack which can reveal information about an attacker. This information is further used to detect and mitigate a DDoS attack. But this technique also contains some drawbacks. For example, the static and passive nature of the honeypots do not guarantee complete concealing from the attacker.

\item Changing IP addresses:
Changing IP address technique \cite{geng2000defeating} allows the computer system changes its IP address to invalidate an old address which may be the potential target of the DDoS attacks. This approach works for an IP address based attack. However, this also incurs a lot of overheads such as updating different entry information. This method works well till the attacker does not get informed about the new IP address. This approach provides the technological futility of addressing the problem solely at the local level.

\item Congestion policing:
The goal of a bandwidth flooding attack is to congest the resource\cite{deng2010implementation}. The impact of this type of attacks can be reduced or eliminated by applying congestion policing mechanism. Re-feedback\cite{briscoe2005policing} and NetFence\cite{liu2010netfence} are two example mechanisms where congestion policing is applied to defend DDoS attacks. Re-feedback ensures metrics in data headers such as time to live and congestion notification will arrive at each relay carrying a truthful prediction of the remainder of their path. NetFence uses a secure congestion policing feedback, to enable robust congestion policing inside the network.

\item Applying security patches:
It is also required to update all security patches regularly to ensure that the system is not affected by bugs or worms. Also, in order to prevent IoT botnet’s generation, it is mandatory to change the default or generic passwords of the IoT devices.86 This is an important awareness from the users side that can fight against the massive IoT–based DDoS attacks. However, the irony is that most of the users of the IoT devices are not aware of the threats or even they do not know the name of the DDoS attacks. When an attacker compromises a device for the DDoS army, the user of the device does not notice any change in their performance or behavior. Thus, stealthily they help in the attack without even knowing the name of the attack that they are participating.

\end{enumerate}

\section{Defenses for Low Rate attacks}
\label{sec:noteffective}
Low-rate DDoS attacks are not handled effectively by the current defense mechanisms. These attacks have ability of concealing their traffic because they are very much like normal traffic. They have the capacity to elude the current anomaly-based detection schemes. Instead of exhausting the network bandwidth, these attacks deny service at the device level. They deplete the resources at the service or the operating system, making it difficult for
the device unable to process legitimate clients' traffic.
These attacks are challenging to detect at the network level since they generate requests at low rate. Examples of such attacks are Flash Crowd Attack, TCP Shrew Attack\cite{kuzmanovic2003low}, Slowloris Attack, Hash-Collision Attack, described in \ref{sec:low}.

Instead of a generic approach, most defenses focus on just a single variant. There is no generic framework
to handle all low-rate DDoS attacks.

Some related work for the slow attacks is as follows:

\begin{enumerate}

\item Reference Architecture: Reference \cite{shtern2014towards} proposes a reference architecture
for defense against low-rate DDoS attacks, enabled by a Software
Defined Infrastructure(SDI). Here, both network
and computational resources are provisioned and modified
on demand to achieve the mitigation goals. The suspicious
traffic is detected and redirected to the Shark Tank system, which is
created on-demand. The Shark Tank system is a separate cluster which possesses full application capabilities designed
to monitor suspicious users. The key feature of the Shark Tank is to act
as a restricted area for the attacker so that they can be
placed under close surveillance. It does absorb
the damage from the DDoS attacks and allows the
system to learn from these attacks for future reference.

\item Collaborative Detection and Filtering: Reference \cite{chen2006collaborative} presents a defense by combining discrete Fourier transform (DFT) and a hypothesis test framework to cope with shrew attacks. By computing the autocorrelation sequence of sampled time series and converting them into frequency-domain spectrum using DFT, the authors find that the power spectrum density (PSD) of a traffic stream containing shrew attacks that has much higher energy in low-frequency band than that appears in the spectrum for legitimate TCP/UDP traffic streams.
Based on this distinction, we develop a distributed collaborative detection and filtering (CDF) scheme to detect and segregate the shrew attack flows from legitimate TCP/UDP traffic flows. In addition to software implementation, the scheme can be implemented by network processor or reconfigurable hardware. The DSP hardware pushes spectral analysis down to the lower packet-processing layer. If the packets are processed by hardware and malicious flows are filtered out timely, the router workload will not increase much in the presence of shrew attacks.

\item Randomization on TCP Retransmission Timeout (RTO): Reference \cite{yang2004defense} proposes randomization on TCP RTO as defense against such attacks. With RTO randomization, an attacker cannot predict the next TCP timeout and consequently cannot
inject the burst at the exact instant. The authors show that randomization can effectively mitigate the impact of DDoS attacks while maintaining fairness and friendliness to other connections.

\item HAWK: Reference \cite{kwok2005hawk} proposes a new stateful adaptive queue management technique called HAWK (Halting Anomaly with Weighted choKing) which works by judiciously identifying malicious shrew packet flows using a small flow table and dropping such packets decisively to halt the attack such that well-behaved TCP sessions can re-gain their bandwidth shares. The authors show that HAWK is highly agile.

\item Shrew Attack Protection: Reference \cite{chang2010taming}  presents a simple priority-tagging filtering
mechanism, called SAP (Shrew Attack Protection), that protects
well-behaved TCP flows against low-rate TCP-targeted
Shrew attacks. In this method, a router maintains a simple set of counters and keeps track of the drop rate for each potential victim. If the monitored drop rates are low, all packets are treated as normal and equally compete to be admitted to the output queue and only dropped based on the AQM (Active Queue Management) policy when the output queue is almost full. However, if the drop rate for a certain victim becomes higher than some dynamically determined threshold (called fair drop rate), the router treats packets for this victim as high-priority, and these high-priority packets are preferentially admitted to the output queue. SAP marks victim packets as high priority until their drop rate is below the fair drop rate. By preferentially dropping normal packets to protect high-priority packets, SAP can prevent low rate TCP-targeted Shrew attacks from causing a well-behaved TCP flow to lose multiple consecutive packets repeatedly. This technique protects well-behaved TCP flows away from near zero throughput when under an attack.

\item Treat as application-level attacks and discern anomalies: A way of detecting low rate DDoS attacks is to treat them as application-level attacks and look for anomalies
in the payload of the incoming service requests. This is the approach taken by DDoS defense providers
and requires costly deep-packet inspection. Reference \cite{jonker2016measuring} presents an investigation in the adoption of cloud-based DDoS Protection Service providers worldwide. The authors focus on nine leading providers according to namely Akamai, CenturyLink, CloudFlare,
DOSarrest, F5 Networks, Incapsula, Level 3, Neustar, and
Verisign. The investigation is done on the basis of long-term, active DNS measurements, which allows for a given domain name, to verify if traffic diversion towards a DDoS Protection Service is in place or not.

\item Feature Feature Score(FFSc): Reference \cite{hoque2016ffsc} introduces a statistical measure called Feature Feature Score
(FFSc) for multivariate data analysis to distinguish DDoS
attack traffic from normal traffic. The authors extract three features of
network traffic, entropy of source IPs, variation of source
IPs and packet rate to analyze the behavior of network traffic
for attack detection. They build profiles of these features during normal operation and detect attack traffic
as the traffic that deviates from these profiles.

\item New Information Metrics: Reference \cite{xiang2011low} proposes that the generalized entropy and information distance
metrics outperform the traditional Shannon entropy
and Kullback–Leibler distance metrics for the low-rate
DDoS attack detection in terms of early detection, lower
false positive rates, and stabilities. It proposes an effective IP traceback scheme based on an
information distance metric that can trace all attacks back
to their own local area networks (LANs) in a short time. The proposed generalized entropy metric can detect attacks several hops earlier than the traditional Shannon metric approach. The proposed information distance metric outperforms the popular Kullback-Leibler divergence approach as it can enlarge the adjudication distance and then obtain the optimal detection sensitivity.

\item Expectation of Packet Size(EPS): Reference \cite{zhou2017low} presents a
measurement, expectation of packet size, that is based on the distribution difference of the packet size to distinguish
low-rate DDoS attacks from legitimate traffic. The authors propose an
expectation of packet size measurement to distinguish
attack traffic from legitimate traffic. The proposed method
is independent of the attack pattern; therefore, it can avoid the inherent shortcomings of the signature-based metric. Moreover, this approach can achieve a large distance gap between the attack traffics and the legitimate traffics, which could further contribute to a low false-positive rate.

\item Traffic Cluster Entropy: Reference \cite{sachdeva2014traffic} proposes traffic cluster entropy as detection metric not only to detect DDoS attacks but also to distinguish DDoS attacks from Flash Events. The authors show that when flash events are triggered, source address entropy increases but the corresponding traffic
cluster entropy does not increase. However, when DDoS attack is launched, traffic cluster entropy also increases along with source address entropy.

\item E--LDAT: Reference \cite{bhuyan2016ldat} proposes E‐LDAT, a lightweight extended‐entropy metric‐based system for both DDoS flooding attack detection and IP traceback. It aims to identify DDoS attacks effectively by measuring the metric difference between legitimate traffic and attack traffic. IP traceback is performed using the metric values for an attack sample detected by the detection scheme. This approach uses a generalized entropy metric with packet intensity computation on the sampled network traffic with respect to time. The authors show that E‐LDAT works well for detecting four classes of DDoS flooding attacks, including constant rate, pulsing rate, increasing rate and subgroup attacks.

\item Machine Learning: Reference \cite{saied2016detection} proposes machine learning to detect and mitigate known and unknown DDoS attacks in real time environments. The authors chose an Artificial Neural Network (ANN) algorithm to detect DDoS attacks based on specific characteristic features (patterns) that separate DDoS attack traffic from genuine traffic.

\item Flow Correlation Coefficient: Reference \cite{yu2012discriminating} proposes a
discrimination algorithm using flow correlation coefficient as a similarity metric among suspicious flows. The authors present theoretical proofs for the feasibility of the proposed discrimination method.

\item Multifractal Detrended Fluctuation Analysis(MF--DFA): Reference \cite{wu2016low} proposes the algorithm of multifractal detrended fluctuation analysis to explore the change in terms of
multifractal characteristics over a small scale of network traffic due to low rate DDoS attacks. Through wavelet analysis, the singularity and
bursty of network traffic under low rate DDoS attacks are estimated by using Holder exponent. The difference values, known as D--value, of Holder
exponent of network traffic between normal and under low rate DDoS attack situations are calculated. The D--value is used as the basis to
determine low rate DDoS attacks. A detection threshold is set based on the statistical results. The presence of low rate DDoS attacks can be confirmed
through comparing D--value with detection threshold. Experiments on detection performance have been performed in the test-bed
network and simulation platform. The extensive experimental results are congruent with the theoretical analysis.

\item FlowTrApp: Reference \cite{buragohain2016flowtrapp} proposes DDoS defense using FlowTrApp. Software Defined Network (SDN)  provides a central
control over the network which helps in getting the global view of the network. FlowTrApp which performs DDoS detection and mitigation using some bounds on two per flow based traffic parameters that is flow rate and flow duration of a flow. It attempts to detect attack traffic ranging from low rate to high rate and long lived to short lived attacks using an SDN engine consisting of sFlow based flow analytics engine sFlow-RT and an OpenFlow controller.

\item K--Nearest
Neighbors traffic classification with correlation analysis (CKNN): Reference \cite{xiao2015detecting} presents a detection approach based on K-nearest
neighbors traffic classification with correlation analysis to detect DDoS attacks. The approach exploits correlation
information of training data to improve the classification accuracy and reduce the overhead caused
by the density of training data. Aiming at solving the huge cost, the authors also present a grid-based method named
r-polling method for reducing training data involved in the calculation.

\item  Matching Pursuit and Orthogonal Matching
Pursuit algorithms: Reference \cite{andrysiak2013ddos} proposes to use Matching Pursuit and Orthogonal Matching
Pursuit algorithms for DDoS detection. The major contribution of the paper is the proposition of 1D
KSVD algorithm as well as its tree based structure representation (clusters), that can be applied to detect low-rate DDoS attacks and network anomalies. The method maintains the features of malicious traffic in form of tree to improve the classification performance and to optimize the implementation in terms of time complexity.

\end{enumerate}

\section{Conclusions}

In this paper, we presented the classification and literature review of DDoS attacks. In spite of scrupulous research over the years to mitigate DDoS attacks, they still exist today, in fact, with larger intensities and impacts as there are different angles to consider. We also described some of the commonly used defense mechanisms. Also, we presented a literature review of the defenses for low-rate DDoS attacks that have not been handled effectively till today.

\begin{acks}

Special thanks to Dr. Jelena Mirkovic and Dr. Genevieve Bartlett for their valuable guidance and advice.

\end{acks}

\bibliographystyle{ACM-Reference-Format}
\input{main.bbl}