Field-controlled quantum anomalous Hall effect in electron-doped CrSiTe$_{ 3 }$ monolayer: a first-principles prediction

We report Chern insulating phases emerging from a single layer of layered chalcogenide CrSiTe$_{3}$, a transition metal trichacogenides (TMTC) material, in the presence of charge doping. Due to strong hybridization with Te $p$ orbitals, the spin-orbit coupling effect opens a finite band gap, leading to a nontrivial topology of the Cr $e_{\mathrm{g}}$ conduction band manifold with higher Chern numbers. Our calculations show that quantum anomalous Hall effects can be realized by adding one electron in a formula unit cell of Cr$_{2}$Si$_{2}$Te$_{6}$, equivalent to electron doping by 2.36$\times$10$^{14}$ cm$^{-2}$ carrier density. Furthermore, the doping-induced anomalous Hall conductivity can be controlled by an external magnetic field via spin-orientation-dependent tuning of the spin-orbit coupling. In addition, we find distinct quantum anomalous Hall phases employing tight-binding model analysis, suggesting that CrSiTe$_{3}$ can be a fascinating new platform to realize Chern insulating systems with higher Chern numbers.


I. INTRODUCTION
Since the discovery of graphene, two-dimensional (2D) materials have been one of the most intensively studied systems owing to their intriguing physical properties and the chance of applications to atomistic low-power devices [1][2][3] . Transition metal chalcogenides (TMC) with layered structures is a representative example of such 2D materials 4,5 . One intensively pursued subject recently in TMC is magnetism originating from the presence of localized magnetic moments in d orbitals at TM ions, where hybridization between TM and chalcogen anions give rise to interesting magnetic behaviours in XPS 3 (X = Mn, Fe, Ni) and CrBTe 3 (B = Si, Ge) [6][7][8][9][10] . Systems like CrBTe 3 or CrI 3 are reported to show Isingtype magnetism and exhibit magnetic orders down to singlelayer limit 11 at finite temperatures 12,13 , defying the Mermin-Wagner theorem that prohibits long-range order in systems with d ≤ 2 dimensions via thermal fluctuations 14 . Successful exfoliations of atomistically-thin CrI 3 and CrGeTe 3 layers presents promising magnetic 2D system with potential device applications 12,13 .
The presence of robust magnetism in two-dimensional system like CrBTe 3 becomes even more interesting in the study of topological properties [15][16][17][18][19][20][21][22][23][24] , where TMC materials has been recently considered as one of good platform of researching topological phases of matters 25,26 and magnetism can be employed as a control knob to tune topological properties in such systems [27][28][29] . Specifically, quantum anomalous Hall effect (QAHE) can be observed in materials with magnetic ordering and sizable band gap, showing quantized integer Hall conductivity in the absence of external magnetic fields. Chern insulator is one of topological phases of matter showing QAHE characterized by Z topological invariant, i.e. Chern numbers, unlike time-reversal-symmetric conventional topological insulators classified by Z 2 -invariants such as Bi 2 Se 3 , Bi 2 Te 3 , Sb 2 Te 3 30 or HgTe 31 . Recent studies suggest Chern insulating phases with high Chern numbers and wide bandgap energy 25,26 can be applied for low-power dissipationless device applications, but exploring suitable candidates is a non-trivial task. CrSiTe 3 , originally reported to be a 2D ferromagnetic semiconductors with trivial band topology 8,[32][33][34] , can be an ideal candidate in the presence of charge doping because of the robust magnetism and the insulating behavior down to the single-layer limit, in addition, due to the presence of strong SOC from Te sites.
In this work, by performing first-principles densityfunctional theory (DFT) calculations, we find Chern insulating phases in CrSiTe 3 monolayer in the presence of electron doping. We find crossing points within the Cr e g conduction band manifolds, where the crossings are removed as spin-orbit-coupling (SOC) is introduced and leading to topologically nontrivial bands with Chern number up to 8. We also find quantized anomalous Hall conductivity (AHC) at one-and two-electron-doping per unit cell, which can be further varied as spin orientation angle changes. Additionally, it is found that the magnetic anisotropy energy (MAE) is suppressed as electron doping is introduced via electrostatic gating, consistently with recent findings [35][36][37] , which makes tuning of AHC via external magnetic field feasible. Finally, we construct a tight-binding (TB) hamiltonian for a deeper understanding of our DFT results and to further pursue the possibility of realizing distinct Chern insulating phases via external stimuli such as strain. Our result reveals that single-layers of magnetic CrSiTe 3 and its structural siblings are promising platforms to realize Chern insulator materials with high Chern numbers.

