The effect of mechanical strain on the Dirac surface states in the (0001) surface and the cohesive energy of the topological insulator Bi₂Se₃

Abstract

The band gap (E_g) engineering and Dirac point tuning of the (0001) surface of 8 QLs (quintuple layers) thick Bi₂Se₃ slab are explored using the first-principles density functional theory calculations by varying the strain. The strain on the Bi₂Se₃ slab primarily varies the bandwidth, modifies the p_z – orbital population of Bi and moves the Dirac point of the (0001) surface of Bi₂Se₃. The Dirac cone feature of the (0001) surface of Bi₂Se₃ is preserved for the entire range of the biaxial strain. However, around 5% tensile uniaxial strain and even lower value of volume conservation strain annihilate the Dirac cone, which causes the loss of topological (0001) surface states of Bi₂Se₃. The biaxial strain provides ease in achieving the Dirac cone at the Fermi energy (E_F) than the uniaxial and volume conservation strains. Interestingly, the transition from direct E_g to indirect E_g state of the (0001) surface of Bi₂Se₃ is observed in the volume conservation strain-dependent E_g. The strain on Bi₂Se₃, significantly modifies the conduction band of Se2 atoms near E_F compared to Bi and Se1, and plays a vital role in the conduction of the (0001) surface of Bi₂Se₃. The atomic cohesive energy of the Bi₂Se₃ slab is very close to that of (0001) oriented nanocrystals extracted from the Raman spectra. The strain-dependent cohesive energy indicates that at a higher value of strain, the uniaxial and volume conservation strain provides better stability than that of the biaxial strain (0001) oriented growth of the Bi₂Se₃ nanocrystals. Our study establishes the relationship between the strained lattice and electronic structures of Bi₂Se₃, and more generally demonstrates the tuning of the Dirac point with the mechanical strain.

---

1. Introduction

The concept of topological insulators (TIs) has changed the conventional way of classifying solids based on symmetry breaking since it possesses symmetry protected topologically active surface states. TIs have created considerable attention among the condensed matter physics community, which behave as an insulator in its bulk but possess topologically active metallic surface states. TIs are characterized by a number known as the Z₂ topological invariant, which can be calculated from the parity of the occupied bands for systems having inversion symmetry.¹ Among the different compounds showing topological properties, Bi₂Se₃ is considered as the prototype of three-dimensional (3D) TIs and has been extensively studied over the last few years due to some of its unique features such as a simple band structure with a direct band gap (E_g) of 0.3 eV, single Dirac cone at the Γ point² and the helical spin texture of opposite spin helicity for the upper and lower Dirac cones.³ The tailoring of the topological surface states has been investigated using different external and internal agents to make the exotic topological phenomena viable for potential applications.^4,5 Some of the widely studied agents are stress, electromagnetic field, chemical substitutions, stacking defects,⁶etc. In particular, strain plays a significant role in tuning the physical properties of TIs. For example, the first-principles density functional theory (DFT) calculations show the shift of the Dirac point of the helical surface states,⁷ enhancement or destruction of surface states,⁸ improvement in thermoelectric properties,⁹ decrease in the bulk carrier density,¹⁰ decrease in the coulombic gap, and increase in the strength of the spin–orbit interaction¹¹ in Bi₂Se₃. Using van der Waals DFT and semi-classical Boltzmann theory, the thermoelectric figure of merit of a single quintuple layer (QL) of Bi₂Se₃ is found to be 0.27, which is more than the bulk value 0.10, and can be further increased to 0.30 by introducing 2.5% compressive strain.¹² DFT calculations have also established that strain plays a crucial role to tune the ground state of the monolayer (ML) In₂Se₃. The strain fully suppresses the polarization in ML-In₂Se₃ and causes out-of-plane polarization in 2D non-polar materials, i.e., ML-Bi₂Se₃, ML-Bi₂Te₃, ML-Sb₂Te₃ and ML-Sb₂Se₃.¹³ Strain engineering also predicts the occurrence of new 3D TI compounds with anti-perovskite structures in ternary cubic centrosymmetric compounds (M₃N)Bi where M = Ca, Sr, Ba.¹⁴ Strain engineering plays a unique role in controlling the band structure and influences the Dirac point energy, Fermi velocity, metallic character, and the topology of the Bi based chalcogenide compounds. Ab initio calculations suggest the modification in the band structures, the absence of bulk-free carrier states, enhancement in the E_g opening in highly strained Bi₂Se₃ films that create the possibilities for a TI-based field-effect transistor.¹⁵

