First-principles study of native point defects in Bi2Se3

Using first-principles method within the framework of the density functional theory, we study the influence of native point defect on the structural and electronic properties of Bi$_2$Se$_3$. Se vacancy in Bi$_2$Se$_3$ is a double donor, and Bi vacancy is a triple acceptor. Se antisite (Se$_{Bi}$) is always an active donor in the system because its donor level ($\varepsilon$(+1/0)) enters into the conduction band. Interestingly, Bi antisite(Bi$_{Se1}$) in Bi$_2$Se$_3$ is an amphoteric dopant, acting as a donor when $\mu$$_e$$<$0.119eV (the material is typical p-type) and as an acceptor when $\mu$$_e$$>$0.251eV (the material is typical n-type). The formation energies under different growth environments (such as Bi-rich or Se-rich) indicate that under Se-rich condition, Se$_{Bi}$ is the most stable native defect independent of electron chemical potential $\mu$$_e$. Under Bi-rich condition, Se vacancy is the most stable native defect except for under the growth window as $\mu$$_e$$>$0.262eV (the material is typical n-type) and $\Delta$$\mu$$_{Se}$$<$-0.459eV(Bi-rich), under such growth windows one negative charged Bi$_{Se1}$ is the most stable one.


I. INTRODUCTION
The narrow-band-gap semiconductor Bi 2 Se 3 (E g ∼0.35eV) 1,2 has been best known for a long time as an excellent thermoelectric material because of their unique near-gap electronic structure and high thermoelectric figure of merit ZT=S 2 σT/K. [3][4][5] Recently, with the development of topological insulators(TIs), [6][7][8] it attracts the research attention again. TI is a new state having an energy gap in its bulk band structure and metallic helical states on its surface, which distinct from simple metal or insulator. 9,10 Bi 2 Se 3 is a strong three dimensional(3D) topological insulator with its surface states consisting of single Dirac cone at the Γ point which is protected by the time-reversal symmetry from any time-reversal perturbation, such as crystal defects and nonmagnetic impurities. 9,11,12 Particularly, its bulk band gap is much larger than the energy scale of room temperature, making it the most suitable candidate for the high-temperature spintronics application. 9,11 Whether as thermoelectric material or topological insulators, the native defects in Bi 2 Se 3 play important role in the material. As thermoelectric material, it needs well-defined electrical and thermal conductivities, high mobility of free current carriers, and thermoelectric power, all of them will be influenced by the presence of native defects. [13][14][15] As a topological insulator, the intrinsic defects such as antisites or vacancies behave as n type dopant and consequently shift the Fermi level above the Dirac point, which makes it difficult to characterize the topological transport properties and to realize topological devices because both of them severely rely on the behavior of surface Dirac fermions. 12,16,17 Hence, it is vital to investigate the nature of native defects in Bi 2 Se 3 . In experiments, lots of efforts 13,18 have been pursued on such issue. Urazhdin et al. 13 have found that Bi 2 Se 3 single-crystal behaves as n type material under excess Bi or Se growth condition. They predicted that Bi Se antisites contribute shallow acceptors when the samples are doped with excess Bi. However, the solubility of Bi in Bi 2 Se 3 is very low and the n type doping behavior of the samples is derived from the compensating defects. While when the samples doped with excess Se, Se Bi antisites introduce shallow donor in the materials. 13 Bludská et al. 18 concluded that the presence of overstoichiometric Bi atoms can induce the following defects: (i) Se vacancy carrying two positive charges (V +2 Se ), which is attributed to the main reason for n-type Bi 2 Se 3 ; (ii) Bi Se antisite carrying one negative charge (Bi −1 Se ); (iii) Bi vacancy supposed to carry three negative charges (V −3 Bi ); (iv) other complex defects. In theoretical calculations, Shuang-Xi Wang et al. 19 have studied the formation energies of native point defects in Bi 2 Se 3 , but the study only gave the results under two extreme growing conditions: Se-rich and Bi-rich. The systematic study of the evolution of the formation energies according to growth conditions has not been reported as yet.
To this end, in our present work, the stability and formation energies of native defects in Bi 2 Se 3 are systematically investigated using first-principles method based on density function theory. The results indicate that Se1 vacancy is a double donor, Bi vacancy is a triple acceptor, and Se Bi is an always active n-type dopant in Bi 2 Se 3 . Interestingly, we find that Bi Se1 is an amphoteric dopant. By analyzing the formation energies of all defects under different growth conditions (such as Bi-rich or Se-rich), we find that under Se-rich growth condition, Se Bi is the most stable one which is independent of electron chemical potential µ e . However,  under Bi-rich growth condition, Se vacancy is the most stable defect, except for a small growth window which is defined by µ e >0.262eV (the material is typical n-type) and ∆µ Se <-0.459eV(Bi-rich). In that growth window, one negative charged Bi Se1 is the most stable native defect.