II. RESULT
A. Electronic and magnetic properties of stoichimetric CrSiTe 3 CrSiTe 3 is one of MAX 3 -type TMTC compounds. It has a R3-type stacking of neighboring CrSiTe 3 monolayers, where neighboring layers are coupled by vdW interaction. Crystal structure of a CrSiTe 3 monolayer is depicted in Fig. 1  Density-of-States (states/eV-f.u. or atom) total Cr e g Cr t 2g Te p Si p

FIG. 2.
Band structure and pDOS plot. (a) Band structure with majority spin (red line) and minority spin channel (blue line), respectively. (b) Total density of states (black line), projected density of states of Cr e g orbital (red line) and t 2g orbital (blue line), Si p orbital (purple line) and Te p orbital (sky blue line), respectively. In the pDOS plot, plus and minus signs correspond to the majority and minority spin components, respectively. Fermi level is set to be zero for both panels.
Si-dimers located at the centers of Cr-hexagons. In undoped CrSiTe 3 , Cr cations are in the d 3 (Cr 3+ ) high-spin configuration (S = 3/2) with fully occupied t 2g orbitals. In addition, strong d pσ hybridization between Cr d-and Te p-orbitals gives rise to additional magnetic moment contributions, so that the magnitude of Cr spin moment is 3.87 µ B 38 . The ferromagnetic ground state can be attributed to the Goodenough-Kanamori-Anderson (GKA) superexchange mechanism 39,40 . Our results are in agreement with previous reports 8,32,41 where choice of Hubbard U parameters is crucial for determining magnetic ground state of CrSiTe 3 38 . When electron doping is introduced, orbital occupation of CrSiTe 3 changes so that magnitude of the magnetic moment of Cr d orbitals, as well as lattice constant, is also affected. With two-electron doping per formula unit cell, Cr e g orbital of spin majority channel becomes half-filled so that the high-spin Cr 2+ (S = 2) configuration is obtained (additional moment of 0.51 µ B further induced due to d pσ-hybridization). With this electron doping, the in-plane lattice parameter a expands to be 7.26 Å, 5.7% larger than that of undoped case 6,32,38 . We check that the ferromagnetic ground state with insulating properties and also topological phases remain unchanged within the Hubbard U parameters and lattice constant range of 0.5 eV ≤ U ≤ 1.5 eV and 6.87 Å ≤ a ≤ 7.26 Å, respectively. Therefore, hereafter we choose the lattice constant and Hubbard U parameters to be 7.0 Å and 1.5 eV, respectively, as a representative case.

B. Band crossings in CrSiTe 3 e g conduction bands
To study magnetic and topological properties of electrondoped monolayer CrSiTe 3 in the FM state, we first focus on the electronic structure of undoped CrSiTe 3 . Fig. 2 describes band structures and projected density-of-states (pDOS) of ferromagnetic single layer CrSiTe 3 . Majority-and minority-spin components show exchange splitting where the minority-spin parts have large band gap energy compared to the majorityspin parts. Right above the Fermi level, four majority-spin bands consisting of Cr e g -and Te p-orbitals exist and are separated from other bands (except Cr t 2g bands in the minorityspin channel). Interestingly, band crossings are found within the e g bands at Γ, K points and on the Γ-M, Γ-K highsymmetry lines. This observation suggests the possibility for the monolayer CrSiTe 3 to host nontrivial band topology in the presence of electron doping and gap opening via SOC, as discussed in the following section.
C. Berry curvature plot, Chern numbers and chiral edge states Fig. 3 shows magnified Cr e g band structures with Chern numbers assigned to each band, plots of Berry curvature in the momentum space, and edge spectra in the presence of integer electron doping, SOC, and out-of-plane ferromagnetic spin orientation. We consider the cases of one and two electron doping per formula unit cell, where both systems are insulating as shown in Fig. 3(a) and (d). For each band, we find unusually high Chern numbers of up to 8, with about 10 meV order of SOC-induced band gap mainly originating from Te p orbital contributions. doping per unit cell ( Fig. 3(b)), peaks of positive Berry curvature for the lowest Cr e g band appear around Γ point and on the Γ -M lines, while negative peaks appear close to K and K points, yielding the Chern number of 2 for the lowest e g band. In the case of half-filled e g (Fig. 3(e)), two peaks of negative Berry curvature are located at K and K points. In addition, six negative peaks are located between six Γ -K lines, resulting in a total Chern number of −4. Fig. 3(c) and (f) show edge spectra from the one and two-electron doping (per formula unit cell) results, respectively. In agreement with the Chern number calculation results, one can find two and four chiral edge states near the Fermi level in Fig. 3(c) and (f), respectively. The emergence of multiple band crossings in Cr e g bands, as shown in Fig. 2(a), gives rise to nontrivial topology with high Chern numbers. Explanation of microscopic origin of multiple Dirac cones in previous theoretical work 25 can also be applied in our case, where hopping integrals calculation results of CrSiTe 3 are listed in Supplementary information section 1. In addition, descriptions for band gap closing at the midpoint of the Γ -K line are discussed. See more details in Supplementary information section 2.

