proposed in Ref. [64] and [65]. The divergence of the atomic displacement
field $|\boldsymbol{\nabla}\cdot\delta\boldsymbol{U}|\lesssim 10^{-5}$ is tiny
and can be safely neglected.

When $S\neq S^{\prime}$, by $C_{3}$ and $m_{y}$ mirror symmetries, the
summation in Eq. 32 produces

$\displaystyle\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}\left[\delta_{\rho\nu}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)+(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|}$
$\displaystyle=$ $\displaystyle
v_{F}\gamma\left(\left(\tau_{3}\right)_{\rho\nu}(\sigma_{1})_{SS^{\prime}}+\left(\tau_{1}\right)_{\rho\nu}(\sigma_{2})_{SS^{\prime}}\right)$
(36) $\displaystyle\mbox{with}\quad v_{F}\gamma=$
$\displaystyle\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})}\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{x}\right]^{2}\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|}$
(37)

To further simplify the notation, we introduce the layer dependent pseudo-
vector field $\mathcal{A}^{(j)}$ as

$\displaystyle\mathcal{A}^{(j)}_{\mu}(\boldsymbol{x})=\pm
R\left(\frac{\theta_{j}}{2}\right)_{\mu\nu}\mathcal{A}_{\nu}\approx\pm(\mathcal{A}_{x}-\frac{\theta_{j}}{2}\mathcal{A}_{y}\
,\ \mathcal{A}_{y}+\frac{\theta_{j}}{2}\mathcal{A}_{x})_{\mu}$ with
$\displaystyle\mathcal{A}=(\partial_{x}\delta U_{x}-\partial_{y}\delta U_{y}\
,\ -(\partial_{x}\delta U_{y}+\partial_{y}\delta
U_{x}))\approx(2\partial_{x}\partial_{y}\varepsilon^{U}(\boldsymbol{x})\ ,\
(\partial_{x}^{2}-\partial_{y}^{2})\varepsilon^{U}(\boldsymbol{x}))$ (38)

where the sign $+$ and $-$ are for the top and bottom layers respectively.

Now, combining with the expansion of the Jacobi factor and keep the terms up
to the second order of derivatives,

$\displaystyle v_{F}\gamma\int{\rm d}^{2}\boldsymbol{x}\
\Psi_{j,S}^{\dagger}(\boldsymbol{x})\left(\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\cdot\mathcal{A}^{(j)}-\left(\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\cdot\mathcal{A}\right)\frac{1}{2}\boldsymbol{\nabla}\cdot\delta\boldsymbol{U}\right)\Psi_{j,S^{\prime}}(\boldsymbol{x})$
(39)

###  A.5 second order derivative terms

In this section, we considered the next order terms, i.e. the 2nd order of
combined derivative of $\delta\boldsymbol{U}$ and gradient of fermion fields.

####  A.5.1 second order gradient of fermion field

First, consider the term

$\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)\frac{1}{2}(-iy_{\mu})(-iy_{\nu})\Psi^{\dagger}_{j,S}(\boldsymbol{x})p_{\mu}p_{\nu}\Psi_{j,S^{\prime}}(\boldsymbol{x})$

Applying the Poisson summation formula,

$\displaystyle-\frac{1}{2}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)y_{\mu}y_{\nu}=-\frac{1}{2}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\left.\frac{\partial^{2}}{\partial
q_{\mu}\partial
q_{\nu}}\tilde{t}(\boldsymbol{q})\right|_{\boldsymbol{q}=\boldsymbol{G}+\boldsymbol{K}}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)$
(40)

For $S=S^{\prime}$, $\delta\boldsymbol{\tau}=0$. We obtain

$\beta_{0}=-\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}a_{\mu}a_{\nu}t(|\boldsymbol{a}|)=-\frac{1}{4}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}|\boldsymbol{a}|^{2}t(|\boldsymbol{a}|)\
.$

For $S\neq S^{\prime}$, by $C_{3}$ and $m_{y}$ symmetry, we obtain

$\displaystyle-\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)=\beta_{1}\left[(\tau_{3})_{\mu\nu}\sigma_{1}+(\tau_{1})_{\mu\nu}\sigma_{2}\right]_{SS^{\prime}}$
(41) with
$\displaystyle\beta_{1}=-\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})}\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{x}\right]^{2}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)\
.$ (42)

Thus, combining together, we obtain

$\displaystyle\int{\rm d}^{2}\boldsymbol{x}\
\sum_{j}\sum_{SS^{\prime}}\Psi_{j,S}^{\dagger}(\boldsymbol{x})\left[\beta_{0}\boldsymbol{p}^{2}\delta_{SS^{\prime}}+\beta_{1}\left((p_{x}^{2}-p_{y}^{2})\sigma_{1}+2p_{x}p_{y}\sigma_{2}\right)_{SS^{\prime}}\right]\Psi_{j,S}(\boldsymbol{x})\
.$ (43)

####  A.5.2 Cross terms between atomic displacement gradients and gradients of
the fermion field

Next, consider the cross term between the first order derivative of
$\delta\boldsymbol{U}$ and the gradient of the fermion field.

$\displaystyle
A_{mlg}^{-1}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)\frac{i}{2}y_{\mu}\partial_{\mu}\delta
U_{\rho}(\boldsymbol{G}+\boldsymbol{K})_{\rho}(-i)y_{\nu}\frac{1}{2}\left((p_{\nu}\Psi_{j,S}(\boldsymbol{x}))^{\dagger}\Psi_{j,S^{\prime}}+\Psi_{j,S}^{\dagger}(\boldsymbol{x})p_{\nu}\Psi_{j,S^{\prime}}(\boldsymbol{x})\right)$
(44)

Applying the Poisson summation formula, we can obtain:

$\displaystyle\frac{1}{2}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)y_{\mu}y_{\nu}(\boldsymbol{G}+\boldsymbol{K})_{\rho}=-\frac{1}{2}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\left.\frac{\partial^{2}}{\partial
q_{\mu}\partial
q_{\nu}}\right|_{\boldsymbol{q}=\boldsymbol{G}+\boldsymbol{K}}(\boldsymbol{G}+\boldsymbol{K})_{\rho}$
$\displaystyle=$
$\displaystyle-\frac{i}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}\left[\delta_{\mu\rho}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)+\delta_{\nu\rho}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)\right.$
$\displaystyle\left.+(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|}\right]$
(45)

For $S=S^{\prime}$, the above formula leads to

$\displaystyle
C_{\mu\nu\rho}=-\frac{i}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}a_{\mu}a_{\nu}a_{\rho}\frac{t^{\prime}(|\boldsymbol{a}|)}{|\boldsymbol{a}|}$
(46)

By $C_{3}$ and $m_{y}$ symmetries, we can prove that the only non-zero
components of the tensor $C_{\mu\nu\rho}$ are

$C_{xxx}=-C_{xyy}=-C_{yyx}=-C_{yxy}=C_{0}\quad\mbox{with}\quad
C_{0}=-\frac{i}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}\left(a_{x}\right)^{3}\frac{t^{\prime}(|\boldsymbol{a}|)}{|\boldsymbol{a}|}$

Then Eq. 44, when $S=S^{\prime}$, can be written as

$\displaystyle C_{0}\int{\rm d}^{2}\boldsymbol{x}\
\Psi_{j,S}^{\dagger}(\boldsymbol{x})\frac{1}{2}\left\\{\boldsymbol{p}\
,\mathcal{A}(\boldsymbol{x})\right\\}\Psi_{j,S}(\boldsymbol{x})$ (47)

Now, for $S\neq S^{\prime}$, the summation in Eq. 45 gives

$\displaystyle\frac{v_{F}}{2}\left[\delta_{\mu\rho}(\bar{\boldsymbol{\sigma}}_{SS^{\prime}})_{\nu}+\delta_{\nu\rho}\left(\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\right)_{\mu}\right]+D_{0}\left[\delta_{\mu\rho}(\bar{\boldsymbol{\sigma}}_{SS^{\prime}})_{\nu}+\delta_{\nu\rho}\left(\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\right)_{\mu}+\delta_{\mu\nu}(\bar{\boldsymbol{\sigma}}_{SS^{\prime}})_{\rho}\right]$
$\displaystyle=$
$\displaystyle\left(\frac{v_{F}}{2}+D_{0}\right)\delta_{\mu\rho}\bar{\sigma}_{\nu,SS^{\prime}}+\left(\frac{v_{F}}{2}+D_{0}\right)\delta_{\nu\rho}\bar{\sigma}_{\mu,SS^{\prime}}+D_{0}\delta_{\mu\nu}\bar{\sigma}_{\rho,SS^{\prime}}$
(48) $\displaystyle\mbox{with}\quad D_{0}$
$\displaystyle=-\frac{i}{6}\sum_{\boldsymbol{a}}e^{i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})}\left((\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{x}\right)^{3}\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|}\
.$ (49)

This leads to the term

$\displaystyle\pm\Psi_{j,S}^{\dagger}(\boldsymbol{x})\frac{1}{2}\left(p_{\mu}\xi_{\mu,SS^{\prime}}(\boldsymbol{x})+\xi_{\mu,SS^{\prime}}(\boldsymbol{x})p_{\mu}+\left(\frac{v_{F}}{2}+D_{0}\right)\frac{1}{2}\left\\{\boldsymbol{\nabla}\cdot\delta\boldsymbol{U},\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\cdot\boldsymbol{p}\right\\}\right)\Psi_{j,S^{\prime}}(\boldsymbol{x})$
(50) with
$\displaystyle\xi_{x,SS^{\prime}}(\boldsymbol{x})=\left[(\frac{v_{F}}{2}+2D_{0})\partial_{x}\delta
U_{x}\sigma_{1}-\left[(\frac{v_{F}}{2}+D_{0})\partial_{y}\delta
U_{x}+D_{0}\partial_{x}\delta U_{y}\right]\sigma_{2}\right]_{SS^{\prime}}\ ,$
$\displaystyle\boldsymbol{\xi}_{y,SS^{\prime}}(\boldsymbol{x})=\left(\left[(\frac{v_{F}}{2}+D_{0})\partial_{x}\delta
U_{y}+D_{0}\partial_{y}\delta
U_{x}\right]\sigma_{1}-(\frac{v_{F}}{2}+2D_{0})(\partial_{y}\delta
U_{y})\sigma_{2}\right)_{SS^{\prime}}\ .$ (51)

####  A.5.3 $(\boldsymbol{\nabla}\delta\boldsymbol{U})^{2}$ terms

Lastly, we consider the terms containing the square of the gradient of
$\delta\boldsymbol{U}$:

$\displaystyle
A_{mlg}^{-1}\frac{1}{2}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)\frac{i}{2}y_{\mu}\frac{i}{2}y_{\nu}\partial_{\mu}\delta
U_{\rho}(\boldsymbol{G}+\boldsymbol{K})_{\rho}\partial_{\nu}\delta
U_{\sigma}(\boldsymbol{G}+\boldsymbol{K})_{\sigma}\Psi_{j,S}(\boldsymbol{x})^{\dagger}\Psi_{j,S^{\prime}}(\boldsymbol{x})$
(52)

Again, we consider the coefficient

$\displaystyle-\frac{1}{8}A_{mlg}^{-1}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(|\boldsymbol{y}|)y_{\mu}y_{\nu}(\boldsymbol{G}+\boldsymbol{K})_{\rho}(\boldsymbol{G}+\boldsymbol{K})_{\sigma}$
$\displaystyle=$
$\displaystyle\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}\left.\frac{\partial^{2}}{\partial
x_{\rho}\partial
x_{\sigma}}\left((\boldsymbol{x}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{x}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}t(|\boldsymbol{x}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)\right)\right|_{\boldsymbol{x}=\boldsymbol{a}}$
$\displaystyle=$
$\displaystyle\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}\left[(\delta_{\mu\sigma}\delta_{\rho\nu}+\delta_{\nu\sigma}\delta_{\rho\mu})t(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)+\right.$
$\displaystyle\left(\delta_{\nu\sigma}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}+\delta_{\mu\sigma}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}+\delta_{\rho\sigma}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}+\right.$
$\displaystyle\left.+\delta_{\nu\rho}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\sigma}+\delta_{\rho\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\sigma}\right)\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|}+$
$\displaystyle\left.(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\sigma}\left(\frac{t^{\prime\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|^{2}}-\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|^{3}}\right)\right]\
.$ (53)

For $S=S^{\prime}$, the above summation can be simplified as

$\displaystyle\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}\left[\left(\delta_{\mu\sigma}\delta_{\rho\nu}+\delta_{\nu\sigma}\delta_{\rho\mu}\right)\left(t(|\boldsymbol{a}|)+\frac{7}{8}|\boldsymbol{a}|t^{\prime}(|\boldsymbol{a}|)+\frac{1}{8}|\boldsymbol{a}|^{2}t^{\prime\prime}(|\boldsymbol{a}|)\right)+\delta_{\rho\sigma}\delta_{\mu\nu}\left(\frac{3}{8}|\boldsymbol{a}|t^{\prime}(|\boldsymbol{a}|)+\frac{1}{8}|\boldsymbol{a}|^{2}t^{\prime\prime}(|\boldsymbol{a}|)\right)\right]$
$\displaystyle=$
$\displaystyle\left(\frac{\mu}{8}+C_{1}\right)\left(\delta_{\mu\sigma}\delta_{\rho\nu}+\delta_{\nu\sigma}\delta_{\rho\mu}\right)+C_{2}\delta_{\rho\sigma}\delta_{\mu\nu}\
,$ (54)

with the expression for the two coefficients $C_{1}^{\prime}$ and $C_{2}$
listed below

