\section*{Supplementary Material}
Details of the computational methods, including the details for DFT calculations, and the DCD model; Results of the electrostatic potential for h-AlN, h-GaN and all monochalcogenides from DFT, DCD model and analytical formula; Potential fluctuations as a function of distance for all 2D materials from DFT, DCD model and analytical formula;
Structures of the 2D materials in crystallographic information file (cif) format. \\

\section*{Acknowledgments}
This work was supported by the Northwestern University MRSEC under National Science Foundation grant No.~DMR-1720139 (Q.Z., P.D.). Work performed at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, was supported by the U.S. DOE, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. We gratefully acknowledge use of the Bebop cluster in the Laboratory Computing Resource Center at Argonne National Laboratory. P.D. would like to acknowledge fruitful discussions with Cristian Cortès, Lincoln Lauhon, and Mark Hersam.


\section{Computational Methods}	\label{SecSI:method}
All first-principles calculations based on Density Functional Theory are performed using the Quantum-Espresso package~\cite{Giannozzi2009QE,Giannozzi2017QE}, with Optimized norm-conserving Vanderbilt pseudopotentials~\cite{Hamann2013} obtained from the PseudoDojo library~\cite{Van2018pseudodojo}. Exchange-correlation potentials use Perdew-Burke-Ernzerhof (PBE) parametrization of the generalized gradient approximation~\cite{Perdew1996}. The plane-wave cutoff energy is 90 Ry for all systems studied in this work. $k$-point sampling of 16$\times$16 is used for h-BN with in-plane size of 2.51$\times$2.51 \AA, or equivalent for other cell sizes. The vacuum regions are all larger than 21 \AA. Convergences of total and electronic energies are 10$^{-5}$ eV/atom, and $10^{-6}$ eV, respectively. 

\section{Analytical Derivation: Periodic 2D array of charges not contained within a plane}

The potential from a periodic charge configuration, where the unit cell is a rectangular lattice with lattice constants of $a, b$ is:
\begin{align}
V(x,y,z=d) = k\sum_{n_1,n_2,i} \frac{q_i}{\sqrt{(n_1a+x_i-x)^2+(n_2b+y_i-y)^2+d^2}} \quad \& \quad \sum_i q_i = 0 \label{eq:1}.
\end{align}

Using Fourier transform, see detailed derivation in our previous work\cite{Zhou2021SAM}, Eq.\ref{eq:1} becomes:
\begin{align} \label{eq0}
	V(x,y,z=d)

	=k_e \sum_i \frac{q_i}{a} \sum_{\vec{k}}{}'
	\frac{e^{-2\pi\left[ i k_1(x-x_i)/a+ i  k_2(y-y_i)/b+d \sqrt{k_1^2/a^2+k_2^2/b^2}\right]}}
	{\sqrt{k_1^2b^2/a^2+k_2^2}} 
\end{align}

We can extend this formulation to consider a neutral periodic 2D array of charges that are not contained within one plane, but have varying distances along the perpendicular direction. In this case Eq.~\ref{eq:1} becomes 
\begin{align}
V(x,y,z=d) = k_e\sum_{n_1,n_2,i} \frac{q_i}{\sqrt{(n_1a+x_i-x)^2+(n_2b+y_i-y)^2+(z_i-d)^2}} \quad \& \quad \sum_i q_i = 0 \label{eq:2}.
\end{align}
$\delta$ in Eq.S12 in the SI of our previous work\cite{Zhou2021SAM}  is now defined as $\delta = (d-z_i)/a$, which leads to a modified form of the final expression, eq.~\ref{eq:final}, as 
\begin{align}
V(x,y,z=d)  
&=k_e\sum_i \frac{q_i}{a} \left[\frac{2\pi z_i}{b} + \sum_{\vec{k}}{}' 
\frac{e^{-2\pi\left[ i k_1(x-x_i)/a+ i  k_2(y-y_i)/b+(d-z_i) \sqrt{k_1^2/a^2+k_2^2/b^2}\right]}}
{\sqrt{k_1^2b^2/a^2+k_2^2}}\right]. \label{eq:final}
\end{align}
We note that in the case where the charge distribution yields a net dipole the first term is the potential  due to a periodic array of dipoles in the far-field limit, which can be equated to an infinite parallel-plate capacitor. Note in order to recover the expected potential difference, we need to calculate the potential different between $z=\pm d$, which picks up a factor of two.

