Shear induced deformation twinning evolution in thermoelectric InSb

Twin boundary (TB) engineering has been widely applied to enhance the strength and plasticity of metals and alloys, but is rarely adopted in thermoelectric (TE) semiconductors. Our previous first-principles results showed that nanotwins can strengthen TE Indium Antimony (InSb) through In–Sb covalent bond rearrangement at the TBs. Herein, we further show that shear-induced deformation twinning enhances plasticity of InSb. We demonstrate this by employing large-scale molecular dynamics (MD) to follow the shear stress response of flawless single-crystal InSb along various slip systems. We observed that the maximum shear strain for the (111)[112¯]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(111)[11\bar 2]$$\end{document} slip system can be up to 0.85 due to shear-induced deformation twinning. We attribute this deformation twinning to the “catching bond” involving breaking and re-formation of In–Sb bond in InSb. This finding opens up a strategy to increase the plasticity of TE InSb by deformation twinning, which is expected to be implemented in other isotypic III–V semiconductors with zinc blende structure.


INTRODUCTION
Indium Antimony (InSb) is found to be one of the most potential high-performance thermoelectric (TE) material in III-V semiconductors [1][2][3][4][5][6] . InSb is also widely used as infrared detectors because of its narrow band gap (0.18 eV at 300 K) and the highest mobility among III-V semiconductors (about 7.8 × 10 4 cm 2 V −1 s −1 at 300 K) [7][8][9][10][11] . However, as for TE materials, the cycling of thermal stress under working conditions causes the crack opening, leading to the deterioration of material performance and even the failure. Although as for infrared detectors, its working temperature (at 77 K) is much lower than its storage temperature (at 300 K), leading to the coefficient of thermal expansion mismatch that causes the brittle cleavage [12][13][14][15][16][17][18] . In recent year, some ductile inorganic semiconductors with special intrinsic structure have been found, such as α-Ag 2 S 19,20 and InSe 21 , which have outstanding mechanical properties. However, for the other semiconductors without excellent intrinsic properties such as InSb, these deficiencies in mechanical properties severely limit its extended engineering applications. Therefore, the mechanical properties of InSb need to be improved for its industrial applications.
One effective way to enhance the material mechanical properties is via structural design. A specific two-dimensional structural design, twin boundary (TB), has beneficial effects on structural mechanical properties, which has attracted attention in recent years. TBs provide low energy grain boundaries that can strengthen the materials by minimizing effective grain size to the nanoscale, while hindering dislocation movement, which is well-known as the Hall-Petch effect 22,23 . The dislocations change from glissile dislocations to sessile dislocations while crossing TBs and thus protect the structure from slip failure, which is called as Basinski mechanism 24 . During deformation twinning, the increasing proportion of twins may change the lattice orientation to a hard direction, which is noted as texture hardening. Because of the above effects of twins, introducing TBs to optimize materials has been widely applied in metals and ceramics [25][26][27][28][29] , yet its effects on TE semiconductors remain less explored. Growth twins can be produced in TE semiconductor InSb during preparation [30][31][32][33] , with a very low twin interfacial energy (20.7 mJ m −2 ) 34 , making introducing TBs appropriate to tune mechanical performance.
In our previous study, we used density functional theory (DFT) to investigate the shear stress response of nanotwin InSb and found that the directional covalent In-Sb bond rearrangements at TBs strengthen structural stiffness, which improves the ideal shear strength of InSb by 11% 34 . However, the nucleation mechanism for nanotwins and the mechanical properties of large-scale nanotwinned InSb remain unexplored. In order to investigate the mechanical properties of InSb under realistic conditions, we employed molecular dynamics (MD) simulations to examine shear stress response along various slip systems for flawless singlecrystal InSb. Interestingly, we found deformation twinning forms easily along the ð111Þ½112 slip system in InSb. Similar phenomena have been observed in other materials that deformation may improve plasticity with strain-induced dislocation 35,36 . To investigate this plastic behavior in TE semiconductors, we examined the evolution of twinning in several stages. By identifying the chemical bond change during twin nucleation, we explained how the shear load along the ð111Þ½112 direction induces the deformation twinning, yet retains structural integrity. By calculating the energy of slipping and twinning, we explained why deformation twinning, rather than slipping, is more plausible under shear loading along the ð111Þ½112 direction. These findings help provide improved understanding of deformation twinning in enhancing the plasticity of semiconductor with zinc blende structure.