II. COMPUTATIONAL METHOD AND PROCEDURES
Bi 2 Se 3 has the rhombohedral crystal structure with space group D 5 3d (R3m). Its primitive unit cell contains three Se atoms and two Bi atoms, and Se atoms have two nonequivalent sites Se1 and Se2, as shown in Fig.1(a). The conventional hexagonal unit cell has a layered structure with three basic unit cells (one unit cell named a quintuple layer(QL)) weakly bound to each other by the van der Waals forces. 20,21 In each QL, five atomic planes with atomic order Se1-Bi-Se2-Bi-Se1, see as fig.1(b), and the chemical bond between Bi and Se atoms is of the covalent-ionic type. 22 Our calculations are performed in a 3×3×1 supercell with 135 atoms by using the plane-wave pseudopotential code Vienna ab initio simulation package (VASP). 23,24 The core-electron interaction is modeled by the projector augmented wave method 25,26 and the generalized gradient approximation of Perdew, Burke, and Ernzerhof, exchange-correlation functional 27 is adopted. The plane wave cutoff energy is chosen as 240eV, and the Brillouin zone is sampled by using 5×5×1 Gammacentered Mon-khorst-Pack grids. All structure parameters of Bi 2 Se 3 are chosen from experimental data as listed in Table I. The energy convergent criterion is 10 −5 eV per unit cell. We relax all atoms without spin-orbit coupling (SOC) because the SOC only slightly affect the relaxation results and is too much time consuming. The forces on all relaxed atoms are less than 0.01eV/Å. The selfconsistent calculations are performed with and without spin-orbit coupling (SOC) to make a comparation. We find that the formation energies with SOC are lower than that without SOC and the difference of transition energy levels is larger than 0.04eV. In this paper, we emphasize on the results with SOC. Our calculation gives a band gap of 0.31eV, which is in good agreement with experimental and other theoretical reports. 2,11,31,32 We consider three types of vacancy point defect in Bi 2 Se 3 , namely vacancy on the Se1, Se2 and Bi sublattices, denoted as V Se1 , V Se2 , and V Bi . We sign their detail locations in Fig.2 (a). As for antisite defects, we consider Se Bi (substitute one Bi atom by Se atom) and Bi Se1 (substitute one Se1 atom by Bi atom) because according to the formation energies discussed below, the formation energy of V Se2 is 0.4eV higher than that of V Se1 , as shown in Fig.2 It is well known that the calculations based on GGA typically underestimate the band gap of a semiconductor and cannot give the absolute position of the defect level. In this paper, we emphasize the analysis of the defect formation energy and the transition energy levels. 33,34 The formation energies of all defects in Bi 2 Se 3 are given by: where H is the formation energy of a defect in and E perf ect totle , are the total energy of the super-cell with defect under charge state q and pure Bi 2 Se 3 in the same super-cell size, respectively. ∆n i is the change in the number of atoms of species i, which represents Se or Bi atom in the present study. µ i is the chemical potential of element i. E def ect V BM is the valence band maximum (VBM) of the super-cell with defect. µ e is the electron chemical potential. The defect transition energy level ε(q/q ′ ) is the electron chemical potential in Eq. (1), at which the formation energy of the defect with charge q is equal to that of another charge q ′ .
To describe formation energy more accurately, we introduce two corrections as follows: i) The total energies of the charged systems should be corrected for the interaction of the charged defect with the compensating background and its periodic images. We use Makov-Payne (M-P) corrections, 35 formulated as q 2 α/2ǫL, where L is the linear dimension of the supercell, ǫ is the static dielectric constant and α is Madelung constant; ii) Because the defects induce significant distortion of the band structure and fluctuation of the bandgap, the E V BM obtained from the first-principles method directly can not be used in Eq.(1). We adopt the correction introduced by previous works. [36][37][38][39] Firstly, we assume that the potentials in the perfect super-cell are similar to those far from a defect in a defective super-cell. Then, the average potential of the plane farthest from the defect in the defective system (V def ect aν ) and the average potential of the corresponding plane in the perfect system (V perf ect aν ) are determined. The difference of the average potentials between the per-fect and defective super-cells is used to determine E V BM of the defective super-cell as follows: The first term on the right-hand side of Eq. (2) is the VBM of perfect super-cell and can be obtained by , where E T (perfect : q) is the total energy of a perfect super-cell under charge state q.

