First Principles Modeling of Topological Insulators: Structural Optimization and Exchange Correlation Functionals

TL;DR

This paper evaluates various density functional theory methods for modeling topological insulators, focusing on structural optimization and electronic properties, and recommends GGA with van der Waals interactions as the most suitable approach.

Contribution

It provides a comprehensive assessment of DFT formalisms for TIs and identifies the most accurate and computationally efficient method for structural and electronic property predictions.

Findings

GGA with vdW interactions best predicts atomic structures

MetaGGA offers slightly higher accuracy in electronic structure

LDA tends to underestimate TI slab thickness

Abstract

Topological insulators (TIs) are materials that are insulating in the bulk but have zero band gap surface states with linear dispersion and are protected by time reversal symmetry. These unique characteristics could pave the way for many promising applications that include spintronic devices and quantum computations. It is important to understand and theoretically describe TIs as accurately as possible in order to predict properties. Quantum mechanical approaches, specifically first principles density functional theory (DFT) based methods, have been used extensively to model electronic properties of TIs. Here, we provide a comprehensive assessment of a variety of DFT formalisms and how these capture the electronic structure of TIs. We concentrate on Bi$_2$Se$_3$ and Bi$_2$Te$_3$ as examples of prototypical TI materials. We find that the generalized gradient (GGA) and kinetic density functional (metaGGA) produce displacements increasing the thickness of the TI slab, whereas we see an opposite behavior in DFT computations using LDA. Accounting for van der Waals (vdW) interactions overcomes the apparent over-relaxations and retraces the atomic positions towards the bulk. Based on an intensive computational study, we show that GGA with vdW treatment is the most appropriate method for structural optimization. Electronic structures derived from GGA or metaGGA employing experimental lattice parameters are also acceptable. In this regard, we express a slight preference for metaGGA in terms of accuracy, but an overall preference for GGA due to compensatory improvements in computability in capturing TI behavior.

Figures

Tables

Table 1. Table 1: Table showing the lattice contants a 𝑎 a and c 𝑐 c , Wyckoff reference of Bi, Se(1) and Se(2) positions, and optical bandgap of bulk Bi 2 Se 3 obtained from various treatments of XCFs in DFT and the corresponding experimental values.

$a$ (Å)	$c$ (Å)	$z$-Bi	$z$-Se(1)	$z$-Se(2)	E${}_{\text{g}}$ (eV)
4.108${}^{\text{a}}$	27.293${}^{\text{a}}$	0.4012${}^{\text{a}}$	0.0${}^{\text{a}}$	0.2093${}^{\text{a}}$	0.042 (0.326)${}^{\text{a}}$
4.188${}^{\text{b}}$	31.429${}^{\text{b}}$	0.3953${}^{\text{b}}$	0.0${}^{\text{b}}$	0.2216${}^{\text{b}}$	0.691 (0.302)${}^{\text{b}}$
4.174${}^{\text{c}}$	28.872${}^{\text{c}}$	0.4001${}^{\text{c}}$	0.0${}^{\text{c}}$	0.2125${}^{\text{c}}$	0.234 (0.198)${}^{\text{c}}$
4.157${}^{\text{d}}$	29.812${}^{\text{d}}$	0.3983${}^{\text{d}}$	0.0${}^{\text{d}}$	0.2160${}^{\text{d}}$	0.353 (0.058)${}^{\text{d}}$
4.140${}^{\text{e}}$	28.657${}^{\text{e}}$	0.4004${}^{\text{e}}$	0.0${}^{\text{e}}$	0.2121${}^{\text{e}}$	0.115 (0.278)${}^{\text{e}}$
4.143${}^{\text{f}}$	28.636${}^{\text{f}}$	0.4008${}^{\text{f}}$	0.0${}^{\text{f}}$	0.2117${}^{\text{f}}$	0.220${}^{\text{g}}$, 0.335${}^{\text{h}}$

Table 2. Table 2: Table showing the lattice contants a 𝑎 a and c 𝑐 c , Wyckoff reference of Bi, Se(1) and Se(2) positions, and optical bandgap of bulk Bi 2 Te 3 obtained from various treatments of XCFs in DFT and the corresponding experimental values.

$a$ (Å)	$c$ (Å)	$z$-Bi	$z$-Se(1)	$z$-Se(2)	E${}_{\text{g}}$ (eV)
4.356${}^{\text{a}}$	29.903${}^{\text{a}}$	0.4016${}^{\text{a}}$	0.0${}^{\text{a}}$	0.2089${}^{\text{a}}$	0.138 (0.130)${}^{\text{a}}$
4.445${}^{\text{b}}$	32.247${}^{\text{b}}$	0.3970${}^{\text{b}}$	0.0${}^{\text{b}}$	0.2163${}^{\text{b}}$	0.555 (0.031)${}^{\text{b}}$
4.432${}^{\text{c}}$	30.446${}^{\text{c}}$	0.3998${}^{\text{c}}$	0.0${}^{\text{c}}$	0.2105${}^{\text{c}}$	0.198 (0.244)${}^{\text{c}}$
4.411${}^{\text{d}}$	31.474${}^{\text{d}}$	0.3981${}^{\text{d}}$	0.0${}^{\text{d}}$	0.2139${}^{\text{d}}$	0.391 (0.079)${}^{\text{d}}$
4.400${}^{\text{e}}$	30.096${}^{\text{e}}$	0.4001${}^{\text{e}}$	0.0${}^{\text{e}}$	0.2096${}^{\text{e}}$	0.179 (0.169)${}^{\text{e}}$
4.386${}^{\text{f}}$	30.497${}^{\text{f}}$	0.4000${}^{\text{f}}$	0.0${}^{\text{f}}$	0.2097${}^{\text{f}}$	0.171${}^{\text{g}}$