Mechanical properties of InSb along different slip systems
First, we investigated the mechanical properties of flawless singlecrystal InSb through calculating the shear response along seven plausible slip systems. The structure of single-crystal InSb is shown in Fig. 1a, it has the F43M space group, with In and Sb located in the 4a (0, 0, 0) and 4c (0.25, 0.25, 0.25) sites, respectively. The lattice parameter from MD simulations is a = 6.479 Å and In-Sb bond length is 2.80 Å at room temperature (300 K), in agreement with experiment (a = 6.476 Å) 37 . The shear stress-strain curves are displayed in Fig. 1b, where the ð111Þ½112 slip system and its opposite slip system ð111Þ½112 have the lowest mechanical strength of 4.6 GPa at shear strain of 0.26, making them the most plausible slip systems. Interestingly, beyond the maximum stress point, we found that the ð111Þ½112 slip system shows an obvious plastic deformation, whereas the other six slip systems show typical brittle failure (as shown in Supplementary Fig. 1). In the following sections, we will focus on explaining the plastic deformation under a shear load along the ð111Þ½112 slip system at three different stages as shown in Fig. 1b:

Stress dropping stage
The structural deformation of the ð111Þ½112 slip system during the stress dropping stage is shown in Fig. 2. The atoms are colored according to their atomic shear strain 38,39 , so that the structure clearly divided into two parts: undeformed blue area and deformed green area. During the rapid stress drop from 4.6 to 1.5 GPa, an obvious three-layered stack gradually forms (stack is judged by structure analysis of dislocation analysis (DXA) 40,41 in open visualization tool (OVITO) 42 ), as shown in Fig. 2a-c. The atomic shear strain in the three-layered stacking increases rapidly with the color changing from blue to green. In the first stack evolution ( Fig. 2d-f), as the shear strain increases from 0.2625 to 0.2690, the 1 6 <112> partial dislocations gradually diffuse (the dislocation type is judged by DXA 40,41 in OVITO 42 ). In the second stack evolution (Fig. 2g-i), the 1 6 <112> partial dislocations nucleate near the edge area of dislocation lines in the first stack (as shown in Fig. 2d). In the third stack evolution ( Fig. 2j-l), the 1 6 <112> partial dislocations nucleate near the edge area of dislocation lines in the second stack (as shown in Fig. 2h). The 1 6 <112> partial dislocations occur from the dislocation reaction in the face centered cubic (fcc) structure 43 : The relationship of these dislocations is also shown in Fig. 2e. Thus, the three-layered stacking evaluation under shear load along the ð111Þ½112 slip system leads to nucleation of deformation twinning, which suddenly reduces the shear stress.

Structural yielding stage
The structural deformation of ð111Þ½112 slip system at structural yielding stage is shown in Fig. 3. With shear strain increasing from 0.269 to 0.30 (from Fig. 2c to Fig. 3a), the stacking ratio (the volume increment of the stacked area within unit strain increment) reduces and the shear stress drop starts to slow down (point B in Fig. 1b).
As the shear strain increases from 0.30 to 0.70 (Fig. 3a-e), stackinginduced deformation twins (green parts) gradually grow with crystal rotation (according to the red and yellow lines), resulting in the texture hardening, corresponding to the stage of C-D in Fig.  1b. Accompanied by the crystal rotation in the deformation twinning, the structure exhibits the yielding stage.