A. The transition energy levels
We investigate the transition energy levels of all defects quantitatively by calculating their formation energies in function of the electron chemical potential µ e , and the chemical potentials µ i are considered to be the total energies of ideal elemental bulk Bi and Se crystal simply(such choice does not influence the results of the transition energy level of native defect), as shown in Fig.3(a). For Se vacancy, we consider both V Se1 and V Se2 , and find that the formation energy of V Se2 is 0.4eV higher than that of V Se1 . The lower formation energy of V Se1 can be explained as a result of weak interaction between QLs. Our results indicate that the stable charge states of V Se1 in Bi 2 Se 3 are V +2 Se1 , V +1

Se1
and V 0 Se1 as the Fermi level moves through the bandgap, which shows that V Se1 acts as a double donor in the systems. The results is in good agreement with the previous report. 18 The donor levels ε(+2/+1)=0.209eV, ε(+1/0)=0.231eV, are 0.101eV and 0.079eV below the conduction band minimum(CBM) respectively. As for Se Bi in Bi 2 Se 3 , it is one positive charged as µ e moves through the bandgap energy window. Consequently, Se Bi in Bi 2 Se 3 behaves as a donor, and the donor level ε(+1/0) enter into the conduction band, indicating that it is an active n-type doping in the system. The anti-site Se protrudes outwards about 0.14Å to the Se1 layer and induce the distance between the Se atoms in the Se1 layer decrease about 10.56%. The stable charge states of isolated Bi vacancy (V Bi ) in Bi 2 Se 3 as the Fermi level moves through the bandgap are V − Bi , V −2 Bi and V −3 Bi , as shown in Fig.3(a), indicating its triple acceptors characteristics. The acceptor levels are ε(0/-1)=-0.07eV (the level is active in the system because it is below the valence band maximum (VBM) ), ε(-1/-2)=0.008eV (shallow acceptor level), and ε(-2/-3)=0.129eV (deep acceptor level), respectively. The results agree with the previous report. 18 The last native defect considered in our present work is Bi Se1 . The results of formation energy of Bi Se1 as shown in Fig.3(a) indicates that the stable charge states of isolated Bi Se1 in Bi 2 Se 3 as the Fermi level moves through the bandgap are Bi + Se1 , Bi 0 Se1 and Bi − Se1 . Its transition energy levels are ε(+1/0)=0.119eV and ε(0/-1)=0.251eV. This means that the Bi Se1 in Bi 2 Se 3 is an amphoteric dopant, acting as a deep donor when the electron chemical potential is close to VBM (the material is typical p-type) and as a deep acceptor when the electron chemical potential is close to CBM (the material is typical n-type). Bi Se1 will behave as a compensation dopant to the electrical type of the material. The result agrees well with the band structure of the Bi 2 Se 3 (111) surface with defect Bi Se1 previously reported 19 with no apparent donor or acceptor feature. Moreover, in experiment, the n-type Bi Se1 was not mentioned and we will discuss the reasons below. As for the geometric configuration around Bi Se1 , we find that the anti-site Bi atom protrudes outward about 0.419Å, 0.425Å and 0.439Å relative to the average height of Se1 layer under neutral, one negative and one positive charge state, respectively. Such results are in good agreement with the experimental observation 13 that there are protrusion point on the surface when Bi 2 Se 3 doped with excess Bi. The protrusion can be attributed to the anti-site Bi on the surface denoted as Se 1 layer in the present work.

