Minimum Viable Model Estimates for Machine Learning Projects

John Hawkins

Transitional AI Research Group, Sydney, Australia

<EMAIL_ADDRESS>
Abstract

Prioritization of machine learning projects requires estimates of both the
potential ROI of the business case and the technical difficulty of building a
model with the required characteristics. In this work we present a technique
for estimating the minimum required performance characteristics of a
predictive model given a set of information about how it will be used. This
technique will result in robust, objective comparisons between potential
projects. The resulting estimates will allow data scientists and managers to
evaluate whether a proposed machine learning project is likely to succeed
before any modelling needs to be done. The technique has been implemented into
the open source application MinViME (Minimum Viable Model Estimator) which can
be installed via the PyPI python package management system, or downloaded
directly from the GitHub repository. Available at https://github.com/john-
hawkins/MinViME

Keywords

Artificial Intelligence, ROI Estimation, Machine Learning Metrics, Cost
Sensitive Learning

## 1\. Introduction

In an ideal world we would priortise all our machine learning projects
according to their expected payoff. This would, in turn, require a reasonable
estimate of both the probability of success and the return given success.
Creating such an estimate is difficult and instead we tend to be guided by
heuristics, implicit and explicit, about the size of the problem and its
difficulty.

Difficulty estimation is often limited to a discussion about what it would
take to feasibly productionise a model. Consideration of the difficulty of
building the model, to the extent that it is done at all, is usually an
informal process relying on the experience of data scientist working on the
project. In this work we present a structured approach to estimating the
difficulty of a machine learning task by estimating the minimum required
performance of a model that would meet the business objectives.

The most common approach to which the details of the business case are
considered in the building of a machine learning solution is through
development of domain specific methods. For example, work on designing loss
functions for specific applications (See [1], [2]). Alternatively, there is a
practice of using cost sensitive learning [3, 4, 5, 6] to permit general
techniques learn a solution that is tuned to the specifics of business problem
(for example [7]). In general, you can achieve the same outcome by learning a
well calibrated probability estimator and adjusting the decision threshold on
the basis of the cost matrix [8]. The cost matrix approach to analysing a
binary classifier can be also be used to estimate its expected ROI [9].

These approaches allow you to improve and customise a model for a given
application, but they do not help you decide if the project is economically
viable in the first place. There have been extensive efforts to characterise
project complexity in terms of data set requirements for given performance
criteria [10], however the task of determining these criteria from business
requirements has not, to our knowledge, been addressed. Fortunately, as we
will demonstrate, some of the machinery for the cost matrix analysis can be
used in reverse to estimate the baseline performance a model would require to
reach minimal ROI expectations.

We use the framework of a cost matrix to develop a method for estimating the
minimum viable binary classification model that satisifies a quantitative
definition of the business requirements. We demonstrate that a lower bound on
precision can be calculated a priori, but that most other metrics require
knowledge of either the expected number of false positives or false negatives.
We devise an additional simplicity metric which allows for direct comparison
between machine learning projects purely on the basis of the business
criteria. Finally, we demonstrate that lower bounds on the AUC and other
metrics can be estimated through numerical simulation of the ROC plot.

We use the minvime ROC simulator to conduct a range of experiments to explore
the manner in which the minimum viable model changes with the dimensions of
the problem landscape. The results of these simulations are presented in a
manner to help develop intuitions about the feasibility of machine learning
projects.

## 2\. Global Criteria

In order to estimate the required performance characteristics of any machine
learning model there are several global criteria that need to be defined.
Firstly, we need to determine the overall ROI required to justify the work to
build and deploy the model. This can take the form of a project cost plus
margin for expected return. Given an overall estimate we then need to amortise
it to the level of analysis we are working at, for example, expected yearly
return.

In addition, we require the expected number of events that will be processed
by the machine learning system for the level of analysis mentioned above. For
example, this may be the number of customers that we need to analyse for
likelihoood of churn every year. Included in this we require to know the base
rate of the event we are predicting. In other words how many people tend to
churn in a year.

Given the above criteria, we next need to evaluate how a model would need to
perform in order to process the events at a sufficient level of compentence.
To do this analysis we need to look at the performance of the model as it
operates on the events we are predicting. We start this process by looking at
the impact a binary classifcation model using a cost matrix analysis.

## 3\. Cost Matrix Analysis

The idea of a cost matrix in machine learning emerged from the literature on
training models on imbalanced datasets. Specifically there are a set of papers
in which the learning algorithm itself is designed such that the prediction
returned is the one that minimises the expected cost [5, 4]. Typically this
means not necessarily returning the class that is most probable, but the one
with the lowest expected cost. In essence, the model is learning a threshold
in the probability space for making class determinations. The same idea can be
applied after training an arbitrary classfier that returns a real-value result
in order to determine a threshold that minimises the error when the model is
applied.

It has been observed that the later approach is not dissimilar from a game
theoretical analysis [11]. In game theory the cost matrix would delineate the
outcomes for all players of a game in which each player’s fortunes are
affected by their own decision and that of the other player. Each axis of the
matrix describes the decisions of a player and by locating the cell at the
nexus of all decisions we can determine the outcome of the game.

In binary classification we replace the players of the game with two sources
of information: the prediction of the model and the ground truth offered up by
reality. For binary classification problems this is a 2x2 matrix in which each
of the cells have a direct corollary in the confusion matrix. In the confusion
matrix we are presented with the number of events that fall into each of the
four potential outcomes. In a cost matrix we capture the economic (or other)
impact of each specific type of event occuring. From here on in we will follow
Elkan [5] and allow the cells of the matrix to be either depicted as costs
(negative quantities) or benefits (positive quantities).

In table 1 we demonstrate the structure of a cost matrix in which we also
depict the standard polarity of the costs and benefits in most applications.

Table 1: Cost/Benefit Matrix for a Binary Classification System. | Actual Negative | Actual Positive
---|---|---
Predict Negative | True Negative (TN) | False Negative (FN)
Predict Positive | False Positive (FP) Cost | True Positive (TP) Benefit

If we assume, without loss of generality, that the positive class corresponds
to the event against which we plan to take action, then as observed by Elkan
[5] the entries of the cost matrix for the predicted negative row should
generally be identical. The impact of taking no action is usually the same. It
is also usually zero because the model is typically being applied to a status
quo in which no action is taken. In most cases you will only need to estimate
the cost/benefits of your true positives and false positives. This is because
the assumed ground state (the status quo) is the one in which the actual costs
of false negatives are already being incurred. Similarly the true negatives
are the cases in which you have taken no action on an event that does not even
represent a loss of opportunity.

