# A Dirac-fermion approach and its application to design high Chern numbers in
magnetic topological insulator multilayers

Dinghui Wang1,#, Huaiqiang Wang1,#, and Haijun Zhang1,2,∗ 1 National
Laboratory of Solid State Microstructures, School of Physics, Nanjing
University, Nanjing 210093, China
2 Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing 210093, China
<EMAIL_ADDRESS>
###### Abstract

Quantum anomalous Hall (QAH) insulators host topologically protected
dissipationless chiral edge states, the number of which is determined by its
Chern number. Up to now, the QAH state has been realized in a few magnetic
topological insulators, but usually with a low Chern number. Here, we develop
a Dirac-fermion approach which is valuable to understand and design high Chern
numbers in various multilayers of layered magnetic topological insulators.
Based on the Dirac-fermion approach, we demonstrate how to understand and tune
high Chern numbers in ferromagentic MnBi2Te4 films through the van der Waals
(vdW) gap modulation. Further, we also employ the Dirac-fermion approach to
understand the experimentally observed high Chern numbers and topological
phase transition from the Chern number $C=2$ to $C=1$ in the
[3QL-(Bi,Sb)1.76Cr0.24Te3]/[4QL-(Bi,Sb)2Te3] multilayers. Our work provides a
powerful tool to design the QAH states with a high Chern number in layered
magnetic topological insulator multilayers.

##  I Introduction

The discovery of topological insulators (TIs) Hasan and Kane (2010); Qi and
Zhang (2011); Bernevig  _et al._ (2006); König  _et al._ (2007); Zhang  _et
al._ (2009); Chen  _et al._ (2009) encourages researchers to study Chern
insulators with a nonzero Chern number. The Chern number, also known as the
Thouless-Kohmoto-Nightingale-Nijs number in condensed matter physics, serves
as a topological invariant to characterize the quantum anomalous Hall (QAH)
state Thouless  _et al._ (1982). In a Chern insulator, there are
topologically-protected gapless chiral edge states Qi  _et al._ (2008); Weng
_et al._ (2015); Liu  _et al._ (2016), the number of which is determined by
its Chern number through the bulk-boundary correspondence. To achieve a
nonzero Chern number, the time-reversal symmetry must be broken. Haldane first
creatively proposed an honeycomb lattice model with a periodic magnetic flux
to realize the QAH state without an external magnetic field Haldane (1988),
though it is still hard to implement in experiments. Interestingly, the QAH
effect has been realized in some recent state-of-the-art experiments, though
the working temperature is still quite low (up to only a few Kelvin). In those
reported experiments, the Chern number $C=1$ is widely observed Liu  _et al._
(2016), for example, in Cr- or V-doped (Bi,Sb)2Te3 films Zhang  _et al._
(2013); Chang  _et al._ (2013); Kou  _et al._ (2014); Checkelsky  _et al._
(2014); Xiao  _et al._ (2018); Grauer  _et al._ (2017); Mogi  _et al._ (2017);
Bestwick  _et al._ (2015), and the recently discovered antiferromagnetic (AFM)
TI MnBi2Te4 films Liu  _et al._ (2020); Deng  _et al._ (2020). Besides, the
QAH effect with the Chern number $C=1$ is also observed in twisted bilayer
graphene, driven by the strong interaction Serlin  _et al._ (2020); Sharpe
_et al._ (2019).

Figure 1: (a) Schematic of a 3D topological insulator (TI) with two Dirac-cone
surface states on the top and bottom surfaces. (b) When the 3D TI is reduced
to a thin film, the two Dirac-cone surface states will gradually couple with
each other, which is defined as a block layer (BL) with the
$\Delta_{\textrm{intra}}$ denoting the intra-layer coupling. (c) Illustration
of a multilayer structure constructed through stacking BLs in (b). The wave
function profile from $k\cdot p$ analysis of the low energy state with
reducing film thickness for (d) the topologically trivial case and (e)
topologically nontrivial case. For the topologically nontrivial (trivial)
case, the wave function indicates that the low energy state is a surface
(bulk) state. (f,g,h) When the thickness of a TI film is gradually reduced,
the Dirac-cone surface states get to couple with each other, seen in (f), (g),
and (h), where the pink peaks schematically indicate the projected Dirac-cone
surface states. The Dirac-cone states can be distinguished even at the ultra-
thin limit (h), which is essentially different from that in the trivial case
(i).

It is well known that the chiral edge states in Chern insulators can carry
dissipationless electric currents, which give rise to a quantized Hall
conductivity $\sigma_{xy}=C\frac{e^{2}}{h}$, where the $C$ denotes the Chern
number. Therefore, the Chen insulator with a high Chern number ($C>1$) is not
only valuable for the fundamental topological physics, but also promising to
improve the performance in next-generation electronic devices. To this end,
Wang et al. Wang  _et al._ (2013) proposed the superlattices of Cr-doped
Bi2(SexTe1-x)3 TI for Chern insulators with $C=2$, and Fang et al. Fang  _et
al._ (2014) proposed to reach Chern numbers between $C=-4$ and $C=4$ in
topological crystalline insulator thin films through the ferromagnetic (FM)
doping. Jiang et al. Jiang  _et al._ (2012) proposed realizing higher Chern
numbers through doping magnetic elements in a multilayer of TI films.
Recently, a Chern number $C=3$ was predicted in double septuple-layer (SL)
MnBi2Te4 thin films under an external electric field Du  _et al._ (2020), and
$C=N$ was predicted in $2N$-SL MnBi2Te4 with alternate stacking orders Zhu
_et al._ (2022). In experiments, a high Chern number (up to $C=5$) has been
realized in [(Bi,Sb)1.76Cr0.24Te3]/[(Bi,Sb)2Te3] multilayers based on the
strategy of direct stacking Zhao  _et al._ (2020). However, Chern insulators
with a high Chern number are still very rare and highly desired.

In this work, we develop an effective Dirac-fermion approach to investigate
the low-energy electronic structure and design high Chern numbers in
multilayers of layered magnetic TIs. In Sec. II, we first show the stable
existence of local Dirac-cone surface states in TI films, even at the ultra-
thin limit. Then, we build the Dirac-fermion approach to describe the low-
energy electronic structures of a large variety of multilayers of layered
(magnetic) TIs. In Sec. III, we demonstrate the topological phase transition
with the Chern number from $C=0$ to $C=1$ in one magnetic block layer (BL),
depending on the competition of the magnetic exchange field and the intra-
layer coupling. We further show how the highest Chern number $C=N$ can be
obtained in a film of $N$ magnetic BLs. As applications of the Dirac-fermion
approach, in Sec. IV, we study the FM MnBi2Te4 films with a high Chern number
predicted. In Sec. V, we further employ the Dirac-fermion approach to
understand the experimentally observed topological phase transitions in [3QL
(Bi,Sb)1.76Cr0.24Te3]/[4QL (Bi,Sb)2Te3] multilayers. Finally, a summary is
given in Sec. IV.

##  II Dirac-fermion approach

TIs are insulating in the bulk but metallic on the surface due to the
existence of gapless Dirac-cone surface states, guaranteed by the topological
invariant $Z_{2}$ Fu  _et al._ (2007); Moore and Balents (2007); Roy (2009),
schematically shown in Fig. 1(a). Here, we take the typical layered TI Bi2Se3
Zhang  _et al._ (2009) as an example to demonstrate how to build the Dirac-
fermion approach for multilayers of layered (magnetic) TIs. It is well known
that Bi2Se3 has a large bulk band gap, in which a gapless helical spin-
momentum locked Dirac-cone surface state arises, as directly observed by angle
resolved photoemission spectroscopy experiments Chen  _et al._ (2009); Xia
_et al._ (2009). The low-energy effective $k\cdot p$ Hamiltonian for the
helical Dirac-cone surface state to linear order in $k$ is given by Liu  _et
al._ (2010a)

$H_{D}(\textbf{k})=\pm v(k_{x}\sigma_{y}-k_{y}\sigma_{x}),$ (1)

where $v$ is the Fermi velocity, the Pauli matrices act on the spin subspace,
and ’$\pm$’ denote the opposite helicities of Dirac states localized on
opposite surfaces. Here, the higher-order hexagonal warping effect Fu (2009)
has been ignored, for we are just interested in the low-energy electronic
structure around the $\Gamma$ point. Because the Dirac-cone surface state
localizes on the surface layers, for example, 2-quintuple-layer (2-QL) region
for Bi2Se3 Zhang  _et al._ (2010), when a TI film is thick enough, the top and
bottom surface states get completely decoupled, as schematically shown in Fig.
1(a). However, along with reducing the film thickness, the top and bottom
surface states will gradually couple with each other with a finite intra-layer
coupling $\Delta_{\textrm{intra}}$, as schematically shown in Fig. 1(b).