Table 3. Table A1: Table showing | ∂ E ∂ k | 𝐸 𝑘 \left|\frac{\partial E}{\partial k}\right| in eV/Å and effective mass m ∗ superscript 𝑚 m^{*} in units of rest mass of electron (refer Eqs. ( 1 ) and ( 3 )) of highest occupied band ( v 𝑣 v ) and lowest unoccupied band ( c 𝑐 c ) along the direction Γ → → Γ absent \Gamma\rightarrow M and Γ → → Γ absent \Gamma\rightarrow K of two-dimensional Brillouin zone for slab models of Bi 2 Se 3 using various XCF treatment. The parity is estimated analyzing the c 𝑐 c and v 𝑣 v band occupancies following the algorithm discussed in Fig.  8 .

LDA	Model	$P$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{c},\Gamma \to M}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{c},\Gamma \to K}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{v},\Gamma \to M}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{v},\Gamma \to K}$	${\left|m^{\ast }\right|}_{\text{c},\Gamma \to M}$	${\left|m^{\ast }\right|}_{\text{c},\Gamma \to K}$	${\left|m^{\ast }\right|}_{\text{v},\Gamma \to M}$	${\left|m^{\ast }\right|}_{\text{v},\Gamma \to K}$
	1QL	$+1$	0.99	1.12	1.12	1.27	0.15	0.16	0.15	0.17
	2QL	$+1$	1.39	1.57	0.41	0.48	0.13	0.14	0.39	0.52
	3QL	$-1$	1.86	2.00	0.06	0.03	0.15	0.16	0.29	0.42
	4QL	$-1$	1.88	2.02	0.09	0.00	0.15	0.16	0.28	0.41
	5QL	$-1$	1.95	2.09	0.16	0.07	0.16	0.16	0.26	0.37
	6QL	$-1$	1.96	2.10	0.18	0.08	0.16	0.17	0.26	0.37
	7QL	$-1$	1.97	2.10	0.19	0.09	0.16	0.17	0.26	0.37
	8QL	$-1$	1.95	2.09	0.18	0.08	0.16	0.17	0.26	0.38
PBE	1QL	$+1$	0.87	0.99	0.95	1.08	0.16	0.17	0.17	0.19
	2QL	$+1$	0.98	1.12	0.24	0.27	0.15	0.16	1.04	1.52
	3QL	$+1$	1.02	1.16	0.13	0.15	0.14	0.15	1.78	2.62
	4QL	$+1$	1.04	1.18	0.09	0.11	0.14	0.15	2.31	3.43
	5QL	$+1$	1.10	1.26	0.04	0.04	0.14	0.15	2.68	2.23
	6QL	$+1$	1.08	1.23	0.03	0.03	0.14	0.15	8.51	33.0
	7QL	$+1$	1.27	1.44	0.29	0.32	0.13	0.14	0.68	0.71
	8QL	$+1$	1.26	1.43	0.26	0.29	0.13	0.14	0.74	0.77
SCAN	1QL	$+1$	0.97	1.10	0.96	1.09	0.15	0.16	0.16	0.17
	2QL	$+1$	1.13	1.29	0.36	0.40	0.13	0.14	0.84	1.47
	3QL	$+1$	1.25	1.42	0.11	0.12	0.13	0.14	4.04	1.64
	4QL	$+1$	1.50	1.68	0.11	0.12	0.13	0.14	1.49	0.91
	5QL	$+1$	1.58	1.76	0.25	0.27	0.13	0.14	1.18	0.80
	6QL	$+1$	1.67	1.84	0.41	0.43	0.13	0.15	1.07	0.81
	7QL	$+1$	1.82	1.98	0.55	0.56	0.14	0.15	1.87	1.01
	8QL	$+1$	1.86	2.02	0.59	0.59	0.15	0.16	2.25	1.05
PBE+vdW	1QL	$+1$	0.87	0.99	0.98	1.11	0.16	0.17	0.17	0.19
	2QL	$+1$	1.11	1.26	0.26	0.28	0.13	0.14	1.79	11.9
	3QL	$+1$	1.45	1.62	0.14	0.14	0.13	0.14	1.87	1.02
	4QL	$+1$	1.75	1.91	0.44	0.44	0.15	0.16	49.3	1.37
	5QL	$+1$	1.88	2.03	0.62	0.61	0.15	0.17	11.6	1.44
	6QL	$+1$	2.11	2.23	0.82	0.79	0.18	0.19	0.92	22.2
	7QL	$+1$	2.08	2.21	0.79	0.76	0.18	0.19	0.96	21.3
	8QL	$+1$	2.18	2.29	0.90	0.86	0.19	0.20	0.76	6.11
SCAN+vdW	1QL	$+1$	0.97	1.10	0.99	1.12	0.15	0.16	0.16	0.17
	2QL	$+1$	1.25	1.42	0.42	0.47	0.12	0.13	0.62	1.04
	3QL	$+1$	1.67	1.84	0.02	0.02	0.13	0.14	0.76	1.91
	4QL	$-1$	1.90	2.05	0.26	0.21	0.14	0.16	0.53	1.35
	5QL	$-1$	2.01	2.15	0.37	0.30	0.16	0.17	0.41	0.96
	6QL	$-1$	2.07	2.20	0.44	0.37	0.16	0.17	0.40	0.94
	7QL	$-1$	2.08	2.21	0.45	0.38	0.17	0.18	0.38	0.85
	8QL	$-1$	2.10	2.23	0.48	0.40	0.17	0.18	0.39	0.85