D. AHC under electron doping
We further examine the AHC of Cr e g bands manifold, the total density-of-states of CrSiTe 3 and MAE illustrated in the total density-of-states vanishes at one-and two-electron doping as depicted in Fig. 4(b). Therefore, Cr e g bands are separated from each other together with nonzero Chern numbers so that Chern insulating phases will be realized under one-or two-electron doping in the formula unit cell.

E. Magnetic anisotropy and its doping dependence
So far, our DFT calculations with SOC effect consider spin aligned to the out-of-plane direction. However, this out-ofplane spin configuration can be suppressed or become unstable as the electron doping is introduced because magnetic anisotropy may depend on e g occupation. Fig. 4(c) shows MAE as a function of electron doping concentrations, which oscillate between easy-plane and easy-axis anisotropies. For example, easy-axis anisotropy occurs at one electron doping, while two-electron doping per formula unit cell favors easyplane anisotropy. It is worth mentioning that MAE may vanish as doping is introduced, which may enable tuning the direction of FM moments and the resulting electronic structure via external magnetic fields. Since spins are magnetic dipoles, magnetic long-range dipole-dipole interactions between local moments may change magnetic anisotropy. Because these inter-dipole interactions are not captured within DFT, we employ Ewald's lattice summation technique to compute dipolar energy and estimate its effect on magnetic anisotropy [42][43][44][45] . Table I lists dipolar interaction and anisotropy energies from DFT calculations as a function of electron doping. It is shown that dipolar interactions favor the easy-plane configuration over the easyaxis one in FM states (as shown in the negative values of D-MAE in Table I). Combining anisotropy energies from DFT (MAE) and dipolar interactions (D-MAE), it can be seen that easy-plane spin configurations are more favored except for undoped and one electron doping conditions. Interestingly, total magnetic anisotropy energy (MAE + D-MAE) at ∆n = 1 is reduced to 0.289 meV/f.u., equivalent to 1.18 Tesla of magnetic field strength. Hence controlling spin alignment and the resulting topological properties of one-electron-doped (per formula unit cell) ferromagnetic monolayer CrSiTe 3 via applying an external magnetic field become achievable, which will be discussed further in the following section (Sec. II F).
To investigate more details about the suppression of magnetic anisotropy energy under electron doping as depicted in Fig. 4(c), we calculate total energy changes induced by varying spin directions as illustrated in Fig. 5(a). Without the change in the electronic structure (i.e. when the local spin moment picture is robust), total energy as a function of spin directions should have a form of (M·z) 2 ∼ cos 2 θ (see Fig. 5(b) for the angle definition). Such trend is maintained until θ is increased up to 30 • . Beyond that angle, the anisotropy energy deviates from the simple local moment picture. The result implies an onset of additional terms that comes into play to favor the easy-plane configuration around θ = 30 • , which turns out to originate from the evolution of the electronic structure for θ (see below the discussion on Fig. 6 for more detailed discussion). As a result, shape of total energy shows local minimum near θ ∼ 50 • and magnitude of magnetic anisotropy energy is suppressed about 1 meV at θ = 90 • in Fig. 5(a).