A natural question is whether the Dirac-cone surface states stably exist, when
the TI film gets to its ultra-thin limit (e.g. 1 QL Bi2Se3 film)? To address
this issue, we begin with the four-band bulk $k\cdot p$ model of three-
dimensional Bi2Se3-family TIs Zhang  _et al._ (2009),

$H_{\mathrm{bulk}}(\textbf{k})=\begin{pmatrix}M(\textbf{k})&A_{1}k_{z}&0&A_{2}k_{-}\\\
A_{1}k_{z}&-M(\textbf{k})&A_{2}k_{-}&0\\\
0&A_{2}k_{+}&M(\textbf{k})&-A_{1}k_{z}\\\
A_{2}k_{+}&0&-A_{1}k_{z}&-M(\textbf{k})\end{pmatrix}+\epsilon_{0}(\textbf{k}).$
(2)

where $k_{\pm}=k_{x}\pm ik_{y}$,
$\epsilon_{0}=C+D_{1}k_{z}^{2}+D_{2}k_{\parallel}^{2}$ and
$M(\textbf{k})=M+B_{1}k_{z}^{2}+B_{2}k_{\parallel}^{2}$. The parameters of
$A_{1,2},$ $C$, $D_{1,2}$, $M$ and $B_{1,2}$ can be taken from the previous
work on Bi2Se3 Zhang  _et al._ (2009); Liu  _et al._ (2010a), in which the
curvature parameters $B_{1,2}$ have the same positive sign ($B_{1,2}$ is
henceforth set to be positive). In the above $k\cdot p$ model, $M<0$ indicates
an inverted band structure at $\Gamma$, and describes a Bi2Se3-like TI, while
$M>0$ describes a topologically trivial insulator. In order to see the
evolution of Dirac-cone surface states along with reducing the thickness of a
TI film, we employ the method based on the quantum well states to calculate
the energy spectrum and wavefunction distributions Liu  _et al._ (2010b).

For the topologically trivial case ($M>0$), the wave functions of the low-
energy state are calculated for the films with reducing the thickness, shown
in Fig. 1(d). It is clear that the low-energy state is not a surface state but
a bulk state, which is not surprising because the $k\cdot p$ model with $M>0$
describes a topologically trivial insulator without topologically protected
surface state. Differently, for the topologically nontrivial case ($M<0$), the
calculated wave functions of the low-energy state are shown in Fig. 1(e). As
expected, surface states are clearly seen to exist locally on the top and
bottom surfaces when the film is thick enough. When reducing the film
thicknessZhou  _et al._ (2008), the top and bottom surface states gradually
merge with a finite intra-layer coupling $\Delta_{\textrm{intra}}$. For the
weak coupling, the top and bottom surface states can be visibly distinguished,
as shown in Fig. 1(f). When the intra-layer coupling further increases, the
two surface states get to significantly merge with each other, as shown in
Fig. 1(g). Especially, when the intra-layer coupling is strong enough, for
example, at the ultra-thin limit of 1-QL Bi2Se3, the wave functions seem to
have one peak, but the surface states can still be decomposed (see Appendix
A), as shown in Fig. 1(h), which is essentially different from the single peak
of wave functions in the topologically trivial case ($M>0$) in Fig. 1(i).
Therefore, the results indicate the stable existence of Dirac-cone surface
states on both surfaces for TI thin films, as long as its bulk Hamiltonian has
$M<0$, which guarantees the validity to build the Dirac-fermion approach for
various multilayers of layered (magnetic) TIs. It is also worth mentioning
that the previous work Shan  _et al._ (2010) has strictly derived the local
Dirac-cone surface states for TI films, which remains valid down to the ultra-
thin limit.

Then, we show how to build a Dirac-fermion model for a multilayer of TIs.
First of all, we take a TI film as one BL, which can consist of one or several
unit layers. For Bi2Se3, the unit layer is a QL (Se-Bi-Se-Bi-Se). For each BL
with two Dirac-cone surface states coming from its top and bottom surfaces,
respectively, the effective low-energy Hamiltonian by the Dirac-fermion model
is written as

$H_{\textrm{BL}}(\textbf{k})=v(k_{x}\sigma_{y}-k_{y}\sigma_{x})\tau_{z}+\Delta_{\textrm{intra}}\tau_{x},$
(3)

where $\tau_{x,z}$ are Pauli matrices acting in the subspaces of top and
bottom surface states. Here, $\tau_{z}$ denotes the opposite helicity of the
top and bottom Dirac-cone states, and $\Delta_{\textrm{intra}}$ denotes the
intra-layer coupling of the two Dirac-cone states within each BL. Strictly
speaking, apart from the constant coupling term, the $\Delta_{\textrm{intra}}$
also contains $k$-dependent terms, such as the quadratic terms derived in
Refs. Lu  _et al._ (2010); Yu  _et al._ (2010). However, since we are mainly
concerned with the low-energy physics around the $\Gamma$ point, we treat the
$\Delta_{\textrm{intra}}$ as a constant term and neglect these higher-order
terms. Moreover, although it has been shown that the $\Delta_{\textrm{intra}}$
may even change its sign as the BL’s thickness reduces Lu  _et al._ (2010), it
is found that its sign does not lead to essential difference for the following
discussion, so the $\Delta_{\textrm{intra}}$ is set to be positive unless we
specifically state.

Furthermore, TI BLs can be stacked into a multilayer, as schematically shown
in Fig. 1(c). Here, we consider the BLs are coupled by a weak van der Waals
(vdW) interaction, which is valid for layered TIs. Naturally, the bottom
Dirac-cone state of each BL couples with the top Dirac-cone state of its below
adjacent BL. We take the parameter $\Delta_{\textrm{inter}}$ to denote the
inter-layer coupling of these two Dirac-cone states. Therefore, the Dirac-
fermion Hamiltonian of a multilayer can be directly written as,

$H(\textbf{k})=\begin{pmatrix}H_{\textrm{BL}}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}&&\\\
\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}&\\\
&\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}(\textbf{k})&\ddots\\\
&&\ddots&\ddots\end{pmatrix}.$ (4)

Here, $\tilde{\Delta}_{\textrm{inter}}=\Delta_{\textrm{inter}}\tau_{-}$, with
$\tau_{-}=(\tau_{x}-i\tau_{y})/2$, where $\tau_{x,y}$ are Pauli matrices
acting on the top and bottom surface states for each BL.

The above Dirac-fermion approach is unified to describe the low-energy
electronic structures of a wide diversity of multilayers of various layered
(magnetic) TIs. It not only offers a large freedom in tuning parameters
($\Delta_{\textrm{intra}}$, $\Delta_{\textrm{inter}}$, electric field $E$,
magnetic exchange field $m$, etc.), but also can often provide a simple and
clear physical picture. For example, as an important regulatory means, the
external electric field can be directly simulated through tuning on-site
potentials in the Dirac-fermion approach, because the Dirac-cone state is
local in real space. When a magnetic TI is considered, a mass term
$m\sigma_{z}$ induced by breaking the time-reversal symmetry is directly
introduced into the Dirac-cone Hamiltonian $H_{D}(\textbf{k})$. A notable
example is MnBi2Te4, known as a layered A-type AFM TI, which can be regarded
as a multilayer of a magnetic TI. Interestingly, its phase diagram can be
investigated with the Dirac-fermion approach as we show in the following
section. It is worth mentioning that, in our previous work, we have employed a
three-Dirac-fermion model Wang  _et al._ (2022) to discover a topologically-
protected gapless surface state for layered magnetic TIs under the vdW
modulation, which provides a unique intrinsic mechanism for understanding the
puzzle of the observed gapless surface state of MnBi2Te4 Chen  _et al._
(2019); Li  _et al._ (2019a); Shikin  _et al._ (2021); Hao  _et al._ (2019);
Nevola  _et al._ (2020). Further, the MnBi2Te4-family materials, such as,
MnBi4Te7 Ding  _et al._ (2020); Vidal  _et al._ (2019); Wu  _et al._ (2020);
Yan  _et al._ (2020) and MnBi6Te10 Shikin  _et al._ (2022); Shi  _et al._
(2019); Vidal  _et al._ (2021); Yan  _et al._ (2020), can be regarded as
multilayers of magnetic TI MnBi2Te4 and nonmagnetic TI Bi2Te3, and can also be
described by the Dirac-fermion approach. In addition, multilayers of TIs and
topologically trivial insulators can also be well investigated through the
Dirac-fermion approach, for example, the multilayer of TI (Bi,Sb)2Te3 and
heavily Cr-doped (Bi,Sb)2Te3.