The overall outcome of using a model can be calculated by multiplying the
cells of the confusion matrix by the corresponding values in the cost matrix
and summing the results. This tells us what the overall impact would have been
had we used the model on the sample of data used to determine the confusion
matrix. The advantage of using a cost/benefit structure in the matrix is that
we can then read off whether the result of using the model was net-positive or
net-negative depending on the polarity of the total impact. This approach can
then be used to optimise a decision threshold by varying the threshold and
recalculating the expected impact.

Before continuing with the analysis of the cost matrix we make some
observations about the process of determining the content of the cost matrix.

###  3.1 Deterministic Outcomes

In some situations we can say in advance what the impact of a specific
conjunction of prediction and action will be. We might know, for example, that
giving a loan to a customer who will default will result in a loss of money.
If we take the action of rejecting that loan then the money is no longer lost.

This determinism in outcomes is typically true of situations where once we
have taken our action there is no dependence on additional decisions by other
people to impact whether our action is effective. This is not always true.
There are certain circumstances in which knowing what will happen does not
guarantee a given outcome. In these cases there is stochastic relationship
between the prediction and the outcome.

###  3.2 Stochastic Outcomes

In the stochastic situation we understand that even in the event of possessing
a perfect oracle for the future events, we would still be left with an element
of randomness in the outcomes of our actions. This is typically true where the
action we take will involve third parties or processes beyond our control to
achieve their result.

The canonical example of a predictive process with a stochastic outcome is a
customer churn model. A churn model is used to intervene and attempt to
prevent a customer from churning. However, our ability to influence a customer
to not churn is not guaranteed by accurately predicting it.

In these situations we need to make additional estimates in order to get the
contents of the cost matrix. For example, if we expect that our interventions
will succeed one in every five attempts, and the value of a successful
intervention is $1,000, then the benefit of a true positive is $200. Defining
a cost/benefit matrix for a situation with stochastic outcomes will require
additional assumptions.

###  3.3 Metrics from the Matrix

If we have an existing model then we can use the cost matrix to estimate the
economic impact by combining it with a confusion matrix. How then can we
estimate the required properties of a model such that it renders a project
viable?

Once the model is deployed it will make predictions for the $N$ cases that
occur in the period of analysis. We can demarcate the content of the unknown
confusion matrix with the symbol $tp$ true positives and $fp$ for false
positives. The benefit of each true positive is captured by the value
$\mathcal{B}$ from our cost matrix, similarly the cost of each false positive
is demarcted $\mathcal{C}$. In order to satisfy that the model meets our ROI
minimum $\mathcal{M}$ then the following will need to hold:

$tp\cdot\mathcal{B}-fp\cdot\mathcal{C}>=\mathcal{M}$ (1)

From Equation 1 we can derive explicit defintions for both the number of true
positive and false positive predictions that would satisfy our minimal
criteria.

$tp>=\frac{fp\cdot\mathcal{C}+\mathcal{M}}{\mathcal{B}}$ (2)

$fp<=\frac{tp\cdot\mathcal{B}-\mathcal{M}}{\mathcal{C}}$ (3)

Note, that these expressions involve a reciprocal relationship between the
true positives and false positives. Moderated by the ratio between the costs
and benefits of taking action in relation to the overall required return. We
would like to be able to estimate the True Positive Rate (TPR) or recall, and
the False Positive Rate (FPR) or fall-out. To define these we need to
introduce one of the other global terms required to estimate model
requirements: the base rate $r$ at which the event we are predicting occurs.
Which allows us to define TPR and FPR as follows:

$TPR=\frac{tp}{N\cdot r}$ (4) $FPR=\frac{fp}{N\cdot(1-r)}$ (5)

One of the most important metrics for evaluating a binary classification model
is precision (or Positive Predictive Value). This determines the proportion of
all predicted positive events that were correct. We can now generate an
expression for what this would be at the exact point of meeting of ROI
requirements.

$precision=\frac{1}{1+\frac{\mathcal{B}}{\mathcal{C}}-\frac{\mathcal{M}}{tp\cdot\mathcal{C}}}$
(6)

Unfortunately, the presence of the number of true positives in this expression
resists additional simplification. However, if we focus our attention on the
break even point we can eliminate this component and derive the following
lower bound on the precision of the model.

$precision>\frac{1}{1+\frac{\mathcal{B}}{\mathcal{C}}}$ (7)

This lower bound corresponds to the line on which the benefits of the true
positives are counteracted by the cost of the false positives. Equation 6
corresponds to the line that lies $\frac{\mathcal{M}}{\mathcal{B}}$ units
above this lower bound. We depict both of these lines schematically in Figure
1.

Figure 1: Schematic representation of potential model outputs depicting the
boundary between viable and inviable models.

This geometric representation suggests an additional metric for the
feasibility of any given project. Note, that the triangle in the upper left
corner represents all possible model outputs that would satisfy the stated
criteria. The dimensions of this triangle are shown in the legend by the lines
of blue and green squares. We can calculate the area of this triangle and
express it as a proportion of the total area of the rectangle. This will
result in a metric which is the proportion of all potential model outputs that
are economically viable. We can define an expression for this ratio as shown
in Equation 8.

$simplicity=\frac{\mathcal{B}}{2\cdot\mathcal{C}\cdot N^{2}\cdot
r(1-r)}\cdot\left(N\cdot r-\frac{\mathcal{M}}{\mathcal{B}}\right)$ (8)

This simplicity measure allows us to rank the difficulty of problems based
purely on their business criteria. As this metric is to be calculated in a
discrete space, and it must deal with edge cases in which the triangle is
truncated, we refer the reader to the source code for its complete
calculation.

Our efforts to estimate additional requirements in terms of standard machine
learning metrics are hindered by the requirement to know either the number of
true positives or false positives in the minimal viable model. To move beyond
this limitation we will employ numerical techniques to estimate the bounds on
the model performance.

###  3.4 ROC Simulation

The ROC plot provides us with a method of examinining the performance of a
binary classification model in a manner that does not depend on a specific
threshold. As shown in the example plot 2 it depicts a range of TPR and FPR
values across all the different potential thresholds of the model.
Additionally, the ROC plot permits the definition of a non-threshold specific
metric to evaluate the overall discriminative power of the model (regardless
of chosen threshold), the Area Under the Curve (AUC) [12].

We can generate synthetic ROC plots by exploiting their key properties.

  1. 1.

They are necessarily monotonically increasing

  2. 2.

They will necessarily pass through the points (0,0) and (1,1)

  3. 3.

They remain above or on the diagonal between these points (unless the model is
inversely calibrated)

We can simulate such curves by designing a parameterised function that
satisfies these criteria. We have designed the function shown in equation 9
for this purpose. It is composed of two parts, both of which independantly
satisfy the criteria.

$y=\alpha\cdot(-(x-1)^{2\beta}+1)+(1-\alpha)\cdot x$ (9)