Table 4. Table A2: Table showing | ∂ E ∂ k | 𝐸 𝑘 \left|\frac{\partial E}{\partial k}\right| in eV/Å and effective mass m ∗ superscript 𝑚 m^{*} in units of rest mass of electron (refer Eqs. ( 1 ) and ( 3 )) of highest occupied band ( v 𝑣 v ) and lowest unoccupied band ( c 𝑐 c ) along the direction Γ → → Γ absent \Gamma\rightarrow M and Γ → → Γ absent \Gamma\rightarrow K of two-dimensional Brillouin zone for slab models of Bi 2 Te 3 using various XCF treatment. The parity is estimated analyzing the c 𝑐 c and v 𝑣 v band occupancies following the algorithm discussed in Fig.  8 .

LDA	Model	$P$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{c},\Gamma \to M}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{c},\Gamma \to K}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{v},\Gamma \to M}$	${\left|\frac{\partial E}{\partial k}\right|}_{\text{v},\Gamma \to K}$	${\left|m^{\ast }\right|}_{\text{c},\Gamma \to M}$	${\left|m^{\ast }\right|}_{\text{c},\Gamma \to K}$	${\left|m^{\ast }\right|}_{\text{v},\Gamma \to M}$	${\left|m^{\ast }\right|}_{\text{v},\Gamma \to K}$
	1QL	$+1$	1.31	1.47	1.09	1.19	0.14	0.16	0.30	0.12
	2QL	$-1$	1.13	1.27	1.54	1.64	0.17	0.18	0.53	8.42
	3QL	$-1$	1.36	1.53	0.72	0.79	0.16	0.17	0.37	1.23
	4QL	$-1$	1.68	1.82	0.41	0.52	0.21	0.21	0.20	0.37
	5QL	$+1$	1.59	1.75	0.52	0.62	0.20	0.20	0.20	0.33
	6QL	$+1$	1.67	1.82	0.44	0.55	0.21	0.21	0.18	0.29
	7QL	$+1$	1.69	1.84	0.42	0.52	0.22	0.22	0.18	0.28
	8QL	$+1$	1.70	1.84	0.40	0.51	0.22	0.22	0.18	0.28
PBE	1QL	$+1$	1.25	1.40	1.07	1.18	0.15	0.17	1.31	0.34
	2QL	$+1$	2.25	2.41	0.44	0.41	0.22	0.26	0.95	2.81
	3QL	$+1$	2.40	2.54	0.63	0.64	0.28	0.31	0.26	0.87
	4QL	$+1$	2.74	2.87	1.02	1.04	0.45	0.42	0.20	0.98
	5QL	$+1$	2.55	2.68	0.68	0.67	0.45	0.37	0.18	0.45
	6QL	$-1$	2.86	2.99	1.08	1.10	1.12	0.51	0.15	0.39
	7QL	$+1$	2.69	2.82	0.83	0.79	0.95	0.43	0.17	0.39
	8QL	$+1$	2.81	2.93	1.01	0.99	3.09	0.50	0.14	0.32
SCAN	1QL	$+1$	2.13	2.29	0.04	0.09	0.19	0.23	0.93	14.7
	2QL	$-1$	1.13	1.29	0.36	0.40	0.13	0.14	0.84	1.47
	3QL	$+1$	1.99	2.17	0.22	0.25	0.19	0.20	0.27	1.34
	4QL	$+1$	2.32	2.47	0.17	0.12	0.28	0.27	0.19	0.43
	5QL	$-1$	2.42	2.56	0.27	0.21	0.35	0.30	0.19	0.44
	6QL	$+1$	2.46	2.60	0.30	0.23	0.38	0.31	0.19	0.43
	7QL	$+1$	2.42	2.62	0.31	0.24	0.42	0.31	0.19	0.43
	8QL	$+1$	2.49	2.63	0.33	0.25	0.47	0.33	0.18	0.39
PBE+vdW	1QL	$+1$	1.21	1.36	1.05	1.15	0.15	0.17	1.24	0.20
	2QL	$+1$	1.36	1.53	0.79	0.83	0.16	0.18	2.79	1.05
	3QL	$+1$	2.05	2.18	0.08	0.01	0.27	0.26	0.22	0.53
	4QL	$+1$	1.83	1.99	0.19	0.25	0.22	0.22	0.22	0.50
	5QL	$+1$	2.03	2.17	0.05	0.03	0.28	0.26	0.19	0.35
	6QL	$+1$	2.09	2.22	0.10	0.01	0.30	0.27	0.19	0.35
	7QL	$+1$	2.11	2.24	0.12	0.03	0.31	0.28	0.18	0.35
	8QL	$+1$	2.05	2.19	0.05	0.03	0.29	0.26	0.20	0.38
SCAN+vdW	1QL	$+1$	1.30	1.45	1.10	1.20	0.14	0.16	0.60	0.16
	2QL	$-1$	1.30	1.47	1.16	1.23	0.15	0.17	0.91	1.97
	3QL	$+1$	1.89	2.05	0.26	0.33	0.20	0.20	0.25	0.65
	4QL	$+1$	1.82	1.99	0.39	0.47	0.19	0.20	0.21	0.45
	5QL	$+1$	2.01	2.16	0.18	0.27	0.23	0.23	0.18	0.32
	6QL	$+1$	2.09	2.23	0.09	0.20	0.25	0.24	0.17	0.29
	7QL	$+1$	2.07	2.21	0.11	0.22	0.25	0.24	0.17	0.29
	8QL	$+1$	2.14	2.28	0.05	0.16	0.27	0.25	0.17	0.32