##  III Dirac-fermion approach for magnetic multilayers

Figure 2: (a) Schematic of a BL with two local Dirac-cone states (D1 and D2).
(b) Band structure from the Dirac model of a BL in the topological phase
transitions from $C=1$ to $C=0$. Here, $|1(2),\uparrow(\downarrow)\rangle$
denote the Dirac-cone states D1 and D2, where $\uparrow(\downarrow)$ indicate
the spin state. $|\pm,\uparrow(\downarrow)\rangle$ denote the bonding ($+$)
and antibonding ($-$) states of $|1(2),\uparrow(\downarrow)\rangle$. For the
condition of $\Delta_{\textrm{intra}}=0$ and $m\neq 0$, both D1 and D2 open an
energy gap and contribute the Chern number $C_{1,2}=1/2$ for a total
$C=C_{1}+C_{2}=1$. A finite intra-layer coupling $\Delta_{\textrm{intra}}$
couples the D1 and D2 to form bonding and antibonding states
$|\pm,\uparrow(\downarrow)\rangle$. If we gradually increase
$\Delta_{\textrm{intra}}$ from 0, when $m>\Delta_{\textrm{intra}}$, the total
Chern number preserves $C=1$; when $m=\Delta_{\textrm{intra}}$, the energy gap
is closed; when $m<\Delta_{\textrm{intra}}$, the total Chern number become to
be $C=0$. (c) A FM bi-BL multilayer with four Dirac-cone states (D1, D2, D3,
D4). $\Delta_{\textrm{inter}}$ denotes the inter-layer coupling. (d) Schematic
band structure of a FM bi-BL multilayer when $m=0$. V1,2 denote the valence
bands, and C1,2 denote the conduction bands. If the magnetic exchange field is
gradually induced, the V1,2 and C1,2 bands will split and lead to two
topological phase transitions from $C=0$ to $C=1$ at $m=m_{1}$, then to $C=2$
at $m=m_{2}$. (e) The phase diagram on the parameters
($m/\Delta_{\textrm{inter}}$,
$\Delta_{\textrm{intra}}/\Delta_{\textrm{inter}}$) of the FM bi-BL multilayer,
calculated by the Dirac-fermion model. (f) A FM 5-BL multilayer and its phase
diagram by the Dirac-fermion model. (g) Schematic magnetic multilayers consist
of magnetic TIs.

### A. A single BL with FM magnetization

We first consider a BL with a FM magnetization. For simplicity, the out-of-
plane magnetization is taken as an example. The FM exchange field would induce
Zeeman-type spin splitting through adding a mass term $m\sigma_{z}$ to Eq.
(3), which reads,

$H_{\textrm{BL}}^{mm}(\textbf{k})=v(k_{x}\sigma_{y}-k_{y}\sigma_{x})\tau_{z}+\Delta_{\textrm{intra}}\tau_{x}+m\sigma_{z},$
(5)

where the two superscripts $m$ and $m$ denote the FM exchange fields for the
top and bottom Dirac-cone surface states, respectively, which are identical
here by considering a uniform FM order. In addition, the sign of $m$ is
determined by the FM order in the $+z$ or $-z$ direction. For simplicity,
unless otherwise specified, $m$ is set to be positive with the FM order in the
$+z$ direction.

For the topologically nontrivial BL ($M<0$), there are two local Dirac-cone
states (D1 and D2) on the top and bottom surfaces, as schematically shown in
Fig. 2(a). The competition of the intra-layer coupling
$\Delta_{\textrm{intra}}$ and the FM exchange field $m$ can drive the two
Dirac-cone surface states into different topological phases. In order to have
a clear physical picture, we start from the limit of
$\Delta_{\textrm{intra}}=0$, where D1 and D2 are decoupled. Under the FM
exchange field $m$, both D1 and D2 will open an energy gap, and each Dirac
cone contributes a half Chern number ($C_{1,2}=1/2$) to the total Chern number
($C=C_{1}+C_{2}=1$) of the BL, and their bands are doubly degenerate, shown in
Fig. 2(b). Then we introduce the intra-layer coupling
$\Delta_{\textrm{intra}}$ to lead to the bonding and anti-bonding states of D1
and D2, so the degeneracy of the bands of D1 and D2 is lifted, but the Chern
number of the BL remains unchanged as $C=1$ as long as the energy gap keeps
open for a small $\Delta_{\textrm{intra}}$. When further increasing the
$\Delta_{\textrm{intra}}$, the energy gap is expected to close at a critical
value of the $\Delta_{\textrm{intra}}$ and reopen for a larger
$\Delta_{\textrm{intra}}$. The gap closing-and-reopening process is
accompanied by a topological phase transition with the Chern number changing
from $C=1$ to $C=0$, shown in Fig. 2(b). Such a process can be confirmed by
the calculation based on the Dirac-fermion model in Eq. (5). The energy
spectra of the four bands are given by $E_{\bf
k}=\pm\sqrt{k^{2}+\left(m\mp\Delta_{\textrm{intra}}\right)^{2}}$, where
$k=\sqrt{k_{x}^{2}+k_{y}^{2}}$. The Fermi level is located between the second
and the third bands. The band gap at $\Gamma$ is
$2\left|m-\Delta_{\textrm{intra}}\right|$ and the Chern number is obtained as

$C=\frac{1}{2}\left[{{\text{sgn}}\left({m+{\Delta_{\textrm{intra}}}}\right)+{\text{sgn}}\left({m-{\Delta_{\textrm{intra}}}}\right)}\right].$
(6)

This confirms that the Chern number is $C=1$ for $\Delta_{\textrm{intra}}<m$,
and a topological phase transition occurs at $\Delta_{\textrm{intra}}=m$ with
a closed band gap, and the Chern number becomes $C=0$ for
$\Delta_{\textrm{intra}}>m$. Therefore, we can conclude that the competition
between the FM exchange field and the intra-layer coupling induces a
topological transition for a single FM BL.

### B. A bi-BL structure with FM magnetization

Then, we investigate a FM bi-BL structure to see the effect of the inter-layer
coupling $\Delta_{\textrm{inter}}$, as schematically in Fig. 2(c). With
respect to Eqs. (4,5), the bi-BL Hamiltonian can be directly written as

$H_{\textrm{FM}}^{bi}(\textbf{k})=\begin{pmatrix}H_{\textrm{BL}}^{mm}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}\\\
\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}^{mm}(\textbf{k})\end{pmatrix}.$
(7)

There are four local Dirac-cone states, namely, D1 and D2 on the upper BL, D3
and D4 on the lower BL, seen in Fig. 2(c). In the first stage, we turn on the
intra-layer coupling $\Delta_{\textrm{intra}}$ and the inter-layer coupling
$\Delta_{\textrm{inter}}$ in sequence but tentatively suppress the FM exchange
field (e.g. $m=0$) for clarity. The bands of the four Dirac-cone states
(D1-D4) should form four doubly degenerate bands (valence bands V1 and V2,
conduction bands C1 and C2), shown in Fig. 2(d), because of the time-reversal
symmetry and the implicit inversion symmetry in this model. At this stage, the
Chern number has to be $C=0$ due to the preserved the time-reversal symmetry
in the absence of the FM exchange field. Then, we turn on the FM exchange
field, and the double degeneracy of the bands will be gradually lifted. As the
FM exchange field increases, two band inversions sequentially happen between
the conduction and valence bands, corresponding to two topological phase
transitions of ($C=0$ to $C=1$) and ($C=1$ to $C=2$), seen in Fig. 2(d). The
first topological phase transition happens along with the gap closing-and-
reopening process between one V1 band and one C1 band. By diagonalizing the
Hamiltonian of Eq. (7), the corresponding band gap can be calculated as
$\left|2m-\left(\bar{\Delta}-\Delta_{\textrm{inter}}\right)\right|$, where
$\bar{\Delta}=\sqrt{4\Delta_{\textrm{intra}}^{2}+\Delta_{\textrm{inter}}^{2}}$.
When the FM exchange field strength reaches
$m_{1}=\left(\bar{\Delta}-\Delta_{\textrm{inter}}\right)/2$, the bi-BL
structure undergoes a topological phase transition from $C=0$ to $C=1$. With
further increasing the FM exchange field, the second topological phase
transition occurs at the gap closing point between one V2 band and one C2
band, at the critical value
$m_{2}=\left(\bar{\Delta}+\Delta_{\textrm{inter}}\right)/2$, corresponding to
a change of the Chern number from $C=1$ to $C=2$.