In this function $\alpha$ is a weighting between zero and one that determines
the mix between the two component functions. The first component function is
an inverted even-powered polynomial function of x, offset so that its origin
lies at coordinates (1,1). The second of function is just $x$ itself, which is
the diagonal line. Both functions necessarily pass through the points $(0,0)$
and $(1,1)$, and their mixture determines the shape of the curve between these
points.

Figure 2 shows some examples of simulated ROC plots with different values of
the parameters $\alpha$ and $\beta$. Note, it is critical to restrict the
space of these parameters such that the search space does not result in
unrealistic ROC plots. The third plot in Figure 2 demonstrates that at extreme
values the curve starts to represent a spline of two linear functions.

Figure 2: Examples of ROC plots generated through simulation.

By running simulations across a variety of values of $\alpha$ and $\beta$ we
can generate multiple synthetic ROC plots. Each plot provides a set of
hypothetical values for the true positive rate and false positive rate across
a range of thresholds. If we combine these with the pre-defined
characteristics of the problem we are solving (the number of cases, baseline
rate, costs and benefits). Then we can identify which of these synthetic ROC
plots would represent acceptable models. We can then rank them such that we
search for the ROC curve with the minimum AUC that satisfies our ROI
constraints.

Given a minimum viable AUC, we can then calculate estimates for the other
metrics discussed above by taking the TPR and FPR values at the minimum viable
threshold.

We have built a Python tool for running such simulations and released it as a
Python package to the PyPI library (Available at $https://github.com/john-
hawkins/MinViME$).

## 4\. Simulations

The MinViME tool allows us to explore the viable model landscape across a
broad range of business criteria in order to understand what drives
feasibility. In order to do so we need to represent the space of business
problems in an informative and comprehensible way. We do so using the
following three ratios.

Figure 3: Relationship between the AUC of the minimum viable model and the
criteria that determine the business problem. Note: We use a log scale on the
X axis to demonstrate how the difficulty varies with each order of magnitude
change in the nature of the problem.

  1. 1.

Benefit to Total ROI. In other words how significant to the total required ROI
is each additional true positive. We explore a range of values between zero
and one.

  2. 2.

Cost to Benefit. In other words how many false positives can we tolerate for
every true positive case we identify. In principle this could take values
greater than one, however this is almost never seen in practice so we explore
a range of values between zero and one

  3. 3.

Base Rate. The rate at which the event being predicted occurs in the
population.

In addition, we use a static value of one million for the number of cases
processed per unit of time This value is realistic for many problems and our
experimentation with values an order of magnitude either side demonstrates
little variation on the results presented below.

###  4.1 Results

In Figure 3 we see the mean minimum required AUC across each of the three
dimensions. The gray bounds indicate the lower and upper quartiles. We see
that problems get easier as either the base rate or the benefit to total ROI
increases. The opposite occurs for the cost to benefit ratio as we would
expect. The base rate plot demonstrates the most variability across the
spectrum, with base rates between $10^{-3}$ and $10^{-2}$ demonstrating the
largest variability. Above this range problems generally have less strenuous
requirements and below it most problems require demanding model performance.

Elkan [5] observes that if we scale the entries in a cost matrix by a positive
constant or if we add a constant to all entries, it will not affect the
optimal decision. However, this is not true for determining the performance
characterists of the minimum viable model. The interaction between the number
of cases, the required ROI and the values in the cost matrix will mean that
these changes can affect the characteristics of the minimum viable model.

When the benefit to total ROI ratio is very small we see that the model
requirements are predominantly for an AUC above $0.8$. To further understand
the relationship between the base rate and the cost to benefit ratio we plot
the surface of required AUC for these values when the benefit to total ROI is
$10^{-4}$.

Figure 4: Surface of the AUC landscape for business problems in which each TP
prediction returns $10^{-4}$ of the total required ROI

In Plot 4 we can see clearly that very low base rates can generate intractable
problems that are likely to extremely difficult to solve. As the base rate
increases we see a general weakening of the effect of the cost/benefit ratio.
However, at very low cost/benefit ratios the impact of the base rate is almost
entirely mitigated.

## 5\. Conclusion

We have discussed the problem of estimating the minimum viable performance
characteristics of a machine learning model. The goal is to be able to do this
before projects are undertaken so that we can make informed decisions about
where our resources are most effectively deployed.

We investigated whether analytical lower bounds of standard machine learning
performance metrics could be calculated from business criteria alone. We
derived this for the model precision, and derived a novel simplicity metric
which permits a priori comparisons of project complexity for binary
classification systems.

We then demonstrated a numerical approach to estimating the minimum viable
model performance through simulation of ROC plots. This in turn allows lower
bound estimation of all other standard binary classification metrics. Using
this method we explored the space of minimum viable models for a range of
business problem characteristics. What we observe are some critical non-linear
relationships between the base rate and the cost/benefit ratio that will
determine whether a project is feasible or likely to be intractable.