Structural strengthening and failure stage
The structural deformation of ð111Þ½112 slip system at structural strengthening and failure stage is shown in Fig. 4. When the shear strain increases to 0.73 (Fig. 4a), the system only contains two-layer matrix. At this moment, the matrix can hardly transform to the deformation twins, but rather it will start to slip along the [112] direction with the increasing shear strain, as shown in Fig. 4b. This leads to a slight stress drop (point D in Fig. 1b). However, the slip direction ([112]) is not parallel to the shear direction (½112) and the slip-induced dislocation movement is hindered at TBs, as shown in Fig. 4c. Thus, the structure retains integrity, resulting in structural strengthening (stage D-E in Fig. 1b). The contribution of the dislocation hindering to the plastic deformation during the stage of structural strengthening (stage D-E in Fig. 1b) is less than that of the twining during the stage of structural yielding (stage B-D in Fig. 1b); thus, the twinning is the dominant mechanism for the plastic deformation in InSb. At 0.88 shear strain (Fig. 4c), the structure starts to slip at TBs, whereas the matrix orientation mainly changes to its symmetry orientation (ð111Þ½112), indicating that the ð111Þ½112 slip system is more likely to slip than the ð111Þ½112 slip system. By changing the crystal orientation, the deformation twinning contributes to the activation of slipping. Then at 0.89 shear strain (Fig. 4d), the slippages rapidly propagate along ½112 direction (origin orientation) and cross the nanotwins, leading to the structural failure (stage E-F in Fig. 1b).

Atomistic explanation of deformation twinning evaluation
To investigate the deformation twinning evaluation of InSb under shear loading along the ð111Þ½112 slip system, we extract the local atomistic configurations in Fig. 2c to illustrate the chemical bonding responses against shear strain, as shown in Fig. 5. As the shear strain increases from 0.2698 to 0.2702 (Fig. 5a-d), In2 and Sb2 move relatively along the ð111Þ½112 direction with In2-Sb2 bond rotating, which gradually changes the matrix orientation to its symmetry orientation (ð111Þ½112). Under shear load along the ð111Þ½112 direction, In2 and Sb1 move along the ð111Þ½112 direction relatively with a decreased distance. Meanwhile, In2 and Sb3 move along the ð111Þ½112 direction relatively with increased distance, until the In2-Sb3 bond breaks at 0.2701 shear strain (Fig.  5c). Beyond breaking the In2-Sb3 bond, Sb1 immediately attracts its nearby atom In2 to form the In2-Sb1 bond at 0.2702 shear strain (Fig. 5d), as they both lack one covalent bond to reach covalent bond saturation. As the In2-Sb1 bond forms, the In2-Sb2 bond rotates to keep the neighbor structure of In2 intact. We refer to this bond-breaking formation process between the neighboring atom as "catching bond," which has been observed in our previous DFT calculations 34 . However, shearing along the symmetry direction of the ð111Þ½112, the ð111Þ½112 direction cannot cause such decreased distance between In2-Sb1 or In3-Sb2, let alone the "catching bond." The "catching bond" can maintain the structural integrity during the extremely crystal rotation in the twin nucleation (stage A-B in Fig. 1b) and further resulting in the twin growth (stage B-C-D in Fig. 1b).
To further investigate the origin of deformation twinning nucleation in InSb, we calculated the bond lengths and angles during the elastic stage, as shown in Fig. 6. Three bond types and three angle types are classified as shown in Fig. 6a. The elongation of the bond 1 is less than that of bond 2 at the elastic stage, but bond 1, rather than bond 2, breaks at the beginning of the twin nucleation (point A), as shown in Fig. 6b. This indicates that the bond stretching is not the main reason for breaking bond 1. The bond angle 1, consisting of bond 1 and 2, is bent significantly at the elastic stage from 109.4°to 124.4°. Similarly, the angle 3 consisting of bond 1 and 3 is remarkably bent as well, from 109.4°t o 100.5°, although the angle 2 consisting of bond 2 and 3 is slightly bent from 109.4°to 106.6°. These angle changes are shown in Fig. 6c. This indicates that the directionality of covalent bond 1 changes significantly against shear stress, which leads to breaking bond 1 at point A. This data analysis explains the breaking of the In2-Sb3 bond in Fig. 5. The breaking of bond 1 leads to twinning on the ð111Þ glide-set plane, which is a nonshear plane, as shown in Fig. 6a. Therefore, by breaking bond 1, the shearing along the ð111Þ½112 direction leads to twinning along the ð111Þ½112 non-shear direction. This atomic explanation explains why the deformation twinning nucleates along the ð111Þ½112 direction, as shown in Fig. 2a.