Equations

$$E(k)=E_0+\left(\frac{\partial E}{\partial k}\right)k+\frac{1}{2}\left(\frac{{\partial }^{2}E}{\partial k^2}\right)k^2+O(k^3),$$

$$\frac{{\partial }^{2}E}{\partial k^2}\to 0.$$

$$m^{\ast }=\frac{1}{{\hslash }^2}\left(\frac{{\partial }^{2}E}{\partial k^2}\right).$$

Full text

I Introduction

Within the last decade, it was reported that the pnictogen chalcogenides Bi2Se3 and Bi2Te3 exhibit the properties of a three-dimensional topological insulator (TI)—small energy gap at the Fermi level, inverted parity of the band edge states leading to a Dirac band dispersion, and quantum oscillations for topological surface states (TSS)  Zhang et al. (2009); Xia et al. (2009); Hsieh et al. (2009a); Hasan and Kane (2010); Moore (2010); Zhang et al. (2010); Taskin et al. (2012); Tsipas et al. (2014). This discovery has resulted in significant interest in these two materials, as their relatively large bulk bandgaps and chemical simplicity make the study of their topological electronic physics readily accessible to theory and experiment work alike. TIs and TSS are projected to be the basis of many cutting-edge device applications, including spintronics and quantum computers Das Sarma et al. (2006); Wang et al. (2016); Fan and Wang (2016). TIs also offer fertile ground for fundamental studies of exotic electronic phenomena arising from wave function topology. Hasan and Kane (2010); Zhang et al. (2013a, b); Kane and Zhang (2015); Metlitski et al. (2015); Teo et al. (2008); Hsieh et al. (2009b); Fu and Kane (2008); Kane and Mele (2005); Wehling et al. (2014)

Structurally, Bi2Se3 and Bi2Te3 are most often described in terms of quintuple layers (QLs), in which atomic layers are arranged in sets of five sublayers along the out-of-plane $c$-axis of Bi2Se3 and Bi2Te3, a configuration called the quintuple layer (QL). The arrangement of atoms in a QL is B(1)-A-B(2)-A-B(1) (for A2B3) along the [0001]-direction of the hexagonal crystal system with space group $\text{R}\mathrm{bar}{3}\text{m}$ (No. 166) Nakajima (1963). The layers in a QL are bounded by stronger chemical bonding, but the inter-QL interaction is comparatively weaker and of van der Waals (vdW) type.

The number of QLs needed to achieve a robust TSS is of primary interest for the accurate description of the TI behavior. Different factors such as quantum effects due to finite size, undercoordinated surface atoms, and hybridization of orbitals from the surface and bulk atoms also contribute to the emergence or suppression of a TSS. Zhang et al. Zhang et al. (2010) have reported the appearance of a Dirac point in samples with thicknesses 6QLs and greater for Bi2Se2, while, for Bi2Te3, Liu et al. Liu et al. (2012) have reported that 4QLs are sufficient. Using a combined model Hamiltonian study and first-principles calculations Liu et al. Liu et al. (2010) have reported an oscillatory crossover from a two-dimensional to a three-dimensional TI as a function of film thickness (number of QLs). Since both intra-QL and inter-QL interactions are crucial to stabilize the TSS Young et al. (2011); Liu et al. (2011), it is important to standardize the settings of density functional theory (DFT) modeling to converge towards a common perspective as applied to this general class of materials.

DFT is one of the most popular theoretical tools to study structural and electronic properties of materials. It relies upon the Kohn-Sham formulation, where the knowledge of electron density is sufficient to calculate several materials properties that depend on the electronic and atomic structure. It is an efficient means of establishing the ground state, even for relatively large systems. Furthermore, because it is a first-principles method, it allows an unambiguous comparison to model experiments. It is important to note that DFT’s robustness is not based only in its capacity for verification and explanation of experimental results. It also allows for predicting properties for previously unexplored systems to guide experimental work. While the validity of principles of DFT is established, the deficiencies of exchange-correlation functionals (XCFs) may limit the quality of its predictive capabilities. A substantial amount of knowledge has been accumulated to quantify these shortcomings, and ways to circumvent them have been suggested with mixed success of transferability Cohen et al. (2008, 2012); Santra and Perdew (2019); Zunger et al. (2010); Adeagbo et al. (2014); Nayak et al. (2015); Dion et al. (2004); Berland et al. (2019); Liechtenstein et al. (1995); Becke and Johnson (2006); Zhao and Truhlar (2006).