The detailed phase diagram, depending on the parameters
($m/\Delta_{\textrm{inter}}$,
$\Delta_{\textrm{intra}}/\Delta_{\textrm{inter}}$), can be obtained from the
Dirac-fermion model in Eq. (7) and is shown in Fig. 2(e). In the limiting case
of $\Delta_{\textrm{intra}}\rightarrow 0$, the two Dirac-cone states (D1,2 or
D3,4) of each BL are completely decoupled, so two Dirac-cone states (D1 and
D4) of the bi-BL structure become isolated, and the other two Dirac-cone
states (D2 and D3) couple together through the inter-layer coupling
$\Delta_{\textrm{inter}}$. The isolated D1 and D4 contribute a Chern number
($C=1$), and the contributed Chern number of the coupled D2 and D3 depends on
the competition of $m$ and $\Delta_{\textrm{inter}}$, for example, the Chern
number $C=1$ for $m>\Delta_{\textrm{inter}}$, and $C=0$ for
$m<\Delta_{\textrm{inter}}$, similar to the consideration in a FM BL.
Therefore, based on the Dirac-fermion model, the highest Chern number of a FM
bi-BL structure can reach but not exceed $|C|=2$. In order to obtain a higher
Chern number, we should stack more FM BLs together. Straightforwardly, at the
limit $\Delta_{\textrm{intra}}\rightarrow 0$, the Chern number of a FM $N$-BL
structure is directly obtained as,

$C_{N}=1+\frac{N-1}{2}[\text{sgn}(m+\Delta_{\textrm{inter}})+\text{sgn}(m-\Delta_{\textrm{inter}})],$
(8)

where the first constant $1$ comes from the two isolated Dirac-cone states on
the top and bottom surfaces of the FM $N$-BL structure, and the second term is
contributed by ($N-1$) pairs of coupled Dirac-cone states between adjacent
BLs. The highest Chern number can be realized as $C=N$ when and only when
$m>\Delta_{\textrm{inter}}$. As an example, five FM BLs are stacked together
to the the highest Chern number $C=5$ at the limit
$\Delta_{\textrm{intra}}\rightarrow 0$, seen in Fig. 2(f). Therefore, we
should stack at least $N$ FM BLs to realize a Chern number $C=N$ ($N\geq 1$),
but we do not need to care for its details of each BL, within the framework of
the Dirac-fermion approach. For example, we know that only the Chern number
$C=1$ is realized in a Cr-doped (Bi,Sb)2Te3 5-QL film in experiments Chang
_et al._ (2013), which can be effectively considered as a FM BL, though it
consists of 5-QL Cr-doped (Bi,Sb)2Te3.

### C. Magnetic TI multilayers

Based on the Dirac-fermion approach, we now consider constructing magnetic TI
multilayers through stacking magnetic/non-magnetic TI layers, for example,
alternate stacking of a FM TI layer and a non-magnetic TI layer ($M<0$), shown
in Fig. 2(g). The FM TI layer can be described by a FM $N$-BL structure, where
$N$ is determined by the highest Chern number the FM TI layer can achieve, and
the non-magnetic TI layer is described by a non-magnetic BL. In this sense, FM
MnBi4Te7 can be effectively considered as a FM TI multilayer of a FM MnBi2Te4
BL and a non-magnetic Bi2Te3 BL. A FM TI multilayer can also consist of a FM
TI layer and a non-magnetic trivial insulator layer, where the trivial
insulator layer acts as a buffer layer to decrease the inter-layer coupling
$\Delta_{\textrm{inter}}$. Besides the FM TI multilayers, AFM TI multilayers
also can be considered. For example, an AFM TI multilayer consists of two FM
TI layers with the opposite magnetization direction, shown in Fig. 2(g). With
the same consideration, each FM TI layer is described by a FM $N$-BL structure
and $N$ is determined by the highest Chern number of the FM TI layer. For
example, AFM MnBi2Te4 can be effectively considered as an AFM TI multilayer,
in which one SL is considered as a FM BL.

##  IV FM $\text{MnBi}_{2}\text{Te}_{4}$ films

Figure 3: (a,b) The phase diagrams of FM 5-BL and 13-BL multilayers. The black
points denote the phase of MnBi2Te4 without the vdW gap modulation. The red
lines denote the phase evolution due to the vdW gap modulation with fixed
$\Delta_{\textrm{intra}}$ and $m$. $\delta d>0$ ($\delta d<0$) indicates the
vdW gap expansion (contraction). (c) The schematic of the vdW gap modulation
($\delta d$). (d) The phase diagram on the parameters (the number of SL $N$,
the vdW gap modulation $\delta$) of AFM MnBi2Te4 films is calculated by the
Dirac-fermion model with $m=40$ meV, $\Delta_{\textrm{intra}}=86.2$ meV, and
$\Delta_{\textrm{inter}}=74.8$ meV from fitting first-principles calculations.
The two red dashed lines in (d) correspond to the phase evolutions (black
solid lines) in (a, b).

In this section, we revisit the realization of the Chern insulator state with
a high Chern number in FM MnBi2Te4 films. MnBi2Te4 is an intrinsic magnetic TI
which provides an ideal platform to study the interplay between magnetism and
topology Gong  _et al._ (2019); Zhang  _et al._ (2019); Li  _et al._ (2019b);
Otrokov  _et al._ (2019a); Sun  _et al._ (2019); Gu  _et al._ (2021); Otrokov
_et al._ (2019b); Wang  _et al._ (2020). MnBi2Te4 has a layered crystal
structure stacked by (Te-Bi-Te-Mn-Te-Bi-Te) SLs with a weak vdW interaction.
Its magnetic ground state is an A-type AFM order with an out-of-plane magnetic
anisotropy along the (111) direction. Although the AFM order breaks the time-
reversal symmetry, the combined symmetry $S=T\tau_{1/2}$ of the time reversal
operation ($T$) and a half-translation operation ($\tau_{1/2}$) is preserved
in MnBi2Te4, which plays the role of the time-reversal symmetry and leads to a
$Z_{2}$ topological classification with $M<0$ in the $k\cdot p$ Hamiltonian
Mong  _et al._ (2010); Zhang  _et al._ (2019). Interestingly, a moderate
external magnetic field ($\sim 10$ T) can flop spins of the AFM phase to a FM
phase MnBi2Te4 Deng  _et al._ (2020); Ge  _et al._ (2020); Ying  _et al._
(2022). The FM MnBi2Te4 was predicted to be an ideal Weyl semimetal in the
bulk Zhang  _et al._ (2019); Li  _et al._ (2019b), which indicates that a
Chern insulator state with a high Chern number should be realized in FM
MnBi2Te4 films.

Next, we investigate the Chern insulator state in FM phase MnBi2Te4 films
through employing the Dirac-fermion approach. First of all, if a MnBi2Te4 SL
is considered as a BL with $M<0$, then one FM MnBi2Te4 $N$-SL film can be seen
as a FM $N$-BL TI multilayer. Secondly, the parameter of
$\Delta_{\textrm{intra}}$ denotes the intra-layer coupling of the two Dirac-
cone states of each MnBi2Te4 SL, the parameter of $\Delta_{\textrm{inter}}$
denotes the inter-layer coupling of the two Dirac-cone states between the
adjacent MnBi2Te4 SLs, and the parameter $m$ denotes the FM exchange field
within each MnBi2Te4 SL.

We then take the FM $5$-SL MnBi2Te4 film as an example. The Dirac-fermion
model can be directly set up for the FM $5$-BL TI multilayer through Eq. (4),
and a phase diagram can be obtained on the parameters
($\Delta_{\textrm{intra}}/\Delta_{\textrm{inter}}$,
$m/\Delta_{\textrm{inter}}$), shown in Fig. 3(a). Though the Chern numbers
from 0 to 5 are expected to be realized through tuning
$\Delta_{\textrm{intra}}/\Delta_{\textrm{inter}}$ and
$m/\Delta_{\textrm{inter}}$, the realistic Chern number of the FM $5$-SL
MnBi2Te4 film is calculated to be $C=1$ from material-dependent parameters
($\Delta_{\textrm{intra}}$, $\Delta_{\textrm{inter}}$ and $m$). Moreover, we
have also calculated the phase diagram for the FM $13$-BL TI multilayer, shown
in Fig. 3(b), where we can see that the highest Chern number $C=13$ can be
obtained at the limiting situation of
($\Delta_{\textrm{intra}}/\Delta_{\textrm{inter}}\rightarrow 0$,
$m/\Delta_{\textrm{inter}}>1$), while the realistic Chern number of the FM
phase $13$-SL MnBi2Te4 film is calculated to be $C=2$. On the whole, the
results of FM phase $5$\- and $13$-SL MnBi2Te4 films are consistent with the
reported experiments Deng  _et al._ (2020); Ge  _et al._ (2020); Ying  _et
al._ (2022).

