Band bending at the surface of Bi$_2$Se$_3$ studied from first principles

The band bending (BB) effect on the surface of the second-generation topological insulators implies a serious challenge to design transport devices. The BB is triggered by the effective electric field generated by charged impurities close to the surface and by the inhomogeneous charge distribution of the occupied surface states. Our self-consistent calculations in the Korringa-Kohn-Rostoker framework showed that in contrast to the bulk bands, the spectrum of the surface states is not bent at the surface. In turn, it is possible to tune the energy level of the Dirac point via the deposited surface dopants. In addition, the electrostatic modifications induced by the charged impurities on the surface induce long range oscillations in the charge density. For dopants located beneath the surface, however, these oscillations become highly suppressed. Our findings are in good agreement with recent experiments, however, our results indicate that the concentration of the surface doping cannot be estimated from the energy shift of the Dirac cone within the scope of the effective continuous model for the protected surface states.


I. INTRODUCTION
The theoretical discovery of the second-generation topological insulators 1 (2GTIs) triggered an intensive experimental effort to observe the predicted surface states 2-10 (SSs) being protected by time-reversal symmetry. 11 It turned out, that the physical properties of the prepared samples are greatly affected by the electron acceptor/donor impurities, that can be found either in the bulk or on the surface. The created Bi 2 Se 3 samples are typically electron doped by the inner point defects. 12 However, it has been shown that the Fermi level can be tuned by the insertion of further bulk dopants into the system. 4,6 The presence of the charged dopants in the system and the inhomogeneous charge distribution close to the surface generates an effective electrostatic field, which can be probed experimentally, for example by second-harmonic generation. 13,14 In addition, the evolved electrostatic field induces a band bending (BB) in the bulk band structure close to the surface, which was successfully observed by recent angle-resolved photomeisson (ARPES) experiments as well. [2][3][4]15 The experimental manipulation of the BB field was also accomplished by the insertion of bulk Cu dopants into the Bi 2 Se 3 matrix. 9 Moreover, ARPES experiments also demonstrated the possibility to shift the Dirac cone by applying charge dopants on the Bi 2 Se 3 surface. 6 An experimental evidence for a large shift of the Dirac cone towards the conduction band was also reported by gated terahertz cyclotron resonance measurements performed on thin Bi 2 Se 3 film. 16 Besides these comprehensive experimental studies, numerous theoretical works were also devoted to the description of the physical properties of the 2GTIs, including first principle calculations, 1,17-19 tight binding 21,22 or effective continuous 11 models. First principle calculations revealed that the structure of the 2GTIs can be described by a sequence of weakly bound quintuple layers (QLs), each consisting of five atomic layers. The ef-fect of the bulk dopants on the electronic behavior of Bi 2 Se 3 crystal was also addressed in recent theoretical studies. 22,23 In addition, Ref. 22 also suggested that the BB profile and the energy of the Dirac point can be controlled by electrostatic effects. The importance of the BB has been demonstrated in other 2GTI based nanostructures as well, including topological/normal insulator interfaces. [24][25][26][27] Motivated by these research findings in this work we study the BB effect and the Dirac cone shift on the surface of 2GTIs theoretically. In particular, we perform screened Korringa-Kohn-Rostoker (SKKR) calculations to examine the role of the charged dopants at Bi 2 Se 3 surface. The slow dynamics of the band bending process suggests that the charge accumulation at the surface is coupled to a much slower surface lattice relaxation. 2 In this work we do not aim to describe the time dependence of the outlined process, but rather to examine the effect of the surface dopants on the band structure. Thus, in our surface calculations we consider a rigid lattice excluding any structural relaxation processes of the surface layers.
The rest of the paper is organized as follows. In Sec. II we present the details of our numerical approach to study the bulk and surface properties of Bi 2 Se 3 . Than in Sec. III we summarize our results organized in three subsections: (i) first we calculate the band structure of a bulk Bi 2 Se 3 system, (ii) then we present our results obtained for the surface band structure, (iii) in the third subsection we examine the BB and the effect of the surface dopants on the dispersion of the SSs. Finally we summarize our work in Sec. IV.