B. Stability of the native defect
We analyze the stability of all the native defects by calculating their formation energies under different growth environments (such as Bi-rich or Se-rich) using the method proposed by Hashibon et al. 40 Using ∆µ Bi (∆µ Se ), the difference of the chemical potential of Bi(Se) between in the bulk Bi 2 Se 3 and in ideal elemental bulk Bi(Se), the method can accurately simulate the relationship between growth condition and the formation energy of native defect. The detail process of this method is as follows: Firstly, the chemical potential of Bi and Se in bulk where µ bulk Bi2Se3 is the chemical potential of one formula unit of Bi 2 Se 3 . The chemical potentials of Bi and Se in ideal elemental bulk are denoted as µ 0 Bi and µ 0 Se , respectively. Then the the difference of the chemical potential of Bi(Se) is defined as: ∆µ Bi =µ Bi -µ 0 Bi (∆µ Se =µ Se -µ 0 Se ). The Gibbs free energy in the process of formation of one formula unit of Bi 2 Se 3 is given as: From Eq. 1, 3 and 4, the formation energy formula of point defects in Bi 2 Se 3 as a function of ∆µ Se can be described as: Similarly, the formation energy formula as a function of ∆µ Bi also can be obtained as above. Eq. 5 and 6 indicate that the formation energy is a function of µ e and ∆µ Se . In addition, Eq.(4) tells us that to form Bi 2 Se 3 , ∆µ Bi and ∆µ Se should meet the following condition: The formation energies under different grown environments are shown in Fig. 3(b) and the typical results for µ e =0eV, µ e =0.31eV under ∆µ Se =0eV (Se-rich) and ∆µ Bi =0eV (Bi-rich) conditions are listed in Tab. II. The results as shown in Fig. 3(b) and (c) indicate that under Se-rich growth condition, the growth window denote as I in Fig. 3(c), Se Bi is the most favorable native defect in the system within the full energy window of the electron chemical potential µ e , it behaves as a single donor. This means that under Se-rich growth condition Se Bi is the predominant native defect in the system and responsible for the n-type characteristics of the system under such growth condition. For example, when ∆µ Se =0eV and µ e =0 eV, the formation energy of Se Bi is 0.175eV, 1.338eV and 1.244eV lower than that of V Se1 , V Bi , and Bi Se1 , respectively. When µ Se =0eV and µ e =0.31eV, the formation energy of Se Bi is 0.305eV, 0.235eV, and 0.994eV lower than that of V Se1 , V Bi , and Bi Se1 , respectively. When the growth window is in the region II as shown in Fig. 3(c), the most favorable native defect is V Se1 , whereas in the growth window III which is confined by µ e >0.262eV (the material is typical n-type) and ∆µ Se <-0.459eV(Bi-rich) one negative charged Bi Se1 is the most favorable one. This means that under Bi-rich growth condition, V Se1 is the predominant native defect and it is the main reason for n-type electrical characteristics of the systems. For example, when ∆µ Bi =0eV (∆µ Se =-0.491eV), the formation energy of V Se1 according to µ e =0 eV is 0.213eV, which is 2.391eV, 1.543eV and 0.333eV lower than that of V Bi , Se Bi and Bi Se1 , respectively. Only when ∆µ Se <-0.459eV (the material is typical n-type), p-type Bi Se1 is the dominant native defect and it behaves as a compensation dopant for the the material. The formation energy of Bi Se1 according to µ e =0.31eV is 0.605eV, which is 0.048eV, 1.206eV and 1.461eV lower than that of V Se1 , V Bi and Se Bi , respectively.
In experiments, the electronic properties of bulk Bi 2 Se 3 crystals are usually dominated by electron donors, resulting in n-type conductivity. Based on above analysis, such n-type conductivity mainly derives from two reasons: First, the formation energy of p-type native defect V Bi is higher than that of n-type native defect, such as Se Bi and V Se1 . Second, one negative charged p-type native defect Bi Se1 is only stable when the material is strong n-type and Bi-rich growth condition (∆µ Se <-0.459eV) behaving as a compensation. Moreover, we also find that one positive charged Bi Se1 only forms when material is p-type and its formation energy is higher than that of V Se1 which is difficult to be found in experiment.