$\displaystyle C_{1}$
$\displaystyle=\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}\left(\frac{7}{8}|\boldsymbol{a}|t^{\prime}(|\boldsymbol{a}|)+\frac{1}{8}|\boldsymbol{a}|^{2}t^{\prime\prime}(|\boldsymbol{a}|)\right)$
(55) $\displaystyle C_{2}$
$\displaystyle=\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}\left(\frac{3}{8}|\boldsymbol{a}|t^{\prime}(|\boldsymbol{a}|)+\frac{1}{8}|\boldsymbol{a}|^{2}t^{\prime\prime}(|\boldsymbol{a}|)\right)$
(56)

This leads to the term

$\displaystyle\left[\left(\frac{\mu}{8}+C_{1}\right)\left((\boldsymbol{\nabla}\cdot\delta
U)^{2}+(\partial_{\mu}\delta U_{\nu})(\partial_{\nu}\delta
U_{\mu})\right)+C_{2}(\partial_{\mu}\delta U_{\nu})(\partial_{\mu}\delta
U_{\nu})\right]\Psi^{\dagger}_{j,S}(\boldsymbol{x})\Psi_{j,S}(\boldsymbol{x})$
(57)

For $S\neq S^{\prime}$, the first term in the square bracket of Eq. 53
vanishes by $C_{3}$ symmetry. Introducing the notation

$\displaystyle R_{\mu\nu\rho\sigma,SS^{\prime}}$
$\displaystyle=\frac{1}{4}v_{F}\gamma\left[\delta_{\mu\sigma}\left((\tau_{3})_{\rho\nu}\sigma_{1}+(\tau_{1})_{\rho\nu}\sigma_{2}\right)+(\mu\leftrightarrow\rho)+(\mu\sigma\leftrightarrow\rho\nu)\right]_{SS^{\prime}}$
(58) $\displaystyle V_{\mu\nu\rho\sigma,SS^{\prime}}$
$\displaystyle=\frac{1}{8}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\mu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\nu}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\rho}(\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}})_{\sigma}\left(\frac{t^{\prime\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|^{2}}-\frac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{SS^{\prime}}|^{3}}\right)$
$\displaystyle=\frac{D_{1}}{6}\left[\delta_{\nu\mu}\left((\tau_{3})_{\rho\sigma}\sigma_{1}+(\tau_{1})_{\rho\sigma}\sigma_{2}\right)+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\rho\leftrightarrow\nu)+(\nu\leftrightarrow\sigma)+(\rho\mu\leftrightarrow\nu\sigma)\right]_{SS^{\prime}}$
$\displaystyle\
+D_{2}\left(\cos(\frac{\pi}{2}n_{2})\sigma_{1}-\sin(\frac{\pi}{2}n_{2})\sigma_{2}\right)_{SS^{\prime}}$
(59) $\displaystyle\zeta_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\left(R_{\mu\nu\rho\sigma,SS^{\prime}}+V_{\mu\nu\rho\sigma,SS^{\prime}}\right)\partial_{\mu}\delta
U_{\rho}\partial_{\nu}\delta U_{\sigma}$ (60)

where $n_{2}$ in Eq. 59 is the number of times that the index $y$ appears in
the subscripts $\mu$, $\nu$, $\rho$, and $\sigma$, i.e.

$n_{2}=\delta_{\mu y}+\delta_{\nu y}+\delta_{\rho y}+\delta_{\sigma y}\ .$

Eq. 59 is derived by $C_{3}$ symmetry, and the constants $D_{1}$ and $D_{2}$
must be real as constrained by $m_{y}$ symmetry. To be more specific,

$\left\\{\begin{array}[]{l}D_{1}=\dfrac{1}{16}\sum\limits_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})}\left(\dfrac{t^{\prime\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|^{2}}-\dfrac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|^{3}}\right)\left(\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{x}\right]^{4}+\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{y}\right]^{4}\right)\\\
D_{2}=\dfrac{1}{16}\sum\limits_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})}\left(\dfrac{t^{\prime\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|^{2}}-\dfrac{t^{\prime}(|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|)}{|\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB}|^{3}}\right)\left(\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{x}\right]^{4}-\left[(\boldsymbol{a}+\delta\boldsymbol{\tau}_{AB})_{y}\right]^{4}\right)\end{array}\right.$

Thus, we have the terms

$\displaystyle\sum_{j}\sum_{SS^{\prime}}\Psi^{\dagger}_{j,S}\left(\zeta_{SS^{\prime}}+\frac{1}{2}v_{F}\gamma\left(\boldsymbol{\nabla}\cdot\delta\boldsymbol{U}\right)\bar{\boldsymbol{\sigma}}_{SS^{\prime}}\cdot\mathcal{A}\right)\Psi_{j,S^{\prime}}$
(61)

intra-layer |  $\alpha$ (eV) |  $C_{1}$ (eV) |  $C_{2}$ (eV) |  $D_{1}$ (eV) |  $D_{2}$ (eV)
---|---|---|---|---|---
Ref. [66] | $-1.0542$ | $-0.1165$ | $0.1448$ | $-0.9829$ | $3.5837$
Ref. [67] | $2.5472$ | $0.4822$ | $-0.1546$ | $-1.0862$ | $-2.1330$
Table S1: Numerical values of the parameters in $\delta H_{intra}$ in Eq.(30)
and (5) for two different microscopic tight binding models.

###  A.6 Deformation Potential

In this subsection, we consider how the strain can couple with the charge
density of the electrons on a single layer. In the previous subsections, we
have calculated the effective continuum Hamiltonian $H_{eff}^{\mathbf{K}}$,
that contains the term

$\displaystyle H_{dp,1}=\alpha\int{\rm d}^{2}\boldsymbol{x}\
\sum_{j,S}\pm\left(\boldsymbol{\nabla}\cdot\delta\boldsymbol{U}\right)\Psi_{j,S}^{\dagger}\Psi_{j,S}=2\alpha\int{\rm
d}^{2}\boldsymbol{x}\
\sum_{j,S}\left(\boldsymbol{\nabla}\cdot\boldsymbol{U}_{j}^{\parallel}\right)\Psi_{j,S}^{\dagger}\Psi_{j,S}$
(62)

with
$\alpha=\frac{1}{4}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}|\boldsymbol{a}|t^{\prime}(|a|)$.
In the above formula, $+$ and $-$ sign are form the top and bottom layers
respectively, and we have used the fact that
$\boldsymbol{\nabla}\cdot\boldsymbol{U}_{t}^{\parallel}=-\boldsymbol{\nabla}\cdot\boldsymbol{U}_{b}^{\parallel}=\frac{1}{2}\boldsymbol{\nabla}\cdot\delta\boldsymbol{U}$.
However, we have assumed that the all the hopping constants depend only on the
displacement, and neglect the possible dependence of the onsite energy on the
nearby atomic configurations. If the onsite energy $\epsilon$ depends on the
lengths of three nearest bonds, it is expected that this dependence should
also induce the coupling between the strain and the fermion density. In this
case, $\epsilon$ is a function of
$\left\\{|\boldsymbol{n}_{j,S}^{(i)}(\boldsymbol{x})|\right\\})$
($i=1,\cdots,3$) where $\boldsymbol{n}^{(i)}_{j,S}(\boldsymbol{x})$ is the
nearest neighbor bond on layer $j$, sublattice $S$, and at the position of
$\boldsymbol{x}$ [62]. As a result, the onsite energy can be expanded to the
linear order of the strain [62]:

$\displaystyle\epsilon\left(\left\\{|\boldsymbol{n}_{j,S}^{(\alpha)}(\boldsymbol{x})|\right\\}\right)\approx\epsilon_{0}+\kappa\boldsymbol{\nabla}\cdot\boldsymbol{U}_{j}^{\parallel}\quad\mbox{with}\quad\kappa=\frac{\sqrt{3}}{2}a\frac{\partial\epsilon}{\partial|\delta_{S}^{\alpha}|}$
(63)

This induces another term that couples the strain with the fermion density,

$\displaystyle H_{dp,2}^{\mathbf{K}}=\sum_{j,S}\int{\rm d}^{2}\boldsymbol{x}\
\left(\epsilon_{0}+\kappa\boldsymbol{\nabla}\cdot\boldsymbol{U}_{j}^{\parallel}\right)\Psi_{j,S}^{\dagger}\Psi_{j,S}\
.$ (64)

Now, combining Eq. 62 and 64, we obtain

$\displaystyle H_{dp}^{\mathbf{K}}=\sum_{j,S}\int{\rm d}^{2}\boldsymbol{x}\
\left(\epsilon_{0}+\alpha_{dp}\boldsymbol{\nabla}\cdot\boldsymbol{U}_{j}^{\parallel}\right)\Psi_{j,S}^{\dagger}\Psi_{j,S}\
.$ (65)

with

$\displaystyle\alpha_{dp}=\kappa+2\alpha=\frac{\sqrt{3}}{2}a\frac{\partial\epsilon}{\partial|\delta_{S}^{\alpha}|}+\frac{1}{2}\sum_{\boldsymbol{a}}e^{-i\boldsymbol{K}\cdot\boldsymbol{a}}|\boldsymbol{a}|t^{\prime}(|\boldsymbol{a}|)$
(66)

Since Ref. [66] and [67] does not provide enough information on how the onsite
energy depends on the nearby atomic configuration, we can calculate only the
value of $\alpha$, but not $\kappa$. The numerical value of $\kappa$ can be
obtained either by ab initio method, such as DFT, or calibrated by the
experimental measurements on $\alpha_{dp}$ with the formula
$\kappa=\alpha_{dp}-2\alpha$.

##  Appendix B Inter-layer Tunneling

In this section, we will present the detailed derivation and the symmetry
analysis of the inter-layer tunneling terms. Starting from the master formula
in Eq. 26, and keep only the first order gradient of the fermion fields, we
obtain

$\displaystyle H_{inter}$
$\displaystyle=A_{mlg}^{-1}\sum_{\boldsymbol{G}}\sum_{SS^{\prime}}e^{i\boldsymbol{G}\cdot\boldsymbol{\tau}_{SS^{\prime}}}e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot
2\boldsymbol{U}_{t,S}^{\parallel}(\boldsymbol{x})}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(\boldsymbol{y}+d_{0}\hat{z},\\{\boldsymbol{n}^{(\alpha)}_{t,S}(\boldsymbol{x}+\boldsymbol{y}/2)\\},\\{\boldsymbol{n}_{b,S^{\prime}}^{(\alpha)}(\boldsymbol{x}-\boldsymbol{y}/2)\\})\times$
$\displaystyle\
\left[\Psi^{\dagger}_{t,S}(\boldsymbol{x})\Psi_{b,S^{\prime}}(\boldsymbol{x})+\frac{\boldsymbol{y}}{2}\cdot\left(\boldsymbol{\nabla}\Psi^{\dagger}_{t,S}(\boldsymbol{x})\Psi_{b,S^{\prime}}(\boldsymbol{x})-\Psi^{\dagger}_{t,S}(\boldsymbol{x})\boldsymbol{\nabla}\Psi_{b,S^{\prime}}(\boldsymbol{x})\right)\right]+h.c.$
(67)

where we have applied the formula that the lattice distortion
$\boldsymbol{U}^{\parallel}_{t}=-\boldsymbol{U}^{\parallel}_{b}$ and are
independent of the sublattice. In addition, the microscopic inter-layer
hopping function are assumed to depend not only on the hopping displacement
$\boldsymbol{y}+d_{0}\hat{z}$, but also on the direction of the nearest
neighbor vectors $\boldsymbol{n}^{(\alpha)}_{t,S}$ and
$\boldsymbol{n}^{(\alpha)}_{b,S^{\prime}}$ on each layer [67]. Note that the
terms in Eq. 26 are expanded only to the first order of the gradients of
fermion fields, and therefore, the inter-layer tunneling can be written as

$\displaystyle H_{inter}=\int{\rm d}^{2}\boldsymbol{x}\
\Psi_{t,S}(\boldsymbol{x})\left[T_{SS^{\prime}}(\boldsymbol{x})+\frac{1}{2}\left\\{\boldsymbol{p},\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})\right\\}\right]\Psi_{b,S^{\prime}}(\boldsymbol{x})+h.c.$
(68)

In the rest of this section, we will study the properties of these two
tunneling fields $T_{SS^{\prime}}(\boldsymbol{x})$ and
$\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})$ and express them in terms
of the microscopic hopping functions.

Different from the intra-layer terms that depends only on the gradient of
$\boldsymbol{U}^{\parallel}$, the inter-layer tunnelings are the function of
$\boldsymbol{U}^{\parallel}$ that is not tiny compared with
$|\boldsymbol{a}|$. As a consequence, we do not expand the inter-layer
tunneling terms in the powers of $\delta\boldsymbol{U}$.

Figure S1: The vectors $\boldsymbol{q}_{\mu,l}$ in the first six shells.

As mentioned in the main text, in general,
$\boldsymbol{U}^{\parallel}_{t}(\boldsymbol{x})=\frac{\theta}{2}\hat{z}\times\boldsymbol{x}+\frac{1}{2}\delta\boldsymbol{U}(\boldsymbol{x})$,
leading to

$e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot
2\boldsymbol{U}^{\parallel}(\boldsymbol{x})}=e^{i\theta(\boldsymbol{G}+\boldsymbol{K})\cdot(\hat{z}\times\boldsymbol{x})}e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot\delta\boldsymbol{U}}=e^{i\boldsymbol{x}\cdot(\boldsymbol{q}_{1,1}+\boldsymbol{g}_{\boldsymbol{G}})}e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot\delta\boldsymbol{U}(\boldsymbol{x})}$

where $\boldsymbol{q}_{1,1}=-\theta z\times\boldsymbol{K}$ and
$\boldsymbol{g}_{\boldsymbol{G}}=-\theta\hat{z}\times\boldsymbol{G}$. For
small twist angle $\theta\ll 1$, it can be shown that
$\boldsymbol{g}_{\boldsymbol{G}}$ is a reciprocal lattice vector of the moire
superlattice. For $\boldsymbol{G}_{i}$ and $\boldsymbol{g}_{i}$ defined in
Fig. 2, it is easy to see that