For TI materials like Bi2Se3 and Bi2Te3, DFT can be utilized in two major ways: The first is the description of the ground-state atomic positions (relaxation). The second is the description of the converged charge density. Both play a significant role in quantifying TSS and the selection of the appropriate XCFs. TI materials have been studied with XCFs including the local density approximation (LDA) and the generalized gradient approximation (GGA) Zhang et al. (2009, 2010); Yazyev et al. (2012, 2010); Li et al. (2014); Tao and Tsymbal (2018); Datzer et al. (2017); Xia et al. (2009); Seixas et al. (2015); Seibel et al. (2016); Rusinov et al. (2013), but the role of the kinetic density functional (metaGGA) has not been hitherto explored. The many-body perturbation theory in the $GW$ approximation, where $G$ is the single particle Green function and $W$ is the screened Coulomb interaction, has been applied to improve description of the band structure Aguilera et al. (2013); Förster et al. (2015, 2016). It is also important to realize that atomic structure is intricately related to electronic properties. For example, the $GW$ approximations are often applied to experimental or DFT-derived atomic structures. For layered TI materials like Bi2Se3 and Bi2Te3, treatment of the vdW forces has proven particularly important, and many studies have demonstrated that doing so greatly affects DFT’s description of electronic structures. This is because, for a slab or film of a fixed number number of QLs, a relationship exists between the thickness of the slab or film (i.e. due to strain) and the possibility of observing a TSS. VdW forces mediate the interactions between QLs across the inter-QL space—depending on the magnitude of their influence, they create smaller or larger vdW gaps. As a result, an XCF treated to explicitly account for the vdW forces might predict a different overall thickness than one that has not been, and, thus, it might predict a different surface state Liu et al. (2013); Luo et al. (2012). However, it must be noted that the vast majority of the thickness of a slab is made up of the QLs themselves. In order to satisfactorily account for a given XCF’s contribution to change in thickness relative to the experimental parameters and those produced by other XCFs, we must assess the role played by the intra-QL space in overall thickness as well, which has not been taken into account by previous studies.

In this study, we carry out a systematic analysis in pursuit of a consensus regarding the role of XCFs in the outcome of structural optimization and its subsequent influence over electronic structure calculations in Bi2Se3 and Bi2Te3. We specifically asses a variety of XCFs including: LDA, GGA, metaGGA, GGA+vdW, and metaGGA+vdW. We find that the structural optimization process is highly sensitive to the type of XCF employed. LDA and GGA produce relaxation trends opposite of each other. The effective role of the vdW treatment is to bring the relaxation toward the experimental bulk positions. While the outcomes of the electronic structure calculations are dependent on the specific XCF, we used the vdW treatments of Grimme et al. Grimme et al. (2010), which do not depend on charge density, making it redundant to apply them to the self-consistent charge density (SCCD) and band structure phases of our calculations. As a result, our vdW-treated XCFs preserve the charge density distribution generated by their root functionals. It is important to note that our structural optimizations with the vdW treatments indicate little deviation (below 2%) from experimental bulk coordinates of Nakajima (Ref. Nakajima (1963)). Therefore, we deduce that the experimental bulk coordinates of Nakajima are sufficient to create structural models of Bi2Se3 and Bi2Te3 in which TSSs that exhibit thickness-dependence in agreement with experiment can be expected, so long as the GGA or metaGGA XCFs (the root functionals of our vdW-corrected XCFs) are used. Lastly, mindful of the range of XCFs, number of QLs (slab thicknesses), and computational strategies pursued in the theoretical literature on TSSs of Bi2Se3 and Bi2Te3, we extend these conclusions to a prescription for a general method for generating models suitable for studying TI physics in binary layered systems of space group $R\mathrm{bar}{3}m$ (No. 166), with DFT.

In this way, we seek to evaluate with new precision the outcome of structural optimizations with a variety of XCFs and benchmark a procedure for reliably modeling TSSs in these systems. We expect this study should act as the groundwork for the production of models that might be used to generate data that would serve the development of novel devices.

II Methodology

Calculations were performed with DFT using the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996a, b) which solves the scalar-relativistic Kohn-Sham Hamiltonian. In all calculations, spin-orbit coupling was treated with the perturbation theory. We used the following three root XCFs: LDA with Ceperley-Alder parametrization Ceperley and Alder (1980), GGA with Perdew-Burke-Ernzerhof parametrization (PBE) Perdew et al. (1996), and metaGGA with Strongly Constrained and Appropriately Normed (SCAN) approach Sun et al. (2015). Slab models were prepared from the unit bulk hexagonal lattice with lattice parameters taken from experiments Nakajima (1963). This set of experimental lattice parameters matches the DFT band structure as reported in Ref. Luo et al. (2012). The schematic picture of the bulk hexagonal unitcell and the slab models are shown in Fig. 1(a–c). The slabs were centered along the vertical direction. A vacuum layer of thickness of $\sim$8.8 Å was added to either side of the surfaces making a total vacuum thickness of $\sim$17.6 Å. This range of vacuum thickness ensures negligible dipole interaction between the surfaces through the vacuum and is regularly used for modeling surfaces and nanomaterials Sahoo et al. (2018, 2010). A systematic study from one quintuple layer (1QL) to eight quintuple layers (8QL) models were performed. The calculations were performed using the projector-augmented wave pseudopotentials, a symmetrized 7$\times$7$\times$1 $k$-point grid to span the two-dimensional Brillouin zone, the kinetic energy cutoff of the plane-waves set to 500 eV and the precision tag set to ‘accurate’. An additional support grid was included in the calculation through the tag ‘addgrid’, which accounts for accurate inter-atomic forces and hence leads to higher quality geometrical optimizations. The total energy convergence for self-consistent field cycles is set to $10^{-7}$ eV and the force convergence for geometrical optimization is set to $0.001$ eV/Å.