Inspired from the theoretical prediction of the possible tuning of bulk band gap and surface states of Bi₂Se₃ by the elastic strain, the strain has been applied on Bi₂Se₃ using different configurations, such as substrate-induced strain, stretching of the flexible substrate, doping, and the intercalation. Chae et al. have grown Bi₂Se₃ films on graphene/SiO₂ and SiO₂ surface.¹⁶ The 3 QLs thick film of Bi₂Se₃ grown on graphene/SiO₂ exhibits strain values of 5.4%, while the film on SiO₂ surface shows 7.6%. Interestingly, the strain on 2 QLs Bi₂Se₃ films is much lower than that of the 3 QLs film, implying that the interfacial strain from the graphene substrate is much more dominant as the thickness decreases within 3 QLs.¹⁶ However, the strain induced by the GaN substrate on Bi₂Se₃ is relatively low between 3% to 2.7%, although the lattice mismatch is 30%.¹⁷ The substrate-induced stress increases with thickness in the case of In₂Se₃ and decreases in the case of Bi₂Se₃ grown on sapphire because of the opposite lattice mismatch and different thermal expansion coefficients between the sample and the substrate.¹⁸ Flötotto et al. have grown ultrathin epitaxial Bi₂Se₃ films bonded onto conductive polyimide foils (Kapton).¹⁹ By using the strain holder, the strain has been applied to Kapton. This technique allows inducing around 2.1% in-plane strain on Bi₂Se₃. In addition, elemental doping in Bi₂Se₃ is a fundamental approach to introduce a mismatch in the effective ionic radii, i.e., strain. Though there are several studies on elemental doping in Bi₂Se₃, the measurement of the induced strain is very scarce. Qi et al. have measured the strain induced by Mn in Bi₂Se₃ after doping. The doping of Mn introduces 6% strain in Bi₂Se₃.²⁰ The other technique to introduce strain in Bi₂Se₃ is through intercalation of various species. The predominantly ionic nature of an intercalant requires either a change in the host lattice oxidation states or the presence of atomic vacancies to maintain charge neutrality, thus limiting the intercalant concentration.²¹ A zero-valent intercalant does not require a change in the oxidation state of the host lattice, thus allowing a high intercalation concentration.²² Using Rietveld refinement of the Bi₂Se₃ host lattice structure, Koski et al. have measured lattice parameters of Bi₂Se₃ and Cu-intercalated Bi₂Se₃.²¹ A 10% Cu intercalation in Bi₂Se₃ introduced strain of 9.66% in the plane and 0.35% along out-of-plane directions. Through an alternative measurement, electron diffraction confirmed, the introduction of 8.69% in-plane and 19.64% out-of-plane strain due to the 60% intercalation of Cu in Bi₂Se₃.²¹

Through low-temperature magneto-transport measurement, it is shown that the TI surface under the compressive strain of ±0.1% experiences a significant Dirac point shift (∓30 meV) as compared to the relaxed surface. The carrier mobility of TI is also increased for the surface under tensile strain.²³ Although numerous reports available in the literature, a systematic study of strain engineering on surface states of Bi₂Se₃ is still lacking. In the present work, we have performed first-principles DFT calculations to study the electronic structure evolution of (0001) surface states of Bi₂Se₃ under the influence of in-plane and out-of-plane strain. The (0001) surface band structures show a substantial variation in the topological properties with the direction of strain (anisotropic behaviour).

2. Computational details

Band structure calculations were performed for (0001) Bi₂Se₃ surface using the Projector Augmented Wave (PAW) pseudopotential and plane-wave basis set as implemented in the Quantum Espresso (QE) package.²⁴ Exchange-correlation potential was approximated through the Perdew–Burke–Ernzerhof general gradient approximation (PBE-GGA) functional. Further, the dispersion corrections were included through the semi-empirical Grimme-D2 van der Waals correction.²⁵ The kinetic energy cut-off to fix the number of plane waves was taken as 50 Ry with the charge density cut-off 200 Ry. The Brillouin zone integration was carried out using a 6 × 6 × 1 Monkhorst pack k-point grid for structural relaxation, and a denser k-mesh 8 × 8 × 1 was used for electronic structure calculations. The convergence criterion for self-consistent energy was taken to be 10⁻⁸ Ry and the atomic positions were optimized until the force on each atom was lower than 10⁻³ eV Å⁻¹. The surface calculations were performed using a slab model. An eight QLs slab of hexagonal Bi₂Se₃ with 40 atoms in the unit cell was generated using VESTA²⁶ and XCrySDen software. A vacuum of 15 Å was added to the top QL to avoid the interaction among the surfaces of neighbouring slabs.

3. Experimental methods

Bi₂Se₃ nanostructures were prepared at 250 °C by the chemical hot injection method. Commercial grade Bi₂O₃ (99.999%) and elemental Se (99.99%) powders were chosen as precursors, while 1-octadecane and oleic acid were used as solvent and reducing agent, respectively. To prepare the Bi precursor, 1 mmol of Bi₂O₃ was mixed with 5 ml of oleic acid and 20 ml of octadecane in a round bottom flask, mixed vigorously using a magnetic stirrer for 60 min, then heated to 100 °C and were kept at that temperature for 30 min for degassing. Then, the Bi precursor was prepared at 250 °C for 360 min under a nitrogen flow. In a typical synthesis, the elemental Se powder was mixed in 15 ml of octadecane in another round bottom flask, heated to 250 °C, after degassing at 100 °C, then the Bi-precursor solution was injected into the flask. Finally, the solution was cooled to room temperature, washed with ethanol and hexane (1 : 3 ratio) several times, centrifuged, and finally, the product was dried at 60 °C.

The phase of these nanoplates was confirmed using a Rigaku Smart lab X-ray diffractometer with Cu-K_α radiation (λ = 1.5405 Å) (see Fig. S2†). The morphological studies were carried out by using a high-resolution transmission electron microscope (HRTEM) (FEI Tecnai-G² T20) with an operating voltage of 200 kV (see Fig. S3†). The Raman spectra were recorded on a Jobin-Yvon LabRAM HR800UV spectrometer instrument equipped with a highly efficient thermo-electrically cooled charge-coupled device (CCD). The spectra were taken at different temperatures in the backscattering configuration using a 632 nm emission line of a He–Ne laser with laser power of 65 μW on the sample surface.

4. Result and discussion

Bulk Bi₂Se₃ belongs to the tetradymite type crystal with a hexagonal structure of the space group R3̄m (D⁵_3d). Bi₂Se₃ possess a layered structure along the hexagonal c-axis [Fig. S1a†]. The conventional unit cell of the hexagonal Bi₂Se₃ contains 15 atoms, with five atomic layers arranged in a particular order called ‘Quintuple Layer’ (QL) along the c-direction. There are three quintuple layers in a conventional hexagonal unit cell. Each QL consists of five atomic layers in the order Se1–Bi–Se2–Bi–Se1 [Fig. 1a]. The top and bottom Se1 atoms have the same Wyckoff positions and are in a similar interaction environment. The Se2 layer acts as the inversion center with two equivalent Bi and Se1 atoms [Fig. S1b†]. Due to the presence of inversion symmetry, the topological invariant can be calculated through the parity of occupied bands at ‘Γ’ point.² Within a QL, the atomic interaction forms stronger covalent bonds, whereas the inter-quintuple layer interaction is weak van der Waals type. The thickness of each QL is about 0.96 nm.²⁷