$\boldsymbol{g}_{\boldsymbol{G}_{1}}=-(\boldsymbol{g}_{1}+\boldsymbol{g}_{2})\quad\mbox{and}\quad\boldsymbol{g}_{\boldsymbol{G}_{2}}=\boldsymbol{g}_{1}\quad\Longrightarrow\quad\boldsymbol{g}_{m_{1}\boldsymbol{G}_{1}+m_{2}\boldsymbol{G}_{2}}=(m_{2}-m_{1})\boldsymbol{g}_{1}-m_{2}\boldsymbol{g}_{2}$

In addition, the lattice relaxation
$\delta\boldsymbol{U}(\boldsymbol{x})=\delta\boldsymbol{U}(\boldsymbol{x}+\boldsymbol{L}_{i})$
($i=1$, $2$) is a periodic function, and so is
$e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot\delta\boldsymbol{U}}$. This suggests
that

$e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot\delta\boldsymbol{U}}=\sum_{\boldsymbol{g}}u_{\boldsymbol{g}}(\boldsymbol{G})e^{i\boldsymbol{g}\cdot\boldsymbol{x}}\quad\Longrightarrow\quad
e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot
2U^{\parallel}_{t}(\boldsymbol{x})}=e^{i\boldsymbol{x}\cdot\boldsymbol{q}_{(1,1)}}\sum_{\boldsymbol{g}}e^{i(\boldsymbol{g}_{\boldsymbol{G}}+\boldsymbol{g})\cdot\boldsymbol{x}}u_{\boldsymbol{g}}(\boldsymbol{G})$

where $\sum_{\boldsymbol{g}}$ sums over all the reciprocal lattice vector of
the moire superlattice. This gives the Fourier transformation of
$e^{i(\boldsymbol{G}+\boldsymbol{K})\cdot
2\boldsymbol{U}^{\parallel}_{t}(\boldsymbol{x})}$.

Now, we obtain the expression of the inter-layer tunneling terms:

$\displaystyle T_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\boldsymbol{\tau}_{SS^{\prime}}}e^{i\boldsymbol{q}_{(1,1)}\cdot\boldsymbol{x}}\sum_{\boldsymbol{g}}e^{i(\boldsymbol{g}_{\boldsymbol{G}}+\boldsymbol{g})\cdot\boldsymbol{x}}u_{\boldsymbol{g}}(\boldsymbol{G})\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(\boldsymbol{y}+d_{0}\hat{z},\\{\boldsymbol{n}^{(\alpha)}_{t,S}(\boldsymbol{x}+\boldsymbol{y}/2)\\},\\{\boldsymbol{n}_{b,S^{\prime}}^{(\alpha)}(\boldsymbol{x}-\boldsymbol{y}/2)\\})$
(69) $\displaystyle\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\boldsymbol{\tau}_{SS^{\prime}}}e^{i\boldsymbol{q}_{(1,1)}\cdot\boldsymbol{x}}\sum_{\boldsymbol{g}}e^{i(\boldsymbol{g}_{\boldsymbol{G}}+\boldsymbol{g})\cdot\boldsymbol{x}}u_{\boldsymbol{g}}(\boldsymbol{G})\int{\rm
d}^{2}\boldsymbol{y}\
(-i\boldsymbol{y})e^{-i(\boldsymbol{G}+\boldsymbol{K})\cdot\boldsymbol{y}}t(\boldsymbol{y}+d_{0}\hat{z},\\{\boldsymbol{n}^{(\alpha)}_{t,S}(\boldsymbol{x}+\boldsymbol{y}/2)\\},\\{\boldsymbol{n}_{b,S^{\prime}}^{(\alpha)}(\boldsymbol{x}-\boldsymbol{y}/2)\\})$
(70)

As shown in the next section, the Fourier transformation of the inter-layer
tunneling in both microscopic tight binding models decay fast as a function of
the momentum. Therefore, it is convenient to express and calculate the Fourier
transformation of the tunneling terms. Since the direction of the nearest bond
$\boldsymbol{n}_{j,S}$ is also periodic with the period of
$\boldsymbol{L}_{i}$, its’ Fourier transformation can be written as

$\displaystyle
T_{SS^{\prime}}(\boldsymbol{x})=\sum_{\mu,l}e^{i\boldsymbol{q}_{\mu,l}\cdot\boldsymbol{x}}T_{SS^{\prime}}^{(\mu,l)}\
,\qquad\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})=\sum_{\mu,l}e^{i\boldsymbol{q}_{\mu,l}\cdot\boldsymbol{x}}\boldsymbol{\Lambda}_{SS^{\prime}}^{(\mu,l)}$
(71)

where the vectors $\boldsymbol{q}_{(\mu,l)}$ are plotted in Fig. S1 for the
first six $\boldsymbol{q}$ shells ($\mu\leq 6$). The subscripts $\mu$ refers
to the shell of the vectors, and $l$ is used to distinguish different vectors
in the same shell. It is easy to see that
$\boldsymbol{q}_{\mu,l}-\boldsymbol{q}_{(1,1)}$ is a reciprocal lattice vector
of the moire superlattice.

Here, we consider the microscopic tight binding models with the inter-layer
hopping depending only on the hopping displacement, as in Ref. [66]. We
introduce the notation for the Fourier transformation of the hopping

$\tilde{t}_{d_{0}}(\boldsymbol{q})=A_{mlg}^{-1}\int{\rm d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{q}\cdot\boldsymbol{y}}t(\boldsymbol{y}+d_{0}\hat{z})\quad\Longrightarrow\quad
A_{mlg}^{-1}\int{\rm d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{q}\cdot\boldsymbol{y}}\boldsymbol{y}t(\boldsymbol{y}+d_{0}\hat{z})=i\boldsymbol{\nabla}_{\boldsymbol{q}}\tilde{t}_{d_{0}}(\boldsymbol{q})$

Therefore,

$\displaystyle T_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=e^{i\boldsymbol{q}_{(1,1)}\cdot\boldsymbol{x}}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\boldsymbol{\tau}_{SS^{\prime}}}\tilde{t}_{d_{0}}(\boldsymbol{G}+\boldsymbol{K})\sum_{\boldsymbol{g}}e^{i(\boldsymbol{g}_{\boldsymbol{G}}+\boldsymbol{g})\cdot\boldsymbol{x}}u_{\boldsymbol{g}}(\boldsymbol{G})$
(72) $\displaystyle\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=e^{i\boldsymbol{q}_{(1,1)}\cdot\boldsymbol{x}}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\boldsymbol{\tau}_{SS^{\prime}}}\left.\boldsymbol{\nabla}_{\boldsymbol{q}}\tilde{t}_{d_{0}}(\boldsymbol{q})\right|_{\boldsymbol{q}=\boldsymbol{G}+\boldsymbol{K}}\sum_{\boldsymbol{g}}e^{i(\boldsymbol{g}_{\boldsymbol{G}}+\boldsymbol{g})\cdot\boldsymbol{x}}u_{\boldsymbol{g}}(\boldsymbol{G})$
(73)

With the rigid twist only, $\delta\boldsymbol{U}=0$, and thus
$u_{\boldsymbol{g}}(\boldsymbol{G})=\delta_{\boldsymbol{g},0}$. If we only
focus on the first $\boldsymbol{q}$ shell ($\mu=1$), we see that