\subsection{Formula Derivation for hexagonal BN, AlN and GaN}
\begin{figure}[H]
	\centering
	\includegraphics[width=0.35\linewidth]{Figures/BNcell.jpg}
	\caption{An Orthorombic Cell of monolayer h-BN, same for AlN and GaN where the N atoms are in gray.} \label{fig:BNcell}
\end{figure}
As shown in Fig.~\ref{fig:BNcell}, each orthorombic cell of h-BN (same for 2D AlN and GaN) with lattice parameters $a, b$, has four atoms (two B atoms and two N atoms), two B atoms with charge $q$ at position $\boldsymbol r(x,y,z=0)$ of  $(0,0,0), (-1/2a,-1/2b,0)$ and two N atoms with charge $-q$ at $(0,1/3b,0), (-1/2a, -1/6b, 0)$ when defining the h-BN plane at $z=0$. The potential at position $(x,y,z=d)$ is then (\emph{hereinafter, $x,y$ are in fractional coordinates, and $d$ has the same unit as $a$ and $b$}):
\begin{align}
	V(x,y,z=d) & = k_e \frac{q}{ab} \sum_{\vec{k}}{}' \frac{e^{-2\pi d M_k}}{M_k} \left\{ \cos[2\pi (k_1x+k_2y)] +  \cos[2\pi(k_1(x-\frac{1}{2})+k_2(y-\frac{1}{2})]\right. \nonumber \\ 
	& \left. {\qquad} - \cos[2\pi(k_1x+k_2(y-\frac{1}{3}))] + \cos[2\pi k_1(x-\frac{1}{2})+2\pi k_2(y-\frac{5}{6})] \right\} \\
	& = k_e\frac{4q}{ab} \sum_{\vec{k}}^{k_1+k_2=2n}{}' \frac{e^{-2\pi d \sqrt{k_1^2/a^2+k_2^2/b^2}}}{\sqrt{k_1^2/a^2+k_2^2/b^2}} \left[- \sin(2\pi k_1x + 2\pi k_2y -\frac{1}{3}\pi k_2)\sin(\frac{\pi k_2}{3})\right]
\end{align}

We could achieve $>$ 99\% accuracy with $k_1=[-1,1], k_2=[-2,2]$, and $b=\sqrt{3}a$. Therefore,
\begin{align}\label{eq:BN}
		V(x,y,d) & =k_e\frac{2A_1q}{a} e^{\frac{-4\pi d}{\sqrt{3}a}}
\end{align}
where $A_1=-2\sqrt{3}\sin(2\pi y-\frac{\pi}{3})\left[\cos(2\pi x)+\cos(2\pi y-\frac{\pi}{3})\right]$

The in-plane potential modulation defined as $\Delta V$ is:
\begin{align}
	\Delta V & = V_{max} - V_{min} = V(0,1/3b,d) - V(0,0,d)  \nonumber \\
	& = k_e\frac{18q}{a} e^{\frac{-4\pi d}{\sqrt{3}a}}
\end{align}

In the case of a rotation of the h-BN plane in angle $\theta$, the potential keeps the same form as Eq.~\ref{eq:BN}, but with a angle dependent $A_1(\theta)$, therefore
\begin{equation}
V(x,y,d,\theta) = k_e\frac{2A_1(\theta)q}{a}e^{\frac{-4\pi d}{\sqrt{3}a}}
\end{equation} 
where
\begin{equation}
	A_1(\theta)=-2\sqrt{3}\sin(2\pi y(\theta)-\frac{\pi}{3})\left[\cos(2\pi x(\theta))+\cos(2\pi y(\theta)-\frac{\pi}{3})\right] \nonumber
\end{equation}
and $x(\theta)= x\cos \theta +\frac{b}{a}y\sin \theta, y(\theta)=-\frac{a}{b}x\sin \theta + y\cos \theta$

\subsection{Formula Derivation for Hexagonal BN, AlN and GaN with Buckling}