Fig. 1 (a) Schematic of a single quintuple layer of Bi₂Se₃ indicating the order of the atomic arrangement. (b) The total density of states of the (0001) surface of Bi₂Se₃ without spin–orbit coupling under ε_a strain. Projected density of states of (c) Bi, (d) Se2 and (e) Se1 atoms on the (0001) surface of Bi₂Se₃ without spin–orbit coupling under tensile ε_c, ε_a and ε_c = −ε_a strain.

The Slab model has been used for the surface states calculations of (112̄1) Bi₂Se₃.²⁸ The calculations show that the energy gap at the Dirac point is closed for slab thickness larger than 6 QLs of Bi₂Se₃. However, the (0001) surface E_g of the Bi₂Se₃ is vanishingly small but finite for 6 QLs thick slab and becomes ideally zero with slab thickness above 6 QLs.²⁹ Therefore, we have chosen an 8 QLs thick slab with 40 atoms to minimize the computational error and acquire reliable results. The unstrained 8 QLs slab constructed using the lattice parameters of the hexagonal-shaped nanocrystals of Bi₂Se₃ prepared by adopting a chemical hot-injection method is used.²⁹ The crystal structures of these nanocrystals are established from the Rietveld refinement analysis, which confirmed the R3̄m (D⁵_3d) crystallographic group of Bi₂Se₃. The refinement of the X-ray diffraction profile with a high degree of precision reveals the cell parameters; a = 4.136 Å and c = 28.59 Å of Bi₂Se₃ [Fig. S2†].²⁹ The hexagonal shape morphology and the (0001) orientation of nanocrystals are further confirmed by the transmission electron microscopy measurement [Fig. S3†]. Thus, the reference lattice parameters of 8 QLs hexagonal Bi₂Se₃ slab with 15 Å vacuum layer are a_ref = 4.136 Å and c_ref = 100.9 Å. These lattice parameters are varied individually or simultaneously to apply strain on the slab and study the band structure of the strained Bi₂Se₃ (0001) surface. The out-of-plane lattice parameter ‘c’ is varied to apply uniaxial strain, which can be defined as . The in-plane lattice parameter ‘a’ is changed for the application of biaxial strain . While the ε_a = −ε_c strain with the conservation of volume is generated by varying both ‘a’ and ‘c’, simultaneously. The strain on the Bi₂Se₃ slab primarily varies the band width, which is reflected in the density of states (DOS).

Before the studies of the (0001) surface band structure of Bi₂Se₃, the bulk band structure calculations are carried out after relaxing the atomic coordinates to the minimum energy. The bulk band structure shows a direct E_g ≈ 0.19 eV [Fig. S4a†], which is very close to the reported experimental^30–32 and theoretical^33,34 values. The kinetic energy in the Hamiltonian without spin–orbit coupling term does not depend on the spin degrees of freedom. However, in real crystals, the spin and orbital motions of the electrons couple together. The incorporation of SOC in the band structure calculation of Bi₂Se₃ removes the spin degeneracy with the metallic edge or surface states are consistent with the high-momentum-resolution angle-resolved photoemission spectroscopy experiment.³⁵ Interestingly, after including the SOC, the band structure shows an expansion of the E_g = 0.38 eV [Fig. S4b†], the bulk E_g consistent with previous reports.^2,11,34,36 The representative total DOS of the (0001) surface of Bi₂Se₃ with ε_a strain are shown in Fig. 1b. The compressive ε_a strain expands the band width, while tensile ε_a strain compresses the band width. The band width of Bi and Se represents the orbital occupancy, and the modification in the band width indicates the change in the electron density, orbital overlap and E_g. The band width of the surface QL atoms due to the tensile ε_a = −ε_c strain expands, while that of the tensile ε_a compresses, as compared to the tensile ε_c strain [Fig. 1c–e]. The Bi atoms in the surface QL remain insulating irrespective of the nature or direction of the strain applied to the slab. However, Se atoms are very sensitive to the nature or directions of the strain applied to the slab. The Se atoms in the surface QL present away from the E_F for tensile ε_c and ε_a = −ε_c strain as compared to that of the tensile ε_a strain [Fig. 1c–e]. The conduction band (CB) band width of Se2 shrinks as well as shifts away from the E_F as compared to that of Bi and Se1, irrespective of the nature or direction of the applied strain. Thus, Se2 atoms play an essential role in the Bi₂Se₃ (0001) surface conduction.³⁶

The Bi atoms in Bi₂Se₃ contribute to the CB, while Se atoms contribute to the valence band (VB) near the E_F. The p_x and p_y orbitals of both Bi and Se atoms are degenerate. The p_x and p_y orbitals of Bi atoms do not contribute at the Γ point and spread over the Γ–M–K path, while p_z orbital of Bi atoms disperse at and around the Γ point for the compressive ε_c strain on the Bi₂Se₃ slab [Fig. 2a and b]. As the compressive ε_c strain on the (0001) surface of Bi₂Se₃ decreases the p_x and p_y orbitals of Bi atoms move toward the E_F along the Γ–M path and the p_z orbital of Bi atoms move away from E_F [Fig. 2g and h]. Similar to Bi atoms, the p_x and p_y orbitals of Se atoms do not contribute at the Γ point, but populate over the Γ–K and Γ–M–K paths and p_z orbital populates around the Γ point due to compressive ε_c strain on the Bi₂Se₃ slab [Fig. 2c–f]. On decreasing the ε_c compression on the (0001) surface of Bi₂Se₃, the population near the Γ point due to the p_x, p_y and p_z orbitals of Se atoms increase, and interestingly the p_z orbitals of Se1 atoms spread over the Γ–M–K path [Fig. 2i–l]. It is interesting to note that the orbital population of Bi and Se after the compressive ε_a strain on the Bi₂Se₃ slab is qualitatively similar to that of the tensile ε_c and the orbitals dispersion of Bi and Se due to the tensile ε_a strain is analogous to that of the compressive ε_c strain. The orbital population because of ε_a = −ε_c strain on the Bi₂Se₃ slab is similar to that of the ε_c strain.