## References

  * [1] J. Johnson and T. Khoshgoftaar, “Medicare fraud detection using neural networks,” _Journal of Big Data_ , vol. 63, no. 6, 2019. [Online]. Available: https://doi.org/10.1186/s40537-019-0225-0
  * [2] C. Hennig and M. Kutlukaya, “Some thoughts about the design of loss functions,” _REVSTAT - Statistical Journal Volume_ , vol. 5, pp. 19–39, 04 2007.
  * [3] P. Domingos, “Metacost: A general method for making classifiers cost-sensitive,” _Proceedings of the Fifth ACM SIGKDD Int’l. Conf. on Knowledge Discovery & Data Mining_, pp. 155–164, 8 1999.
  * [4] D. D. Margineantu, “On class-probability estimates and cost-sensitive evaluation of classifiers,” in _In Workshop on Cost-Sensitive Learning at the Seventeenth International Conference on Machine Learning (WCSL at ICML2000_ , 2000.
  * [5] C. Elkan, “The foundations of cost-sensitive learning,” _Proceedings of the Seventeenth International Conference on Artificial Intelligence: 4-10 August 2001; Seattle_ , vol. 1, 05 2001.
  * [6] Y. Tian and W. Zhang, “Thors: An efficient approach for making classifiers cost-sensitive,” _IEEE Access_ , vol. 7, pp. 97 704–97 718, 2019.
  * [7] H. Fatlawi, “Enhanced classification model for cervical cancer dataset based on cost sensitive classifier,” 08 2017.
  * [8] N. Nikolaou, N. Edakunni, M. Kull, P. Flach, and G. Brown, “Cost-sensitive boosting algorithms: Do we really need them?” _Machine Learning_ , vol. 104, 08 2016.
  * [9] O. Ylijoki, “Guidelines for assessing the value of a predictive algorithm: a case study,” _Journal of Marketing Analytics_ , vol. 6, 01 2018.
  * [10] S. Raudys and A. Jain, “Small sample size effects in statistical pattern recognition: Recommendations for practitioners,” _Pattern Analysis and Machine Intelligence, IEEE Transactions on_ , vol. 13, pp. 252–264, 04 1991.
  * [11] I. E. Sanchez, “Optimal threshold estimation for binary classifiers using game theory,” _ISCB Comm J_ , p. 2762, 2017. [Online]. Available: https://doi.org/10.12688/f1000research.10114.3
  * [12] A. P. Bradley, “The use of the area under the roc curve in the evaluation of machine learning algorithms,” _Pattern Recognition_ , vol. 30(7), pp. 1145–1159, 1997.

© 2020 By AIRCC Publishing Corporation. This article is published under the
Creative Commons Attribution(CC BY) license.


negligible for interaction strength in regime II with
$0.04E_{0}<U(a^{M}_{1})<0.09E_{0}$ but a strong band mixing is found in regime
III with larger Coulomb interaction $0.09E_{0}<U(a^{M}_{1})$. Fig. S12 shows
the self-consistent solutions for the initial density matrix
$\rho^{\text{R}}_{ini}=\tau_{0}\sigma_{y}$ as an example. Here Fig. S12(a)-(c)
are for $U(a_{1}^{\text{M}})=0.13E_{0}$ (regime III) while (d)-(f) are for
$U(a_{1}^{\text{M}})=0.08E_{0}$ (regime II). The filled band
$|\tilde{\psi}^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$ in Fig. S12(a) for
regime III has the Chern number $C=0$ and that in Fig. S12(d) for regime II
has $C=+1$. Fig. S12(b) and (e) show the distribution of Berry curvature
$\Omega_{xy}$ in the moiré BZ for regime III and regime II, respectively. Fig.
S12(c) and (f) show the projection of
$|\tilde{\psi}^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$ into non-interacting
states $|u^{\text{CB1}}_{\pm i}({\textbf{{k}}})\rangle$ and
$|u^{\text{CB2}}_{\pm i}({\textbf{{k}}})\rangle$ in regime III and regime II,
respectively. One can see that the interacting ground state
$|\tilde{\psi}^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$ has a strong
component from $|u^{\text{CB2}}_{\pm i}({\textbf{{k}}})\rangle$, in addition
to $|u^{\text{CB1}}_{\pm i}({\textbf{{k}}})\rangle$, due to the strong band
mixing in regime III, while only the $|u^{\text{CB1}}_{\pm
i}({\textbf{{k}}})\rangle$ part dominates the interacting ground state in
regime II. Thus, from Fig. S12, we show that the Coulomb interaction can drive
the interacting ground state into a trivial Mott insulator phase Sorella and
Tosatti (1992); Neto et al. (2009) via band mixing between CB1 and CB2 when
there is inversion symmetry.

Figure S13:  (a) The energies of the self-consistent Hartree Fock states with
different initial density matrices $\rho^{\text{R}}_{ini}({\textbf{{k}}})$.
(b) The spectra (orange) of $H^{\text{HF}}({\textbf{{k}}})$ with
$\rho^{\text{R}}_{ini}({\textbf{{k}}})=\sigma_{y}\tau_{0}$ and
$U(a_{1}^{\text{M}})=0.13E_{0}$. The cyan lines are non-interacting spectrum.
(c) Berry curvature for the lowest band of $H^{\text{HF}}({\textbf{{k}}})$ in
(b). (d) The overlap between the ground states
$|\tilde{\psi}^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$ for four-band
Hartree-Fock calculations and the ground states
$|\psi^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$ for two-band Hartree-Fock
calculations. Here the calculation is for the four-band model with the
parameters $\phi=0,\alpha=0,V_{0}/E_{0}=1.2$.

For the case $\alpha=0,V_{0}/E_{0}=1.2,\phi=0$ without inversion symmetry, the
self-consistent solutions are summarized in Fig. S13. As shown in the phase
diagram of Fig. 4(b) in the main text, the system stays in the QAH state with
$C=+1$ for the realistic Coulomb interaction strength
$U(a_{1}^{\text{M}})=U_{0}\approx 0.13E_{0}$, which is quite different from
the inversion symmetric case. Here we show more details of this calculation in
Fig. S13 for $U(a_{1}^{\text{M}})=0.13E_{0}$. Fig. S13(a) shows that the self-
consistent solutions with the initial complex density matrices that break the
$C_{2z}\mathcal{T}$ still have lower energy. We consider $\sigma_{0}\tau_{y}$
as an example and show the energy dispersion of HF bands in Fig. S13(b). The
distribution of the Berry curvature $\Omega_{xy}$ in the moiré BZ is shown in
Fig. S13(c). We further project the interacting ground state of the four-band
model, denoted as $|\tilde{\psi}^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$,
into that of the two-band model, denoted as
$|\psi^{\text{HF}}_{j=1}({\textbf{{k}}})\rangle$, and find their overlap is
almost 1 in the whole moiré BZ, as shown in Fig. S13(d). Thus, the inter-band
mixing is negligible in the inversion asymmetric case for
$U(a_{1}^{\text{M}})=U_{0}\approx 0.13E_{0}$.

Figure S14:  (a) Spectrum for the moiré system for smaller moiré unit cells.
(b) The spectra (orange) of $H^{\text{HF}}$ with
$\rho^{\text{R}}_{y}({\textbf{{k}}})$ for mirror polarized states. Here the
calculation is for the two-band model with the parameters
$\phi=1/3,\alpha=0.16,V_{0}/E_{0}^{\prime}=1.1$.

###  II.5 Coulomb interaction for smaller moiré unit cells

In this section, we discuss the moiré systems with a smaller moiré lattice
constant, $|{\textbf{{a}}}_{1}^{\text{M}}|=14$nm, for a twist angle
$\theta=1.0^{\circ}$. The Coulomb interaction scales inversely with the moiré
unit cell length and its strength can be estimated as
$U({\textbf{{a}}}_{1}^{\text{M}})=10$ meV.

The spectrum is shown in Fig. S14(a). The parameters for the spectrum are
$\Delta_{1}=14$meV, $\phi=1/3$, and $m_{0}=30$meV. The energy scale is
$E_{0}^{\prime}=v|{\textbf{{b}}}_{1}^{\text{M}}|=77$meV. The bandwidth of CB1
is $7.4$meV and the direct gap between CB1 and CB2 is $4$meV. The ratio
between $U({\textbf{{a}}}_{1}^{\text{M}})$ and bandwidth is smaller for the
smaller $|{\textbf{{a}}}_{1}^{\text{M}}|$. When the Coulomb interaction is
considered for CB1 with the density matrix
$\rho^{\text{R}}_{y}({\textbf{{k}}})$ as shown in Fig. S14(b), the mirror
polarized states with $C=+1$ can be induced.

##  III Computational Methods for DFT calculations and moiré lattice

The DFT calculations were performed with the Vienna Ab initio Simulation
Package (VASP)${}^{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{kresse1996efficient}{\@@citephrase{(}}{\@@citephrase{)}}}}$.
The exchange-correlation functional was chosen as the Perdew-Burke-Ernzerhof
type generalized-gradient approximation${}^{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{perdew1996generalized}{\@@citephrase{(}}{\@@citephrase{)}}}}$
and the projector-augmented-wave method was used for the core-electron
potentials${}^{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{blochl1994projector,kresse1999ultrasoft}{\@@citephrase{(}}{\@@citephrase{)}}}}$.
The energy cutoff was set as 340 eV for all calculations. The convergence
criterion was set as 10-5 eV for self-consistent electronic calculations and
the $k$-point meshes were set as $13\times 13\times 1$ to sample the Brillouin
zone. We used the DFT-D3 method${}^{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{grimme2010consistent}{\@@citephrase{(}}{\@@citephrase{)}}}}$
to correctly describe the van der Waals interactions.