For monolayer hexagonal structures, we consider here the case when the atoms are not planar, i.e.two point charges $q$ at position $\boldsymbol r(x,y,z)$ of  $(-1/2a,-1/b,z_1), (0,0,z_2)$ and two N atoms with charge $-q$ at $(-1/2a, -1/6b,z_3), (0,1/3b, z_4)$. when defining the average position of all the atoms in $z$ direction as $z=0$, we have:
\begin{align}\label{eq:BNz}
V(x,y,-d) 
& = k_e \frac{q}{ab}\sum_{\vec{k}}{}' cos[2\pi(k_1x+k_2y)] \left((-1)^{k_1+k_2} e^{-2\pi z_1 M_k} +  e^{-2\pi z_2M_k}\right) \frac{e^{-2\pi dM_k}}{M_k} \nonumber \\
& - k_e \frac{q}{ ab}\sum_{\vec{k}}{}' cos[2\pi k_1x + 2\pi k_2(y-\frac{1}{3})] \left( (-1)^{k_1+k_2} e^{-2\pi z_3 M_k} +e^{-2\pi z_4M_k}\right) \frac{e^{-2\pi dM_k}}{M_k} \nonumber \\ 
& + k_e\frac{q}{ab}(z_1+z_2-z_3-z_4)
\end{align}
where the last term in Eq.~\ref{eq:BNz} is correcting for the dipole effects. We then can derive the analytical formula for two different cases here: 

\noindent \textbf{Case \Romannum{1}: $z_1=z_2=-z_3=-z_4=\Delta z$,}
\begin{align}
V(x,y,z=-d) &=  k_e \frac{q}{ab}\sum_{\vec{k}}{}' \frac{e^{-2\pi d M_k}}{M_k}
\{ 
\cos\left[2\pi (k_1x+k_2y)\right]e^{-2\pi \Delta zM_k} -  
\cos\left[2\pi (k_1x+k_2(y-\frac{1}{3}))\right] e^{2\pi \Delta zM_k} 
\}
\end{align}
Using Taylor expansion, we can achieve $>95$\% accuracy to the first order, and $>$99\% accuracy to the second order (Eq.~\ref{eq:z195}), therefore
\begin{align}
V(x,y,z=-d)=& k_e \frac{q}{a}e^{-\frac{4\pi d}{\sqrt{3}a}} 
\{
(A_1+A_1^{\prime}) e^{-\frac{4\pi \Delta z}{\sqrt{3}a}} 
- (A_1^{\prime}-A_1) e^{\frac{4\pi \Delta z}{\sqrt{3}a}}
\} \\
\approxeq & 
\begin{cases}
k_e\frac{2q}{a}\left[ A_1(1+\frac{8\pi^2 \Delta z^2}{3a^2}) - A_1^{\prime}\frac{4\pi \Delta z}{\sqrt{3}a} \right] e^{-\frac{4\pi d}{\sqrt{3}a}} \\%  \label{eq:z199} \\
k_e\frac{2q}{a}\left[ A_1- A_1^{\prime}\frac{4\pi \Delta z}{\sqrt{3}a} \right] e^{-\frac{4\pi d}{\sqrt{3}a}}  \label{eq:z195} 
\end{cases} 
\end{align}
where $A_1 = -2\sqrt{3}\sin(2\pi y-\frac{\pi}{3})\left[\cos(2\pi x)+\cos(2\pi y-\frac{\pi}{3}) \right]$, and\\
$A_1^{\prime} = 2\cos(2\pi y-\frac{\pi}{3})\left[\cos(2\pi x)-\cos(2\pi y-\frac{\pi}{3}) \right]+1$. 


\vspace{12pt}
\noindent \textbf{Case \Romannum{2}: $z_1=-z_2=-z_3=z_4=\Delta z$,}
\begin{align}
	V(x,y,z=-d) = &k_e\frac{q}{ab} \sum_{\vec{k}}{}' \frac{e^{-2\pi dM_k}}{M_k} \left[e^{2\pi \Delta zM_k} + (-1)^{k_1+k_2}e^{-2\pi \Delta zM_k} \right] \nonumber \\
	& \left\{ \cos \left[2\pi (k_1x+ k_2 y)\right]-(-1)^{k_1+k_2} \cos \left[2\pi(k_1x+k_2(y-\frac{1}{3})) \right] \right\}