$\displaystyle T_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=e^{i\boldsymbol{q}_{(1,1)}\cdot
x}\sum_{\boldsymbol{G}}e^{i\boldsymbol{G}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\tilde{t}_{d_{0}}(\boldsymbol{G}+\boldsymbol{K})e^{i\boldsymbol{g}_{\boldsymbol{G}}\cdot\boldsymbol{x}}$
$\displaystyle\approx
e^{i\boldsymbol{q}_{(1,1)}\cdot\boldsymbol{x}}\tilde{t}(\boldsymbol{K})+e^{i\boldsymbol{q}_{1,2}\cdot\boldsymbol{x}}e^{-i\boldsymbol{G}_{1}\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\tilde{t}_{d_{0}}(-\boldsymbol{G}_{1}+\boldsymbol{K})+e^{i\boldsymbol{q}_{1,2}\cdot\boldsymbol{x}}e^{-i(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})\cdot\delta\boldsymbol{\tau}_{SS^{\prime}}}\tilde{t}_{d_{0}}(-\boldsymbol{G}_{1}-\boldsymbol{G}_{2}+\boldsymbol{K})$
(74)

Due to $C_{3}$ symmetry,
$\tilde{t}_{d_{0}}(\boldsymbol{K})=\tilde{t}_{d_{0}}(-\boldsymbol{G}_{1}+\boldsymbol{K})=\tilde{t}_{d_{0}}(-\boldsymbol{G}_{1}-\boldsymbol{G}_{2}+\boldsymbol{K})$
and are real by $C_{2x}$ symmetry. Thus, by calculating the component of
$w_{i}^{(1,l)}$, we see that

$w_{0}^{(1,l)}=w_{0}\ ,\quad w_{2}^{(1,1)}=0\ ,\quad w_{3}^{(1,l)}=0\ ,\quad
w_{1}^{(1,2)}=w_{1}^{(1,3)}=-\frac{1}{2}w_{1}^{(1,1)}\ ,\quad\quad
w_{2}^{(1,2)}=-w_{2}^{(1,3)}=-\frac{\sqrt{3}}{2}w_{0}\ .$

where $w_{0}=\tilde{t}_{d_{0}}(\boldsymbol{K})$. Thus, if we further neglect
the gradient coupling fields $\boldsymbol{\Lambda}_{SS^{\prime}}$, we fully
recover the BM model [61] in which $w_{0}=w_{1}$.

The values of $w_{i}^{(\mu,l)}$ and $\boldsymbol{\lambda}_{i}^{(\mu,l)}$ have
been listed in Table. 3 and S2 for the two microscopic tight binding models
proposed in Ref. [66] and [67]. Notice that we only listed the values of $w$s
and $\boldsymbol{\lambda}$s only for the first several $\boldsymbol{g}$
shells. The contributions from other shells are negligible.

In Sec. D.2, we will discuss the inter-layer tunnelings terms for another more
complicated microscopic model in which the inter-layer hopping depends not
only on the hopping displacement, but also on the direction of the nearest
neighbor bonds on the two layers.

| $w_{0}^{(4,1)}$ | $w_{1}^{(4,1)}$ | $w_{2}^{(4,1)}$ | $w_{3}^{(4,1)}$ | $w_{0}^{(5,1)}$ | $w_{1}^{(5,1)}$ | $w_{2}^{(5,1)}$ | $w_{3}^{(5,1)}$ | $w_{0}^{(6,1)}$ | $w_{1}^{(6,1)}$ | $w_{2}^{(6,1)}$ | $w_{3}^{(6,1)}$
---|---|---|---|---|---|---|---|---|---|---|---|---
| Ref. [66]
---
relaxed
$-0.39$ | $1.13$ | $-0.20$ | $0$ | $5.79$ | $6.36$ | $0$ | $0$ | $5.18$ | $-2.78$ | $-4.98$ | $0$
| Ref. [67]
---
relaxed
$0.16$ | $0.34$ | $0.10$ | $-0.07$ | negligible
| $\boldsymbol{\lambda}_{0}^{(4,1)}/a$ | $\boldsymbol{\lambda}_{1}^{(4,1)}/a$ | $\boldsymbol{\lambda}_{2}^{(4,1)}/a$ | $\boldsymbol{\lambda}_{3}^{(4,1)}/a$ | $\boldsymbol{\lambda}_{0}^{(5,1)}/a$ | $\boldsymbol{\lambda}_{1}^{(5,1)}/a$ | $\boldsymbol{\lambda}_{2}^{(5,1)}/a$ | $\boldsymbol{\lambda}_{3}^{(5,1)}/a$ | $\boldsymbol{\lambda}_{0}^{(6,1)}/a$ | $\boldsymbol{\lambda}_{1}^{(6,1)}/a$ | $\boldsymbol{\lambda}_{2}^{(6,1)}/a$ | $\boldsymbol{\lambda}_{3}^{(6,1)}/a$
| Ref. [66]
---
relaxed
$(-0.1,-0.2)$ | $(0.5,0.1)$ | $(-0.1,-0.7)$ | $(0,0)$ | $(-5.2,0)$ | $(-5.1,0)$ | $(0,-0.2)$ | $(0,0)$ | $(2.2,-4.1)$ | $(-1.1,2.1)$ | $(-2.1,3.3)$ | $(0,0)$
| Ref. [67]
---
relaxed
$(0.3,0.1)$ | $(0.3,-0.1)$ | $(-0.1,0)$ | $(0,0)$ | negligible
Table S2: Parameters of the inter-layer tunneling terms for two microscopic
tight binding models, in the presence of the lattice relaxation. $a$ is the
magnitude of the primitive lattice vector, and all numbers are in the unit of
meV.

##  Appendix C Symmetry Analysis

In this section, we consider the constraints on them by various symmetries.
Here, we focus on three different symmetry transformations,
$C_{2}\mathcal{T}$, $C_{3}$, and $C_{2x}$, and the fermion fields transform as

$\displaystyle C_{2}\mathcal{T}:\quad$
$\displaystyle\Psi_{j,S}(\boldsymbol{x})\longrightarrow\mathcal{K}\Psi_{j,\bar{S}}(-\boldsymbol{x})$
(75) $\displaystyle C_{3}:\quad$
$\displaystyle\Psi_{j,S}(\boldsymbol{x})\longrightarrow
e^{-i\frac{2\pi}{3}(\sigma_{z})_{SS}}\Psi_{j,S}\left(R\left(-\frac{2\pi}{3}\right)\boldsymbol{x}\right)$
(76) $\displaystyle C_{2x}:\quad$
$\displaystyle\Psi_{j,S}(\boldsymbol{x})\longrightarrow\Psi_{\bar{j},\bar{S}}(m_{y}\boldsymbol{x})\
,$ (77)

where $\bar{j}$ and $\bar{S}$ are the layer and sublattice index different
from $j$ and $S$ respectively. $R\left(-\frac{2\pi}{3}\right)\boldsymbol{x}$
is the vector $\boldsymbol{x}$ rotated clockwisely by the angle of $2\pi/3$,
i.e.

$R\left(-\frac{2\pi}{3}\right)(x,y)^{T}=\left(-\frac{x}{2}+\frac{\sqrt{3}}{2}y\
,-\frac{y}{2}-\frac{\sqrt{3}}{2}x\right)^{T}$

and $m_{y}$ is the reflection symmetry through $xz$ plane, i.e.
$m_{y}(x,y)^{T}=(x,-y)^{T}$.

The inter-layer tunneling terms are invariant under $C_{2}\mathcal{T}$
transformation. It leads to the constraints

$\displaystyle\left(\sigma_{1}T(\boldsymbol{x})\sigma_{1}\right)_{SS^{\prime}}=\left(T_{SS^{\prime}}(-\boldsymbol{x})\right)^{*}\quad\left(\sigma_{1}\boldsymbol{\Lambda}(\boldsymbol{x})\sigma_{1}\right)_{SS^{\prime}}=\left(\boldsymbol{\Lambda}_{SS^{\prime}}(-\boldsymbol{x})\right)^{*}$
(78) $\displaystyle\Longrightarrow$
$\displaystyle\left(\sigma_{1}T^{(\mu,l)}\sigma_{1}\right)_{SS^{\prime}}=\left(T_{SS^{\prime}}^{(\mu,l)}\right)^{*}\
,\quad\left(\sigma_{1}\boldsymbol{\Lambda}^{(\mu,l)}\sigma_{1}\right)_{SS^{\prime}}=\left(\boldsymbol{\Lambda}_{SS^{\prime}}^{(\mu,l)}\right)^{*}$
(79)

This suggests that we can write

$\displaystyle
T_{SS^{\prime}}^{(\mu,l)}=\left(w_{0}^{(\mu,l)}\sigma_{0}+w_{1}^{(\mu,l)}\sigma_{1}+w_{2}^{(\mu,l)}\sigma_{2}+iw_{3}^{(\mu,l)}\sigma_{3}\right)_{SS^{\prime}}$
(80)
$\displaystyle\boldsymbol{\Lambda}_{SS^{\prime}}^{(\mu,l)}=\left(\boldsymbol{\lambda}_{0}^{(\mu,l)}\sigma_{0}+\boldsymbol{\lambda}_{1}^{(\mu,l)}\sigma_{1}+\boldsymbol{\lambda}_{2}^{(\mu,l)}\sigma_{2}+i\boldsymbol{\lambda}_{3}^{(\mu,l)}\sigma_{3}\right)_{SS^{\prime}}$
(81)

where all $w_{i}^{(\mu,l)}$ ($i=0$, $\cdots$, $4$) are real numbers and
$\boldsymbol{\lambda}^{(\mu,l)}_{i}$ are two-component real vectors.

The invariance under $C_{3}$ transformation leads to the following constraints
on $T_{SS^{\prime}}$ and $\boldsymbol{\Lambda}_{SS^{\prime}}$:

$\displaystyle T_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\left(e^{i\frac{2\pi}{3}\sigma_{z}}T\left(R\left(\frac{2\pi}{3}\right)\boldsymbol{x}\right)e^{-i\frac{2\pi}{3}\sigma_{z}}\right)_{SS^{\prime}}$
$\displaystyle\Lambda_{\alpha,SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\left(e^{i\frac{2\pi}{3}\sigma_{z}}R_{\alpha\beta}\left(-\frac{2\pi}{3}\right)\Lambda_{\beta}\left(R\left(\frac{2\pi}{3}\right)\boldsymbol{x}\right)e^{-i\frac{2\pi}{3}\sigma_{z}}\right)_{SS^{\prime}}\
.$ (82)

For notational convenience, we introduce
$\boldsymbol{q}_{(\mu,l^{\prime})}=R(2\pi/3)\boldsymbol{q}_{(\mu,l)}$. Then,

$\displaystyle
T^{(\mu,l)}_{SS^{\prime}}=\left(e^{i\frac{2\pi}{3}\sigma_{z}}T^{(\mu,l^{\prime})}e^{-i\frac{2\pi}{3}\sigma_{z}}\right)_{SS^{\prime}}\quad,\quad\Lambda_{\alpha,SS^{\prime}}^{(\mu,l)}=R_{\alpha\beta}\left(-\frac{2\pi}{3}\right)\left(e^{i\frac{2\pi}{3}\sigma_{z}}\Lambda_{\beta}^{(\mu,l^{\prime})}e^{-i\frac{2\pi}{3}\sigma_{z}}\right)_{SS^{\prime}}\
.$ (83)

Eventually, the contact and gradient coupling constants should satisfy the
following relations:

$\displaystyle w_{0}^{(\mu,l^{\prime})}=w_{0}^{(\mu,l)}\ ,\qquad
w_{3}^{(\mu,l^{\prime})}=w_{3}^{(\mu,l)}\ ,\qquad
w_{1}^{(\mu,l^{\prime})}-iw_{2}^{(\mu,l^{\prime})}=e^{i\frac{2\pi}{3}}\left(w_{1}^{(\mu,l)}-iw_{2}^{(\mu,l)}\right)$
$\displaystyle(\boldsymbol{\lambda}_{0}^{(\mu,l^{\prime})})_{\alpha}=\left(R\left(\frac{2\pi}{3}\right)\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{0}^{(\mu,l)}\right)_{\beta}\
,\qquad(\boldsymbol{\lambda}_{3}^{(\mu,l^{\prime})})_{\alpha}=\left(R\left(\frac{2\pi}{3}\right)\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{3}^{(\mu,l)}\right)_{\beta}$
$\displaystyle(\boldsymbol{\lambda}_{1}^{(\mu,l^{\prime})}-i\boldsymbol{\lambda}_{2}^{(\mu,l^{\prime})})_{\alpha}=e^{i\frac{2\pi}{3}}\left(R\left(\frac{2\pi}{3}\right)\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{1}^{(\mu,l)}-i\boldsymbol{\lambda}_{2}^{(\mu,l^{\prime})}\right)_{\beta}\
.$ (84)

Lastly, we consider $C_{2x}$ symmetry. It impose the constraints:

$\displaystyle T_{SS^{\prime}}(\boldsymbol{x})$
$\displaystyle=\left(\sigma_{x}T(m_{y}\boldsymbol{x})\sigma_{x}\right)_{S^{\prime}S}^{*}\
,\qquad\boldsymbol{\Lambda}_{\alpha,SS^{\prime}}(\boldsymbol{x})=(\tau_{3})_{\alpha\beta}\left(\sigma_{x}\boldsymbol{\Lambda}_{\beta}(m_{y}\boldsymbol{x})\sigma_{x}\right)_{S^{\prime}S}^{*}\
.$ (85)

Again, we introduce the notation
$\boldsymbol{q}_{(\mu,n)}=-m_{y}\boldsymbol{q}_{(\mu,l)}$. It leads to

$\displaystyle w_{0}^{(\mu,n)}=w_{0}^{(\mu,l)}\ ,\qquad
w_{1}^{(\mu,n)}=w_{1}^{(\mu,l)}\ ,\qquad w_{2}^{(\mu,n)}=-w_{2}^{(\mu,l)}\
,\qquad w_{3}^{(\mu,n)}=w_{3}^{(\mu,l)}$
$\displaystyle\left(\boldsymbol{\lambda}_{0}^{(\mu,n)}\right)_{\alpha}=\left(\tau_{3}\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{0}^{(\mu,l)}\right)_{\beta}\
,\quad\left(\boldsymbol{\lambda}_{1}^{(\mu,n)}\right)_{\alpha}=\left(\tau_{3}\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{1}^{(\mu,l)}\right)_{\beta}\
,\quad\left(\boldsymbol{\lambda}_{2}^{(\mu,n)}\right)_{\alpha}=-\left(\tau_{3}\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{2}^{(\mu,l)}\right)_{\beta}$
$\displaystyle\left(\boldsymbol{\lambda}_{3}^{(\mu,n)}\right)_{\alpha}=\left(\tau_{3}\right)_{\alpha\beta}\left(\boldsymbol{\lambda}_{3}^{(\mu,l)}\right)_{\beta}\
.$ (86)

The formulas above are the general constraints for the inter-layer tunneling
fields $T$ and $\boldsymbol{\Lambda}$. As an example, we explicitly write down
the formula here for $\boldsymbol{q}$s in the innermost shell $\mu=1$. It is
obvious that $n=1$ for $l=1$. Combining Eq. 84 and 86,

$\displaystyle w_{0}^{(1,l)}=w_{0}\ ,$ $\displaystyle w_{3}^{(1,l)}$
$\displaystyle=w_{3}\ ,$ $\displaystyle w_{2}^{(1,1)}$ $\displaystyle=0$ (87)
$\displaystyle w_{1}=w_{1}^{(1,1)}\ ,$ $\displaystyle w_{1}^{(1,2)}$
$\displaystyle=w_{1}^{(1,3)}=-\frac{w_{1}}{2}\ ,$ $\displaystyle
w_{2}^{(1,2)}$ $\displaystyle=-w_{2}^{(1,3)}=-\frac{\sqrt{3}}{2}w_{1}$ (88)
$\displaystyle\boldsymbol{\lambda}_{0}^{(1,1)}=(\lambda_{0,x},0)\ ,$
$\displaystyle\boldsymbol{\lambda}_{1}^{(1,1)}$
$\displaystyle=(\lambda_{1,x},0)\ ,$
$\displaystyle\boldsymbol{\lambda}_{3}^{(1,1)}$
$\displaystyle=(\lambda_{3,x},0)\ ,$
$\displaystyle\boldsymbol{\lambda}_{2}^{(1,1)}=(0,\lambda_{2,y})\ .$ (89)

##  Appendix D Microscopic Hopping Function

###  D.1 Slater-Koster like hopping parameterization

First, we consider the microscopic tight binding model proposed in Ref. [66],
in which

$\displaystyle t(\boldsymbol{d})$
$\displaystyle=V_{pp\pi}^{0}e^{-\frac{|\boldsymbol{d}|-a_{0}}{\delta}}\left[1-\left(\frac{\boldsymbol{d}\cdot\hat{z}}{|\boldsymbol{d}|}\right)^{2}\right]+V_{pp\sigma}^{0}e^{-\frac{|\boldsymbol{d}|-d_{0}}{\delta}}\left(\frac{\boldsymbol{d}\cdot\hat{z}}{|\boldsymbol{d}|}\right)^{2}$
(90)

where $V^{0}_{pp\pi}=-2.7$eV, $V^{0}_{pp\sigma}=0.48$eV. $a_{0}=0.142$nm is
the distance between the two nearest-neighbor carbon atoms on the same layer;
$d_{0}=0.335$nm is the inter-layer distance. The decay length for the hopping
is $\delta=0.319a_{0}$. The intra-layer hopping thus can be expressed as

$\displaystyle
t_{intra}(\boldsymbol{d})=V_{pp\pi}^{0}e^{-\frac{|\boldsymbol{d}|-a_{0}}{\delta}}$
(91)

where $\boldsymbol{d}$ is the in-plane hopping displacement.It is easy to show
that the Fourier transformation of the intra-layer hopping is

$\displaystyle\tilde{t}_{intra}(\boldsymbol{q})=A_{mlg}^{-1}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{y}\cdot\boldsymbol{q}}t(\boldsymbol{y})=V_{pp\pi}^{0}\frac{2\pi\delta^{2}}{A_{mlg}}e^{a_{0}/\delta}\left[1+(q\delta)^{2}\right]^{-\frac{3}{2}}$
(92)

It is now clear that the $\tilde{t}_{intra}$ decays in $q^{-3}$ in the
momentum space, but the intra-layer hopping $t$ exponentially decays in real
space. The Fourier transformations $\tilde{t}$ of the inter-layer hoppings are

$\displaystyle\tilde{t}_{d_{0}}(|\boldsymbol{q}|)=\frac{1}{\mathcal{A}_{mlg}}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{y}\cdot\boldsymbol{q}}t(\boldsymbol{y}+d_{0}\hat{z})$
$\displaystyle=$ $\displaystyle\frac{2\pi
d_{0}^{2}}{\mathcal{A}_{mlg}}\left[V_{pp\pi}^{0}\int_{0}^{\infty}{\rm d}y\
yJ_{0}(qd_{0}y)e^{-\frac{d_{0}}{\delta}\left(\sqrt{y^{2}+1}-a_{0}/d_{0}\right)}\frac{y^{2}}{y^{2}+1}+V_{pp\sigma}^{0}\int_{0}^{\infty}{\rm
d}y\
yJ_{0}(qd_{0}y)e^{-\frac{d_{0}}{\delta}\left(\sqrt{y^{2}+1}-1\right)}\frac{1}{y^{2}+1}\right]$
(93)

With the Fourier transformation in Eq. 93, we are able to calculate all the
components $w_{i}^{(\mu,l)}$ and $\boldsymbol{\lambda}_{i}^{(\mu,l)}$ of the
inter-layer tunneling fields $T_{SS^{\prime}}(\boldsymbol{x})$ and
$\boldsymbol{\Lambda}_{SS^{\prime}}(\boldsymbol{x})$, with their values listed
in Tab. S2.

###  D.2 Wannier based hopping parameterization

In this subsection, we consider the model proposed in Ref. [67]. Note that
Ref. [67] does not provide a general formula for the intra-layer hoppings, but
lists its magnitude for a set of discrete hopping distance. Here, we fit the
intra-layer hopping with the following formula:

$\displaystyle
t_{intra}(\boldsymbol{r})=t_{0}e^{-\alpha_{0}\bar{r}^{2}}\cos(\beta_{0}\bar{r})+t_{1}\bar{r}^{2}e^{-\alpha_{1}(\bar{r}-r_{1})^{2}}$
(94)