Energy explanation of competition between twin and slip
The double fcc lattice structure has two different kinds of closepacked (111) planes: glide-set plane (narrowly spaced atomic layer) and shuffle-set plane (widely spaced atomic layer) 44 . Slipping on these two different planes is induced by different dislocations. On the shuffle-set plane, slipping can be induced by a perfect dislocation along the f111g<110> direction, resulting in cleavage, as observed in InSb under shear along different direction ( Supplementary Fig. 1). However, this slipping cannot be induced by leading-trailing partial dislocations, because partial dislocations on the shuffle-set plane are unstable. On the glide-set plane, the slipping cannot be induced directly by a perfect dislocation along the f111g<110> direction, because the density of covalent bonds on the glide-set plane is much higher than that on the shuffle-set plane, making it more likely to be induced by the leading-trailing partial dislocations. The twinning can be formed only on the glide-set plane, because the partial dislocation exits only on the glide-set plane. These three kinds of shear deformation modes (the slipping on the shuffle-set plane and the glide-set plane, and the twinning on the glide-set plane) are illustrated in Supplementary Fig. 2.
To discuss the competitive mechanism between these three shear deformation modes, we employed DFT to calculate the cleavage energy on the shuffle-set plane and the generalized stacking fault energy (GSFE) on the glide-set plane, as shown in Fig. 7a. The details of the DFT calculations are described in "Methods." We observed that the cleavage energy γ c (1440.0 mJ m −2 ) is higher than the unstable stacking energy γ ust (778.9 mJ m −2 ) (Fig. 7a), indicating that the slipping-induced cleavage on the shuffle-set plane is more difficult to occur than the stacking on the glide-set plane, as well as this stacking-induced slipping and twinning. By further observing the structures of γ ust points, the atoms near the stack plane and its neighboring atoms form approximate equilateral triangles, whereas the distance of In1-Sb1 (2.76 Å) is only 4.2% shorter than the In-Sb bond in single crystal (2.88 Å), as shown in Fig. 7b. This indicates that the structure retains its integrity at this point. On the other hand, the unstable twinning energy γ ut (798.5 mJ m −2 ) is very close to the unstable stacking energy γ ust (778.9 mJ m −2 ) (Fig. 7a), indicating that, after the stacking on the glide-set plane, twinning with a second layer stack is easy to substitute the slipping with a trailing partial dislocation 45 . In addition, both the stacking fault energy γ sf (36.9 mJ m −2 ) and the extrinsic stacking fault energy 2γ t * (37.5 mJ m −2 ) are very low, indicating that the twinning transition structures, both one-layer stack and two-layer microtwin, are very Fig. 6 Bond lengths and angles with increasing shear strain. a Schematic diagram of the bonds and angles, and the glide-set and shuffle-set planes. b Bond lengths with increasing shear strain. c Bond angles with increasing shear strain. The double fcc lattice has two different kinds of slip planes on {111}: glide-set plane (narrowly spaced atomic layer) and shuffle-set plane (widely spaced atomic layer) 44 . In a, the glide-set plane and shuffle-set plane are represented by the green line and the yellow line, respectively. In b and c, point A in Fig. 1b is represented by black dashed lines to distinguish the elastic and plastic stage. Fig. 7 The cleavage energy and GSFE curves. a The cleavage energy and GSFE curves. b Optimized model at γ ust points. In a, the maximum energy of the stacking energy curve (black square), the twinning energy curve (red circle), and the cleavage energy curve (blue triangle) are the unstable stacking energy γ ust , the unstable twinning energy γ ut , and the cleavage energy γ c , respectively. The minimum energy in the middle of stacking energy curve (black square) is the stacking fault energy γ sf and the minimum energy at the end of twinning energy curve (red circle) is the extrinsic stacking fault energy 2γ t *. These definitions are consistent with ref. 46 . In b, the section of the atoms in the black dashed rectangle is displayed on the right. stable 46 . Therefore, the twinning on the glide-set plane can be easily induced by suitable shear load, such as shear load along the ð111Þ½112 direction. This shear can lead to the "catching bond" in case of the breakages, as we have explained before.
By analyzing the energy barrier of different shear processes, under the shear load along the ð111Þ½112 direction, the reason for generating twinning, rather than slipping, can be summarized as: at first, the slipping-induced cleavage on the shuffle-set plane is difficult to occur due to the high γ c , whereas the stacking on the glide-set plane may occur due to the low γ ust . Then, due to the twinning energy is close to the stacking energy, on the glide-set plane, the second stacking along the ½112 shear direction is more likely to be activated than the trailing dislocation along the ½211 direction or ½121 direction that causes the slipping. This finally results in the deformation twinning on the glide-set plane.