\end{align}
We can achieve $>99$\% accuracy with $-1\leq k_1 \leq 1, -2 \leq k_2 \leq 2$ and to second order Taylor expansion, therefore
\begin{align}
V(x,y,z=-d) = k_e \frac{2q}{a}\left[A_1 \left(1+\frac{8\pi^2 \Delta z^2}{3a^2}\right)e^{\frac{-4\pi d}{\sqrt{3}a}} + A_2 \frac{\Delta z}{a} e^{\frac{-2\pi d}{\sqrt{3}a}}  + A_3 \frac{\Delta z}{a}e^{\frac{-2\pi d}{a}}  +A_4 \frac{\Delta z}{a} e^{\frac{-2\sqrt{7}\pi d}{\sqrt{3}a}} \right]
\end{align}

and $>98$\% accuracy to the first order in Taylor expansion, where
\begin{align}
V(x,y,z=-d) = k_e \frac{2q}{a}\left[A_1e^{\frac{-4\pi d}{\sqrt{3}a}} + A_2 \frac{\Delta z}{a} e^{\frac{-2\pi d}{\sqrt{3}a}}  + A_3 \frac{\Delta z}{a}e^{\frac{-2\pi d}{a}}  +A_4 \frac{\Delta z}{a} e^{\frac{-2\sqrt{7}\pi d}{\sqrt{3}a}} \right]
\end{align}

where $A_2=4\pi \cos(2\pi y-\frac{\pi}{3})/\sqrt{3}, A_3 = 8\pi \cos(2\pi x)/\sqrt{3}, A_4=8\pi \cos(2\pi x)\cos(4\pi y-\frac{2\pi}{3})$.


\subsection{Formula For Metal MonoChalcogenides}

\begin{figure}[H]
	\centering
	\includegraphics[width=0.4\linewidth]{Figures/MMCq.jpg}
	\caption{(a) top view and (b) side view of the unit cell and the positions of the point charges for metal monochalcogenides shown as black dots.} \label{fig:MMCq}
\end{figure}

Monolayer of metal monochalcogenides (MMC) MX where M=Ge, Sn, Pb and X=S,Se,Te has four atoms per unit cell as shown in Fig.~\ref{fig:MMCq} (b). We can use two point charges $\hat{r}_1 (-x_0a, -y_0b, -\Delta z)$ and $\hat{r}_2 (x_0a, y_0b, \Delta z)$ with charges of $q$ and $-q$, respectively, to represent each MX as shown in Fig.~\ref{fig:MMCq}.
Since the total dipole moment is zero, the potential is\\
\begin{align}
	V(x,y,z=d) &=k_e\sum_i \frac{q_i}{ab}  \sum_{\vec{k}}{}' \frac{e^{-2\pi (d-z_i)M_k }}{M_k} cos\left\{2\pi \left[k_1(x-x_i)+k_2(y-y_i)\right]\right\} \\ 
	& = k_e\frac{q}{ab} \sum_{\vec{k}}{}' \frac{e^{-2\pi dM_k}}{M_k} \left\{  \cos\{2\pi\left[k_1(x+x_0)+k_2(y+y_0)\right]\}e^{2\pi \Delta z M_k} \right.  \nonumber \\
	& \qquad \qquad \qquad \qquad  \left. {} - \cos\{2\pi\left[k_1(x-x_0)+k_2(y-y_0)\right]\}e^{-2\pi \Delta z M_k}
	\right\}	
\end{align}