F. Switching AHC via external magnetic fields
Because the magnetic anisotropy energy of Cr e g bands in the presence of one electron doping conditions is small enough to tune the FM via external fields, we investigate the behavior of AHC as the spin orientation is changed from the out-of-plane to in-plane direction. Fig. 5(b) shows the magnitude of AHC, averaged over the azimuthal angle φ of the net magnetization, as a function of polar angle θ of the FM spin orientation to the layer-normal direction. From θ = 0 • to θ = 45 • , AHC remains quantized at 2e 2 /h where there are deviations with respect to azimuthal spin angle φ. Surprisingly, as θ increases, AHC begins to reduce and goes to almost zero at θ = 90 • . We also confirm the anti-symmetric behavior of AHC as a function of spin polar angle θ. One can imagine that spin directions for θ = 180 • is opposite to θ = 0 • so that rotation of corresponding chiral edge modes are reversed, i.e. clockwise to counterclockwise. Because the direction of the FM moment can be switched between the out-of-plane and in-plane direction via external magnetic fields at one-electron doping per formula unit cell, the quantum anomalous Hall phase at this doping can also be switched on-off via external magnetic fields of about 1.18 Tesla estimated in Section II E.
To understand the changes of band features and AHC at one-electron doping as the net magnetization is tilted, we plot Berry curvatures in the BZ and along with the e g band disper-sion with tilting the spin orientation direction as summarized in Fig. 6. Here we set spin azimuthal angle φ = 0 to plot Berry curvatures and band structures. As θ increases, electron and hole pockets start to develop, close to θ = 45 • , in the middle of the Γ-M 1,2,3 lines and around M 1,3 points (see Fig. 6(e) and (f)). Especially, the presence of hole pockets and their expansion contribute to the reduction of the AHC as θ increases beyond 30∼45 • (compare with Fig. 5(b)), while the AHC contributions from electron pockets around M 1,3 points are almost vanishing.
As θ is further increased beyond 60 • , sign of Berry curvature distribution around M 1 and M 3 points are flipped (compare Fig. 6(g) with (i) and (k)). This behavior is attributed to the band touching below the Fermi level and the resulting sign reversal of Berry curvature of occupied bands. Comparing panel Fig. 6(h), (j), and (l), it can be noticed that band crossings occur on Γ-M 1 and Γ-M 3 lines, slightly below the Fermi level, around θ = 75 • . The crossings give rise to sign flipping of the Berry curvature of involved bands in the vicinity of the crossing point, which leads to cancellation of net Berry curvature and vanishing AHC at θ = 90 • .
In addition, note that the band gap between highestoccupied and lowest-unoccupied bands at Γ is suppressed as θ is increased, so that at θ = 90 • the quadratic band touching is restored. This occurs because i) the bands close to the Γ point consist mostly of Te p x,y -orbitals, and ii) SOC in the presence of FM behaves as an orbital Zeeman fields λ SO m ·L (λ SO and m being SOC strength and net magnetization, respectively). When m = m z z the λ SO m zLz splits the p x,y doublet intoL z eigenstates (p x ± ip y ), but when m is in-plane the splitting cannot occur due to the absence of p z character. Hence tilting the spin direction leads to the quenching of SOC in the vicinity of Γ point.