As discussed in the main text and the next section IV, the moiré potential
$\Delta({\textbf{{r}}})$ in the twisted Sb2/Sb2Te3 hetero-structures could be
obtained from the uniform potential
$\tilde{\Delta}({\textbf{{d}}}_{{\textbf{{R}}}})$ for the hetero-structure
with a uniform shift ${\textbf{{d}}}_{{\textbf{{R}}}}$ between the Sb2 layer
and Sb2Te3 layer. Therefore, we will describe below our first principles
calculations of the band structures with different stacking configurations
that correspond to different shifting vectors
${\textbf{{d}}}_{{\textbf{{R}}}}$. As shown in Fig. S15(a), the
heterostructure model with the specific stacking contains $1\times 1$ Sb2
monolayer and $1\times 1$ 2QL Sb2Te3 thin film. A vacuum layer with 20 Å was
added along the z direction to avoid the interaction between adjacent slabs.
Because the uniform potential
$\tilde{\Delta}({\textbf{{d}}}_{{\textbf{{R}}}})$ is induced by coupling
between Sb2 monolayer and Sb2Te3 thin films, we fixed the value of the van der
Waals gap inside 2QL Sb2Te3 thin film as the bulk value (2.708 Å), and then
let Sb2 monolayer and its neighboring atoms fully relax until the calculated
forces are smaller than 0.01 eV/Å. In Fig. 5 of the main text, we plot the
relaxed lattices and corresponding band structures for heterostructures with
AA, AB, and BA stacking, whose interlayer distances between Sb2 monolayer and
2QL Sb2Te3 thin film are 3.92 Å, 2.81 Å, and 2.91 Å respectively.

Besides the stacking configurations shown in the main text, other stacked
heterostructures can be also found in the moiré pattern (see Fig. 5 in the
main text). In order to calculate moiré potential $\Delta({\textbf{{r}}})$
accurately (see Section IV below), we consider extra nine stacked
configurations, as shown in Fig. S15(b-d), which are located in the
intermediate regions among AA, AB, and BA stackings, named AAmAA-X (Fig.
S15(b)), AAmAB-X (Fig. S15(c)), and AAmBA-X (Fig. S15(d)) with X = I, II, III.
We take the heterostructure models with AAmAA stackings as an example to show
their lattices in details and the other intermediate stacked configurations
could be obtained by using the same method. There are three types of AAmAA
stacked structures in the moiré pattern, corresponding to the shifting vector
${\textbf{{d}}}_{R}$ as
$\tilde{{\textbf{{a}}}}_{1}/2+0\tilde{{\textbf{{a}}}}_{2}$,
$0\tilde{{\textbf{{a}}}}_{1}-\tilde{{\textbf{{a}}}}_{2}/2$, and
$\tilde{{\textbf{{a}}}}_{1}/2+\tilde{{\textbf{{a}}}}_{2}/2$ for AAmAA-I,
AAmAA-II, and AAmAA-III respectively. These three configurations are related
by C3z rotation, and thus we only need to calculate the electronic structure
for one of them. On the other hand, because these intermediate stacked
structures are not in the local minimum of the potential energy surface, we
only relax the $z$ direction coordinate while fix the $x$ and $y$ coordinates
of the Sb2 monolayer (see Fig. S15(a)). The corrugation effect, which is
crucial for predicting the correct band structure in twisted bilayer
grapheneLucignano et al. (2019); Uchida et al. (2014); Koshino et al. (2018),
is taken into account after this lattice relaxation process. Using the same
method, we obtained the relaxed intermediate stacked heterostructures for
AAmAB-X and AAmBA-X (X = I, II, III). The interlayer distances between Sb2
monolayer and 2QL Sb2Te3 thin film are 2.94 Å, 3.38 Å, and 3.42 Å for AAmAA,
AAmAB, and AAmBA stacked heterostructures respectively and the related band
structures are shown in Fig. S15(e) correspondingly.

Figure S15:  (a) The side view of the Sb2/Sb2Te3 heterostructure with AA
stacking. The solid black lines mark the unitcell used in DFT calculations. Sb
atoms in Sb2 monolayer are marked as gray, Sb atoms and Te atoms in Sb2Te3
films are marked as red and blue. The atoms in the region with the green
background are frozen when we relax the lattice structures, the $x$ and $y$
coordinates of Sb atoms in the region with the yellow background are fixed.
The value of the van der Waals gap inside 2QL Sb2Te3 films is 2.078Å. (b-d)
The top view for heterostructures with the stacking of AAmAA-X, AAmAB-X, and
AAmBA-X (X = I, II, III). Corresponding ${\textbf{{d}}}_{R}$ is shown by the
green arrows. The black arrows show the lattice vector of Sb2 monolayer.  (e)
Calculated band structures of the heterostructure with the stacking of AAmAA,
AAmAB, and AAmBA, respectively. The Fermi levels are set as zero. Figure S16:
(a) The lattice structures for the twisted Sb2 (gray) on top of Sb2Te3 (blue
and orange) at a location same as Fig.4(c) of the main text.
$\tilde{{\textbf{{a}}}}_{1,2}$ are primitive lattice vectors for the Sb2Te3
layer and $\tilde{{\textbf{{a}}}}_{1,2}^{\prime}$ are primitive lattice
vectors for Sb2 layer. ${\textbf{{d}}}_{{\textbf{{R}}}_{1,2}}$ are the local
lattice shifts at sites ${\textbf{{R}}}_{1,2}$. (b) Commensurate Sb2 (gray) on
top of Sb2Te3 with a constant shift ${\textbf{{d}}}_{{\textbf{{R}}}_{2}}$ as
an approximation for the local lattice structures in (a). (c) BZs for Sb2Te3
(orange), Sb2 (gray), and moiré superlattice (blue).
$\tilde{{\textbf{{b}}}}_{1,2}$, $\tilde{{\textbf{{b}}}}_{1,2}^{\prime}$,
$\tilde{{\textbf{{b}}}}^{\text{M}}_{1,2}$ are reciprocal lattice vectors for
for Sb2Te3, Sb2, moiré superlattice, respectively.