As we can achieve accuracy $> 98$\% with $-1 \leq k_1, k_2  \leq 1$, therefore
\begin{align}\label{eq:mmc}
V(x,y,z=d) = & k_e \frac{q}{a} \left\{ \frac{4\left[C_x C_y e^{2\pi \Delta z \sqrt{1/a^2+1/b^2}} - D_x D_y e^{-2\pi \Delta z\sqrt{1/a^2+1/b^2}} \right]}{\sqrt{b^2/a^2+1}} e^{-2\pi d\sqrt{1/a^2+1/b^2}} \nonumber \right. \\ 
& \left. {\qquad} + 2  \left[C_y e^{2\pi \Delta z/b} - D_y e^{-2\pi \Delta z/b} \right]e^{\frac{-2\pi d}{b}} \right. \nonumber \\
& \left. {\qquad} + \frac{2a }{b} \left[C_x e^{2\pi \Delta z/a} - D_x e^{-2\pi \Delta z/a} \right]e^{\frac{-2\pi d}{a}}  \right\}
\end{align}
where $C_x=\cos[2\pi (x+x_0)], D_x=\cos[2\pi (x-x_0)],C_y=\cos[2\pi (y+y_0)],D_y=\cos[2\pi (y-y_0)]$.
We can achieve $>98$\% accuracy using Taylor expansion to the second order (Eq.~\ref{eq:mmc2}), and $>95$\% to the first order(Eq.~\ref{eq:mmc1}), we have
\begin{align} \label{eq:mmc2}
V(x,y,d) & = k_e \frac{q}{a}  \left\{ B_1(x,y,\Delta z) e^{-2\pi d\sqrt{1/a^2+1/b^2}} \nonumber \right. \\
& \left. + B_2(x,y,\Delta z) e^{\frac{-2\pi d}{b}} + B_3(x,y,\Delta z) e^{\frac{-2\pi d}{a}}  \right\}
\end{align}
\begin{align} \label{eq:mmc1}
V(x,y,z=d) = & k_e \frac{q}{a} \left\{ \frac{B_1+ B_1^{\prime} 2\pi \sqrt{\frac{b^2}{a^2}+1}\frac{\Delta z}{b} }{\sqrt{b^2/a^2+1}} e^{-2\pi d\sqrt{1/a^2+1/b^2}} \nonumber \right. \\ 
& \left. {\quad} + \left[B_2+B_2^{\prime}2\pi \frac{\Delta z}{b}\right]e^{\frac{-2\pi d}{b}} +  \frac{a}{b}\left[B_3+B_3^{\prime} 2\pi\frac{\Delta z}{a}\right] e^{\frac{-2\pi d}{a}} \right\} 
\end{align}
where $B_1 = 4(C_xC_y-D_xD_y), B_1^{\prime}=4(C_xC_y+D_xD_y), B_2=2(C_y-D_y), B_2^{\prime}=2(C_y+D_y), B_3 = 2(C_x-D_x), B_3^{\prime}= 2(C_x+D_x)$.


\subsection{Analytical Formula for Moir\'{e} Structures}
For a unit cell rotated by angle $\theta$, the potential can be computed by replacing $x,y$ in Eq.\ref{eq:BN} and Eq.\ref{eq:mmc2} by (note that $x,y$ are in fractional coordinates of the cell vectors):
\begin{align}
    & x(\theta)=x\cos \theta + \frac{b}{a}y\sin \theta \nonumber \\
    & y(\theta) = -\frac{a}{b}x sin \theta + y\sin \theta
\end{align}
Therefore, by adding the potential with the above $x(\theta), y(\theta)$, we can obtain angle-specific analytical formula for moir\'{e} structures. For example, for two layers of h-BN in the configuration as shown in Fig. 4(a), the potential at the center 2D material can be computed as:
\begin{align}
    V(x,y,d_1,d_2,\theta)=k_e\frac{2q}{a}\left[A(x,y)e^{\frac{-4\pi d_1}{\sqrt{3}a}}+A(x(\theta),y(\theta))e^{\frac{-4\pi d_2}{\sqrt{3}a}}\right]
\end{align}
where $d_1, d_2$ are the distance of the two layers of h-BN, which are considered as equal in Fig.4. Same can be obtained for transition metal monochalcogenides using Eq.\ref{eq:mmc2}.

Meanwhile, using the analytical formula for the potential, we can derive their electric field by:
\begin{align}
    \vec{F}(x,y,z,\theta) = -\left[\frac{\partial V(x,y,d,\theta)}{\partial x} \vec{x} + \frac{\partial V(x,y,d,\theta)}{\partial y} \vec{y} +\frac{\partial V(x,y,d,\theta)}{\partial z} \vec{z}\right]
\end{align}

\section{Lattice Parameters and Charges}

\begin{table}[H]
	\centering
	\captionsetup{margin={1.0cm,1.0cm}}
	\caption{\ lattice parameters and the absolute magnitude of the point charges $q$ for each 2D material studied in this work. The positions for the positive and negative point charges are shown in Fig.3(d), Fig.\ref{figs:AlNz}(d), Fig.\ref{figs:GaNz}(d) for the hexagonal monolayer, and Fig.\ref{figs:MMCall} for MMCs.} 
	\label{tbl:qstr}
	\begin{tabular*}{0.7\textwidth}{@{\extracolsep{\fill}}cccc}
		\hline
		2D & a (\AA) & b (\AA) & $|q|$ \\ \hline
		h-BN	& 2.51	& 2.51	& 0.89
 \\
		h-AlN	& 3.13	& 3.13	& 1.47
\\
		h-GaN	& 3.21	& 3.21	& 1.18
\\
	    GeS	& 3.65	& 4.50	& 0.23
\\
		GeSe	& 3.97	& 4.30	& 0.25
\\
		GeTe	& 4.39	& 4.24	& 0.18
\\
		SnS	& 4.09	& 4.30	& 0.26
\\
		SnSe	& 4.29	& 4.39	& 0.30
\\
		SnTe	& 4.55	& 4.58	& 0.25
\\

		PbS	& 4.25	& 4.25	& 0.50
\\
		PbSe	& 4.41	& 4.41	& 0.48
\\
		PbTe	& 4.64	& 4.64	& 0.36\\
		\hline
	\end{tabular*} 
\end{table}	

\section{Potential Results for h-AlN and h-GaN}
\begin{figure}[H]
	\centering
	\includegraphics[width=0.9\linewidth, center]{Figures/FigAlNz.pdf}
	\caption{In-plane electrostatic potential at $d=3.03$ \AA from single layer of h-AlN computed from (a) DFT, (d) the DCD model, (g) formula. The left two columns show the potential for h-AlN with Raman-active (second column) and infrared-active (third column) phonon mode, see atomic configurations on the top, computed from DFT, DCD model and formula, respectively.} \label{figs:AlNz}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=0.9\linewidth, center]{Figures/FigGaNz.pdf}
	\caption{In-plane electrostatic potential at $d=3.03$ \AA from single layer of h-GaN computed from (a) DFT, (d) the DCD model, (g) formula. The left two columns show the potential for h-GaN with Raman-active (second column) and infrared-active (third column) phonon mode, see atomic configurations on the top, computed from DFT, DCD model and formula, respectively.} \label{figs:GaNz}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=0.9\linewidth, center]{Figures/FigBNVdline.pdf}
	\caption{$\Delta V_{pot}=V_{max}-V_{min}$ as a function of $d$ for (a) h-BN, (b) h-AlN, (c) h-GaN computed from DFT, model and formula.} \label{figs:BNdline}