Each QL of a multi-QL system is weakly bounded to the others by the vdW force. Past works have demonstrated the importance of accounting for the influence of the vdW forces in layered systems like Bi2Se3 and Bi2Te3 in DFT calculations by using specially adapted treatments Luo et al. (2012); Liu et al. (2013). We incorporated such treatments into another set of calculations with the GGA functional using the zero-damping DFT-D3 method of Grimme et al. Grimme et al. (2010). Although the metaGGA-SCAN functional accounts for vdW interactions to some extent, we wanted to cross-check the outputs including the vdW interactions, so a set of calculations using the SCAN+vdW treatment was carried out.

III Structural Properties

The experimental bulk lattice parameters and corresponding results from the DFT calculations for Bi2Se3 and Bi2Te3 are shown in Tab. 1 and Tab. 2, respectively. Our calculations are delineated by rows by XCF. The data shows that SCAN+vdW predictions of lattice parameters most closely reflect experimental results, followed in this regard by PBE+vdW. LDA and PBE XCFs deviate more from the experimental values, with LDA predicting structural parameters lower than experiment and PBE predicting structural parameters higher than experiment. Relaxation along $a$ resulted in far less change overall compared with relaxations along $c$.

The slab models generated from the bulk experimental coordinates of Ref. Nakajima (1963) show variation in the direction and extent of relaxation with different XCF treatments. The relative change of thickness of the slab $\Delta t$ = $t_{\text{XCF}}-t_{\text{Expt.}}$, where $t_{\text{XCF}}$ and $t_{\text{Expt.}}$ are the optimized thickness and ideal bulk-cut thickness, respectively, as a function of number of QLs for Bi2Se3 and Bi2Te3 is shown in Fig. 2(a) and Fig. 2(b), respectively. The horizontal line passing through zero for $\Delta t$ represents the thickness of the structure with experimental lattice parameters. For all calculations conducted with the LDA functional, a negative value of $\Delta t$ is observed. Naturally, this implies the LDA functional induces a reduction in the thickness of all slabs. An opposite trend is observed for the PBE functional, implying that PBE predicts the thickness of any slab should increase. The SCAN functional, for both materials, predicts a less severe increase of the slab thickness than PBE. When we apply the vdW treatments to both PBE and SCAN, both functionals predict that the thickness of the slab should be bulk-like. These trends, on a functional-by-functional basis, are reflected in the resultant bulk coordinates of each XCF (Tab. 1 and Tab. 2). By introducing SOC, the bandgap shrinks. This is consistent with the expectation that these materials are composed of high-$Z$ elements that have sizable SOC interactions. It is notable, however, that, for the Bi2Se3 LDA model, enabling SOC increases the bandgap. This result is opposite those for all other bulk models of both Bi2Se3 and Bi2Te3. We surmise that this is due to the large reduction in the $a$ lattice parameter.

While it is obvious each functional affects the general thickness of the slab, a focus on the thickness alone obscures the control each XCF exerts over the role of more specific mechanisms in the outcome of a given structural optimization. In the left panels of Fig. 3, we illustrate the change in the position of each atom along $c$ relative to its original index position for 3QL and 6QL models for Bi2Se3 (Fig. 3(a) and (b)) and Bi2Te3 (Fig. 3 (c) and (d)), respectively, taking into account the five discussed XCFs. The $\Delta z$ (= $z_{\text{relaxed}}$ $-$ $z_{\text{initial}}$) is the amount of change incurred by an atom in comparison to its position in the initial model generated from the experimental bulk lattice parameters. Since the zeroth index in the left panel of Figs. 3(a)–(d) is set to the center layer of the slab for models with an odd number of QLs and one of the two center layers for models with an even number of QLs, an atom’s position relative to the horizontal at zero is indicative of the direction of its relaxation. Presented such a way, the positive value of $\Delta z$ for the negative layer indices and negative value of $\Delta z$ for the positive layer indices imply inward relaxation of the two ends of the slab. Similarly, the negative value of $\Delta z$ for the negative layer indices and positive value of $\Delta z$ for the positive layer indices imply outward relaxation. The trend of the change in the layer positions within a QL—identifiable as the five-point clusters demarcated by the grey boxes—is generally linear, inward or outward, for all calculations. A step between each QL represents the vdW gap (not to be confused with the vdW treatments for XCFs). PBE predicts the largest outward change in the size of the vdW gap relative to the experimental observations. PBE+vdW and SCAN+vdW, on the other hand, predict a far smaller change in the vdW gap, slightly inward or outward depending on the system and the thickness of the slab (refer to the left panel). The remaining functionals predict absolute changes in the vdW gap larger than the vdW-corrected functionals, but smaller than pure PBE. Comparing to the bulk, which is the horizontal line passing through $\Delta z=0$ in the left panel of Figs. 3(a)–(d), it is observed that SCAN+vdW and PBE+vdW deviate minimally from the corresponding bulk atomic positions. Roy et al. Roy et al. (2014a, b) have performed surface X-ray diffraction studies on (0001) Bi2Se3 grown as an ultrathin film on Si(111) using MBE. They find outward relaxation of the top Se-Bi layer by $\sim$2%-4% as compared to the bulk coordinates, and $\sim$3% contraction of the top inter-QL separation compared to the same. Our results suggest that PBE+vdW relaxation match more closely to the experimental results for the outer Se-Bi layers. However, the contraction of the top inter-QL layer is not represented in our calculations.

It is clear from the right panel of Figs. 3(a)–(d) that the majority of the thickness of each relaxed slab is occupied by the intra-QL space and a minority by the inter-QL space. That the thickness of the vdW-corrected functionals maintains approximately the same ratio of intra-QL space to inter-QL space as the other functionals, but exhibits far lower thickness in general, suggests two possibilities: 1. By merit of the large absolute differences between the thicknesses of the non-vdW-corrected and vdW-corrected models, that the vdW force is the primary mediating force during structural optimization. 2. That the vdW force controls the inter-QL distances as well as the the intra-QL distances.