##  IV Moiré potentials from the fitting to the DFT band structure

In this section, we discuss the method to obtain Moire potential
$\Delta({\textbf{{r}}})$ in Eq. 1 of the main text from the above DFT
calculation in Sec. III Jung et al. (2014).

We first consider the effective Hamiltonian $H^{\text{DFT}}$ in Eq. (3) of the
main text. Compared to the original Hamiltonian $H_{0}$ in Eq. 1, the
spatially dependent moiré potential term $\Delta({\textbf{{r}}})$ is changed
to a uniform potential term $\tilde{\Delta}({\textbf{{d}}}_{R})$ for a fixed
${\textbf{{d}}}_{R}$ that describes the relative shift between Sb2 and Sb2Te3
layers. $H^{\text{DFT}}$ describes the effective model for the hetero-
structure with a uniform shift between two atomic layers, and different values
of the shifting vector ${\textbf{{d}}}_{R}$ describe different stacking
configurations. Thus, we can use the energy dispersion of $H^{\text{DFT}}$ to
fit to that from the DFT calculations.

The effective Hamiltonian $\hat{H}^{\text{DFT}}({\textbf{{d}}}_{R})$ is then
given by

$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}}_{R})|{\textbf{{k}}}_{2},\beta_{2}\rangle=\delta({\textbf{{k}}}_{1}-{\textbf{{k}}}_{2})(\begin{pmatrix}h_{D}^{t}({\textbf{{k}}}_{1})&ms_{0}\\\
ms_{0}&h_{D}^{b}({\textbf{{k}}}_{1})\end{pmatrix}+\begin{pmatrix}\tilde{\Delta}({\textbf{{d}}}_{R})s_{0}&0\\\
0&\alpha\tilde{\Delta}({\textbf{{d}}}_{R})s_{0}\end{pmatrix})=\delta({\textbf{{k}}}_{1}-{\textbf{{k}}}_{2})H^{\text{DFT}}({\textbf{{k}}}_{1},{\textbf{{d}}}_{\textbf{{R}}}),$
(S81)

where $H^{\text{DFT}}({\textbf{{k}}},{\textbf{{d}}}_{\textbf{{R}}})$ is just
Eq.(4) in the main text, $h_{D}^{t/b}({\textbf{{k}}})$ are the top/bottom
Dirac surface states same as Eq. (1) of the main text, $s_{0}$ are the
identical matrix in spin space, $m$ is the tunnelling between two surfaces,
and $\alpha$ captures the difference in the potentials on two surfaces created
by the Sb2 layer. $|{\textbf{{k}}},\beta\rangle$ is the atomic Bloch states
for the Sb2Te3 and Sb2 lattice with a constant shift. $\beta_{1,2}=1,...,4$
represents both the spin and layer degrees of freedom. The spectra of this
model are given by

$E^{\text{DFT}}_{\eta,\xi}({\textbf{{k}}},{\textbf{{d}}})=\frac{1+\alpha}{2}\tilde{\Delta}({\textbf{{d}}})+\eta\sqrt{m^{2}+\left(\frac{1-\alpha}{2}\tilde{\Delta}({\textbf{{d}}})+\xi
v^{2}k^{2})\right)^{2}}.$ (S82)

with $\eta=\pm,\xi=\pm$. By fitting
$E^{\text{DFT}}_{\eta,\xi}({\textbf{{k}}},{\textbf{{d}}})$ to the spectrum
calculated from DFT in Fig.4(b) of the main text and Fig. S15(f), the
parameters $\alpha$, $\tilde{\Delta}({\textbf{{d}}})$, $v$ and $m$ in the
model Hamiltonian can be obtained. From the three-fold rotation symmetry
$C_{3}$ of the underlying lattice, one can obtain

$\langle{\textbf{{k}}},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}})|{\textbf{{k}}},\beta_{2}\rangle=\langle
C_{3}{\textbf{{k}}},\beta_{1}|\hat{H}^{\text{DFT}}(C_{3}{\textbf{{d}}})|C_{3}{\textbf{{k}}},\beta_{2}\rangle$
(S83)

and

$E^{\text{DFT}}_{\eta,\xi}({\textbf{{k}}},{\textbf{{d}}})=E^{\text{DFT}}_{\eta,\xi}(C_{3}{\textbf{{k}}},C_{3}{\textbf{{d}}})=E^{\text{DFT}}_{\eta,\xi}({\textbf{{k}}},C_{3}{\textbf{{d}}}),$
(S84)

where the last line following from $E^{\text{DFT}}_{\eta,\xi}$ depends on
$|{\textbf{{k}}}|$ in Eq. (S82). So, only one of three stacking related by
$C_{3}$ shown in Fig. S15(b)-(d) needs to be fitted.

Next we will establish the relation between the moiré Hamiltonian $H_{0}$ in
Eq.1 of the main text and the Hamiltonian $H^{\text{DFT}}$ constructed from
the fitting to the DFT energy spectrum. We start from the momentum-space moiré
Hamiltonian
$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{k}}}_{2},\beta_{2}\rangle$
in Eq.1 of the main text, where $\beta_{1,2}=1,\cdots,4$ label both layer and
spin indices and $|{\textbf{{k}}},\beta\rangle$ are atomic Bloch states for
each layer underlying the moiré superlattice. As the moiré Hamiltonian $H_{0}$
does not preserve atomic lattice translation, the crystal momentum k is not a
good quantum number and $H_{0}$ can mix different k states. The atomic Bloch
wave function $|{\textbf{{k}}},\beta\rangle$ is related to atomic Wannier
function $|{\textbf{{R}}},\beta\rangle$ by

$\displaystyle|{\textbf{{k}}},\beta\rangle=\sum_{{\textbf{{R}}}}e^{i{\textbf{{k}}}\cdot{\textbf{{R}}}}|{\textbf{{R}}},\beta\rangle,$
(S85)

so the moiré Hamiltonian is transformed into the form on the atomic Wannier
function basis as

$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{k}}}_{2},\beta_{2}\rangle=\sum_{{\textbf{{R}}}_{1}{\textbf{{R}}}_{2}}e^{-i{\textbf{{k}}}_{1}\cdot{\textbf{{R}}}_{1}}\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{R}}}_{2},\beta_{2}\rangle
e^{i{\textbf{{k}}}_{2}\cdot{\textbf{{R}}}_{2}}.$ (S86)