The values of the fitted parameters are listed in Table. S3. It is obvious
that the intra-layer hopping decays exponentially as a function of
$|\boldsymbol{r}|$.

According to Ref. [67], the general form of the inter-layer hopping can be
written as

$\displaystyle t_{inter}(\boldsymbol{r})=$ $\displaystyle
V_{0}(r)+V_{3}(r)\left(\cos(3\theta_{12})+\cos(3\theta_{21})\right)+V_{6}(r)\left(\cos(6\theta_{12})+\cos(6\theta_{21})\right)$
(95) $\displaystyle V(r)=$
$\displaystyle\lambda_{0}e^{-\xi_{0}\bar{r}^{2}}\cos(\kappa_{0}\bar{r})$ (96)
$\displaystyle V_{3}(r)=$
$\displaystyle\lambda_{3}\bar{r}^{2}e^{-\xi_{3}(\bar{r}-x_{3})^{2}}$ (97)
$\displaystyle V_{6}(r)=$
$\displaystyle\lambda_{6}e^{-\xi_{6}(\bar{r}-x_{3})^{2}}\sin(\kappa_{6}\bar{r})$
(98)

where the vector $\boldsymbol{r}$ is the in-plane projected vector of the
hopping displacement and $\bar{r}=r/a$. The variables $\theta_{12}$ and
$\theta_{21}$ are the angles between $\boldsymbol{r}$ and the nearest neighbor
bond vectors $\boldsymbol{n}_{j,S}$ on two layers, i.e.

$\displaystyle\theta_{12}=\cos^{-1}\left(-\frac{\boldsymbol{r}\cdot\boldsymbol{n}_{j,S}}{r|\boldsymbol{n}_{j,S}|}\right)=\theta_{\boldsymbol{r}}-\theta_{j,S}+\pi,$
(99)
$\displaystyle\theta_{21}=\cos^{-1}\left(\frac{\boldsymbol{r}\cdot\boldsymbol{n}_{j^{\prime},S^{\prime}}}{r|\boldsymbol{n}_{j^{\prime},S^{\prime}}|}\right)=\theta_{\boldsymbol{r}}-\theta_{j^{\prime},S^{\prime}}\
.$ (100)

In the above formula, we define $\theta_{\boldsymbol{r}}$ to be the angle
between the vector $\boldsymbol{r}$ and the $x$ axis, and $\theta_{j,S}$
($\theta_{j^{\prime},S^{\prime}}$) to be the angle between the bond vector
$\boldsymbol{n}_{j,S}$ ($\boldsymbol{n}_{j^{\prime},S^{\prime}}$) and the $x$
axis. Note that each carbon atom has three different bond vectors. We label
the angles of these three bonds by $\theta_{j,S}^{(\alpha)}$
($\theta_{j^{\prime},S^{\prime}}^{(\alpha)}$), where the superscript $\alpha$
is the index of the bond vectors. Without the lattice distortion (e.g. as for
a rigid twist), the three in-plane nearest neighbors of a carbon atom are
$C_{3}$ symmetric about the carbon atom, and
$\theta^{(\alpha)}_{j,S}=\theta_{j,S}^{(1)}+2\pi(\alpha-1)/3$. Therefore, the
angles $\theta_{12}$ and $\theta_{21}$ could differ by $2\pi/3$ if choosing a
different nearest neighbor bond, leading to the same $\cos(3m\theta_{12})$ and
$\cos(3m\theta_{21})$ with $m$ being an integer. As a consequence, Eq. 95 is
independent of the choice of the bond vectors. The values of all microscopic
hopping parameters are listed in Table. S3.

Inter |  $\lambda_{i}$ (eV) | $\xi_{i}$ | $x_{i}$ | $\kappa_{i}$ |  |
---|---|---|---|---|---|---
$V_{0}$ | $0.3155$ | $1.7543$ |  | $2.0010$ |  |
$V_{3}$ | $-0.0688$ | $3.4692$ | $0.5212$ |  |  |
$V_{6}$ | $-0.0083$ | $2.8764$ | $1.5206$ | $1.5731$ |  |
Intra |  $t_{0}$ (eV) | $\alpha_{0}$ | $\beta_{0}$ |  $t_{1}$ (eV) | $\alpha_{1}$ | $r_{1}$
| $-18.4295$ | $1.2771$ | $2.3934$ | $-3.7183$ | $6.2194$ | $0.9071$
Table S3: Parameters for the inter-layer and intra-layer hoppings in the model
proposed in Ref. [67].

In the presence of the lattice distortion, however, the local $C_{3}$ symmetry
is broken, and thus Eq. 95 depends on the choice of the bond vectors. In this
case, we set the inter-layer hopping as

$\displaystyle t_{inter}(\boldsymbol{r})=$ $\displaystyle
V_{0}(r)+V_{3}(r)\left(\frac{1}{3}\sum_{\alpha=1}^{3}\cos(3\theta_{12}^{(\alpha)})+\frac{1}{3}\sum_{\alpha=1}^{3}\cos(3\theta_{21}^{(\alpha)})\right)+V_{6}(r)\left(\frac{1}{3}\sum_{\alpha=1}^{3}\cos(6\theta_{12}^{(\alpha)})+\frac{1}{3}\sum_{\alpha=1}^{3}\cos(6\theta_{21}^{(\alpha)})\right)$
(101)

where
$\theta_{12}^{(\alpha)}=\theta_{\boldsymbol{r}}-\theta^{(\alpha)}_{j,S}+\pi$
and
$\theta_{21}^{(\alpha)}=\theta_{\boldsymbol{r}}-\theta^{(\alpha)}_{j^{\prime},S^{\prime}}$.

To obtain the Fourier transformation of the hopping function, consider

$\displaystyle\int{\rm d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{q}\cdot\boldsymbol{y}}t_{n}(|\boldsymbol{y}|)e^{in(\theta_{\boldsymbol{y}}-\theta_{j})}=\int_{0}^{\infty}{\rm
d}y\ yt_{n}(y)\int_{0}^{2\pi}{\rm d}\theta
e^{-iqy\cos(\theta-\theta_{q})}e^{in(\theta-\theta_{j})}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}{\rm d}y\ yt_{n}(y)\int_{0}^{2\pi}{\rm
d}\theta\
\sum_{m=-\infty}^{\infty}(-ie^{i(\theta-\theta_{q})})^{m}J_{m}(qy)e^{in(\theta-\theta_{j})}=2\pi\int_{0}^{\infty}{\rm
d}y\ yt_{n}(y)(-i)^{-n}J_{-n}(qy)e^{in(\theta_{q}-\theta_{j})}$
$\displaystyle=$ $\displaystyle
2\pi(-i)^{n}e^{in(\theta_{q}-\theta_{j})}\int_{0}^{\infty}{\rm d}y\
yt_{n}(y)J_{n}(qy)$ (102)

where $\theta_{q}$ is the angle between the vector $\boldsymbol{q}$ and the
$\hat{x}$ axis. In the derivation above, we have applied the formula
$e^{\frac{z}{2}(t-t^{-1})}=\sum_{m=-\infty}^{\infty}t^{n}J_{n}(z)$ and
$J_{-n}(z)=(-1)^{n}J_{n}(z)$. Starting from Eq. 102, it is easy to obtain

$\displaystyle\int{\rm d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{q}\cdot\boldsymbol{y}}t_{n}(|\boldsymbol{y}|)\cos(n(\theta_{\boldsymbol{y}}-\theta_{j}))=2\pi(-i)^{n}\cos(n(\theta_{q}-\theta_{j}))\int_{0}^{\infty}{\rm
d}y\ yJ_{n}(qy)t_{n}(y)$ (103)

For notation convenience, we introduce

$\displaystyle\tilde{V}_{i}(q)=\frac{2\pi}{A_{mlg}}\int_{0}^{\infty}{\rm d}y\
yV_{i}(y)J_{i}(qy)$ (104)

where $i=0$, $3$, and $6$.

We first consider the Fourier transformation when the lattice is locally
$C_{3}$ symmetric, i.e. the lattice relaxation is absent. In this case,

$\displaystyle\tilde{t}_{C_{3}}(\boldsymbol{q},\theta_{j,S},\theta_{j^{\prime},S^{\prime}})=\tilde{t}(\boldsymbol{q},\theta_{j,S},\theta_{j^{\prime},S^{\prime}})=A_{mlg}^{-1}\int{\rm
d}^{2}\boldsymbol{y}\
e^{-i\boldsymbol{q}\cdot\boldsymbol{y}}t(\boldsymbol{y}+d_{0}\hat{z},\theta_{j,S},\theta_{j^{\prime},S^{\prime}})$
$\displaystyle=$
$\displaystyle\tilde{V}_{0}(q)+i\left(-\cos(3(\theta_{q}-\theta_{j,S}))+\cos(3(\theta_{q}-\theta^{(\alpha)}_{j^{\prime},S^{\prime}}))\right)\tilde{V}_{3}(q)-\left(\cos(6(\theta_{q}-\theta_{j,S}))+\cos(6(\theta_{q}-\theta_{j^{\prime},S^{\prime}}))\right)\tilde{V}_{6}(q)$
(105)

where $\tilde{t}_{C_{3}}$ is introduced to refer to the Fourier transformation
of the inter-layer hopping function when the lattice is locally $C_{3}$
symmetric.

In the presence of the lattice distortion that varies slowly in the real
space, the angle $\theta_{j,S}^{(\alpha)}$ can be approximated as

$\displaystyle\theta_{j,S}^{(\alpha)}=\theta_{\boldsymbol{\delta}_{S}^{(\alpha)}}+\delta\theta^{(\alpha)}_{j,S},$
(106)
$\displaystyle\delta\theta_{j,S}^{(\alpha)}=\frac{(\hat{z}\times\boldsymbol{\delta}_{S}^{(\alpha)})\cdot\delta\boldsymbol{n}_{j,S}^{(\alpha)}}{|\boldsymbol{\delta}_{S}^{(\alpha)}|^{2}}=\frac{\epsilon_{\mu\nu}}{|\boldsymbol{\delta}_{S}^{(\alpha)}|^{2}}\delta_{S,\mu}^{(\alpha)}\frac{\partial
U_{j,\nu}^{\parallel}}{\partial x_{\rho}}\delta_{S,\rho}^{(\alpha)}.$ (107)

where
$\theta_{\boldsymbol{\delta}_{B}^{(\alpha)}}=\theta_{\boldsymbol{\delta}_{A}^{(\alpha)}}+\pi$,
and for our choice of the coordinate system,
$\theta_{\boldsymbol{\delta}_{A}^{(1)}}=\pi/6$,
$\theta_{\boldsymbol{\delta}_{A}^{(2)}}=\pi/6+2\pi/3$,
$\theta_{\boldsymbol{\delta}_{A}^{(3)}}=\pi/6-2\pi/3$. Since
$\boldsymbol{U}^{\parallel}_{t,S}=-\boldsymbol{U}^{\parallel}_{b,S}=\frac{1}{2}(\theta\hat{z}\times\boldsymbol{y}+\delta\boldsymbol{U})$,
we obtain

$\delta\theta_{t,S}^{(\alpha)}=-\delta\theta_{b,S}^{(\alpha)}=\frac{\theta}{2}+\frac{\epsilon_{\mu\nu}}{2|\boldsymbol{\delta}_{S}^{(\alpha)}|^{2}}\delta_{S,\mu}^{(\alpha)}\frac{\partial\delta\boldsymbol{U}_{\nu}}{\partial
x_{\rho}}\delta_{S,\rho}^{(\alpha)}\ .$

Consequently, the hopping function is approximated as

$\displaystyle
t(\boldsymbol{y}+d_{0}\hat{z},\\{\boldsymbol{n}^{(\alpha)}_{j,S}\\},\\{\boldsymbol{n}^{(\alpha)}_{j^{\prime},S^{\prime}}\\})\approx
t(\boldsymbol{y}+d_{0}\hat{z},\\{\theta_{\boldsymbol{\delta}_{S}^{(\alpha)}}\\},\\{\theta_{\boldsymbol{\delta}_{S^{\prime}}^{(\alpha)}}\\})+\sum_{\alpha}\frac{\partial
t}{\partial\theta^{(\alpha)}_{j,S}}\delta\theta^{(\alpha)}_{j,S}+\sum_{\alpha}\frac{\partial
t}{\partial\theta^{(\alpha)}_{j^{\prime},S^{\prime}}}\delta\theta^{(\alpha)}_{j^{\prime},S^{\prime}}$
$\displaystyle=$ $\displaystyle V_{0}(y)+V_{3}(y)\left(-\cos
3(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S}})+\cos
3(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S^{\prime}}})\right)+V_{6}(y)\left(\cos
6(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S}})+\cos
6(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S^{\prime}}})\right)$
$\displaystyle+V_{3}(y)\left(\left(\sum_{\alpha}\delta\theta_{j,S}^{(\alpha)}\right)\sin
3(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S}})-\left(\sum_{\alpha}\delta\theta_{j^{\prime},S^{\prime}}^{(\alpha)}\right)\sin
3(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S^{\prime}}})\right)$
$\displaystyle-2V_{6}(y)\left(\left(\sum_{\alpha}\delta\theta_{j,S}^{(\alpha)}\right)\sin
6(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S}})+\left(\sum_{\alpha}\delta\theta_{j^{\prime},S^{\prime}}^{(\alpha)}\right)\sin
6(\theta_{\boldsymbol{y}}-\theta_{\boldsymbol{\delta}_{S^{\prime}}})\right)$
(108)

Correspondingly, its Fourier transformation can be written as