II. DETAILS OF THE NUMERICAL CALCULATIONS
In order to describe the electron structure of Bi 2 Se 3 , we used the relativistic spin-polarized SKKR method. 28 The first-principles calculations were performed by density functional theory using the local spin-density approximation and the Ceperley-Alder parametrization of the exchange correlation function 29 within the the atomicsphere approximation (ASA). In our calculations we used an angular momentum cutoff l max = 2.
The 2GTIs posses a rhombohedral lattice structure were the atoms are located in parallel layers forming a triangular lattice perpendicularly to the crystallographic (111) axis. This lattice structure can be described by a periodic sequence of QLs, whereas each of the QL consists of five strongly bound atomic layers. The QLs, on the other hand, are weakly bound to each other by van der Waals couplings. In particular, the QLs in the Bi 2 Se 3 crystal consist of atomic layers Se1-Bi-Se2-Bi-Se1, whereas Se1 and Se2 are selenium atoms at inequivalent geometrical positions. Fig. 1(a) shows the structure of one QL in the lattice. Since the QLs are weakly bound to each other, the crystal surface in the (111) crystallographic direction is favored to be terminated by Se atoms. Thus, in our calculations we assumed a flat surface formed by the last Se atomic plane of a QL. The (111) surface of Bi 2 Se 3 has a hexagonal structure, with lattice constant a = 4.138Å and lattice vectors given by a 1 = (a, 0, 0) and a 2 = (− 1 2 a, √ 3 2 a, 0), where the axis z is parallel to the direction (111). The shortest period in the lattice structure along the z direction is given by three successive QLs 1 , since the nearest neighbor QLs are shifted against each other by a lattice vector a 3 = (5a, 5 √ 3 a, z QL ), where z QL = 9.547Å is the height of one QL.
According to the SKKR method 28 , we modeled the bulk system by a single QL surrounded by semi-infinite bulk regions from both the left and right sides. The surface Greens functions of the semi-infinite regions were calculated by the iterative method discussed in Ref. 30 and 31. The charge density distribution was calculated by means of the energy-dependent Greens function, integrated up to the Fermi energy. The position of the Fermi energy was determined to satisfy the total charge neutrality condition. In particular, because of the ASA used in our calculations, however, the value of the Fermi level is obtained by a numerical error. Though, the correct position of the Fermi energy is essential for insulator systems. Refs. 32-34 proposed a procedure to correct the value of the Fermi energy in a self-consistent way based on the Lloyd's formula. 35 Following this procedure we obtained the proper Fermi energy by re-normalizing the wave functions in order to obtain the correct spaceintegrated charge distribution.
For the surface calculations we considered an interface region surrounded by a semi-infinite bulk system from the left and by a vacuum from the right [see Fig. 1 The interface region was constructed from six QLs to support a smooth transition of the atomic potentials as we proceed from the bulk layers to the vacuum. We also included additional (in total eight) vacuum layers between the surface of the interface region and the semi-infinite vacuum side. From an electrostatic point of view, the atomic potentials (as well as the charge density distribution) were obtained by imposing zero electrostatic field far away from the surface of the crystal.
The described numerical method is sufficient to obtain plausible results for the studied system, 19 however, the inclusion of empty spheres between the atomic layers further stabilized our numerical approach. In our calcu- lations we used identical geometrical parameters for the lattice structure as in Ref. 20, however, we optimized the positions and the Wigner-Seits radii of the empty spheres to reproduce the main features of the experimentally observed band structure of the SSs. 4 Focusing on the direct band gap at the Γ point and on the slope of the Dirac cone, the obtained numerical parameters are summarized in Table I. Finally, the band structure can be obtained from the Bloch spectral function (BSF). For surface calculations the layer-resolved BSF depends on the energy E and on the parallel momentum k : where G + (r n , r n , E, k ) is the retarded Greens function in the nth atomic sphere at position r n . The integral is taken over the volume of the nth atomic sphere Ω n , and the trace is taken over the quantum numbers of the total angular momentum. Thus, the BSF is ideal to study the surface states, whereas the three-dimensional bulk band structure is projected onto the two-dimensional Brillouin zone (BZ) corresponding to the (111) surface [see Fig.1(b)].

III. RESULTS
A. The bulk band structure of Bi2Se3 compound First we determine the self-consistent atomic potentials of the bulk Bi 2 Se 3 crystal. Then we calculated the layer resolved BSF given by Eq. (1) for the bulk system at complex energies with a small imaginary part, and summed the contributions of the individual atomic and empty layers within one QL. The result of this procedure Figure 2. Projection of the BSF onto the two-dimensional BZ summed over the layers of one QL inside the bulk plotted along the KΓM cross section of the two-dimensional BZ. The blue and red colors represents areas of low and high density of states, respectively. The Fermi level is at the top of the valence band. The BSF was calculated at complex energies with small imaginary part of ∼ 0.7 meV.
can be interpreted as the projection of the bulk band structure onto the two-dimensional BZ corresponding to the surface perpendicular to the (111) crystallographic direction. The calculated BSF along the KΓM cross section of the BZ is shown in Fig.2. The high density areas (forming narrow lines close to the extremal points of the band structure) indicate electron states of low dispersion along the (111) direction. The Fermi level is at the top of the valence band.
It should be noted that the Fermi level is sensitive to the presence of impurities, even at small concentration. Several ARPES measurements indicated that due to the bulk dopants the Fermi energy was tuned into the bulk conduction band. 4 According to our calculations, the direct bulk band gap at the Γ point is ∼ 447 meV. Our results are in good agreement with other theoretical calculations. 1,36,37 B. The surface band structure of Bi2Se3 We now turn our attention to the SSs formed on the surface of Bi 2 Se 3 crystal. Making use of the atomic potentials determined by the self-consistent calculations described in Sec.II we calculated the layer-resolved BSF in the QLs beneath the surface.   3 shows the calculated BSF summed over the layers of the individual QLs. The narrow lines correspond to the bands with low dispersion in the (111) direction, including the SSs. Within the bulk band gap the dispersion of the protected SSs forms a Dirac cone, which is anisotropic in the k plane. The Dirac point is located ∼ 250 meV below the bulk conduction minimum, in good agreement with the ARPES measurements. 6 The SSs penetrates below the surface up to the third QL, where the intensity of the SSs eventually vanishes. The signatures of the bulk band structure, however, can be observed in all the QLs up to the topmost one. Thus, the SSs spatially overlap with the bulk states, hence they can hybridize. Consequently, one can observe an increased width of the SS bands in the vicinity of the bulk bands.
The SSs that are not hybridized with the bulk states, on the other hand, can be described by a 2 × 2 effective Hamiltonian proposed by Fu et al. 11 up to the third order of the momentum k = ( k x , k y ): Here v k = v 0 + αk 2 , k ± = k x ± ik y and σ x,y,z stand for the Pauli matrices acting in the spin space. Parameters m, v 0 , α and λ are to be determined either from first principles or from fits to the experimental data. The energy eigenvalues of the electron states are given by the expression where ± labels the bands above/below the Dirac point. The effective mass m introduces an asymmetry between the upper and lower side of the Dirac cone that is, indeed, significant according to the calculated band structure shown in Fig. 3. As the main objective of our work we discuss our results on the BB effect induced by the electric charge accumulation and/or inhomogeneous charge distribution close to the surface. We also show that the Dirac cone can be shifted in energy due to the effective electric field generated by the deposited surface dopants. Our findings are in good agreement with recent ARPES experiments. 6 Comparing the band structure shown in Fig. 3.(a) to the ones of Figs. 3.(b)-(d), one can observe the BB in the topmost QL even without any extra charge dopants added to the system. Due to the effective electrostatic field induced by the inhomogeneous charge distribution close to the surface, the bulk conduction band is repelled upward by about 100 meV in the first QL compared to the conduction band minimum in the other QLs (see Fig. 3). We expect that charge dopants deposited on the Bi 2 Se 3 surface further modifies the electrostatic configuration of the top layers. Indeed, in Ref. 6 the surface of the Bi 2 Se 3 sample was dosed by NO 2 molecules and accomplished to shift the Dirac point towards the conduction band. The presence of charged impurities in the sample, like electron donor Se vacancies generated during the sample preparation process, also plays a crucial role in the electrostatic properties of the surface. The slow migration of the Se vacancies 3,12 results in an increased concentration of positively charged impurities close to the surface. Thus, a time dependent BB was observed in ARPES experiments 2,3,6 , where the conduction band minimum at the surface was gradually bent under the Fermi level.
In our calculations we simulated the effect of the surface dopants by a planar capacitor situated on the surface of the Bi 2 Se 3 crystal, and charged by δq s per twodimensional unit cell. For example, the presence of the electron donor Se vacancies can be described by a positively charged capacitor. On the other hand, the presence of the electron acceptor NO 2 molecules used in the ARPES experiments of Ref.6 can be modeled by a negatively charged capacitor.
The electric field of the charged capacitor is accounted for by means of the Poisson equation, that is solved selfconsistently within the SKKR code. In the vacuum the electric field of the capacitor is canceled due to the attracted (or repelled) electrons from (to) the bulk reservoir. Fig. 4 shows the induced charge excess (∆Q) on the individual layers compared to the undoped system. The results were obtained for a capacitor located close to the surface, and being charged to δq = 0.01q e , where q e ≈ 1.6 × 10 −19 is the elementary charge unit. In particular, Fig. 4 shows the charge excess obtained for two different positions of the capacitor. As expected, the charge of the capacitor is compensated by the accumulated electrons within the first two QLs, since the induced charge excess can be hosted only by the unsaturated SSs. It can also be noticed that the induced effective electric field  might generate further charge transfer between the layers in the low-lying QLs as well. For example, in Fig. 4(a) one can clearly observe oscillations of ∆Q up to the 5th QLs. However, since the oscillations are centered around ∆Q = 0, the net charge of the corresponding QLs remains zero within the numerical precision of our calculations. In the following we present our result obtained for the planar capacitor located at position corresponding to Fig. 4(a). Recently the energy shift of the Dirac cone was successfully demonstrated by ARPES measurements as well. 6 As it is shown in Fig. 5, the surface bands are shifted upward or downward equally in each QL. This behavior differs from the standard BB mechanism of the bulk bands. As the bend of the bulk bands varies continuously with the distance measured from the surface. Indeed, in Fig. 3(a)-(d) one can observe a gradual decrease of the conduction band minimum in the successive QLs. The layer resolved BSF (not shown in the manuscript) also confirms that the energy shift of the conduction band minimum varies smoothly from layer-to-layer.
We now examine the relation between the energy shift of the Dirac cone and the concentration of the charge dopants. In case the Fermi energy is located inside the bulk band gap the charge excess induced by the surface doping can be hosted only by the unsaturated SSs, forming a Dirac cone at the center of the two-dimensional BZ. One can then expect, that the energy shift of the Dirac  cone is in correspondence with the increment or reduction of the occupied electron states required to host the induced charge excess. The number of the electron states on the Dirac cone per two-dimensional unit cell (Ω) between energies E 1 and E 2 can be calculated as where is the density of states and v g (k) = |grad k E ± (k)| is the group velocity. The integral path is taken over the constant energy contour Γ E of the SS spectrum (3) at energy E. Fig. 6 shows the calculated energy shift of the Dirac cone as a function of the deposited charge excess. Surprisingly, we have found a significant difference between the predictions of the effective model and the SKKR results. The effective model highly overestimates the energy shift of the Dirac cone compared to the results obtained by the SKKR model. One can explain this inconsistency by the following arguments. Within the effective model we assumed that the surface bands are shifted without any notable changes of the dispersion. However, this approximation is valid only at energies close to the Dirac point. Indeed, in Fig. 5 the shape of the Dirac cones are similar in the doped and undoped cases. One can, however, observe remarkable changes in the surface band structure around the crossing points with the other bands when surface dopants are present in the system (see Fig. 5). Consequently, the charge density distribution on these parts of the band structure becomes undergoes to a significant change as well. Thus, the low-energy segments of the unsaturated surface bands also play an important role in the screening of the surface dopants.
Our results indicate that the energy shift of the Dirac cone cannot be estimated within the effective model given by Hamiltonian (2), since a significant portion of the induced charge excess is hosted by the low-energy segments of the surface bands. Besides a planar capacitor we also considered a scenario of spatially distributed charge dopants below the surface, that might be closer to the realistic case. In such calculations we used a model where the charge dopants were concentrated into the first two QLs below the surface. However, we did not find any qualitative differences compared to the case of a planar capacitor.

IV. SUMMARY
In summary, we have calculated the bulk and surface band structure of Bi 2 Se 3 topological insulator by using the SKKR method. Our results are in good agreement with other theoretical and experimental results. In order to examine the effect of the charged dopants on the properties of the Bi 2 Se 3 surface, we also calculated the surface band structure in the presence of a charged planar capacitor situated close to the surface. We have found, that for a Fermi energy located in the bulk band gap the induced charge excess is hosted by the unsaturated SSs. Thus, the charge of the surface dopants is also screened within the first two QLs below the surface. In addition, due to the charge excess and the inhomogeneous charge distribution close to the surface, the bulk bands undergo a BB effect. Consequently, the bulk bands becomes bent in the atomic layers close to the surface, but one can already recover the properties of the bulk crystal starting from the third QL. Our results also indicate, that the BB profile can be tuned via the deposited surface dopants.
In contrast to the bulk bands, the Dirac cone (being formed by the SS bands inside the bulk band gap) becomes shifted in energy due to the deposited surface dopants without a BB. The magnitude and the direction of this energy shift depends on the concentration and on the sign of the deposited surface dopants. However, this effect cannot be described within the scope of the effective continuous model of the SSs. Our self-consistent numerical results showed that besides the Dirac cone the low-energy segments of the surface bands also play an important role in the electrostatic properties of the surface.
In order to check experimentally our findings, one need to independently measure the doping concentration on the surface and the energy shift of the Dirac cone. We believe that the combination of the scanning tunneling microscope and ARPES techniques can serve this purpose. Moreover, the surface doping of the 2GTIs can be used to cancel the BB effect on the surface, which is essential to take an advantage of the protected SSs in transport devices. Finally, the possibility of tuning the position of the Dirac cone might also be of great importance for future experimental applications of these materials.