For LDA, PBE, and SCAN models, even when they are as large as 8QL, convergence to the experimental bulk characteristics is not obtained. On the other hand, it should again be noted that the relaxations of the vdW-corrected models, with a slight bias in favor of the SCAN+vdW models in particular, show strong convergence to the corresponding experimental bulk for any number of QLs. This suggests that slabs prepared from the experimental coordinates, and not subject to structural optimization, adequately capture the physics accounted for by the vdW treatments, and are reasonable for studying the electronic structure. But, a conclusion can not be drawn without evaluating the electronic structure directly. This is done in the next section.

In Fig. 4, we plot the relationship between the ratio of the optimized inter-QL space (the most outward intra-QL), $d_{\text{QL}}$, and inter-QL space (the top Bi-Te or Bi-Se layer separation), $d_{\text{iQL}}$, to their respective ultrathin film experimental values. The experimental values are taken from Ref. Roy et al. (2014a) for Bi2Se3 and Ref. Fukui et al. (2012) for Bi2Te3. In both panels, the horizontal and vertical dashed lines intersect at the experimental data point, and the diagonal dashed lines represent an extrapolation of the ratio of the portions of the total thickness occupied by the inter-QL space to the intra-QL space in the experimental films for slabs of variable thickness. The advantage of this set of plots is that the contributions of the changes in the $d_{\text{iQL}}$ and the $d_{\text{QL}}$ to the overall change in thickness after relaxation can be accounted for explicitly. For Bi2Se3, PBE induces a large change in the $d_{\text{iQL}}$; enabling the vdW treatment greatly reduces the predicted growth in the $d_{\text{iQL}}$, while also inducing a comparable reduction in the $d_{\text{QL}}$. By contrast, for Bi2Te3, PBE greatly increases both $d_{\text{QL}}$ and $d_{\text{iQL}}$ inter-QL separation, and in a fashion that deviates strongly from the experimental ratio (the dashed diagonal line). Enabling the vdW treatment leads to a convergence of the data to the experimental ratio, yet not to the experimental values reported in Ref. Fukui et al. (2012). This again implies that the vdW treatment controls the overall thickness by changing both the $d_{\text{iQL}}$ and $d_{\text{QL}}$. For both of the materials, a similar trend can be seen with the SCAN functional and its vdW-treated counterpart. For both materials, LDA predicts a smaller $d_{\text{iQL}}$ and $d_{\text{QL}}$.

IV Electronic Properties

In this section, we seek to relate the selection of the XCF to the relationships between sample thickness and the TSS observed in experiment using the optimized structures generated in the previous section.

First, however, a 6QL model generated from experimental bulk lattice parameters is chosen to show the $xy$-integrated charge density as a function of $z$. Further, the charge density is compared with the LDA-obtained charge density as the reference by the construct $\left(\rho_{z}-\rho_{z}\text{(LDA)}\right)/\rho_{z}\text{(LDA)}$. The comparison of charge density is done for the outer QL and one of the inner QLs (third QL from the outer QL) of Bi2Se3 in Fig. 5(a) and Fig. 5(b), respectively, and similarly for Bi2Te3 in Fig. 5(c) and Fig. 5(d). The results demonstrate that the charge densities of PBE and SCAN show different trends as compared to LDA, and do not change with respect to where the QL is located relative to the center of the model. These 6QL unrelaxed models give band gap values of 0.0002 eV, 0.0027 eV, 0.0027 eV, 0.0050 eV, and 0.0050 eV for LDA, PBE, PBE+vdW, SCAN and SCAN+vdW XCFs for Bi2Se3. Similarly for Bi2Te2, the bandgaps are 0.0007 eV, 0.0012 eV, 0.0012 eV, 0.0008 eV, and 0.0008 eV for LDA, PBE, PBE+vdW, SCAN and SCAN+vdW, respectively. The data suggests that the bandgaps are controlled by the root functionals, with no change appearing with the addition of vdW treatments.

Fig. 6 shows the evolution of the bandgap as a function of slab thickness for Bi2Se3 and Bi2Te3 with all XCFs. It is evident that the best agreement to experimental trends is achieved for Bi2Se3 by the SCAN functional, followed closely by the vdW-corrected functionals and the $GW$ method. For Bi2Te3, best agreement to experimental trends is achieved by the SCAN+vdW functional, followed closely by LDA and PBE+vdW. Given issues of computational efficiency and convergence we confronted with the SCAN functional, and the known tendency of LDA to underestimate bandgap, we conclude PBE+vdW strikes the best balance between computability and accuracy of all the functionals tested for both systems.

We now attempt to theoretically evaluate the TSSs of our models by examining other properties characteristic of the TSS, besides the bandgap: band parity, energy gradient with respect to $k$, and effective mass of the highest occupied band and lowest unoccupied band. The ideal TSS is characterized by linear band dispersion and a bandgap tending to zero. For $k$ close to $\Gamma$, the band energy can be expanded in $k$ using the Taylor’s expansion.

[TABLE]

where $E_{0}$ is band energy at the $\Gamma$-point. In the limit of the ideal TSS;

[TABLE]