G. Tight binding model
To explore possibilities of further tuning the topological nature of this compound and obtain a deeper understanding of high Chern numbers, we construct a model Hamiltonian that correctly captures characteristics of the conduction bands. As seen in the Fig. 2, the four conduction bands are composed in Cr-e g and Te-p orbitals. We can integrate out the Te-p orbitals and construct the e g effective Hamiltonian as follow: where the first and second term represent the spin-independent and and spin-dependent hopping, respectively. The detailed derivation can be found in the Supplementary information section 3. Fig. 7 shows the TB band structure without SOC effect. The hopping parameters are obtained from the Wannier function analysis. Since the TB model is constructed in the perfect octahedron system without trigonal distortion for simplicity, there exists a slight difference in the energy dispersion compared to the DFT conduction band in Fig. 4, but the overall features agree well. When we include the t soc parameters in the model, the bandgap opens, and each band shows [-2, 6, 8, -4] of Chern number, which also verifies the validity of the model Hamiltonian.
In this model, the SOC-dependent Hamiltonian has quiet an interesting form. For instance, the NN SOC Hamiltonian (H soc 1 ) simply looks as whereM denotes the unit vector for magnetization direction.
In this equation, each vectorx,ŷ, andẑ represents the local coordinates denoted in Fig. 1. This magnetization-dependent hopping gives a different effect on each k point. Since the gap closing is observed at some high symmetry k points, let us examine how these parameter changes on each k point. At Γ point, the parameter is proportional toM · c, where c denotes the [001] direction. Thus, the SOC effect becomes maximum when the spins point out-of-plane direction, while it vanishes when the spin lies on the plane. However, we have slightly different behavior on K = (2/3, 1/3) point. In this case, the matrix element vanishes when the magnetization points out-of-plane direction (t soc 1 (K) ∼ĉ · (ẑ +xe −i4π/3 + ye −i2π/3 ) = 0. For the general k points, this term does not vanishes. However, due to the dependence on the local axis, these terms cancel each other and become relatively small. Thus, we can observe that the band structure from the full SOC Hamiltonian when the magnetization lies in-plane looks very similar to the FM band structure without SOC effect.
To examine the relation between the parameter strength and topological character, we also examined the change of Chern number of each band by varying some parameters. Fig. 8 shows the evolution of the Chern number while varying two parameters t 3 (3 rd NN hopping) and t soc onsite (onsite SOC term). Other parameters such as NN and 2 nd NN terms are extracted from the Wannier function analysis. In this figure, the four integer pair (i 1 , i 2 , i 3 , i 4 ) represent the Chern number of each band calculated with model Hamiltonian. We can observe that the Chern number set of (2,−6,8,−4), the actual Chern number of the CrSiTe 3 system, is observed on the lower-left parameter region. With two-electron doping, the total Chern number equals −4, which does not change even if we change t 3 . It change from 4 to 2 with t soc onsite being 0 ∼ −4 meV. Interestingly, this figure shows that structure distortion can achieve in each phase. The fully optimized structure of CrSiTe 3 suffers trigonal distortion. The angle between Cr-Te-Cr is about 86 • . In this case, t 3 and t soc onsite becomes about 180 and -1 meV as marked with red triangle in the figure. If we perform the DFT calculation in the cubic structure where the local octahedron has no distortion(Cr-Te-Cr angle being 90 • ), the Chern number changes from (2,−6,8,−4) to (2,0,2,−4). The hopping parameters corresponding to this case is about 240 and −3 meV for t 3 and t soc onsite , respectively. To determine the effect of structure distortion on the change of the Chern number, we performed the DFT calculations on the various structures where Cr-Te-Cr angles were set to be between 86 (fully optimized structure) and 90 • (no distortion). On each step, we extracted t 3 and t soc onsite parameters. Other parameters also vary, but their changes were minimal compared to these two parameters, so we fixed them. Each parameter increased linearly from the fully optimized structure marked with a dotted line. During this transition, the total Chern number changes from −4 to 2. Even though it is not an optimized structure, we can further distort the local octahedron to achieve the Cr-Te-Cr angle to be 93 • and monitor the change of hopping parameters. In this case, the Chern number becomes (−4,6,2,−4), and the corresponding 3 rd NN hopping parameter becomes about 300 meV. During this process, the Chern number changes when Cr-Te-Cr angles become about 91.8 • (t 3 being about 270 meV). These observation implies that the local distortion near metal site (Cr) on CrSiTe 3 induces large modification of t 3 and t soc onsite parameters and phase transition. The significant modification of t 3 comes from the modification of bond length between Te p orbitals. Table II summarizes the structure information on each calculation. We can see that the Te-Te bond length decreases from 4.23 (fully optimized) to 3.90 Å (trigonally compressed distortion). As discussed, the 3 rd NN hopping contains Te → Te hopping channel. This considerable reduction of bond length between Te sites increases the t 3 parameter. We can also observe that the Cr-Te distance decreases. This parameter is related to t soc onsite , and we also observe that this value slightly increases. Other parameters such as t 2 or t 5 remains almost the same. The local distortion of the octahedron can be achieved by applying external strain, and it induces a change in the hopping parameters. Thus, we can effectively map the effect of distortion to the change of hopping parameters, although this Hamiltonian assumes a cubic structure. Since the Chern number depends on the hopping parameters and the optimized structure lies near the phase boundary, we expect it to be possible to tune the Chern number phase by tuning external strain or change of ligand hopping strength.

III. DISCUSSION
We study the electronic and topological properties of the electron-doped single-layer structure of CrSiTe 3 by performing first-principle calculations based on density functional theory (DFT). The lattice constant and Hubbard U parameters are fixed by confirming that FM ground state and topological characters are invariant under certain conditions. We construct MLWF, the Fourier transform of Kohn-Sham wavefunctions of each band, and calculate Berry curvatures to analyze topological characters using converged MLWF. Nontrivial topology appears within Cr e g band manifold which can lead to Chern insulating phases in CrSiTe 3 . Chiral edge modes are also calculated, consistent with total Chern number calculations. Moreover, we find spin angle-dependent AHC under one-electron doping in the formula unit cell and suppression of magnetic anisotropy. For our calculation, oneelectron doping in the formula unit cell with 7.0Å lattice constant corresponds to carrier density as 2.36×10 14 cm −2 . Ionic liquid gating is representative method of charge doping where order of doping concentration is known as 10 14 cm −2 in 2D systems [46][47][48] . In recent study, about 2×10 14 cm −2 of electron doping for one layer is accomplished in CrGeTe 3 35 almost adjacent to our theoretical goal. Because CrGeTe 3 has the same crystal structure and similar lattice constant as CrSiTe 3 , we expect our theoretical suggestion of Chern insulating phases and magnetic field dependent QAHE can be realized as shown throughout section II. By performing tightbinding (TB) model analysis, we confirm that DFT calculation results of Chern number for Cr e g orbitals are reproduced. We also find distinct Chern insulating phases by adjusting hopping parameters in our constructed TB model and drawing rel-evant phase diagrams. We expect that CrSiTe 3 can be a Chern insulator controlled by the external magnetic field, which can be realized by electron doping.

Density Functional Theory Calculations
To obtain band structures and projected density-of-states (pDOS), we perform ab initio electronic structure calculations based on density-functional theory (DFT) using OpenMX 49,50 code, which employs linear-combination-of-pseudo-atomicorbital basis with norm-conserving pseudopotentials. We choose generalized gradient approximation (GGA) exchangecorrelation functional in the parameterization of Perdew, Burke and Enzerhof (PBE) 51 with Hubbard U parameters chosen to be 1.5 eV for Cr d orbital. SOC effects are included in the calculations via fully-relativistic pseudopotentials implemented in OpenMX 49,50 . Pseudo-atomic orbitals are set to be s 3 p 2 d 2 for Cr, s 2 p 2 d 1 for Si, s 3 p 3 d 2 for Te, respectively. To simulate two-dimensionality, we insert 20 Å of vacuum in the unit cell. 10 × 10 × 1 of k-space mesh is adopted for the momentum space integration. Energy cutoff for choosing realspace grid is set to be 700 Ry (96 × 96 × 360 real space grid). 10 −6 Hartree/Bohr of force criterion is chosen for the optimization of internal coordinates while keeping C 3 rotation, inversion and mirror symmetries. SOC effects are excluded in the process of structural relaxation. We fix Bravais lattice as hexagonal and determine lattice constant by performing total energy minimization calculation as a function of unit cell size.

Analysis of topological characteristics
Berry curvature is given by where integrating Berry curvature in BZ gives Chern number Generalized tight-binding model for e g manifold In addition to the MLWF tight-binding model, we further introduce a Slater-Koster-type tight-binding model for a deeper understanding of SOC and lattice distortion effects.
We first start from including both Te p orbitals and Cr d orbitals, then project out Te p orbitals to obtain an effective four-band Hamiltonian.

DATA AVAILABILITY
All data are available from the corresponding author upon reasonable request.

SUPPLEMENTARY INFORMATION: GRAPHIC OF MLWF AND SIMPLE TIGHT BINDING MODEL APPROACH
We first start with a graphene-like case for a simple tightbinding approach as a toy model. If we consider s band like orbital with NN hopping in a honeycomb lattice, the Dirac cone appears at the K point. When 2 nd and 3 rd NN hopping are taken into account, the Hamiltonian is given by following expressions which becomes more complicated Component of above Hamiltonian f (k) and g(k) is written by where k = (k x , k y ), a is lattice constant and t i is i th NN hopping parameters. Diagonalizing the Hamiltonian expressed in Eq. 7, one can obtain equation of eigenenergy function as follows.
In this case, additional crossing point appears near the midpoint of Γ -K line under the condition of large 3 rd NN hopping parameter compared to 1 st and 2 nd NN hopping terms. This features are depicted in fig. 9. i j = i, 0 |Ĥ | j, r n , where i, j index indicates pair of MLWFs and r n is cell displacement vector which corresponds to n th NN site. Unit of hopping parameters is meV. Table III shows the overlap matrix from our MLWF calculation result, which corresponds to hopping integrals in TB model analysis. One can find a significant 3 rd NN term compared to other terms which satisfy the condition of the above simple model approach. It can be confirmed qualitatively in Fig. 10 which shows converged MLWF graphically. Hopping between NN sites is small due to the orthogonality of mediating Te p orbitals. For 2 nd and 3 rd NN hopping, mediating two Te p orbitals are located in different layers for 2 nd  NN case while in the same layer for 3 rd NN case, which induces much significant 3 rd NN hopping parameters compared to 1 st and 2 nd ones. This is a simple two-band model with slike orbital to understand the shape of band structure briefly; thus, we considered a more realistic tight-binding model considering up to 5 th NN hopping with four band model in the tight-binding analysis part.