Here
$\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{R}}}_{2},\beta_{2}\rangle$
describes the Hamiltonian matrix element between atomic Wannier functions
located at ${\textbf{{R}}}_{1}$ and ${\textbf{{R}}}_{2}$ in the superlattice
shown in Fig. S16(a). As the overlap between atomic Wannier functions decays
quickly as the distance increases, we only consider the local hopping within
the length scale
$|{\textbf{{R}}}_{2}-{\textbf{{R}}}_{1}|\sim\mathcal{O}(|\tilde{{\textbf{{a}}}}_{1}|)$,
where $\tilde{{\textbf{{a}}}}_{1}$ is the atomic primitive lattice vector for
the Sb2Te3 layer. In this atomic length scale, the Hamiltonian matrix element
between two Wannier orbitals near R on the superlattice structure with twist
angle $\theta$ in Fig. S16(a) can be approximated locally by the Hamiltonian
matrix element for two atomic layers with a constant shift
${\textbf{{d}}}_{\textbf{{R}}}$ in Fig. S16(b)Jung et al. (2014), where

${\textbf{{d}}}_{\textbf{{R}}}=\mathcal{R}(\theta){\textbf{{R}}}-{\textbf{{R}}}$
(S87)

and $\mathcal{R}(\theta)$ as the rotation operator for the Sb2 layer with the
rotating angle $\theta$. This approximation is valid for a small twist angle
$\theta$ because the local shift vector ${\textbf{{d}}}_{\textbf{{R}}}$ is
almost uniform at the atomic length scale,

${\textbf{{d}}}_{{\textbf{{R}}}_{2}}\approx{\textbf{{d}}}_{{\textbf{{R}}}_{1}}$
(S88)

for
$|{\textbf{{R}}}_{2}-{\textbf{{R}}}_{1}|\sim\mathcal{O}(|\tilde{{\textbf{{a}}}}_{1}|)$.
The Hamiltonian matrix element between two atomic Wannier orbitals for the
commensurate lattice is captured by
$\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}}_{{\textbf{{R}}}_{2}})|{\textbf{{R}}}_{2},\beta_{2}\rangle$,
so we make the approximation

$\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{R}}}_{2},\beta_{2}\rangle\approx\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}}_{{\textbf{{R}}}_{2}})|{\textbf{{R}}}_{2},\beta_{2}\rangle$
(S89)

and

$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{k}}}_{2},\beta_{2}\rangle=\sum_{{\textbf{{R}}}_{1}{\textbf{{R}}}_{2}}e^{-i{\textbf{{k}}}_{1}\cdot{\textbf{{R}}}_{1}}\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}}_{{\textbf{{R}}}_{2}})|{\textbf{{R}}}_{2},\beta_{2}\rangle
e^{i{\textbf{{k}}}_{2}\cdot{\textbf{{R}}}_{2}}.$ (S90)

To extract ${\textbf{{R}}}_{2}$ in
$\hat{H}^{\text{DFT}}({\textbf{{d}}}_{{\textbf{{R}}}_{2}})$ for the summation,
we transform $\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}})$ to the
momentum-space by

$\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}})=\sum_{\tilde{{\textbf{{G}}}}}e^{-i\tilde{{\textbf{{G}}}}\cdot{\textbf{{d}}}_{\textbf{{R}}}}\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}})$
(S91)

as
$\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}}+x\tilde{{\textbf{{a}}}}_{1}+y\tilde{{\textbf{{a}}}}_{2})=\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}})$
is periodic for atomic lattice vectors ($x,y$ are integers here), as shown in
Fig. S16(b). We also denote the atomic reciprocal lattice vector
$\tilde{{\textbf{{G}}}}=\tilde{{\textbf{{G}}}}_{wz}=w\tilde{{\textbf{{b}}}}_{1}+z\tilde{{\textbf{{b}}}}_{2}$
with integers $w,z$, so the summation over $\tilde{{\textbf{{G}}}}$ is
equivalent to the summation over $w,z$. $\tilde{{\textbf{{b}}}}_{1,2}$ are
atomic reciprocal lattice vectors satisfying
$\tilde{{\textbf{{b}}}}_{i}\cdot\tilde{{\textbf{{a}}}}_{j}=\delta_{ij}$ for
$i,j=1,2$ and shown in Fig. S16(c). Since

$\tilde{{\textbf{{G}}}}_{wz}\cdot{\textbf{{d}}}_{\textbf{{R}}}=\tilde{{\textbf{{G}}}}_{wz}\cdot(\mathcal{R}(\theta){\textbf{{R}}}-{\textbf{{R}}})=(\tilde{{\textbf{{G}}}}_{wz}-\mathcal{R}(\theta)\tilde{{\textbf{{G}}}}_{wz})\cdot(\mathcal{R}(\theta){\textbf{{R}}})={\textbf{{G}}}_{wz}\cdot({\textbf{{R}}}+{\textbf{{d}}}_{\textbf{{R}}})\approx{\textbf{{G}}}_{wz}\cdot{\textbf{{R}}}$
(S92)

with the moiré reciprocal lattice vectors ${\textbf{{G}}}_{wz}$ given by

$\begin{split}{\textbf{{G}}}_{wz}&=\tilde{{\textbf{{G}}}}_{wz}-\mathcal{R}(\theta)\tilde{{\textbf{{G}}}}_{wz}\\\
&=(w\tilde{{\textbf{{b}}}}_{1}+z\tilde{{\textbf{{b}}}}_{2})-\mathcal{R}(\theta)(w\tilde{{\textbf{{b}}}}_{1}+z\tilde{{\textbf{{b}}}}_{2})\\\
&=w(\tilde{{\textbf{{b}}}}_{1}-\mathcal{R}(\theta)\tilde{{\textbf{{b}}}}_{1})+z(\tilde{{\textbf{{b}}}}_{2}-\mathcal{R}(\theta)\tilde{{\textbf{{b}}}}_{2})\\\
&=w(\tilde{{\textbf{{b}}}}_{1}-\tilde{{\textbf{{b}}}}_{1}^{\prime})+z(\tilde{{\textbf{{b}}}}_{2}-\tilde{{\textbf{{b}}}}_{2}^{\prime})\\\
&=w{{\textbf{{b}}}_{1}^{\text{M}}}+z{{\textbf{{b}}}_{2}^{\text{M}}},\end{split}$
(S93)

we have

$\displaystyle
e^{-i\tilde{{\textbf{{G}}}}\cdot{\textbf{{d}}}_{\textbf{{R}}}}\approx
e^{-i{\textbf{{G}}}\cdot{\textbf{{R}}}},\quad\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}})\approx\sum_{w,z}e^{-i{\textbf{{G}}}_{wz}\cdot{\textbf{{R}}}}\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}}_{wz}).$
(S94)