Further, we consider the effect of the vdW gap modulation ($d_{0}+\delta d$)
of FM MnBi2Te4 films through external pressures, where $d_{0}$ denotes the
pristine vdW gap. The in-plane pressure is expected to induce a vdW gap
expansion ($\delta d>0$), and the out-of-plane pressure is to induce a vdW gap
contraction ($\delta d<0$). We assume the inter-layer coupling
$\Delta_{\textrm{inter}}$ changes exponentially to the deviation ($\delta d$)
from the pristine vdW gap, $\Delta_{\textrm{inter}}(\delta
d)=\Delta_{\textrm{inter}}\text{exp}{(-\lambda\delta d)}$ where $\lambda$
denotes the changing rate. The effect of the vdW gap modulation in the phase
diagram of FM $5$\- and $13$-SL MnBi2Te4 films is shown in Figs. 3(a) and
3(b), respectively. For the FM $5$-SL MnBi2Te4 film, we can see that the vdW
gap expansion ($\delta d>0$) can induce a topological phase transition from
$C=1$ to $C=0$, and that the vdW gap contraction ($\delta d<0$) prefers to
keep the ($C=1$) phase. But it is quantitatively different for the FM $13$-SL
MnBi2Te4 film. The vdW gap expansion ($\delta d>0$) can induce two topological
phase transitions from $C=2$ to $C=1$ and then $C=1$ to $C=0$, and the vdW gap
contraction ($\delta d<0$) induces a topological phase transition from $C=2$
to $C=1$. Based on the Dirac-fermion model, we carry out systematic
calculations on the pressure-induced phase diagrams for FM $N$-SL MnBi2Te4
films with $N$ from 1 to 20, shown in Fig. 3(d). We can notice that a large
vdW gap expansion ($\delta d>0$) always leads to a topologically trivial phase
with $C=0$, which is well understood because it is at the decoupled limit
situation. Interestingly, a large vdW gap contraction ($\delta d<0$) also
prefers to induce a low Chern number $C=1$, not a high Chern number. The high
Chern number only arises for a moderate vdW gap modulation, for example, $C=2$
for $N=13$ and $C=3$ for $N=20$, which is an important consequence from the
Dirac-fermion approach.

##  V nonmagnetic TI/magnetic insulator multilayers

Recently, Zhao et al. Zhao  _et al._ (2022, 2020) fabricated magnetic TI
multilayers consisting of alternate [3QL-(Bi,Sb)1.76Cr0.24Te3
(BCT)]/[4QL-(Bi,Sb)2Te3 (BT)] bilayers. Notably, a high Chern number of $C=5$
has been experimentally realized in these multilayers Zhao  _et al._ (2022).
Moreover, a topological phase transition from $C=2$ to $C=1$ has been observed
in a [3QL-BCT/4QL-BT/3QL-BCT/4QL-BT/3QL-BCT] multilayer through changing Cr
doping of the middle 3QL-BCT layer Zhao  _et al._ (2020). In the following, we
will present a understanding of the two phenomena in these multilayers from
the viewpoint of the Dirac-fermion approach.

### A. A high Chern number $C=N$ case

For a TI multilayer structure constructed by $N$ alternating [3QL-BCT]/[4QL-
BT] layers, shown in Fig. 4(a), each 4QL-BT layer can be chosen as a
nonmagnetic [4QL-BT] BL. Two Dirac-cone states on the top and bottom surfaces
of each [4QL-BT] BL are weakly coupled through a weak intra-layer coupling
$\Delta_{\textrm{intra}}$ because of the thickness of [4QL-BT] BL.
Differently, because the spin-orbit coupling (SOC) of Cr atom is negligible,
the effective SOC in the 3QL-BCT layer is expected to be much reduced due to
the heavily doped Cr through replacing Bi by Cr, so the bulk BCT becomes a
topologically trivial FM insulator. On one side, the presence of the 3QL-BCT
layer significantly reduces the inter-layer coupling $\Delta_{\textrm{inter}}$
between neighboring [4QL-BT] BLs, which enables us to take the limit of
$\Delta_{\textrm{inter}}\rightarrow 0$ to have a clear physical picture. On
the other side, the 3QL-BCT layer provides a FM proximity effect to its
adjacent [4QL-BT] BLs, which brings a magnetic exchange field $m$ to the
Dirac-cone states of the [4QL-BT] BL. As discussed in Sec. III, for each [4QL-
BT] BL, the competition between $\Delta_{\textrm{intra}}$ and $m$ should lead
to different topological phases, namely, $C=1$ when
$m>\Delta_{\textrm{intra}}$, and $C=0$ when $m<\Delta_{\textrm{intra}}$.

In this decoupled limit, the total Hamiltonian from the Dirac-fermion model
can be simply written as a direct sum of $N$ [4QL-BT] BLs

$H_{N}(\textbf{k})=\bigoplus_{i=1}^{N}H_{i,\textrm{BL}}^{mm}(\textbf{k}),$ (9)

where $i$ is the [4QL-BT] BL index and $m$ denotes the magnetic exchange field
on each Dirac-cone state. According to the discussion in Sec. III, the total
Chern number of the multilayer is then given by

$C=\frac{N}{2}\left[\text{sgn}\left(m+\Delta_{\textrm{intra}}\right)+\text{sgn}\left(m-\Delta_{\textrm{intra}}\right)\right].$
(10)

Therefore, if $m>\Delta_{\textrm{intra}}$, a high Chern number of $C=N$ can be
obtained in $N$-[3QL-BCT]/[4QL-BT] multilayers, as illustrated in Fig. 4(b),
which is consistent with the experimental observations Zhao  _et al._ (2020).

Figure 4: (a) The schematic structure of [3QL-BCT]/[4QL-BT] multilayers. (b)
The linearly increased Chern number with the number of [3QL-BCT/4QL-BT]
layers, in the case of weak inter-layer and intra-layer couplings. (c) Chern
number evolution from the 2-BL (red line) and 3-BL (blue line) Dirac model due
to decreasing the Cr doping in the middle 3QL-BCT layer. We take the
parameters $m=0.7$, $\Delta_{\textrm{intra}}=0.2$, and
$\Delta_{\textrm{inter}}=0.1$ for the 2-BL model, and
$\Delta_{\textrm{inter}}=1$, and $m_{0}=2m^{\prime}$ for the 3-BL model.
(d,e,f) Phase diagram for the 2-BL model with $m^{\prime}=m$, and the inter-
layer coupling $\Delta_{\textrm{inter}}=0$ (d) and $\Delta_{\textrm{inter}}=1$
(e). A finite $\Delta_{\textrm{inter}}$ induces a $C=1$ phase (e). Topological
phase transition from $C=2$ to $C=1$ arises when the inter-layer coupling
$\Delta_{\textrm{inter}}$ is continuously increased from 0 to 1 (f).

### B. Topological phase transition from $C=2$ to $C=1$

Then, we focus on the [3QL-BCT/4QL-BT/3QL-BCT/4QL-BT/3QL-BCT] multilayer to
understand the experimentally observed topological phase transitions from
$C=2$ to $C=1$ induced by reducing the Cr doping Zhao  _et al._ (2022) of the
middle 3QL-BCT layer. Generally speaking, there are two effects if the Cr
doping amount is decreased from the heavily doping limit. First of all, the FM
exchange field produced by moments of Cr atoms should be accordingly
decreased. Secondly, the SOC in the 3QL-BCT layer should be enhanced, and the
energy gap of the 3QL-BCT layer is accordingly reduced, which would lead to an
enlarged effective inter-layer coupling $\Delta_{\textrm{inter}}$ between two
neighboring [4QL-BT] BLs. It should also be noted that with decreasing the Cr
doping amount, the BCT could possibly becomes topologically non-trivial with a
band inversion ($M<0$) hosting extra Dirac-cone states. In the following,
depending on whether the middle 3QL-BCT layer is topologically trivial or not,
we will set up a 2-BL or 3-BL model to discuss the phase transition from $C=2$
to $C=1$.