Fig. 2 Projected band structures of the p_x and p_z orbitals of the Bi, Se2 and Se1 atoms on the (0001) surface of Bi₂Se₃ under uniaxial strain with −6% (a–f) and 6% (g–l).

The calculated E_g of the (0001) surface of unstrained Bi₂Se₃ slab without SOC (WSOC), i.e., a = 4.136 Å, and c = 100.9 Å is ∼0.3 eV [Fig. 3b], which is closed to the reported value 0.26 eV for six QLs slabs.^12,37 The corresponding band structures show an insulating direct E_g with the position of the conduction band minimum (CBM) near the Fermi energy (E_F), but after including SOC, the surface shows metallic behaviour with the single Dirac cone at the Γ point [Fig. 3e], which is the hallmark for TIs.^2,38,39 The occurrence of the Dirac point (DP) below E_F is consistent with previous reports of Bi₂Se₃.^40,35

Fig. 3 The (0001) surface band structures of Bi₂Se₃ without (a–c) and with SOC (d–f) under the ε_c strain. (g) The ε_c strain-dependent band gap of the (0001) surface of Bi₂Se₃ without and with SOC. (h) The position of the Dirac point energy with the variation of ε_c.

As the 8 QLs of the Bi₂Se₃ slab are compressed uniaxially by decreasing the ‘c’ value keeping ‘a’ and ‘b’ fixed, the direct E_g decreases gradually and becomes indirect at 4.9% compressive strain due to the formation of “M” shape at the VB maximum [Fig. 3a]. The E_g vanishes at a very large 10% compression [Fig. 3g]. This metallic state formation in the (0001) surface of the Bi₂Se₃ slab has been explained by the increase in the band width of p_z orbitals of Bi and Se in the valence band region.³⁶ The indirect E_g in the (0001) surface of the Bi₂Se₃ has been observed previously for different combinations of lattice parameters under uniaxial strain.⁴¹ The compressive uniaxial strain has been imposed on the Sb₂Se₃ for the realization of its metallic state.²⁸ On the other hand, as the tensile stress increases, i.e., the value of ‘c’ increases above 100.9 Å, the E_g monotonically increases and moves the p_z orbital of Bi and Se atoms away from the E_F, preserving the direct E_g [Fig. 3c and g]. Although achieving such a large value of strain is practically challenging, a recent experimental study demonstrates the occurrence of 9.6% strain along ‘a’ and 19.6% along ‘c’ through the intercalation of 60% zero-valent Cu into Bi₂Se₃ nanoribbons without disrupting the host lattice. The ‘a’ and ‘c’ lattice constants of Bi₂Se₃ increase from 4.14 ± 0.01 Å and 28.5 ± 0.01 Å to 4.54 ± 0.74 Å and 34.1 ± 0.76 Å, after the intercalation of Cu, respectively.²¹

The SOC makes the (0001) surface of Bi₂Se₃ conducting for an ample range of ‘c’ values. The Dirac cone appears deep inside the VB region at −0.16 eV due to the −4.9% ε_c strain, and p_x and p_y orbitals of Bi and the p_z orbital of Se1 atoms disperse along the Γ–M–K direction [Fig. 3d]. As the ε_c strain decreases, DP moves up towards E_F and appear very close to E_F (at −0.02 eV) for 4.9% ε_c strain [Fig. 3f], and the flat band feature in the VB is suppressed. The observed shift of DP energy in the band structure of the (0001) surface of Bi₂Se₃ by the variation of strain is consistent with previous studies on the strain-dependent DP energy of the (112̄1) surface of Bi₂Se₃.⁸ On further increasing the tensile ε_c strain, i.e., increasing the ‘c’ value, the band structure opens a finite energy gap (0.05 eV) with the annihilation of the Dirac point [Fig. 3g]. The reports on the strain-induced annihilation of DP of the (0001) surface of the Bi₂Se₃ and the transition to a topologically trivial insulating phase is consistent with previous reports.^8,28

The non-degenerate helical surface states of topological insulators provide many interesting topological phenomena and are the starting point for the Majorana excitations.⁴² The Dirac point should be at the E_F so that the conventionally confined states within the vertex are well separated from the Majorana-type excitations. The observed (0001) surface band structures of Bi₂Se₃ with the variation of ε_c strain confirmed that the uniaxial strain could be an effective tool to tune DP to be present at E_F. The surface band structures show that the compressive strain moves DP away from E_F, whereas the tensile strain moves the DP towards E_F [Fig. 3h]. In addition, the variation of ε_c strain provides the change of (0001) surface of Bi₂Se₃ from metallic-to-topological-to-normal insulator.^28,36 The normal-to-topological insulator transition has been observed in β-As₂Te₃ by the application ε_c strain.⁴³