Here
$\tilde{{\textbf{{a}}}}_{1,2}^{\prime}=\mathcal{R}(\theta)\tilde{{\textbf{{a}}}}_{1,2}\
(\tilde{{\textbf{{b}}}}_{1,2}^{\prime}=\mathcal{R}(\theta)\tilde{{\textbf{{b}}}}_{1,2})$
are primitive (reciprocal) lattice vectors for the twisted Sb2 layer as shown
in Fig. S16(c) and
${\textbf{{b}}}_{1,2}^{\text{M}}=\tilde{{\textbf{{b}}}}_{1,2}-\tilde{{\textbf{{b}}}}_{1,2}^{\prime}$
are the moiré reciprocal lattice vectors. The approximation in Eq. (S92) is
valid as
$|{\textbf{{d}}}_{\textbf{{R}}}|\sim\mathcal{O}(|\tilde{{\textbf{{a}}}}_{1}|)\ll|{\textbf{{a}}}_{1}^{\text{M}}|$
and
${\textbf{{G}}}\cdot{\textbf{{d}}}_{\textbf{{R}}}\ll{\textbf{{G}}}\cdot{\textbf{{a}}}_{1}^{\text{M}}\sim\mathcal{O}(1)$.

Substituting Eq. (S89) and (S94) into Eq. (S86) leads to

$\begin{split}\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{k}}}_{2},\beta_{2}\rangle&\approx\sum_{{\textbf{{R}}}_{1}{\textbf{{R}}}_{2},w,z}e^{-i{\textbf{{k}}}_{1}\cdot{\textbf{{R}}}_{1}}e^{-i{\textbf{{G}}}_{wz}\cdot{\textbf{{R}}}_{2}}\langle{\textbf{{R}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}}_{wz})|{\textbf{{R}}}_{2},\beta_{2}\rangle
e^{i{\textbf{{k}}}_{2}\cdot{\textbf{{R}}}_{2}}\\\
&=\sum_{w,z}\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}}_{wz})|{\textbf{{k}}}_{2}-{\textbf{{G}}}_{wz},\beta_{2}\rangle\\\
&=\sum_{w,z}\delta({\textbf{{k}}}_{2}-{\textbf{{k}}}_{1}-{\textbf{{G}}}_{wz})\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}}_{wz})|{\textbf{{k}}}_{1},\beta_{2}\rangle.\end{split}$
(S95)

The last line comes from the conservation of crystal momenta of
$H^{\text{DFT}}({\textbf{{d}}})$,

$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}})|{\textbf{{k}}}_{2},\beta_{2}\rangle=\int\mathrm{d}^{2}{\textbf{{d}}}e^{i\tilde{{\textbf{{G}}}}\cdot{\textbf{{d}}}}\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}({\textbf{{d}}})|{\textbf{{k}}}_{2},\beta_{2}\rangle=\delta({\textbf{{k}}}_{1}-{\textbf{{k}}}_{2})\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}})|{\textbf{{k}}}_{1},\beta_{2}\rangle.$
(S96)

Eq. (S95) connects $H^{\text{DFT}}$ and $H_{0}$ in atomic Bloch states in
general, which is applied to the moiré potential in our model next.

We next show the relation between the potential
$\tilde{\Delta}({\textbf{{d}}}_{\textbf{{R}}})$ from DFT and the moiré
potential $\Delta({\textbf{{r}}})$ in our model by Eq. (S95). The Fourier
transform of
$\langle{\textbf{{k}}},\alpha|\hat{H}^{\text{DFT}}({\textbf{{d}}}_{\textbf{{R}}})|{\textbf{{k}}},\alpha\rangle$
in Eq. (S81) by Eq. (S96) is

$\langle{\textbf{{k}}},\beta_{1}|\hat{H}^{\text{DFT}}(\tilde{{\textbf{{G}}}})|{\textbf{{k}}},\beta_{2}\rangle=\begin{pmatrix}h_{D}^{t}({\textbf{{k}}})&ms_{0}\\\
ms_{0}&h_{D}^{b}({\textbf{{k}}})\end{pmatrix}\delta_{\tilde{{\textbf{{G}}}}=0}+\begin{pmatrix}\tilde{\Delta}(\tilde{{\textbf{{G}}}})s_{0}&0\\\
0&\alpha\tilde{\Delta}(\tilde{{\textbf{{G}}}})s_{0}\end{pmatrix}.$ (S97)

In atomic Bloch basis, the moiré Hamiltonian from Eq.1 of the main text
without external electrical field is

$\langle{\textbf{{k}}}_{1},\beta_{1}|\hat{H}_{0}|{\textbf{{k}}}_{2},\beta_{2}\rangle=\begin{pmatrix}h_{D}^{t}({\textbf{{k}}}_{1})&ms_{0}\\\
ms_{0}&h_{D}^{b}({\textbf{{k}}}_{1})\end{pmatrix}\delta({\textbf{{k}}}_{2}-{\textbf{{k}}}_{1})+\begin{pmatrix}\Delta({\textbf{{G}}})s_{0}&0\\\
0&\alpha\Delta({\textbf{{G}}})s_{0}\end{pmatrix}\delta({\textbf{{k}}}_{2}-{\textbf{{k}}}_{1}-{\textbf{{G}}}).$
(S98)

By comparison of two Hamiltonian following Eq. (S95), one obtains

$\Delta({\textbf{{G}}}_{wz})=\tilde{\Delta}(\tilde{{\textbf{{G}}}}_{wz})$
(S99)

for $w,z$ as integers,
${\textbf{{G}}}_{wz}=w{\textbf{{b}}}_{1}^{\text{M}}+z{\textbf{{b}}}_{2}^{\text{M}}$,
and
$\tilde{{\textbf{{G}}}}_{wz}=w\tilde{{\textbf{{b}}}}_{1}+z\tilde{{\textbf{{b}}}}_{2}$.
In real space, this leads to

$\Delta({\textbf{{R}}})=\sum_{w,z}e^{-i{\textbf{{G}}}_{w,z}\cdot{\textbf{{R}}}}\Delta({\textbf{{G}}}_{wz})\approx\sum_{w,z}e^{-i\tilde{{\textbf{{G}}}}_{w,z}\cdot(\mathcal{R}(\theta){\textbf{{R}}}-{\textbf{{R}}})}\tilde{\Delta}(\tilde{{\textbf{{G}}}}_{wz})=\tilde{\Delta}(\mathcal{R}(\theta){\textbf{{R}}}-{\textbf{{R}}})=\tilde{\Delta}({\textbf{{d}}}_{\textbf{{R}}})$
(S100)

following Eq. (S92) and Eq. (S99), reproducing Eq. (4) in the main text. It
then can be interpolated to the whole real space by replacing the atomic
lattice vectors R by the continuous variable r,

$\Delta({\textbf{{r}}})\approx\tilde{\Delta}(\mathcal{R}(\theta){\textbf{{r}}}-{\textbf{{r}}}),$
(S101)

because the atomic length scale is much smaller than the moiré length scale
for small twist angles so that it is a good approximation to take the
continuous limit for the atomic length scale.