Firstly, when the middle 3QL-BCT layer is topologically trivial, the whole
multilayer system can be effectively described by the Dirac-fermion model
consisting of two [4QL-BT] BLs, with the Hamiltonian given by

$H_{2\textrm{BL}}(\textbf{k})=\begin{pmatrix}H_{\textrm{BL}}^{mm^{\prime}}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}\\\
\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}^{m^{\prime}m}(\textbf{k})\end{pmatrix},$
(11)

where

$H_{\textrm{BL}}^{mm^{\prime}}=\begin{pmatrix}H_{D}^{m}(\textbf{k})&\Delta_{\textrm{intra}}\\\
\Delta_{\textrm{intra}}&H_{D}^{m^{\prime}}(\textbf{k})\end{pmatrix}.$ (12)

Here, $H_{\textrm{BL}}^{mm^{\prime}}$ indicates the Dirac-fermion model for a
[4QL-BT] BL, and $m$ ($m^{\prime}$) denotes the magnetic proximity exchange
filed acting on top (bottom) Dirac-cone states, respectively. For simplicity,
we first consider the case with $\Delta_{\textrm{inter}}=0$, where the two BLs
are decoupled, which enables us to independently consider
$H_{\textrm{BL}}^{mm^{\prime}}$ or $H_{\textrm{BL}}^{m^{\prime}m}$. Since only
the Cr-doping of the middle 3QL-BCT layer is decreased, $m^{\prime}<m$ should
be expected. Under this condition, it can be shown that, for each [4QL-BT] BL,
when $m^{\prime}>\Delta_{\textrm{intra}}^{2}/m$
($m^{\prime}<\Delta_{\textrm{intra}}^{2}/m$) is satisfied, the Chern number is
obtained as $C_{\textrm{BL}}=1$ ($C_{\textrm{BL}}=0$). Therefore, when
decreasing $m^{\prime}$ from $m$ ($m>\Delta_{\textrm{intra}}$ is assumed)
along with the gradually reduced Cr-doping, the multilayer should have a
topological phase transition directly from $C=1+1=2$ to $C=0+0=0$, without an
intermediate $C=1$ phase, shown in Fig. 4(c). The phase diagram on parameters
($\Delta_{\textrm{intra}}$, $m$) is calculated and shown in Fig. 4(d), and we
also see that there are only two regions of $C=2$ and $C=0$. Therefore, though
the $C=2$ Chen insulator state is well captured, the results do not reproduce
the topological phase transition from $C=2$ to $C=1$ observed in experiments
Zhao  _et al._ (2022). Then, we will show that the missing $C=1$ phase emerges
when taking the inter-layer coupling $\Delta_{\textrm{inter}}$ into account.
For simplicity, we fix $m^{\prime}=m$. When we consider a finite inter-layer
coupling (e.g. $\Delta_{\textrm{inter}}=1$), a $C=1$ phase region will emerge
in the phase diagram between the $C=0$ and $C=2$ phases, as shown in Fig.
4(e). This can be easily understood based on the 2-BL Dirac-fermion model,
schematically shown in Fig. 2(d). A finite inter-layer coupling pushes up the
conduction bands C2 and pulls down the valence bands V2, and then causes a
topological phase transition from $C=2$ to $C=1$, shown in Fig. 4(f).
Therefore, the inter-layer coupling $\Delta_{\textrm{inter}}$ is indispensable
for the phase transition from $C=2$ to $C=1$, and the enhanced inter-layer
coupling with decreasing the Cr doping further facilitates the experimental
observation of this transition in [3QL-BCT/4QL-BT/3QL-BCT/4QL-BT/3QL-BCT]
multilayer Zhao  _et al._ (2022).

Secondly, when the middle 3QL-BCT layer becomes topologically nontrivial in
the low Cr-doping regime, we need to build a 3-BL Dirac-fermion model
including two [4QL-BT] BLs and one [3QL-BCT] BL, and the Hamiltonian is
written as,

$H_{3\textrm{BL}}(\textbf{k})=\begin{pmatrix}H_{\textrm{BL}}^{mm^{\prime}}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}&\\\
\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}^{m_{0}m_{0}}(\textbf{k})&\tilde{\Delta}_{\textrm{inter}}\\\
&\tilde{\Delta}_{\textrm{inter}}^{\dagger}&H_{\textrm{BL}}^{m^{\prime}m}(\textbf{k})\end{pmatrix},$
(13)

where $m_{0}$ is the magnetic exchange filed acting on the Dirac states of the
middle [3QL-BCT] BL. Because the inter-layer coupling
$\Delta_{\textrm{inter}}$ is much enhanced with decreased Cr doping, the three
BLs are strongly bound together, compared with the 2-BL Dirac-fermion model.
It is found that within a large parameter region, the whole multilayer system
retains $C=1$, seen in Fig. 4(c).

We should stress that the 3-BL Dirac-fermion model is valid only for the
topologically nontrivial middle 3QL-BCT layer, while the 2-BL Dirac-fermion
model is applicable for the topologically trivial middle 3QL-BCT layer.
Therefore, if the BCT changes from a topologically trivial phase ($M>0$) to a
topologically nontrivial phase ($M<0$) with decreasing the Cr doping, the 3-BL
and 2-BL Dirac-fermion model should be combined to describe the observed
topological phase transition from $C=2$ to $C=1$ for the [3QL-BCT/4QL-BT/3QL-
BCT/4QL-BT/3QL-BCT] multilayer in the experiments Zhao  _et al._ (2022).

##  VI Summary

To summarize, we build a Dirac-fermion approach to effectively describe the
low-energy electronic structure of various multilayers of layered (magnetic)
TIs. Based on the Dirac-fermion model, we demonstrate how to design Chern
insulator states with a high Chen number in multilayers of layered magnetic
TIs. As an example, we employ the Dirac-fermion model to investigate the Chen
insulator state of FM MnBi2Te4 films under the vdW gap modulation. We
calculated the phase diagram and found that a moderate vdW gap modulation is
necessary to realize a high Chen number for FM MnBi2Te4 thick films. We also
employ the Dirac-fermion approach to understand the topological phase
transition in the [3QL-BCT/4QL-BT] multilayers. Weak intra-layer and inter-
layer couplings and a strong magnetic exchange proximate field are the key to
realizing the high Chern number of $C=N$ in $N$-[3QL-BCT]/[4QL-BT]
multilayers. Further, the Dirac-fermion model reproduces the topological phase
transition from $C=2$ to $C=1$ for the [3QL-BCT/4QL-BT/3QL-BCT/4QL-BT/3QL-BCT]
multilayer when the Cr doping is decreased, which reveals that either a
topologically trivial middle 3QL-BCT layer with an enhanced inter-layer
coupling or a topologically nontrivial middle 3QL-BCT layer can induce the
topological phase transition. Therefore, the developed Dirac-fermion approach
provides a powerful tool to understand the topological phase transition in
various of multilayers of layered magnetic TIs.

###### Acknowledgements.

This work is supported by National Key Projects for Research and Development
of China (Grant No.2021YFA1400400 and No.2017YFA0303203), the Fundamental
Research Funds for the Central Universities (Grant No. 020414380185), Natural
Science Foundation of Jiangsu Province (No. BK20200007), the Natural Science
Foundation of China (No. 12074181, No. 12104217, and No. 11834006) and the Fok
Ying-Tong Education Foundation of China (Grant No. 161006). Dinghui Wang is
supported by the program A/B for Outstanding PhD candidate of Nanjing
University.

##  Appendix A Surface states wavefunction

For a thin film of a three-dimensional TI located within $|z|<L/2$ in the
$z$-direction, we use the quantum-well-based methods Liu  _et al._ (2010a),
where the basis of the quantum well states is given by

$\langle{z}\ket{n\alpha}=\sqrt{\frac{2}{L}}\sin\Big{(}\frac{z}{L}n\pi+\frac{n\pi}{2}\Big{)}\ket{\Lambda_{\alpha}}.$
(14)

where $n=1,2,3...$, and $\ket{\Lambda_{\alpha}}$ ($\alpha=1,2,3,4$) denotes
the basis of the Hamiltonian of the TI in Eq. (2). By the substitution
$k_{z}\rightarrow-i\partial_{z}$ in the original four-band Hamiltonian, we can
obtain
$H_{nm}=\bra{n}H_{3D}(k_{x},k_{y},k_{z}\rightarrow-i\partial_{z})\ket{m}$,
with the help of the relations