$\displaystyle\tilde{t}(\boldsymbol{q},\\{\theta_{j,S}^{(\alpha)}\\},\\{\theta_{j^{\prime},S^{\prime}}^{(\alpha)}\\})\approx\tilde{t}_{C_{3}}(\boldsymbol{q},\theta_{\boldsymbol{\delta}_{S}},\theta_{\boldsymbol{\delta}_{S^{\prime}}})+\frac{1}{3}\left(\sum_{\alpha}\delta\theta_{j,S}^{(\alpha)}\right)\frac{\partial\tilde{t}_{C_{3}}(\boldsymbol{q},\theta_{\boldsymbol{\delta}_{S}},\theta_{\boldsymbol{\delta}_{S^{\prime}}})}{\partial\theta_{\boldsymbol{\delta}_{S}}}+\frac{1}{3}\left(\sum_{\alpha}\delta\theta_{j^{\prime},S^{\prime}}^{(\alpha)}\right)\frac{\partial\tilde{t}_{C_{3}}(\boldsymbol{q},\theta_{\boldsymbol{\delta}_{S}},\theta_{\boldsymbol{\delta}_{S^{\prime}}})}{\partial\theta_{\boldsymbol{\delta}_{S^{\prime}}}}$
(109)

where $\tilde{t}_{C_{3}}$, defined in Eq. 105, is the Fourier transformation
of the inter-layer hopping for locally $C_{3}$ symmetric lattice.

(a) (b) (c) (d)

Figure S2: Spectrum after truncating to a different number of
${\boldsymbol{q}}$ shells in the inter-layer tunneling terms. We have
considered the effective continuum Hamiltonian $H_{eff}^{\mathbf{K}}$
constructed for two microscopic models in (above) Ref. [66] and (below) Ref.
[67]. For both models with the rigid twist only (left), the inclusion of two
${\boldsymbol{q}}$ shells leads to an almost perfect agreement between the
spectrum produced by $H_{eff}^{\mathbf{K}}$ and the microscopic model. In the
presence of the lattice relaxation (right), however, more ${\boldsymbol{q}}$
shells are needed to achieve comparable accuracy. Note that while the tight
binding spectra automatically contain both valleys, for the continuum model we
show the spectra only for one valley in order to avoid clutter.

(a) (b) (c) (d)

Figure S3: Impact of the higher order terms in the effective Hamiltonian on
the energy spectrum near the CNP. The spectra of $H_{eff}^{\mathbf{K}}$
constructed for two microscopic models in (above) Ref. [66] and (below) [67]
produced by (blue) keeping or (green) dropping the second order gradient terms
in the intra-layer continuum Hamiltonian and gradient couplings in the inter-
layer tunneling terms. For comparison, the spectra of the microscopic tight
binding models are also plotted and marked as red. For both models, dropping
the higher order derivative terms leads to a mismatch of $\sim 5$-$10$meV,
that is consistent with the estimate of this energy scale in the main text.
Note that while the tight binding spectra automatically contain both valleys,
for the continuum model we show the spectra only for one valley in order to
avoid clutter.

(a) (b) (c)

Figure S4: The two singular values, $\lambda_{1}(\boldsymbol{k})$ (left
column) and $\lambda_{2}(\boldsymbol{k})$ (right column), of the projected p-h
operator $\hat{\mathcal{P}}$ defined in Sec. IV by turning off different
inter-layer terms in $H_{eff}^{\mathbf{K}}$ for the microscopic tight binding
model in Ref. [67], and keeping all terms in the intra-layer part of
$H^{\mathbf{K}}_{eff}$ in Eq.31. (a) all the terms in Eq. 38 are kept in the
inter-layer tunnelings, (b) the contact couplings $w_{3}^{(\mu,l)}$ for all
the shells are set to $0$, and (c) both the contact $w_{3}^{(\mu,l)}$ and the
gradient couplings $\boldsymbol{\Lambda}_{SS^{\prime}}^{(\mu,l)}$ are set to
$0$. It is clear that the p-h asymmetry is dominated by the contribution from
$w_{3}^{(\mu,l)}$.

##  Appendix E Approximate Formula of $v_{Dirac}$

In this section, we derive the approximate formula of the $v_{Dirac}$ in the
presence of the pseudo-vector field that is induced by the lattice relaxation
proposed in Ref. [65]. Motivated by Fig. 5, we neglect
$\boldsymbol{\nabla}\varphi^{\mathcal{A}}$, the irrotational part of the
pseudo-vector field $\mathcal{A}$, and thus
$\mathcal{A}\approx\boldsymbol{\nabla}\times(\hat{z}\varepsilon^{\mathcal{A}})$.
Furthermore, as demonstrated in Table 1, the lattice relaxation
$\delta\boldsymbol{U}$ in Ref. [65] is dominated by the lowest harmonics.
Because $\varepsilon_{\mathcal{A}}$ is even and real, its Fourier components
must be pure imaginary and also odd. Considering $C_{3}$ symmetry, we can
assume that

$\displaystyle\varepsilon^{\mathcal{A}}(\boldsymbol{x})\approx
i\tilde{\varepsilon}^{\mathcal{A}}_{1}\sum_{j=1}^{3}\left(e^{i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}-e^{-i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}\right)$
(110)

where $\boldsymbol{g}_{3}=-(\boldsymbol{g}_{1}+\boldsymbol{g}_{2})$. Under
these approximation, we can obtain

$\displaystyle\mathcal{A}_{+}=\frac{1}{\tilde{\varepsilon}^{\mathcal{A}}_{1}}\left(\mathcal{A}_{1}+i\mathcal{A}_{2}\right)\approx
i|\boldsymbol{g}_{1}|\sum_{j=1}^{3}\omega^{j-1}\left(e^{i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}+e^{-i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}\right)$
(111)

where $\boldsymbol{g}_{3}=-(\boldsymbol{g}_{1}+\boldsymbol{g}_{2})$ and
$\omega=e^{i2\pi/3}$. Now, for the equation

$\displaystyle
v_{F}\begin{pmatrix}p_{+}+\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}\mathcal{A}_{+}&\alpha
U(\boldsymbol{x})\\\ \alpha
U(-\boldsymbol{x})&p_{+}-\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}\mathcal{A}_{+}\end{pmatrix}\begin{pmatrix}\Phi_{1}(\boldsymbol{x})\\\
\Phi_{2}(\boldsymbol{x})\end{pmatrix}=0$ (112)

where $\alpha=w_{1}/(v_{F}k_{\theta})$, and
$U(\boldsymbol{x})=\sum_{j}\omega^{j-1}e^{i\boldsymbol{q}_{j}\cdot\boldsymbol{x}}$.
We further assume that $\alpha$ and
$\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}$ are small and expand the
wavefunction in terms of the powers of $\alpha$ and
$\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}$:

$\displaystyle\Phi_{1}(\boldsymbol{x})$
$\displaystyle=1+\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}\Phi_{1}^{(0,1)}(\boldsymbol{x})+\alpha^{2}\Phi_{1}^{(2.,0)}(\boldsymbol{x})+(\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1})^{2}\Phi_{1}^{(0,2)}(\boldsymbol{x})+\cdots$
(113) $\displaystyle\Phi_{2}(\boldsymbol{x})$
$\displaystyle=\alpha\Phi_{2}^{(1,0)}(\boldsymbol{x})+\alpha\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}\Phi_{2}^{(1,1)}(\boldsymbol{x})+\cdots$
(114)

where we have used the fact that $\Phi_{1}(\boldsymbol{x})$ is even in
$\alpha$ and $\Phi_{2}(\boldsymbol{x})$ is odd in $\alpha$. Substituting Eq.
113 and 114 into Eq. 112 and comparing the powers of $\alpha$ and
$\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}$, we obtain the equations

$\displaystyle p_{+}\Phi_{1}^{(0,1)}+\mathcal{A}_{+}=0\ ,$ $\displaystyle
p_{+}\Phi_{1}^{(2,0)}+U(\boldsymbol{x})\Phi_{2}^{(1,0)}=0\ ,$ $\displaystyle
p_{+}\Phi_{1}^{(0,2)}+\mathcal{A}_{+}\Phi_{1}^{(0,1)}=0$ (115) $\displaystyle
p_{+}\Phi_{2}^{(1,0)}+U(-\boldsymbol{x})=0\ ,$ $\displaystyle
p_{+}\Phi_{2}^{(1,1)}-\mathcal{A}_{+}\Phi_{2}^{(1,0)}+U(-\boldsymbol{x})\Phi_{1}^{(0,1)}=0\
.$ (116)

We consider the correction of $v_{Dirac}$ by the inclusion of $\mathcal{A}$ by
expanding it to the powers of
$O(\alpha^{2}\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1})$. This allows us to
focus only on the lowest harmonics of the Fourier expansion of $\Phi_{1}$ and
$\Phi_{2}$. After some calculations, we found

$\displaystyle\Phi_{1}^{(0,1)}$
$\displaystyle=-i\sum_{j=1}^{3}\left(e^{i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}-e^{-i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}\right)\
,$ $\displaystyle\Phi_{2}^{(1,0)}$
$\displaystyle=i\sum_{j}e^{-i\boldsymbol{q}_{j}\cdot\boldsymbol{x}}$ (117)
$\displaystyle\Phi_{1}^{(2,0)}$
$\displaystyle=\frac{-i}{\sqrt{3}}\sum_{j}\left(\omega
e^{i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}-\omega^{*}e^{-i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}\right)\
,$ $\displaystyle\Phi_{2}^{(1,1)}$
$\displaystyle\approx-2i\sqrt{3}\sum_{j}e^{-i\boldsymbol{q}_{j}\cdot\boldsymbol{x}}$
(118) $\displaystyle\Phi_{1}^{(0,2)}$
$\displaystyle=-\sum_{j}\left(e^{i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}+e^{-i\boldsymbol{g}_{j}\cdot\boldsymbol{x}}\right)\
.$ (119)

Now, Eq. 72 gives

$\displaystyle
v_{Dirac}\approx\frac{1-6(\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1})^{2}-3\alpha^{2}(1-\frac{14}{3}\sqrt{3}\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1})}{1+3\alpha^{2}+6(\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1})^{2}}\
.$ (120)

Based on Table. 1 and 2, $\gamma\tilde{\varepsilon}^{\mathcal{A}}_{1}\approx
0.06$. Eq. 120 gives $\alpha\approx 0.79$ when $v_{Dirac}=0$. This value is
very close to $\alpha=0.7857$, the numerical result obtained in Sec. V.


* Wiedermann et al. (2021) Wiedermann, M., J. F. Siegmund, J. F. Donges, and R. V. Donner, 2021: Differential Imprints of Distinct ENSO Flavors in Global Patterns of Very Low and High Seasonal Precipitation. Frontiers in Climate, 3.
  * Williams and Patricola (2018) Williams, I. N., and C. M. Patricola, 2018: Diversity of ENSO Events Unified by Convective Threshold Sea Surface Temperature: A Nonlinear ENSO Index. Geophysical Research Letters, 45 (17), 9236–9244, 10.1029/2018GL079203, URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2018GL079203, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2018GL079203.
  * Wills et al. (2018) Wills, R. C., T. Schneider, J. M. Wallace, D. S. Battisti, and D. L. Hartmann, 2018: Disentangling Global Warming, Multidecadal Variability, and El Niño in Pacific Temperatures. Geophysical Research Letters, 45 (5), 2487–2496, https://doi.org/10.1002/2017GL076327, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017GL076327, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017GL076327.
  * Wills et al. (2020) Wills, R. C. J., D. S. Battisti, K. C. Armour, T. Schneider, and C. Deser, 2020: Pattern Recognition Methods to Separate Forced Responses from Internal Variability in Climate Model Ensembles and Observations. Journal of Climate, 33 (20), 8693–8719, 10.1175/JCLI-D-19-0855.1, URL https://journals.ametsoc.org/view/journals/clim/33/20/jcliD190855.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Wolter and Timlin (2011) Wolter, K., and M. Timlin, 2011: El Niño/Southern Oscillation behaviour since 1871 as diagnosed in an extended multivariate ENSO index. Int. J. Climatol., 31, 1074–1087.
  * Wyrtki (1975) Wyrtki, K., 1975: El niño—the dynamic response of the equatorial pacific ocean to atmospheric forcing. 5 (4), 572–584, 10.1175/1520-0485(1975)005¡0572:ENTDRO¿2.0.CO;2, URL https://journals.ametsoc.org/view/journals/phoc/5/4/1520-0485˙1975˙005˙0572˙entdro˙2˙0˙co˙2.xml, publisher: American Meteorological Society Section: Journal of Physical Oceanography.
  * Yu et al. (2010) Yu, J.-Y., H.-Y. Kao, and T. Lee, 2010: Subtropics-Related Interannual Sea Surface Temperature Variability in the Central Equatorial Pacific. Journal of Climate, 23 (11), 2869–2884, 10.1175/2010JCLI3171.1, URL https://journals.ametsoc.org/view/journals/clim/23/11/2010jcli3171.1.xml, publisher: American Meteorological Society Section: Journal of Climate.
  * Zhang et al. (2014) Zhang, H., A. Clement, and P. D. Nezio, 2014: The South Pacific Meridional Mode: A Mechanism for ENSO-like Variability. Journal of Climate, 27 (2), 769–783, 10.1175/JCLI-D-13-00082.1, URL https://journals.ametsoc.org/view/journals/clim/27/2/jcli-d-13-00082.1.xml, publisher: American Meteorological Society Section: Journal of Climate.


# Analysing the Data-Driven Approach of Dynamically Estimating Positioning
Accuracy

Grigorios G. Anagnostopoulos, Alexandros Kalousis  Geneva School of Business
Administration, HES-SO Geneva, Switzerland
Email:
<EMAIL_ADDRESS>
###### Abstract