The (0001) surface band structure of Bi₂Se₃ is investigated by applying the biaxial strain in the range of ±10%. The (0001) surface band structure of Bi₂Se₃ with −8.1% of ε_a strain exhibits the indirect E_g and the p_x and p_y orbitals of Bi atoms and p_z orbitals of Se atoms are dispersed near E_F [Fig. 4a]. As the compressive ε_a strain decreases, p_x and p_y orbitals of the Bi atoms disperse away from E_F, while p_z orbitals of Bi atoms disperse towards E_F, thus, the E_g becomes direct at −5.7% of ε_a strain [Fig. 4b]. In addition to the re-dispersion of orbitals, E_g increases from 0.28 eV to 0.38 eV, and the indirect-to-direct E_g transition occurs by the decrease of ε_a from −8.1% to −5.7%. On further decreasing the ε_a strain, i.e., increasing the ‘a’ value, p_z orbitals of Bi and Se approaches towards E_F, thus, the direct E_g decreases monotonically [Fig. 4c, d and i]. Interestingly, the biaxial strain is very sensitive to tune E_g, the change of in-plane lattice parameter by 0.1 Å corresponds to approximately 2% variation in strain, which changes E_g by 0.02 eV. The change in E_g with ε_a is consistent with the variation of E_g in N-doped Sb₂Te₃, where an appropriate in-plane strain can enlarge the bulk E_g.⁴⁴ It is to be noted that the Bi₂Se₃ slab shows the direct E_g from −5.7% to 8.8% biaxial strain. The observation of direct E_g in a similar range of strain (−6% to +6% along the x-direction and −6% to 10% along the y-direction) is reported previously in Bi₂Se₃.⁴⁵ On incorporating the SOC effect, the Dirac like dispersion is observed for the entire range of ε_a strain [Fig. 4e–h]. Interestingly, the DP is exactly at the Γ point at −8.1% of ε_a strain. As the ε_a strain decreases from −8.1%, the DP energy moves away from E_F [Fig. 4j]. It is interesting to note that both compressive and tensile ε_a strain could not open the E_g at the DP. The Dirac cone is preserved for the entire range of ε_a strain, which indicates that the topological properties remain preserved. A similar variation of tensile biaxial strain on the band structure of the (112̄1) Bi₂Se₃ slab was observed for one to six QLs thickness with the SOC effect.³⁶ However, the result is limited to the tensile ε_a strain study and the effect of ε_a strain without SOC on the E_g of Bi₂Se₃ is very scarce.³⁶

Fig. 4 The (0001) surface band structures of the Bi₂Se₃ without (a–d) and with SOC (e–h) under the ε_a strain. (i) The ε_a strain-dependent band gap of the (0001) surface of the Bi₂Se₃ without and with SOC. (j) The position of the Dirac point energy with the variation of ε_a.

The influence of strength of strain while keeping the volume of the Bi₂Se₃ slab constant, on the (0001) surface band structure of Bi₂Se₃ is also investigated. At −7.0% of compressive volume conservation strain, the (0001) surface of Bi₂Se₃ behaves like a degenerate semiconductor as p_z orbitals of Bi atoms cross the E_F [Fig. 5a]. However, the (0001) surface of Bi₂Se₃ behaves like a direct E_g semiconductor by reducing the compressive volume conservation strain to −6.0%, which moves p_z orbitals of Bi atoms away from E_F [Fig. 5b]. On further reducing the volume conservation strain from the compressive state to the tensile state, the (0001) surface gap increases up to 5.8% [Fig. 5c]. In addition, p_x and p_y orbital dispersions of Bi atoms spread towards E_F along the Γ–M direction. The (0001) surface of Bi₂Se₃ shows that the Dirac surface state feature with the Dirac point below the Fermi level for a smaller value of tensile strain [Fig. 5d]. Above 5.8% volume conservation strain, the (0001) surface gap decreases and becomes indirect as the p_z orbital of the Se atoms come closer to the p_x and p_y orbitals of the Bi atoms [Fig. 5g]. The decrease in the (0001) surface E_g at higher tensile strain could be due to the compressive ε_a, which increases the band width. Note that the study of the band structure of Bi₂Se₃ with the variation of ε_c = −ε_a strain is very rare to compare our results. In the presence of SOC, the (0001) surface of Bi₂Se₃ is conducting if the strength of the applied strain with the conservation of the volume of the slab is between −7.9% to 4.9% [Fig. 5g]. The (0001) surface of Bi₂Se₃ is a normal semiconductor if the strain on the slab is larger than 4.9%. For ε_c = −ε_a > 4.9%, the band structures of the (0001) surface show an opening of the energy gap, which increases slowly with the increase in the strain [Fig. 5e and f]. For ε_c = −ε_a = −9.9%, DP is deep inside VB at −0.345 eV below E_F [Fig. 5h]. As the volume conservation strain decreases from −9.9%, DP moves towards E_F and appear at E_F for ε_c = −ε_a = 5.9%. On further increasing volume conservation strain, the Dirac cone feature of Bi₂Se₃ gets annihilated, and a finite E_g emerges at the Γ point. The position of the DP energy with the variation of ε_c = −ε_a strain is plotted in Fig. 5h, which is qualitatively very similar to the case of ε_c strain variation.

Fig. 5 The (0001) surface band structures of Bi₂Se₃ without (a–c) and with SOC (d–f) under the volume conservation strain. (g) The volume conservation strain-dependent band gap of the (0001) surface of Bi₂Se₃ without and with SOC. (h) The position of the Dirac point energy with the variation of volume conservation strain.