$\displaystyle\braket{n}{m}=\delta_{nm},$ (15)
$\displaystyle\bra{n}-i\partial_{z}\ket{m}=\frac{2}{L}\frac{(-1+(-1)^{m+n})nm}{n^{2}-m^{2}},$
(16)
$\displaystyle\bra{n}-\partial^{2}_{z}\ket{m}=\left(\frac{n\pi}{L}\right)^{2}\delta_{nm}.$
(17)

By choosing a cut-off number of $N$ of the quantum well states and
diagonalizing the $4N\times 4N$ Hamiltonian $H_{nm}$, the $i$-th eigenstates
of the quantum well is then obtained as

$\ket{\varphi^{i}(z)}=\sqrt{\frac{2}{L}}\sum_{\alpha=1}^{4}\sum_{n=1}^{N}c^{i}_{\alpha,n}\sin\Big{(}\frac{z}{L}n\pi+\frac{n\pi}{2}\Big{)}\ket{\Lambda_{\alpha}},$
(18)

where $c^{i}_{\alpha,n}$ is the superposition coefficient. Then the
probabilitiy density for the wavefunction reads

$|\varphi^{i}(z)|^{2}=\frac{2}{L}\sum_{\alpha=1}^{4}\left|\sum_{n=1}^{N}c^{i}_{\alpha,n}\sin\Big{(}\frac{z}{L}n\pi+\frac{n\pi}{2}\Big{)}\right|^{2}.$
(19)

For the topologically nontrivial case, we focus on the bands located inside
the bulk gap, while for topologically trivial case, the highest conduction
bands are chosen. In the numerical calculations, we choose the same parameters
as those in Ref. Liu  _et al._ (2010a).

Following Eq. (2.32) in Ref. Shen (2012), we construct a exponentially decay
function to decomposed the surface wave function, which reads

$f(z)=A\left(e^{\tilde{z}/\xi_{+}}-e^{\tilde{z}/\xi_{-}}\right),$ (20)