The primary expectation from positioning systems is for them to provide the
users with reliable estimates of their position. An additional information
that can greatly help the users utilize the position estimate is the level of
uncertainty that a positioning system assigns to the position estimate it
produced. The concept of dynamically estimating the accuracy of position
estimates of fingerprinting positioning systems has been sporadically
discussed over the last decade in the literature of the field, where mainly
handcrafted rules based on domain knowledge have been proposed. The emergence
of IoT devices and the proliferation of data from Low Power Wide Area Networks
(LPWANs) has facilitated the conceptualization of data-driven methods of
determining the estimated certainty over position estimates. In this work, we
analyze the data-driven based determination of the Dynamic Accuracy Estimation
(DAE), considering it in the broader context of a positioning system. More
specifically, with the use of a public LoRaWAN dataset, the current work
analyses: the repartition of the available training set between the tasks of
determining the location estimates and the DAE, the concept of selecting a
subset of the most reliable estimates, and the impact that the spatial
distribution of the data has to the accuracy of the DAE. The work provides a
wide overview of the data-driven approach of DAE determination in the context
of the overall design of a positioning system.

###### Keywords:

IoT, Fingerprinting, Error Estimation, LoRaWAN, LPWAN, Localization,
Positioning, Reproducibility, Machine Learning

†† Code: 10.5281/zenodo.4099079 / Data: 10.5281/zenodo.4117818
Preprint v1. Please check for updates

##  I Introduction

The interest over the utilization of IoT networks has been increasing over the
recent years. Moreover, the massive amounts of data that the IoT connectivity
produces, invite for intelligent ways of utilizing this massive volume of
information. In this direction, the emergence of Low Power Wide Area Networks
(LPWANs) has offered the grounds for energy efficient localization
applications. The data-driven methods of localization, which are a common
solution for indoor positioning systems, called fingerprinting methods, have
been recently applied in outdoor IoT settings [1, 2, 3, 4, 5]. Apart from the
production of location estimates, the positioning systems may also produce a
Dynamic Accuracy Estimation (DAE), which quantifies the certainty of the
system over its produced location estimate [6, 7].

In this work, we elaborate on the concept of Dynamic Accuracy Estimation, in a
data-driven approach. DAE can be perceived as the claimed accuracy of a
position estimate, according to the positioning system that produced said
estimate. As the name states, DAE is produced dynamically, in an on-line
manner, along with the production of the location estimates, without having
access to the ground truth location.

In order to calculate the DAE in a data-driven manner, a model is trained,
receiving the same feature set as the model that produces the position
estimates, while the targets of the two models are different. Instead of
training the model having as target the ground truth location, as the
positioning model does, the target this time is the distance between the
ground truth location and the location estimate produced by the positioning
model. Thus, this second model learns to predict the amount of the
localization error of the location estimates of the first model. In this way,
when the system processes a new reception of signals at an unknown location,
the two models will be able to produce: (i) a location estimate, provided by
the first (positioning) model, and (ii) an estimate of the localization error,
given by the second (DAE related) model.

The potential of the use of DAE depends on the particularities of the business
case in which a positioning system is meant to be applied. The main idea
behind the use of DAE is that it can be used either for offering to the user a
confidence level along with each position estimate, or for facilitating a
higher-level system in taking related decisions. Two examples of such decision
making could be the use of DAE in hand-off algorithms which could use the DAE
to switch among different technologies, or the option to present to the user
only a percentage of the most accurate position estimates, should that suit a
particular business case. For instance, in business cases where a very
frequent update of the location estimates is not required, the system could
select to return only a subset of the produced position estimates based on
their accuracy as estimated by the DAE. This concept is exemplified in one of
the tests presented in this work.

##  II Related Work

The concept of on-line error estimation of the location estimates that a
positioning system produces has been sporadically studied in the literature of
the field, through more than a decade. This concept of error prediction is far
from constituting a usually advertised feature of indoor positioning systems
and relevant publications, or a subject of comparison in the competitions of
the field [8] which mainly compare the accuracy of the location predictions.
Nevertheless, a framework of evaluating positioning systems, proposed in 2016
[9], lists the DAE as one of the metrics that the framework can evaluate.
Different approaches have been proposed regarding the way in which the
certainty over a produced location estimate can be estimated. The categorical
division of these approaches can be described by two main types of methods.

In the first category of the rule-based methods, which monopolized the
literature of the field for a long time, the aim is to hand-craft analytic or
heuristic methods of estimating the quality of location estimates. A variety
of rule-based methods has been proposed, such as using: the average geographic
distance between the nearest neighbors [10], improved by factoring in as a
weight the proximity of the nearest neighbors and also introducing the
location estimate [11], or by introducing a weighted average of likelihoods
instead of a simple average of euclidean distances [12]; a mechanism to
calculate a Dilution-of-Precision-like value [13]; the geographic distance
from the furthest neighbor to the location estimate [14]; the spatial
distribution of the latest position estimates [15, 16]; and an offline method
based on the Cramer-Rao Lower Bound ratio [17].

On the other hand, the data-driven methods, which constitute the second
category, that recently emerged, aim at learning to predict the quality of
location estimates based on a dataset used to train a machine learning model.
The current work focuses on this second, data-driven approach. A surprisingly
early proposal of a data-driven method of DAE determination dates back in
2007, when Dearman et al. [18] proposed a multiple linear regression approach.
That work uses designer-defined features based on the RSS fingerprints to
train a model which predicts the localization error of a location estimate
produced by a new signal reception.

A pair of recent works, by Lemic et al. [6] and Lemic and Famaey [7] has
revived the data-driven approach. The availability of large datasets of IoT-
based, LPWAN networks, such as the ones presented by Aernouts et al. [1], lays
a proper ground for such approaches. In 2019, Lemic et al. [6] presented an
extensive analysis of the performance of several well-known regression
algorithms (linear regression with regularizers, kNN, SVN, Random Forests) on
predicting the positioning error of location estimates. In a natural
continuation of that work [6], Lemic and Famaey [7] investigated the
capabilities of Neural Networks (NN) in addressing the same regression
problem. By utilizing the LPWAN datasets of [1], the authors show that the NN
approach outperforms kNN, which was the best performing method of [6].

In the current work, we further analyze the data-driven based determination of
Dynamic Accuracy Estimation (DAE), considering it in the context of a broader
positioning system. In this context, we focus on points such as the
repartition of the available training data between the models producing
location estimates and DAE estimates, the location estimate selection based on
the DAE, as well as the impact of the spatial distribution of the data to the
DAE.

##  III The Dataset Used

Aenrouts et al. [1] have made publicly available 3 fingerprinting datasets of
Low Power Wide Area Networks. Two of these datasets were collected in the
urban area of Antwerp, one using Sigfox and another using LoRaWAN. The third
dataset was collected in the rural area between the towns of Antwerp and
Ghent, using Sigfox. Since the initial publication, two major updates have
been published. In this work, the latest version (1.3) of the LoRaWAN dataset
has been utilized, which contains 130430 LoRaWAN messages in the urban area of
Antwerp. The fingerprints are collected in an area of approximately 53 square
kilometers, though the majority of them lay in the central area of Antwerp
which is approximately half the size.

Fingerprinting techniques are often compared to their counterpart, the ranging
techniques such as multilateration, which require a minimum of 3 receiving
gateways to produce a unique position estimate. Even if satisfying results can
be obtained with fingerprinting methods when using messages with fewer than 3
receiving gateways, in this work we limit the dataset to the messages
containing at least 3 receiving gateways. A total of 75054 messages with fewer
than 3 receiving gateways are dropped, while 55375 messages with at least 3
receiving gateways are kept to be used. A common train, validation, test set
split will be used for the experiments of this work. We will use 70% of the
dataset for training purposes, 15% for validation, and 15% as a test set for
reporting the models’ performance. In the following sections we will elaborate
on the way the training set is repartitioned among two distinct training
tasks: the production of position estimates and the production of dynamic
accuracy estimates.

##  IV Methodology and Terminology

In the setting used in the current work, each message reception obtained by
the network of LoRa gateways will be utilized to produce a position estimate,
accompanied with a degree of certainty that the system assigns to the position
estimate it produced. In this work, this second estimate will be referred to
as the Dynamic Accuracy Estimation (DAE), since it concerns an estimation over
the accuracy of produced location estimates, which is produced dynamically, as
the system operates. The relevant literature refers to DAE either as an
’accuracy’ [13, 14, 17] or as ’error’ [18, 10, 11, 6, 12] estimation, or even
as ’confidence’ [15]. The actual, accurate measurement of the error/accuracy
of a position estimate, can only become available when the ground truth
location is known and its distance from the location estimate produced by the
positioning system is measured.

For the positioning system to produce this twofold estimation, two data-driven
models are trained. The first model, $M_{1}$, is trained to predict the
locations from which the message transmissions were made by the mobile device,
based on the signals received by the LoRaWAN gateways. In training time, the
received signals (in this case the RSS values received at each gateway) are
the input features fed to the model, while the target variables are the
latitude and longitude of the ground truth location $Loc_{GT}$. In testing
time, a set of previously unseen message receptions are fed to the trained
model $M_{1}$ which produces location estimates $Loc_{Est}$. The error of each
produced estimate is defined as the geographic distance between $Loc_{GT}$ and
$Loc_{est}$. In this work, the Haversine distance formula has been used to
define this distance.

$Error_{pos}=Error_{M1}=D_{haversine}(Loc_{GT},Loc_{est})$ (1)

The predictive performance of model $M_{1}$ is evaluated by calculating the
errors of all produced location estimates, as by Equation 1, and drawing
statistics such as the mean, the median or specific percentile values of the
error distribution.

The second model, $M_{2}$, aims to predict the error $Error_{pos}$ that the
location estimate $Loc_{est}$ produced by $M_{1}$ will have, given only the
signals received by the gateways. Therefore, in training time, $M_{2}$
utilizes the same input features as $M_{1}$: the received signals of message
transmissions (in this case the RSS values received at each gateway), while
this time the target variable is the positioning error $Error_{pos}$ that
$M_{1}$’s position estimates have for the same message receptions. In testing
time, $M_{2}$ produces the Dynamic Accuracy Estimation, or $DAE$, that is a
prediction regarding the value of $Error_{pos}$ that $M_{1}$ can achieve for a
certain message reception.

In order to quantitatively characterize and evaluate the prediction of
$M_{2}$, the absolute difference between the target variable $Error_{pos}$ and
the predicted value $DAE$ is calculated for each prediction, and cumulative
statistical metrics such as the mean or the median are produced. Thus, the
following metric is commonly used:

$Error_{DAE}=Error_{M2}=\lvert Error_{pos}-DAE\rvert$ (2)

##  V Experimentation and Results

In this section, the performance of model $M_{2}$, and more particularly, the
repercussion that such a model can have on a positioning system as a whole, is
studied. Firstly, since a finite set of fingerprints will be available for the
training of the two models ($M_{1}$, $M_{2}$), we experiment with the effect
that the proportion of training data which are assigned to each of the two
models has on their performance (Subsection V-A). Following, in Subsection
V-B, we exemplify the models’ capabilities in a use-case where a percentage of
the most trustworthy position estimates are selected to be used. Lastly, in
Subsection V-C, we analyze the performance of $M_{2}$ in local spatial areas.

Throughout the experiments of this work, the ExtraTrees regression method has
been used for the training of both $M_{1}$ and $M_{2}$ models. Based on the
results of [3] and on preliminary tests performed for this work, the
ExtraTrees was preferred from the kNN used in other relevant works.
Nevertheless, a detailed comparison among regression methods is out of the
scope of this work.

###  V-A Training set repartition detween $M_{1}$ and $M_{2}$

In this first test, we study the performance of the two models for different
repartitions of the training set between the two training tasks. The overall
dataset has been divided into 70%, 15%, 15% for train, validation and test
purposes respectively. At this point, the further repartition of the training
set into the two training tasks is studied. The positioning task, served by
model $M_{1}$, is the principal task of the positioning system, while the
production of a DAE by a model $M_{2}$ is naturally a secondary task. Given
this fact, it is reasonable for a designer of such a system to assign more
data to the training of $M_{1}$. On the other hand, one may consider
particular use-cases in which it might be requested to prioritize a very high
certainty regarding the quality of position estimates over the average overall
performance of position estimate production. The overview of the performance
trade-off between $M_{1}$ and $M_{2}$ for different repartitions of the
training data is a suitable way of facilitating the designer of such systems
to take informed decisions. Figure 1 presents the performance of $M_{1}$ and
$M_{2}$ on the validation set, for different portions of the training set
assigned to the secondary task, dealt by $M_{2}$.

Figure 1: Performance of models $M_{1}$ and $M_{2}$ on the validation set, for
various portions of the training set assigned to the training of $M_{2}$. The
performance of $M_{1}$ and $M_{2}$ corresponds to the mean values of the
errors described in Equations 1 and 2 respectively, calculated on the
validation set.

In Figure 1, it can be observed that the repercussion of the volume of
training data is greater on the primary task of positioning dealt by $M_{1}$
than on the secondary task of DAE dealt by $M_{2}$. In model $M_{1}$, the
difference in the mean validation error between the case of using 90% of the
training data and the case of only using 10% is 53 meters, which constitutes a
20% increase of the error. On the other hand, regarding $M_{2}$, the reduction
of the mean error between the cases of using 90% and 10% of the training data
for its training is only 16 meters, or 11% of relative error increase. As
discussed previously, the selection of the way that the training data will be
allocated to the two tasks depends on the requirements of the application. For
the rest of the tests of this work we will proceed with 50%-50% sharing of
training data between $M_{1}$ and $M_{2}$, corresponding to 35%-35% of the
overall dataset. With this repartition, the validation error of $M_{1}$ has a
mean of 275 meters and a median of 171 meters, while $M_{2}$ has a 149-meter
mean error and an 85-meter median error. The respective test set error for
$M_{1}$ has a mean of 276 meters and a median 176 meters, while $M_{2}$ has a
147-meter mean error and an 83-meter median error.

###  V-B Location estimate selection based on DAE