There are several DFT studies on the DP energy of Bi₂Se₃ surface,^7,8,28,40 which interestingly indicates that the position of the DP energy depends on the orientation of the Bi₂Se₃ surface. For example, the DP exactly appears at E_F in the unstrained condition of (112̄1) oriented rhombohedral Bi₂Se₃.² In contrast, the DP energy is located below and above E_F in (0001) and (11̄00) surfaces of the unstrained hexagonal Bi₂Se₃, respectively.^40,46 The electrons are naturally accumulated on the surface of the real Bi₂Se₃ single crystals because of Se deficiency, E_F appears in the conduction band, which makes the detection of topological surface states difficult. Thus, the tuning of E_F position in the band gap is necessary for the realization of the topological surface states of Bi₂Se₃.⁴⁷ Interestingly, our band structure calculations of the (0001) surface of Bi₂Se₃ show the tuning of the energy position of DP by the application of strain on Bi₂Se₃ slab. The DP energy can be tuned from −0.318 eV to −0.02 eV by changing uniaxial/volume conservation strain from −9.1% to 6.9%. As the uniaxial/volume conservation strain switches from compressive to the tensile state, the DP energy in CB moves towards E_F. However, DP energy shifts from −0.009 eV to −0.145 eV by changing biaxial strain on the Bi₂Se₃ slab from −8.12% to 8.8%, the opposite effect of uniaxial/volume conservation strain. Moreover, the (0001) surface of Bi₂Se₃ becomes insulating, even in the presence of SOC for large tensile volume conservation strain, which indicates the reduction in the SOC strength. However, it is reported that the volume conservation strain on Bi₂Se₃ with the (11̄00) oriented surface lifts the spin degeneracy and shifts the DP energy away from the Γ point.⁴⁶ For (112̄0) surface Bi₂Se₃, the DP energy appears above and below E_F for the tensile and compressive strain on the Bi₂Se₃ slab, respectively.⁴⁶

Fig. 3–5 show that the occurrence of the redistribution of the orbital population of the (0001) surface due to the variable strain on Bi₂Se₃. The compressive ε_a strain increases the p_x and p_y orbitals overlap, which moves the orbital population towards the E_F along the Γ–M–K in the CB. Thus, the band gap becomes indirect and shows a peak in the band gap vs. strain graph [Fig. 4i]. A similar situation arises while applying tensile volume conservation strain, and a peak appears in the band gap vs. strain graph [Fig. 5g]. However, under compressive ε_c strain, the dispersion of p_z orbitals is dominated at the Γ point, and a direct band gap is observed in the +4.9% to +10% ε_c strain region. The band gap linearly varies with the ε_c strain without the application of SOC.

The result indicates that the (0001) surface of hexagonal Bi₂Se₃ is very sensitive to uniaxial strain, which can drive the (0001) surface state from metallic-to-topological-to-semiconducting nature with a maximum surface gap of 0.6 eV. The topological phase is observed in the range of ±6% variation of ‘c’. In contrast, the biaxial strain controls the dispersion of p_x and p_y bands, which toggles the (0001) surface of Bi₂Se₃ between the metallic and topological states. The topological state is observed for the variation of ‘a’ from −8% to 8%. Thus, the probability of achieving the topological surface states in the (0001) oriented hexagonal Bi₂Se₃ is higher for biaxial strain than the case of the uniaxial and volume conservation strain. Interestingly, biaxial strain could not create a finite energy gap at the Dirac point; hence, it does not play a significant role in destroying topological surface states. Therefore, the annihilation of the Dirac point largely depends on the out-of-plane interaction among the atomic orbitals of Bi and Se. The variation of p_x, p_y and p_z orbitals with ε_c = −ε_a alleviate the (0001) surface of Bi₂Se₃ similar to the case of uniaxial strain. However, the tensile volume conservation strain distinguishes a transition from a direct E_g state to an indirect E_g state from the strain-dependent E_g plot without SOC. Overall, the Dirac cone can be tuned to E_F by the application of tensile uniaxial, compressive biaxial, and tensile volume conservation strain.

Further, the atomic cohesive energy of the (0001) Bi₂Se₃ slab is calculated for different values of uniaxial, biaxial, and volume conservation strain. Fig. 6 shows the variation of cohesive energy per atom in the Bi₂Se₃ slab with the strain. The cohesive energy (E_coh) is calculated using the expression;

where n_Bi and n_Se are the number of atoms in the slab and E_Bi and E_Se are single isolated atom energies of Bi and Se, respectively. The E_coh of the reference structure is 2.57 eV, which is close to the theoretically calculated value and 2.97 eV.⁴⁸

Fig. 6 Atomic cohesive energy of the 8 QLs Bi₂Se₃ slab under the uniaxial, biaxial, and volume conservation strain.

To compare the theoretically calculated cohesive energy with that of Bi₂Se₃ nanocrystals, Raman spectra of Bi₂Se₃ nanocrystals were recorded at various temperatures above the Debye temperature (∼185 K)⁴⁹ of the Bi₂Se₃ [Fig. 7]. The Raman spectra of these nanocrystals in the assessable Raman shift range of the Raman spectrometer exhibits a Raman line appearing at ∼131.2 and ∼173.4 cm⁻¹, which correspond to the E_g and A_1g modes, respectively. The E_g mode appears because of the in-plane vibrations of Se and Bi atoms and A_1g mode appears from the out-of-plane vibrations of Se and Bi atoms.⁵⁰ The peak position of different modes in Raman spectra of Bi₂Se₃ are extracted using the Lorentz fit and plotted in Fig. 8. The Raman shift (ω(z_b,T)) of both modes varies linearly with temperature consistent with the following expression;⁴⁹

where z_b is the bulk coordination number, C_V is the three-dimensional specific heat per bond, and ω(1) is the vibrational frequency of an isolated dimer, which is the reference point for the optical redshift upon nanosolid and bulk. The above expression of ω(z_b,T) is valid only for the material with negligible thermal expansion at high T, especially when the temperature is larger than the Debye temperature. The linear fit to the ω(z_b,T) of the E_g and A_1g modes provides the ω(1) and E_coh of Bi₂Se₃ nanocrystals. The ω(1) is found to be 36.9 cm⁻¹, which is very close to the value reported for Bi₂Se₃.⁴⁹ The E_coh of Bi₂Se₃ nanocrystals calculated from the ω(z_b,T) expression is 2.58 eV, which is larger than the values reported for Bi₂Se₃ nanostructures as 1.24 eV (ref. 49) and 1.37 eV.⁵¹ However, E_coh extracted from Raman scattering is very close to the value calculated from the first principles density functional theory. The cohesive energy remains close to 2.57 eV for ± 4% the variation of strain, which could be the foundation of the (0001) oriented growth of Bi₂Se₃ nanocrystals. Note that, the variation of cohesive energy is higher in biaxial strain as compared to the uniaxial and volume conservation strain state for the strain larger than ± 4%. The strain-dependent cohesive energy indicates that at a higher value of strain, uniaxial and volume conservation strain provides better stability than that of the biaxial strain (0001) oriented growth of Bi₂Se₃ nanocrystals. The strain (±4%) independent cohesive energy and the tuning of the Dirac point position with strain will shed light on achieving the dissipationless (0001) surface transport in Bi₂Se₃, which is the essential requirement for the application of topological insulators in topological quantum computation and low-power spintronic devices.