\end{figure}

\section{Potential Results for all Transition Metal MonoChalcogenides}
\begin{figure}[H]

	\subfloat{\includegraphics[width=0.99\linewidth,center]{Figures/FigGeSnX1d_top.pdf}} \\%\hspace{0.5pt}
	\subfloat{\includegraphics[width=0.5\linewidth,center]{Figures/FigPbX1d_top.pdf}}
	\caption{Potential from transition metal monochalcogenides at distances of $\sim$ 3.1 \AA above computed using DFT, model and analytical formula in Eq.\ref{eq:mmc2}. The black dots show the positions of the bottom two atoms. The plus and minus symbols mark the positions of the postive and negative charges. } \label{figs:MMCall}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=0.9\linewidth, center]{Figures/FigMMCVdline.pdf}
	\caption{$\Delta V_{pot}=V_{max}-V_{min}$ as a function of $d$ for MMCs computed from DFT, model and formula.} \label{figs:MMCdline}
\end{figure}

\section{Impacts of Stacking and Twist Angle on the Potential and Electric Field for MMCs}
\begin{figure}[H]
	\centering
	\includegraphics[width=0.9\linewidth, center]{Figures/FigSVF0.pdf}
	\caption{(a-d) Potential at the center 2D materials from two layers of SnTe $d=3.18$\AA\ above and below, see the configuration in Fig.4(a), with twist angle $\theta=0$ but different stackings where $\textbf{R}$ is the in-plane translation of the cell in fractional coordinates. (e-h) in-plane electric field where the length and direction of the arrows show the magnitude and direction of the electric field, respectively.} \label{figs:MMCt0}
\end{figure}

\begin{figure}[H]
	\centering
	\includegraphics[width=0.6\linewidth, center]{Figures/FigSVF20.pdf}
	\caption{In-plane electric field distribution for a moir\'{e} structure. For a heterostructure as shown in Fig.4 (a) for SnTe with a twist angle $\theta=20^o$, the (a) potential and (b) in-plane electric field.} \label{figs:MMCt20}
\end{figure}




