Unraveling the dislocation core structure at a van der Waals gap in bismuth telluride

Tetradymite-structured chalcogenides such as bismuth telluride (Bi2Te3) are of significant interest for thermoelectric energy conversion and as topological insulators. Dislocations play a critical role during synthesis and processing of such materials and can strongly affect their functional properties. The dislocations between quintuple layers present special interest since their core structure is controlled by the van der Waals interactions between the layers. In this work, using atomic-resolution electron microscopy, we resolve the basal dislocation core structure in Bi2Te3, quantifying the disregistry of the atomic planes across the core. We show that, despite the existence of a stable stacking fault in the basal plane gamma surface, the dislocation core spreading is mainly due to the weak bonding between the layers, which leads to a small energy penalty for layer sliding parallel to the van der Waals gap. Calculations within a semidiscrete variational Peierls-Nabarro model informed by first-principles calculations support our experimental findings.


Supplementary Note 1: Semidiscrete variational Peierls-Nabarro model
The dislocation disregistry vector δ i (index i labels the Cartesian components) is defined as the difference in displacement between the material just above and just below the slip plane containing the dislocation. Here, we consider a perfectly straight dislocation with the line parallel to the z-direction and the slip plane normal to the y-direction. If the dislocation is planar, i.e., only spreads along the slip plane, then δ i is only a function of the x-coordinate. Using this coordinate system, δ 1 corresponds to the edge component of the disregistry, δ 3 to the screw component, and δ 2 to the out-of-plane component.
The classical Peierls-Nabarro model [1,2] represents the dislocation core by a continuous spread of infinitesimal Volterra dislocations. Thus δ i (x) is a continuous function satisfying the boundary conditions δ i (∞) − δ i (−∞) = b i , where b i is the dislocation Burgers vector. By assuming a particular functional form for δ i (x), the dislocation energy can be evaluated by integrating the interactions of the infinitesimal dislocations with each other and with the crystal over all x and adding the misfit energy.
In the semidiscrete variational Peierls-Nabarro (SDVPN) model [3][4][5], the disregistry δ i is only evaluated at discrete values of x that are equally spaced by ∆x. Typically, ∆x is taken to correspond to a spacing of atomic columns along the x-direction, giving the resulting model an atomic character. The full dislocation is represented by a number of equally spaced partial dislocations. Their density ρ i (x) is defined as the numerical derivative of the disregistry, which is assumed to be constant between the evaluated x coordinates. By its nature and definition, The total energy of the SDVPN dislocation at zero applied stress is represented as the sum of the energy arising from the dislocation's elastic strain field and the crystal misfit due to the dislocation spreading: Both terms are expressed as functions of δ i (x), which allows for the disregistry profile to be solved for by minimizing the total energy. The elastic energy of the dislocation accounts for the interaction energy between all α-β pairs of partial dislocations positioned at the N evaluated x-coordinates: where K ij is the energy coefficient tensor and χ is a multiplicative factor that depends on the distance between the dislocation positions x α and x β : The leading coefficient in Eq. (4) depends on how K ij is defined. Here, K ij is taken such that in the isotropic elasticity case K 11 = K 22 = µ/(1 − ν) and K 33 = µ, where µ and ν are the shear modulus and Poisson's ratio, respectively. The misfit energy is obtained by mapping the disregistry profile onto the gamma-surface (i.e., the generalized stacking fault energy γ as a function of the translation vector): The γ function can either be 2D, in which case it only depends on the in-plane components of the disregistry, δ 1 and δ 3 , or 3D, in which case it includes the out-of-plane disregistry component δ 2 . Some SDVPN models include an additional correction term to account for either the surface effect of the slip plane [6,7] or non-local interactions between the evaluated positions [8,9]. Neither correction was included here as they require assumptions of the dislocation symmetry along the slip plane or are numerically fitted using atomistic data. The anisotropic energy coefficient tensor K ij was obtained by solving the Eshelby elastic model for a straight dislocation [10] using the Stroh method [2,11,12] with elastic constants calculated with each of the DFT functionals. The resulting K ij has the form The four non-zero elements represent a screw component K 33 , an edge component in the slip plane K 11 , an edge component perpendicular to the slip plane K 22 , and a edge coupling component K 12 . For all sets of the elastic constants C ij explored here, the general trends in the K ij components are that K 11 and K 22 are comparable but not equal, K 33 is slightly smaller than the edge components, and K 12 is roughly an order of magnitude smaller than the other components. Solving for the dislocation disregistry requires inputs for γ(δ i ) and K ij at the relaxed lattice parameter a. These values were obtained using four different exchange-correlation functionals: LDA [13], optPBE-vdW and optB88-vdW [14][15][16][17], and DFT-D2 [18]. Table  1 summarizes the lattice and elastic constants of Bi 2 Te 3 computed by DFT in comparison with experiment. (Note that for rhombohedral crystals such as Bi 2 Te 3 , 2C 66 = C 11 − C 12 .) While the lattice constant displays decent agreement across the DFT functionals, the elastic constants show more variability, including C 13 and C 33 . The sensitivity of C 13 and C 33 to the choice of the DFT functional is due to the component 3 direction being normal to the van der Waals-bonded layers in the Bi 2 Te 3 crystal structure.
An initial disregistry guess was made using an arctangent expression similar to the one used in the classic Peierls-Nabarro model, where ξ is the dislocation halfwidth and S is a scaling factor introduced to satisfy the boundary conditions despite the fact that x(0) and x(N ) are finite. The ξ value used in the initial disregistry was the one that minimized Eq. (3) using the disregistry of the form of Eq. (9). For the calculations performed here, we chose N = 300 and ∆x = a √ 3/3, which resulted in x values spanning the range ±38 nm. This wide range was necessary to ensure that the boundary conditions did not affect the dislocation core with its large spreading. As the gamma-surface obtained by the DFT calculations only considered the in-plane translations, the out-of-plane disregistry component δ 2 was excluded from the calculations. As a result, the energy coefficient tensor only included the K 11 and K 33 components. Starting with the initial guess (9), the disregistry profile was optimized using Powell's minimization method [23] subject to the boundary conditions (10). For each calculation, the minimization algorithm was restarted ten times for better convergence.