where $A$ is the normalized coeffecient, and
$\xi_{\pm}^{-1}=\lambda_{0}-\lambda$. This function can well capture
properties of density of the surface states. In a narrow quantum well, this
function is used to decompose the combined surface states in Figs. 1(f)-(h).

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Bernevig  _et al._ (2006) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * König  _et al._ (2007) M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Zhang  _et al._ (2009) H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nat. Phys. 5, 438 (2009).
  * Chen  _et al._ (2009) Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Science 325, 178 (2009).
  * Thouless  _et al._ (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Qi  _et al._ (2008) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008).
  * Weng  _et al._ (2015) H. Weng, R. Yu, X. Hu, X. Dai, and Z. Fang, Adv. Phys. 64, 227 (2015).
  * Liu  _et al._ (2016) C.-X. Liu, S.-C. Zhang, and X.-L. Qi, Annu. Rev. Condens. Matter Phys. 7, 301 (2016).
  * Haldane (1988) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Zhang  _et al._ (2013) J. Zhang, C.-Z. Chang, P. Tang, Z. Zhang, X. Feng, K. Li, L. li Wang, X. Chen, C. Liu, W. Duan, K. He, Q.-K. Xue, X. Ma, and Y. Wang, Science 339, 1582 (2013).
  * Chang  _et al._ (2013) C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, _et al._ , Science 340, 167 (2013).
  * Kou  _et al._ (2014) X. Kou, S.-T. Guo, Y. Fan, L. Pan, M. Lang, Y. Jiang, Q. Shao, T. Nie, K. Murata, J. Tang, Y. Wang, L. He, T.-K. Lee, W.-L. Lee, and K. L. Wang, Phys. Rev. Lett. 113, 137201 (2014).
  * Checkelsky  _et al._ (2014) J. G. Checkelsky, R. Yoshimi, A. Tsukazaki, K. S. Takahashi, Y. Kozuka, J. Falson, M. Kawasaki, and Y. Tokura, Nat. Phys. 10, 731 (2014).
  * Xiao  _et al._ (2018) D. Xiao, J. Jiang, J.-H. Shin, W. Wang, F. Wang, Y.-F. Zhao, C. Liu, W. Wu, M. H. W. Chan, N. Samarth, and C.-Z. Chang, Phys. Rev. Lett. 120, 056801 (2018).
  * Grauer  _et al._ (2017) S. Grauer, K. M. Fijalkowski, S. Schreyeck, M. Winnerlein, K. Brunner, R. Thomale, C. Gould, and L. W. Molenkamp, Phys. Rev. Lett. 118, 246801 (2017).
  * Mogi  _et al._ (2017) M. Mogi, M. Kawamura, A. Tsukazaki, R. Yoshimi, K. S. Takahashi, M. Kawasaki, and Y. Tokura, Sci. Adv. 3, eaao1669 (2017).
  * Bestwick  _et al._ (2015) A. J. Bestwick, E. J. Fox, X. Kou, L. Pan, K. L. Wang, and D. Goldhaber-Gordon, Phys. Rev. Lett. 114, 187201 (2015).
  * Liu  _et al._ (2020) C. Liu, Y. Wang, H. Li, Y. Wu, Y. Li, J. Li, K. He, Y. Xu, J. Zhang, and Y. Wang, Nat. Mater. 19, 522 (2020).
  * Deng  _et al._ (2020) Y. Deng, Y. Yu, M. Z. Shi, Z. Guo, Z. Xu, J. Wang, X. H. Chen, and Y. Zhang, Science 367, 895 (2020).
  * Serlin  _et al._ (2020) M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Sharpe  _et al._ (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Wang  _et al._ (2013) J. Wang, B. Lian, H. Zhang, Y. Xu, and S.-C. Zhang, Phys. Rev. Lett. 111, 136801 (2013).
  * Fang  _et al._ (2014) C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. Lett. 112, 046801 (2014).
  * Jiang  _et al._ (2012) H. Jiang, Z. Qiao, H. Liu, and Q. Niu, Phys. Rev. B 85, 045445 (2012).
  * Du  _et al._ (2020) S. Du, P. Tang, J. Li, Z. Lin, Y. Xu, W. Duan, and A. Rubio, Phys. Rev. Research 2, 022025 (2020).
  * Zhu  _et al._ (2022) W. Zhu, C. Song, H. Bai, L. Liao, and F. Pan, Phys. Rev. B 105, 155122 (2022).
  * Zhao  _et al._ (2020) Y.-F. Zhao, R. Zhang, R. Mei, L.-J. Zhou, H. Yi, Y.-Q. Zhang, J. Yu, R. Xiao, K. Wang, N. Samarth, _et al._ , Nature 588, 419 (2020).
  * Fu  _et al._ (2007) L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Moore and Balents (2007) J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (2007).
  * Roy (2009) R. Roy, Phys. Rev. B 79, 195321 (2009).
  * Xia  _et al._ (2009) Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nat. Nat. 5, 398 (2009).
  * Liu  _et al._ (2010a) C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and S.-C. Zhang, Phys. Rev. B 82, 045122 (2010a).
  * Fu (2009) L. Fu, Phys. Rev. Lett. 103, 266801 (2009).
  * Zhang  _et al._ (2010) W. Zhang, R. Yu, H.-J. Zhang, X. Dai, and Z. Fang, New J. Phys. 12, 065013 (2010).
  * Liu  _et al._ (2010b) C.-X. Liu, H. Zhang, B. Yan, X.-L. Qi, T. Frauenheim, X. Dai, Z. Fang, and S.-C. Zhang, Phys. Rev. B 81, 041307 (2010b).
  * Zhou  _et al._ (2008) B. Zhou, H.-Z. Lu, R.-L. Chu, S.-Q. Shen, and Q. Niu, Phys. Rev. Lett. 101, 246807 (2008).
  * Shan  _et al._ (2010) W.-Y. Shan, H.-Z. Lu, and S.-Q. Shen, New J. Phys. 12, 043048 (2010).
  * Lu  _et al._ (2010) H.-Z. Lu, W.-Y. Shan, W. Yao, Q. Niu, and S.-Q. Shen, Phys. Rev. B 81, 115407 (2010).
  * Yu  _et al._ (2010) R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, Science 329, 61 (2010).
  * Wang  _et al._ (2022) D. Wang, H. Wang, D. Xing, and H. Zhang, arXiv:2205.08204 (2022).
  * Chen  _et al._ (2019) Y. J. Chen, L. X. Xu, J. H. Li, Y. W. Li, H. Y. Wang, C. F. Zhang, H. Li, Y. Wu, A. J. Liang, C. Chen, S. W. Jung, C. Cacho, Y. H. Mao, S. Liu, M. X. Wang, Y. F. Guo, Y. Xu, Z. K. Liu, L. X. Yang, and Y. L. Chen, Phys. Rev. X 9, 041040 (2019).
  * Li  _et al._ (2019a) H. Li, S.-Y. Gao, S.-F. Duan, Y.-F. Xu, K.-J. Zhu, S.-J. Tian, J.-C. Gao, W.-H. Fan, Z.-C. Rao, J.-R. Huang, J.-J. Li, D.-Y. Yan, Z.-T. Liu, W.-L. Liu, Y.-B. Huang, Y.-L. Li, Y. Liu, G.-B. Zhang, P. Zhang, T. Kondo, S. Shin, H.-C. Lei, Y.-G. Shi, W.-T. Zhang, H.-M. Weng, T. Qian, and H. Ding, Phys. Rev. X 9, 041039 (2019a).
  * Shikin  _et al._ (2021) A. M. Shikin, D. A. Estyunin, N. L. Zaitsev, D. Glazkova, I. I. Klimovskikh, S. O. Filnov, A. G. Rybkin, E. F. Schwier, S. Kumar, A. Kimura, N. Mamedov, Z. Aliev, M. B. Babanly, K. Kokh, O. E. Tereshchenko, M. M. Otrokov, E. V. Chulkov, K. A. Zvezdin, and A. K. Zvezdin, Phys. Rev. B 104, 115168 (2021).
  * Hao  _et al._ (2019) Y.-J. Hao, P. Liu, Y. Feng, X.-M. Ma, E. F. Schwier, M. Arita, S. Kumar, C. Hu, R. Lu, M. Zeng, Y. Wang, Z. Hao, H.-Y. Sun, K. Zhang, J. Mei, N. Ni, L. Wu, K. Shimada, C. Chen, Q. Liu, and C. Liu, Phys. Rev. X 9, 041038 (2019).
  * Nevola  _et al._ (2020) D. Nevola, H. X. Li, J.-Q. Yan, R. G. Moore, H.-N. Lee, H. Miao, and P. D. Johnson, Phys. Rev. Lett. 125, 117205 (2020).
  * Ding  _et al._ (2020) L. Ding, C. Hu, F. Ye, E. Feng, N. Ni, and H. Cao, Phys. Rev. B 101, 020412 (2020).
  * Vidal  _et al._ (2019) R. C. Vidal, A. Zeugner, J. I. Facio, R. Ray, M. H. Haghighi, A. U. B. Wolter, L. T. Corredor Bohorquez, F. Caglieris, S. Moser, T. Figgemeier, T. R. F. Peixoto, H. B. Vasili, M. Valvidares, S. Jung, C. Cacho, A. Alfonsov, K. Mehlawat, V. Kataev, C. Hess, M. Richter, B. Büchner, J. van den Brink, M. Ruck, F. Reinert, H. Bentmann, and A. Isaeva, Phys. Rev. X 9, 041065 (2019).
  * Wu  _et al._ (2020) X. Wu, J. Li, X.-M. Ma, Y. Zhang, Y. Liu, C.-S. Zhou, J. Shao, Q. Wang, Y.-J. Hao, Y. Feng, E. F. Schwier, S. Kumar, H. Sun, P. Liu, K. Shimada, K. Miyamoto, T. Okuda, K. Wang, M. Xie, C. Chen, Q. Liu, C. Liu, and Y. Zhao, Phys. Rev. X 10, 031013 (2020).
  * Yan  _et al._ (2020) J.-Q. Yan, Y. H. Liu, D. S. Parker, Y. Wu, A. A. Aczel, M. Matsuda, M. A. McGuire, and B. C. Sales, Phys. Rev. Mater. 4, 054202 (2020).
  * Shikin  _et al._ (2022) A. M. Shikin, N. L. Zaitsev, A. V. Tarasov, T. P. Makarova, D. A. Glazkova, D. A. Estyunin, and I. I. Klimovskikh, JETP Lett. 116, 556 (2022).
  * Shi  _et al._ (2019) M. Z. Shi, B. Lei, C. S. Zhu, D. H. Ma, J. H. Cui, Z. L. Sun, J. J. Ying, and X. H. Chen, Phys. Rev. B 100, 155144 (2019).
  * Vidal  _et al._ (2021) R. C. Vidal, H. Bentmann, J. I. Facio, T. Heider, P. Kagerer, C. I. Fornari, T. R. F. Peixoto, T. Figgemeier, S. Jung, C. Cacho, B. Büchner, J. van den Brink, C. M. Schneider, L. Plucinski, E. F. Schwier, K. Shimada, M. Richter, A. Isaeva, and F. Reinert, Phys. Rev. Lett. 126, 176403 (2021).
  * Gong  _et al._ (2019) Y. Gong, J. Guo, J. Li, K. Zhu, M. Liao, X. Liu, Q. Zhang, L. Gu, L. Tang, X. Feng, D. Zhang, W. Li, C. Song, L. Wang, P. Yu, X. Chen, Y. Wang, H. Yao, W. Duan, Y. Xu, S.-C. Zhang, X. Ma, Q.-K. Xue, and K. He, Chin. Phys. Lett. 36, 076801 (2019).
  * Zhang  _et al._ (2019) D. Zhang, M. Shi, T. Zhu, D. Xing, H. Zhang, and J. Wang, Phys. Rev. Lett. 122, 206401 (2019).
  * Li  _et al._ (2019b) J. Li, Y. Li, S. Du, Z. Wang, B.-L. Gu, S.-C. Zhang, K. He, W. Duan, and Y. Xu, Sci. Adv. 5, eaaw5685 (2019b).
  * Otrokov  _et al._ (2019a) M. M. Otrokov, I. P. Rusinov, M. Blanco-Rey, M. Hoffmann, A. Y. Vyazovskaya, S. V. Eremeev, A. Ernst, P. M. Echenique, A. Arnau, and E. V. Chulkov, Phys. Rev. Lett. 122, 107202 (2019a).
  * Sun  _et al._ (2019) H. Sun, B. Xia, Z. Chen, Y. Zhang, P. Liu, Q. Yao, H. Tang, Y. Zhao, H. Xu, and Q. Liu, Phys. Rev. Lett. 123, 096401 (2019).
  * Gu  _et al._ (2021) M. Gu, J. Li, H. Sun, Y. Zhao, C. Liu, J. Liu, H. Lu, and Q. Liu, Nat. Commun. 12, 3524 (2021).
  * Otrokov  _et al._ (2019b) M. M. Otrokov, I. I. Klimovskikh, H. Bentmann, D. Estyunin, A. Zeugner, Z. S. Aliev, S. Gaß, A. U. B. Wolter, A. V. Koroleva, A. M. Shikin, M. Blanco-Rey, M. Hoffmann, I. P. Rusinov, A. Y. Vyazovskaya, S. V. Eremeev, Y. M. Koroteev, V. M. Kuznetsov, F. Freyse, J. Sánchez-Barriga, I. R. Amiraslanov, M. B. Babanly, N. T. Mamedov, N. A. Abdullayev, V. N. Zverev, A. Alfonsov, V. Kataev, B. Büchner, E. F. Schwier, S. Kumar, A. Kimura, L. Petaccia, G. Di Santo, R. C. Vidal, S. Schatz, K. Kißner, M. Ünzelmann, C. H. Min, S. Moser, T. R. F. Peixoto, F. Reinert, A. Ernst, P. M. Echenique, A. Isaeva, and E. V. Chulkov, Nature 576, 416 (2019b).
  * Wang  _et al._ (2020) H. Wang, D. Wang, Z. Yang, M. Shi, J. Ruan, D. Xing, J. Wang, and H. Zhang, Phys. Rev. B 101, 081109 (2020).
  * Mong  _et al._ (2010) R. S. K. Mong, A. M. Essin, and J. E. Moore, Phys. Rev. B 81, 245209 (2010).
  * Ge  _et al._ (2020) J. Ge, Y. Liu, J. Li, H. Li, T. Luo, Y. Wu, Y. Xu, and J. Wang, Natl. Sci. Rev. 7, 1280 (2020).
  * Ying  _et al._ (2022) Z. Ying, S. Zhang, B. Chen, B. Jia, F. Fei, M. Zhang, H. Zhang, X. Wang, and F. Song, Phys. Rev. B 105, 085412 (2022).
  * Zhao  _et al._ (2022) Y.-F. Zhao, R. Zhang, L.-J. Zhou, R. Mei, Z.-J. Yan, M. H. W. Chan, C.-X. Liu, and C.-Z. Chang, Phys. Rev. Lett. 128, 216801 (2022).
  * Shen (2012) S. Shen,  _Topological insulators_ (Springer-Verlag Berlin Heidelberg, 2012).