Fig. 7 Raman spectra of the Bi₂Se₃ hexagon plates synthesized at 250 °C, recorded at various temperatures.

Fig. 8 Temperature-dependent Raman shift of the E_g and A_1g modes of the Bi₂Se₃ hexagon plates synthesized at 250 °C. The solid line is the fit to the data.

5. Conclusion

In conclusion, the (0001) oriented Bi₂Se₃ nanocrystals were prepared by a hot injection method using a nontoxic solvent. The atomic cohesive energy of these nanocrystals extracted from the temperature-dependent Raman spectra is 2.58 eV. By utilizing the lattice parameters of these nanocrystals, an 8 QLs thick Bi₂Se₃ slab is constructed to explore the effects of strain on the (0001) surface states and the cohesive energy calculated using the first-principles DFT calculations without and with the presence of SOC effects. The CB band width of the Se2 shrinks as well as shifts away from E_F compared to that of Bi and Se1 irrespective of the nature or directions of the applied strain. The band structures of the (0001) surface of Bi₂Se₃ show that Se2 is very sensitive to mechanical strain, and strain can tune the Dirac point energy. The Dirac cone feature of the (0001) surface of Bi₂Se₃ is preserved for the entire range of the biaxial strain. However, around 5% tensile uniaxial strain and even lower value of volume conservation strains annihilate the Dirac cone, which causes the loss of topological (0001) surface states of Bi₂Se₃. The biaxial strain provide ease in achieving the Dirac cone at E_F than uniaxial and volume conservation strains. The volume conservation strain distinguishes a transition from a direct E_g state to an indirect E_g state in the strain-dependent E_g plot without the SOC. The strain-dependent cohesive energy indicates that at a higher value of strain, uniaxial and volume conservation strain provides better stability compared to that of the biaxial strain (0001) oriented growth of Bi₂Se₃ nanocrystals. The application of mechanical strain on Bi₂Se₃ slab has a significant influence on the band width of ‘p’ orbitals of Bi and Se, E_g and Dirac point energy, which may provide a new pathway to control many physical properties of Bi₂Se₃ with strain and shade light on designing artificial topological materials for technological applications.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

We acknowledge the use of the computing resources at HPCE, IIT Madras.

References