Supplementary Note 2: Burgers circuit construction
The Burgers vectors of the experimentally observed dislocations were determined by circuit mapping. This approach consists of constructing a closed path of crystallographic translation vectors around the atomically resolved image of the defect. Using the FSRH convention [2], with the line direction out of the board, the Burgers vector b is given by: where C i are the individual segments of the circuit path.
LDA [13] optPBE-vdW [14][15][16][17] DFT-D2 [18] optB88-vdW [14][15][16][17] Experim. [19][20][21][22] a (nm) In the example shown in Figure 1c, To demonstrate the possible effect of the dislocation core structure in Bi 2 Te 3 on electronic properties, DFT calculations (using the optB88-vdW exchange-correlation functional) have been performed for a Bi 2 Te 3 slab containing six quintuplet layers (we employ a hexagonal unit cell with 15 atoms). We considered two systems: one without stacking faults and another with a stacking fault in the middle of the slab. The stacking fault structure is representative of the dislocation core structure. The Supplementary Figure 1(a,b) shows the electronic band structure of the two systems, presented in the form of a spectral function for an energy range near the Fermi level (set to zero) and for wave-vectors near the Γ point along the Γ-M direction of the Brillouin zone of the corresponding hexagonal reciprocal space. In the system without the stacking fault (panel a), the overall bands tructure is in agreement with that obtained in previous DFT studies (compare to Fig. 4a in [24]). We note the presence of topological surface states intersecting at the Dirac point located near zero energy at the Γ point. These states are not affected significantly by the presence of the stacking fault in the middle of the slab. To disentangle the surface states from the rest we show in Supplementary Figure 1(c-d) the electronic band structure of the two systems projected on one of the middle quintuplets. In the absence of a stacking fault (panel c) the bulk-like states show a gap at the Γ point of about 0.25 eV. These states are impacted by the stacking fault, which leads to a valence-conduction gap separation becoming significantly smaller, namely about 0.18 eV. We believe that this local band gap renormalization is due to the change in local strain caused by the stacking fault. Indeed, previous ab initio calculations have shown that the band gap of bulk Bi 2 Te 3 at the Γ point tends to decrease when subject to uniaxial tensile strain [25]. We find that the quintuplet separation at the stacking fault is about 10% larger than in the absence of the stacking fault. This expansion can be deduced from a plot of the electron charge density distribution function, as shown in Supplementary Figure 1 for a particular isovalue surface (about a quarter of the maximum value). Comparing the two cases it is apparent that the structural relaxation in the stacking fault region increases the gap between the quintuplet layers. The local renormalization (decrease) of the band gap is likely to impact (increase) the concentration of free charge carriers near the stacking fault. A more detailed analysis (outside the scope of this work) is needed to better understand the impact of dislocation cores on the electronic transport properties of Bi 2 Te 3 .