In real materials, we see a deviation from this ideal condition. Knowledge of $\left(\frac{\partial E}{\partial k}\right)$ and $\left(\frac{\partial^{2}E}{\partial k^{2}}\right)$ at the $\Gamma$-point provides an idea of the degree of deviation from the ideal TSS. The quantity $\left(\frac{\partial^{2}E}{\partial k^{2}}\right)$ is an important property for carrier charges as it is related to the effective mass via the relation

[TABLE]

The calculated value of $\left(\frac{\partial E}{\partial k}\right)$ and $m^{*}$ in the direction $\Gamma\rightarrow\text{M}$ and $\Gamma\rightarrow\text{K}$ obtained from finite difference method is shown in Tab. A1 and Tab. A2 for Bi2Se3 and Bi2Te3, respectively, for models varying in slab thickness (1QL–8QL) and XCFs. We find that the gradient $\partial E/\partial k$ of valence ($v$) and conduction bands ($c$) obtained from PBE have smaller values as compared to other XCFs for Bi2Se3. However, a similar tendency is not reflected for Bi2Te3. The effective masss of conduction band $m^{*}_{c}$ obtained from different XCF treatments are similar for Bi2Se3, while $m^{*}_{c}$ and $m^{*}_{v}$ have similar values (within reasonable error) for Bi2Te3. The parity is calculated with the algorithm taken from Ref. Geilhufe et al. (2015) (see Fig. 8, stated for Bi2Se3 as an example). The description of parity using this algorithm indicates the band inversion tendency for $P=-1$. Since the TSS is a limiting case, the parity of TSS can turn out to be $+1$ or $-1$ for such an algorithm.

Using a slightly different flavor of GGA (optPBE) and vdW treatments, the results in the band gap obtained for Bi2Se3 goes to zero smoothly as a function of thickness Liu et al. (2013), which is similar to our results. However, the band gap of Bi2Te3 appears to converge towards zero from 2QL onwards, which in our case (PBE+vdW) goes to zero for 3QL and larger, with all XCFs converged to zero from 4QLs. This is consistent with the experimental report of Ref. Liu et al. (2012). Comparing the trends that can be observed in Fig. 6, in which bandgap as a function of the number of QLs for both experimental and vdW-treated XCFs data is illustrated, we project that six QLs for Bi2Se3 and four QLs for Bi2Te3 are minimum models that capture a stable TSS.

Finally, in Fig. 7, we compare the band structures produced from models optimized with PBE and PBE+vdW for Bi2Se3 and Bi2Te3. The color coding of the band structure is based on the algorithm described in Fig. 8, as in Ref. Geilhufe et al. (2015), where $w$ is the average contribution of a type of atom for a given orbital to the band. For example, $w(\text{Bi}_{p})$ is the weight of $p$-orbital of Bi atom in a specific band. It must be emphasized that for Bi2Se3 and Bi2Te3 the $p$ orbitals constitute the valence band ($v$) and the conduction band ($c$) states near the Fermi level. As can be seen, the vdW-corrected models produce the gapless $\Gamma$-point states and linear dispersion expected at 6QLs for both systems. In contrast, the PBE-optimized structures do not. Given that both models are derived from the same charge density predictions, this figure demonstrates the importance of the vdW treatments in producing results that agree with experiment purely as a matter of atomic structure.

V Conclusions

Using prototypical TIs Bi2Se3 and Bi2Te3, we have here systematically assessed the role of several XCFs and relaxation effects, particularly as they relate to inter- and intra-QL separations, for prediction of TSS. We show here that there is an inextricable relationship between the number of quintuple layers, effects of structural optimization, and the electronic structure when different versions of XCFs are employed in first principles calculations carried out via DFT. A detailed analysis of the inter- and intra-layer relaxations reveal that GGA overestimates inter- and intra-QL separation, contradicting experimental findings Roy et al. (2014a). On the other hand, the relatively simple LDA functional produces an opposite trend in relaxation effects, also in contradiction of experimental results. When more complicated XCFs such as metaGGA are employed, the relaxed structural parameters are found to be between the values obtained from experimental work and GGA. With vdW treatments applied, GGA and SCAN predict bulk-like crystal structures. In addition to structural parameters, we also investigate the thickness dependence of the bandgap for Bi2Se3 and Bi2Te3 for all XCFs. We ultimately determine that GGA+vdW offers the most reliable results, for both accuracy when compared to experimental studies and computability. We also determine that the GGA+vdW functional produces a structure that does not deviate much from the bulk, and generates a thickness-bandgap dependence that agrees well with experimental results. We thus conclude that structural optimization may not be necessary to model TI physics in binary layered systems of space group $R\mathrm{bar}{3}m$ (No. 166) within a tolerable accuracy. It must be noted that modeling TI systems of this type with complex non-stoichiometric chemistry and doping effects would require structural optimization, where the GGA+vdW is most appropriate. Comparing the trends in bandgap as a function of number of quintuple layers for vdW-treated XCFs and experimental results, we conclude that the six QL model for Bi2Se3 and the four QL model for Bi2Te3 represent the minimum necessary thicknesses capable of capturing a stable TSS.

Acknowledgements

The authors would like to acknowledge useful discussions with Matthias Geilhufe and Gayanath Fernando. AVB acknowledges the support of VILLUM FONDEN via the Center of Excellence for Dirac Materials (Grant No.11744) and the European Research Council under the European Unions Seventh Framework Program (FP/2207-2013)/ERC Grant Agreement No. DM-321031 ERC.

The computational resources employed in this study were provided by High Performance Computing facility of University of Connecticut.

Appendix A Appendix A

Appendix B Appendix B