1. M. Z. Hasan and C. L. Kane, Rev. Mod. Phys., 2010, 82, 3045.
2. H. Zhang, C.-X. Liu, X. L. Qi, X. Dai, Z. Fang and S. C. Zhang, Nat. Phys., 2009, 5, 438.
3. Y. Zhao, Y. Hu, L. Liu, Y. Zhu and H. Guo, Nano Lett., 2011, 11, 2088.
4. J.-M. Zhang, W. Zhu, Y. Zhang, D. Xiao and Y. Yao, Phys. Rev. Lett., 2012, 109, 266405.
5. D. Pesin and A. H. MacDonald, Nat. Mater., 2012, 11, 409.
6. H. Aramberri, J. I. Cerdá and M. C. Muñoz, Nano Lett., 2015, 15, 3840.
7. L. Zhao, J. Liu, P. Tang and W. Duan, Appl. Phys. Lett., 2012, 100, 131602.
8. Y. Liu, Y. Y. Li, S. Rajput, D. Gilks, L. Lari, P. L. Galindo, M. Weinert, V. K. Lazarov and L. Li, Nat. Phys., 2014, 10, 294.
9. Y. Saeed, N. Singh and U. Schwingenschlogl, Appl. Phys. Lett., 2014, 105, 031915.
10. S. H. Park, J. Chae, K. S. Jeong, T.-H. Kim, H. Choi, M.-H. Cho, I. Hwang, M.-H. Bae and C. Kang, Nano Lett., 2015, 15, 3820.
11. S. M. Young, S. Chowdhury, E. J. Walter, E. J. Mele, C. L. Kane and A. M. Rappe, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 085106.
12. Y. Saeed, N. Singh and U. Schwingenschlogl, Appl. Phys. Lett., 2014, 104, 033105.
13. J. Z. Zhao, L. C. Chen, B. Xu, B. B. Zheng, J. Fan and H. Xu, Phys. Rev. B, 2020, 101, 121407(R).
14. Y. Sun, X.-Q. Chen, S. Yunoki, D. Li and Y. Li, Phys. Rev. Lett., 2010, 105, 216406.
15. T.-H. Kim, K. S. Jeong, B. C. Park, H. Choi, S. H. Park, S. Jung, J. Park, K.-H. Jeong, J. W. Kim, J. H. Kim and M.-H. Cho, Nanoscale, 2016, 8, 741.
16. J. Chae, S.-H. Kang, S. H. Park, H. Park, K. Keong, T. H. Kim, S.-B. Hong, K. S. Kim, Y. K. Kwon, J. W. Kim and M.-H. Cho, ACS Nano, 2019, 13, 3931.
17. G.-C. Wang and T.-M. Lu, J. Appl. Phys., 2019, 125, 082401.
18. Q. Li, Y. Wang, T. Li, W. Li, F. Wang, A. Janotti, S. Law and T. Gu, ACS Omega, 2020, 5, 8090.
19. D. Flötotto, Y. Bai, Y.-H. Chan, P. Chen, X. Wang, P. Rossi, C.-Z. Xu, C. Zhang, J. A. Hlevyack, J. D. Denlinger, H. Hong, M.-Y. Chou, E. J. Mittemeijer, J. N. Eckstein and T.-C. Chiang, Nano Lett., 2018, 18, 5628.
20. S. Qi, H. Yang, J. Chen, X. Zhang, Y. Yang and X. Xu, Sci. Rep., 2016, 6, 29161.
21. K. J. Koski, J. J. Cha, B. W. Reed, C. D. Wessells, D. Kong and Y. Cui, J. Am. Chem. Soc., 2012, 134, 7584.
22. M. V. Sidorov, M. J. McKelvy, J. M. Cowley and W. S. Glaunsinger, Chem. Mater., 1998, 10, 3290.
23. C. Schindler, C. Wiegand, J. Sichau, L. Tiemann, K. Nielsch, R. Zierold and R. H. Blick, Appl. Phys. Lett., 2017, 111, 171601.
24. P. Giannozzi, et al. , J. Phys.: Condens. Matter, 2009, 21, 395502.
25. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
26. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272.
27. Y. Sun, H. Cheng, S. Gao, Q. Liu, Z. Sun, C. Xiao, C. Wu, S. Wei and Y. Xie, J. Am. Chem. Soc., 2012, 134, 20294.
28. W. Liu, X. Peng, C. Tang, L. Sun, K. Zhang and J. Zhong, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 245105.
29. S. K. Das and P. Padhan, ACS Appl. Nano Mater., 2020, 3, 274–282.
30. G. Martinez, B. A. Piot, M. Hakl, M. Potemski, Y. S. Hor, A. Materna, S. G. Strzelecka, A. Hruban, O. Caha, J. Novák, A. Dubroka, Č. Drašar and M. Orlita, Sci. Rep., 2017, 7, 6891.
31. M. Orlita, B. A. Piot, G. Martinez, N. K. Sampath Kumar, C. Faugeras, M. Potemski, C. Michel, E. M. Hankiewicz, T. Brauner, Č. Drašar, S. Schreyeck, S. Grauer, K. Brunner, C. Gould, C. Brüne and L. W. Molenkamp, Phys. Rev. Lett., 2015, 114, 186401.
32. A. Dubroka, O. Caha, M. Hroncek, P. Friš, M. Orlita, V. Holý, H. Steiner, G. Bauer, G. Springholz and J. Humlícek, Phys. Rev. B, 2017, 96, 235202.
33. S. K. Mishra, S. Satpathy and O. Jepsen, J. Phys.: Condens. Matter, 1997, 9, 461.
34. M. Jurczyszyn, M. Sikora, M. Chrobak and L. Jurczyszyn, Appl. Surf. Sci., 2020, 528, 146978.
35. Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava and M. Z. Hasan, Nat. Phys., 2009, 5, 398.
36. H. Aramberri and M. C. Muñoz, Phys. Rev. B, 2017, 95, 205422.
37. S. Chege, P. Ning'I, J. Sifuna and G. O. Amolo, AIP Adv., 2020, 10, 095018.
38. O. V. Yazyev, J. E. Moore and S. G. Louie, Phys. Rev. Lett., 2010, 105, 266806.
39. W. Zhang, R. Yu, H.-J. Zhang, X. Dai and Z. Fang, New J. Phys., 2010, 12, 065013.
40. A. Bera, K. Pal, D. V. S. Muthu, U. V. Waghmare and A. K. Sood, J. Phys.: Condens. Matter, 2016, 28, 105401.
41. H. Aramberri and M. C. Muñoz, J. Phys.: Condens. Matter, 2018, 1, 015009.
42. J. Kellner, M. Eschbach, J. Kampmeier, M. Lanius, E. Mły'nczak, G. Mussler, B. Holländer, L. Plucinski, M. Liebmann, D. Grützmacher, C. M. Schneider and M. Morgenstern, Appl. Phys. Lett., 2015, 107, 251603.
43. K. Pal and U. V. Waghmare, Appl. Phys. Lett., 2014, 105, 062105.
44. S. Qi, R. Gao, M. Chang, T. Hou, Y. Han and Z. Qiao, Phys. Rev. B, 2020, 102, 085419.
45. Z. Lu, Y. Wu, Y. Xu, C. Ma, Y. Chen, K. Xu, H. Zhang, H. Zhub and Z. Fang, Nanoscale, 2019, 11, 20620.
46. Y. Takagaki, B. Jenichen and J. Tominaga, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 245302.
47. E. Uesugi, T. Uchiyama, H. Goto, H. Ota, T. Ueno, H. Fujiwara, K. Terashima, T. Yokoya, F. Matsui, J. Akimitsu, K. Kobayashi and Y. Kubozono, Sci. Rep., 2019, 9, 5376.
48. Y. Liu, Y. Xu, Y. Jia and H. Zhang, Phys. Chem. Chem. Phys., 2020, 22, 9685.
49. X. X. Yang, Z. F. Zhou, Y. Wang, R. Jiang, W. T. Zheng and C. Q. Sun, J. Appl. Phys., 2012, 112, 083508.
50. J. Zhang, Z. Peng, A. Soni, Y. Zhao, Y. Xiong, B. Peng, J. Wang, M. S. Dresselhaus and Q. Xiong, Nano Lett., 2011, 11, 2407.
51. L. B. Abdalla, L. Seixas, T. M. Schmidt, R. H. Miwa and A. Fazzio, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 045312.

---

Footnote

† Electronic supplementary information (ESI) available: The schematic of crystal symmetry, X-ray diffraction spectra, high-resolution transmission electron microscopy image, and bulk band structure of Bi₂Se₃ nanocrystals. See DOI: 10.1039/d1na00139f

---