Having trained the two models, we proceed to the exemplification of their
combined usage in a process of selecting a subset of position estimates, based
on the confidence that the system claims over them. To do so, the following
procedure is undertaken. Initially, both models provide their predictions on
the previously unseen validation set: $M_{1}$ provides the position estimates,
and $M_{2}$ the corresponding DAE. The outcomes for all data points of the
validation set are sorted with respect to the estimated DAE values, thus
starting from the position estimates in which $M_{2}$ assigns high certainty
and moving towards the position estimates which $M_{2}$ predicts as being very
erroneous. Figure 2 offers an overview of the positioning system’s performance
with respect to different portions of the validation set, selected based on
the DAE.

Figure 2:  The mean and median positioning error of portions (subsets) of the
validation set, selected according to their DAE value. The data points have
been sorted based on the DAE estimated over them.

Depending on how selective a higher-level system needs to be regarding the
accuracy that it requires from the positioning system, the threshold of
acceptable DAE values can be set accordingly. The designer of the system may
decide regarding the selection method to be used. One approach is to select a
certain percentage of the most accurate position estimates, according to the
DAE. Alternatively, a hard threshold of the acceptable DAE values can be set,
irrespectively of the distribution of the DAE value in the studied dataset.

It is important to note that if it is assumed that the selection process is
not done on-line for each message individually as the messages arrive
sequentially, but cumulatively upon the collection of a dataset of a certain
volume, then the distribution of the DAE values in this set can be considered
available. Consequently, in that case, the median or any other percentile
value of DAE of such an available set (as for instance, the test set of the
current study) can be used as the threshold for selecting the required
percentage of position estimates on which $M_{2}$ has assigned the highest
certainty. On the other hand, if the use-case dictates an on-line decision,
the corresponding distribution of the DAE values can be obtained by the
validation set, which is known.

The results depicted in Figure 2 suggest that $M_{2}$ manages to estimate the
DAE values in a relatively reliable way. In the validation set, in which the
positioning model achieves a 275-meter mean error and a 171-meter median
error, it is possible to learn to select signal readings of good quality,
which produce position estimates whose accuracy is significantly better that
the rest. For instance, selecting 50% of the estimates (those with the lowest
estimated error according to DAE), results in a mean positioning error of 108
meters with a median positioning of 35 meters. These values correspond to a
61% improvement of the mean error and an 80% improvement of the median error.

It is noteworthy that, for use-cases such as this one where the goal is the
distinction between more accurate and less accurate location estimates, the
absolute values of DAE, and consequently its absolute error values, as by
Equation 2, are not as crucial as the relative ordering of the location
estimates based on their DAE. In other words, assuming location estimates $i$
and $j$, with $Error_{M1}^{i}<Error_{M1}^{j}$, there are applications where
the correct ordering of the $M_{2}$ estimates, that is $DAE^{i}<DAE^{j}$, may
be more important than their absolute error values, $Error_{M2}^{i}$ and
$Error_{M2}^{j}$.

Figure 3: The spatial distribution of the data points of the dataset, color
coded into 20 clusters. The cluster ID is indicated in red at the cluster
centres. Gateway locations are depicted in black.

###  V-C Location estimate selection in local clusters

The spatial distribution of the location estimates and its relation to the DAE
assigned to them by the $M_{2}$ model is a very interesting subject. The real-
world dataset of Antwerp used in this work, as most relevant real-world
datasets, does not have a regular, fixed density of collected data throughout
different spatial areas. Using machine learning models implies that their
performance is dependent on the quantity of available training data. Moreover,
the particular setting of geo-localization further implies that the quantity
of data at different areas (the spatial density of training data) determines
the models’ performance at these areas. This is a consequence of the fact that
fingerprinting techniques learn the local particularities of the environment
through a mesh of fingerprints. A scarcity of data in a certain zone cannot be
compensated by a profusion of data at a distant area, where most likely a
distinct set of gateways will be involved in the signal reception, and thus
other features will be of importance in that area.

To study the impact of locality in the two models discussed in this work, the
dataset has been clustered into 20 areas, using the k-means algorithm,
operated on the ground truth locations of the data. In Figure 3, the 20
clusters are presented in a color code, having their ID number indicated at
the cluster centres.

It can be observed that the majority of data points has been collected at the
central area of Antwerp, at the right-hand side of the map. On the other hand,
clusters 5, 9 and 15 are characterized by a rather sparse collection of data.
As mentioned previously, the density of data collection at different areas may
affect the performance of the positioning system in those areas. Nevertheless,
there are other factors that affect the performance of the positioning system
as well. The complexity of the environment over which the signals are
propagated is a major factor. Factors such as the presence of obstacles (high
buildings, hills, etc.) as well as the density and the locations of the
gateways may greatly affect a system’s performance. For instance, in the city
center it is reasonable to expect to have a higher impact of the multipath
effect, than in suburban or rural areas.

Figure 4 highlights in blue the mean positioning error, for each of the 20
clusters. Furthermore, as a metric indicating the spatial density of data
within a cluster, the mean distance of all location pairs in a cluster is
indicated in red. The correlation coefficient of the two metrics of Figure 4
is 0.72.

Figure 4:  The $M_{1}$ error per cluster, along with the mean distance of all
location pairs per cluster.

Figure 5 presents the distribution of positioning error of each cluster.
Moreover, Figure 6 presents the distribution of positioning error after having
selected the 50% most accurate location estimates of each cluster, based on
their DAE, as discussed in Subsection V-B. In both Figures, the number of data
points of each cluster contributing to the statistics is depicted in blue. The
amount of points per cluster has an anticorrelation of coefficient $-0.63$
with the mean positioning error of the clusters in Figure 5 and of $-0.5$ in
Figure 6. These figures provide an insight on how the location estimate
selection process may perform when operated locally.

This is of great interest, since a selection of estimates from the whole
dataset may favour the overall error statistics, but under-represent certain
regions of high error, and over-represent regions of low error. Specifically,
when selecting the 50% of the most accurate (according to DAE) location
estimates of each cluster, the mean positioning error of all selected
locations is 186 meters and the median 111 meters, in comparison to the
respective 108-meter mean and 35-meter median error reported previously for
the case of selecting from the whole dataset. The significant difference,
especially in terms of the median error, is related to the uneven selection of
location estimates from different regions.

Figure 5: The distribution of the positioning error of $M_{1}$, for each
cluster. The mean value is indicated by the green triangle, while in blue the
amount of contributing data points per cluster are depicted. Figure 6:  The
distribution of the positioning error of $M_{1}$, for each cluster, for all
location estimates that have a bellow median DAE value of their cluster.

##  VI Conclusions and Future Work

In this work, we analyzed various aspects of the data-driven based
determination of Dynamic Accuracy Estimation (DAE), in the context of its
potential use in a broader system. In the previous interesting works that have
only recently elaborated on the data-driven approach of DAE determination [6,
7], the authors have utilized the great majority of the training set (87.6%)
for training the DAE model. Though this option allows the presentation of the
performance capabilities of the data-driven determination of the DAE in its
full extend, we considered imperative to investigate the impact that the
training data repartition between the two training tasks has to the overall
system’s performance. In the general case, it is natural to prioritize the
positioning task as the primary task of the system, assigning to it more
training data, rather than the DAE determination task. Moreover, the
experimentation of the current work indicated that when receiving more
training data, the positioning model has a great gain in error reduction (both
in absolute and percentage terms), than what the DAE determining model in the
same case. This work has exemplified, in Subsection V-A, the way to obtain an
overview of this trade-off, which allows informed decisions to be taken.

In Subsection V-B, we have elaborated on the fact that DAE, apart from its
evident usefulness as an additional information to the user, can have other
interesting applications when utilized by the overall system. A utilization of
DAE for the selection of the most accurate location estimates has been
exemplified, returning to the user a subset of estimates with a significantly
increased location accuracy. Such a process may be used in use-cases where
such an option would be appropriate. Moreover, Subsection V-C has studied the
performance of the two models in local, spatial areas, via clustering. In the
particular setting under study, it was observed that certain areas (cluster 3)
may have a high concentration of data, yielding very accurate predictions from
both models, in an unbalanced manner comparing to the rest of the dataset.
Secondly, it was observed that the suburban areas (clusters 5, 9 and 15) had
the highest error. These three clusters were also among the most sparsely
collected areas. Overall, we identified that the mean positioning error per
cluster had a significant correlation with the average distance among (or the
density of) the data points per cluster, and a significant anticorrelation
with the amount of data points per cluster. Lastly, we have highlighted the
fact that in datasets with such diverse density of data collected among
different spatial zones, it is interesting to look into performance statistics
from spatial zones along with the statistics of the overall dataset.

In the spirit of repeatability, reproducibility, verifiability and
comparability of results, it is important not only to report performance
metrics over a public dataset, but to also provide the means for reproducing
the same train/validation/test set split. Ideally, the code which
unambiguously presents the exact way the reported results were produced may be
published as well. In this way, the community is able to repeat, reproduce,
verify and consistently compare results, and is eventually enabled to take
informed decisions over the way positioning systems are designed, implemented
and deployed. For these reasons, the full code implementation of the current
work is openly available at the Zenodo repository.

## References

  * [1] M. Aernouts, R. Berkvens, K. Van Vlaenderen, and M. Weyn, “Sigfox and lorawan datasets for fingerprint localization in large urban and rural areas,” _Data_ , vol. 3, no. 2, 2018. [Online]. Available: http://www.mdpi.com/2306-5729/3/2/13
  * [2] T. Janssen, M. Aernouts, R. Berkvens, and M. Weyn, “Outdoor fingerprinting localization using sigfox,” in _2018 International Conference on Indoor Positioning and Indoor Navigation (IPIN)_ , Sep. 2018, pp. 1–6.
  * [3] G. G. Anagnostopoulos and A. Kalousis, “A reproducible comparison of rssi fingerprinting localization methods using lorawan,” in _2019 16th Workshop on Positioning, Navigation and Communications (WPNC)_ , 2019, pp. 1–6.
  * [4] ——, “A reproducible analysis of rssi fingerprinting for outdoor localization using sigfox: Preprocessing and hyperparameter tuning,” in _2019 International Conference on Indoor Positioning and Indoor Navigation (IPIN)_ , 2019, pp. 1–8.
  * [5] T. Janssen, R. Berkvens, and M. Weyn, “Benchmarking rss-based localization algorithms with lorawan,” _Internet of Things_ , vol. 11, p. 100235, 2020\. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S2542660520300688
  * [6] F. Lemic, V. Handziski, M. Aernouts, T. Janssen, R. Berkvens, A. Wolisz, and J. Famaey, “Regression-based estimation of individual errors in fingerprinting localization,” _IEEE Access_ , vol. 7, pp. 33 652–33 664, 2019.
  * [7] F. Lemic and J. Famaey, “Artificial neural network-based estimation of individual localization errors in fingerprinting,” in _2020 IEEE 17th Annual Consumer Communications Networking Conference (CCNC)_ , 2020, pp. 1–6.
  * [8] J. Torres-Sospedra, A. Jiménez, A. Moreira, T. Lungenstrass, W.-C. Lu, S. Knauth, G. Mendoza-Silva, F. Seco, A. Pérez-Navarro, M. Nicolau, and et al., “Off-line evaluation of mobile-centric indoor positioning systems: The experiences from the 2017 ipin competition,” _Sensors_ , vol. 18, no. 2, p. 487, Feb 2018. [Online]. Available: http://dx.doi.org/10.3390/s18020487
  * [9] C. M. de la Osa, G. G. Anagnostopoulos, M. Togneri, M. Deriaz, and D. Konstantas, “Positioning evaluation and ground truth definition for real life use cases,” in _2016 International Conference on Indoor Positioning and Indoor Navigation (IPIN)_ , Oct 2016, pp. 1–7.
  * [10] H. Lemelson, M. B. Kjærgaard, R. Hansen, and T. King, “Error estimation for indoor 802.11 location fingerprinting,” in _Location and Context Awareness_ , T. Choudhury, A. Quigley, T. Strang, and K. Suginuma, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 138–155.
  * [11] P. Marcus, M. Kessel, and M. Werner, “Dynamic nearest neighbors and online error estimation for smartpos,” _International Journal On Advances in Internet Technology_ , vol. 6, pp. 1–11, 01 2013.
  * [12] Y. Li, Z. He, Z. Gao, Y. Zhuang, C. Shi, and N. El-Sheimy, “Toward robust crowdsourcing-based localization: A fingerprinting accuracy indicator enhanced wireless/magnetic/inertial integration approach,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 3585–3600, April 2019.
  * [13] V. Moghtadaiee, A. G. Dempster, and B. Li, “Accuracy indicator for fingerprinting localization systems,” in _Proceedings of the 2012 IEEE/ION Position, Location and Navigation Symposium_ , 2012, pp. 1204–1208.
  * [14] D. Zou, W. Meng, and S. Han, “An accuracy estimation algorithm for fingerprint positioning system,” in _2014 Fourth International Conference on Instrumentation and Measurement, Computer, Communication and Control_ , 2014, pp. 573–577.
  * [15] R. Elbakly and M. Youssef, “CONE: zero-calibration accurate confidence estimation for indoor localization systems,” _CoRR_ , vol. abs/1610.02274, 2016. [Online]. Available: http://arxiv.org/abs/1610.02274
  * [16] S. Khandker, R. Mondal, and T. Ristaniemi, “Positioning error prediction and training data evaluation in rf fingerprinting method,” in _2019 International Conference on Indoor Positioning and Indoor Navigation (IPIN)_ , 2019, pp. 1–7.
  * [17] A. Nikitin, C. Laoudias, G. Chatzimilioudis, P. Karras, and D. Zeinalipour-Yazti, “Indoor localization accuracy estimation from fingerprint data,” in _2017 18th IEEE International Conference on Mobile Data Management (MDM)_ , 2017, pp. 196–205.
  * [18] D. Dearman, A. Varshavsky, E. de Lara, and K. N. Truong, “An exploration of location error estimation,” in _UbiComp 2007: Ubiquitous Computing_ , J. Krumm, G. D. Abowd, A. Seneviratne, and T. Strang, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 181–198.