IV. ELECTRONIC STRUCTURES
The density of states (DOS) of perfect Bi 2 Se 3 system and three systems with native defect of V Se1 , Se Bi or Bi Se1 are shown in Fig. 4. In comparison with that of perfect system, it is clearly seen that each defect induces two defect states into the band gap of the material. For V Se1 , the defect states just below the Fermi level mainly come from p states of three nearest neighbor (NN) Bi atoms of V Se1 and are completely filled by two electrons. The electrons filled the defect states are easily released leading the defect a double donor. The defect states of Se Bi mainly come from p states of four Se atoms(the Se atom substituted the Bi and its three NN Se atoms), which are half filled. The electron filled the defect state is easily released leading the defect a single donor. The defect states of Bi Se1 mainly come from p states of four Bi atoms (the Bi atom substituted Se1 and its three NN Bi atoms), which are half filled. Depending on the conductive type of the material, either the electron filled the defect states is released leading the defect a donor or the unfilled state capture one electron leading the defect an acceptor. Bi Se1 behaves as an amphoteric dopant in the system.
Comparing the DOS and the formation energy results, we find that the defect levels in DOS are lower than the results of formation energy. The phenomenon may come from two possible reasons as follows: (i) the different relaxed geometries between neutral and charged system. Taking V Se for example, the distance between Bi atoms (nearest to the Se1 layer with vacancy) changes obviously from 3.834Å for q=0 to 4.370Å for q=+2; (ii)the correction of VBM for each charged system tend to make both acceptor and donor level higher. Moreover, we also note that although both Se Bi and Bi Se1 induce half filled defect states around the Fermi level, Se Bi acts as a single donor whereas Bi Se1 behaves as amphoteric dopant. The reason is that the un-filled defect level above the Fermi level of Se Bi is higher than that of the Bi Se1 , which make it difficult to capture electron acting as p-type dopant.

V. CONCLUSION
In present work, we find that V Se in Bi 2 Se 3 is a double donor; V Bi is a triple acceptor; Se Bi is an active donor, whereas Bi Se1 is an amphoteric dopant acting as a donor when µ e <0.119eV(the material is typical p-type) and as an acceptor when µ e >0.251eV(the material is typical n-type). The analysis of the stabilities of all defects show that under Se-rich growth condition, Se Bi is most stable according to all electron chemical potential µ e . Under Bi-rich growth condition, V Se is the main defect. One negative charged Bi Se1 is most stable in a small growth window which defined by µ e >0.262eV (the material is typical n-type) and ∆µ Se <-0.459eV(Bi-rich).