DISCUSSION
In summary, we identified the deformation twinning formation in TE InSb under shear loading along the ð111Þ½112 direction at 300 K, which is rarely observed in TE semiconductors, but a typical plastic deformation behavior in metal and alloy. By analyzing the chemical bond response, we found that the shear load along the ð111Þ½112 direction can induce deformation twinning by breaking the directionality of the In-Sb covalent bond to change its neighboring structure. Meanwhile, to achieve saturation of covalent bond under the new neighboring structure, the "catching bond" occurs, which further maintains the structural integrity. By analyzing the cleavage energy on the shuffle-set plane and GSFE curves on the glide-set plane, we found that the energy barrier of the slipping-induced cleavage on the shuffle-set plane is higher than that of the twinning on the glide-set plane, making deformation twinning is preferable. The deformation twinning with the crystal re-orientation leads to a large yielding stage, resulting in a maximum shear strain of~0.85. These findings should provide guidelines helpful to improve the plasticity of zinc blende structure semiconductors through TBs engineering.

MD simulations
For the MD simulations, we chose the Tersoff potential 47 to describe the interactions of atoms in InSb, whose format is shown below: where f R (r ij ) is the attractive function, f A (r ij ) is the repulse function, b ij represents the many-body term, and f c (r ij ) is the cutoff function. All parameters used here are from ref. 48 52 . The average lattice parameter a is~6.48 Å, which also coincides with experimental results a = 6.476 Å 37 . These verify the feasibility of this force field. Four structural models of InSb were built with size of~140 Å along all three dimensions as shown in Supplementary Fig. 3 and seven kinds of slip systems were considered for shear, which are listed in Supplementary  Table 1. All the shear MD simulations were performed using the Largescale Atomic Molecular Massively Parallel Simulator 53,54 . The equation of motion was integrated with a time step of 1 × 10 −3 ps. Periodic boundary conditions were applied in all three directions to eliminate the surface effects. All systems were relaxed under isothermal-isobaric ensemble with isotropic barostat set to 0 Pa and Nose-Hoover thermostat 55 at 300 K for DFT calculations DFT calculations were run on the Vienna Ab Initio Simulation Package 56-58 , using the generalized gradient approximation with the Perdew-Burke-Ernzerhof 59 exchange-correlation function. The energy tolerance and the force tolerance applied for the convergence were 1 × 10 −6 eV and 1 × 10 −2 eV Å −1 , respectively. An energy cutoff of 500 eV was applied for the geometry optimization. A 7 × 7 × 7 Monkhorst-Pack 60 k-point was used for the single-crystal InSb geometry optimization, whereas an 8 × 1 × 14 Monkhorst-Pack grid was used for cleavage energy and GSFE calculations. The lattice parameter from the DFT geometry optimization is a = 6.648 Å. For the other energy calculations, a single cell was created along the direction of ½112 ½111 ½110 with 12 atoms and expanded to a 1 × 4 × 1 supercell with 12 layers along the [111] direction. Periodic boundary conditions were applied along all three dimensions. For cleavage energy calculation, a vacuum buffer of 20 Å was inserted in a shuffle-set plane near the y-direction periodic boundary. The vacuum buffer is thick enough that the relative energy change between the flawless structure and the cleavage structure is independent of further increasing vacuum thickness. This relative energy change divided by the separated plane area is the cleavage energy γ c . For GSFE calculations, the model with vacuum buffer of 20 Å was further displaced rigidly on different glide-set plane. For stacking process, the model was displaced rigidly on a glide-set plane in the middle of the model along the ð111Þ½112 direction. For twinning process, the one-layer stacked model was further displaced rigidly on another nearby glide-set plane along the ð111Þ½112 direction. The structures for the GSFE calculations were optimized only along the [111] direction. The relative energy change per unit area on shear plane under different displacements constitute the GSFE curves.

DATA AVAILABILITY
All necessary data generated or analyzed during this study are included in this published article and its Supplementary Information files. Extra data are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
Code sharing not applicable to this article, as no custom computer code or algorithm were applied during the current study.