Dislocation structure and mobility in the layered semiconductor InSe: A first-principles study

The structure and mobility of dislocations in the layered semiconductor InSe is studied within a multiscale approach based on generalized Peierls--Nabarro model with material-specific parametrization derived from first principles. The plasticity of InSe turns out to be attributed to peculiarities of the generalized stacking fault relief for the interlayer dislocation slips such as existence of the stacking fault with a very low energy and low energy barriers. Our results give a consistent microscopic explanation of recently observed [Science {\bf 369}, 542 (2020)] exceptional plasticity of InSe.


I. INTRODUCTION
Two-dimensional (2D) semiconductors have attracted the great interest due to their electrical and optical properties that make them promising candidates for future electronic and opto-electronic applications [1][2][3].
Indium selenide (InSe) is a typical representative of the layered van der Waals (vdW) semiconductors, which has received considerable attention in recent years because of the possibility to fabricate InSe structures down to monolayer thickness.Atomic-thick InSe crystals demonstrate high environmental stability and attractive electronic properties [4][5][6][7][8], appealing for practical applications.
The presence of a vdW gap between different layers makes the material easy to cleave, thus making it difficult to retain its structural integrity through plastic deformations.As a consequence, plastic deformability is highly unexpected for bulk vdW inorganic semiconductors at room temperature.However, recent studies demonstrate a number of counterexamples [9][10][11][12].In fact, the ability to sustain large deformations without fracture is an exceptional property of layered semiconductors [9,10], which makes them advantageous for applications compared to conventional organic semiconductors.Apart from vdW inorganic semiconductors, high plasticity is typical to layered carbides and nitrides, constituting a broad class of materials, the so-called MAXphases [13,14].
Mechanical properties of InSe are characterized by superior in-plane flexibility of individual layers [9].As it has been recently shown [10], layered phase of InSe demonstrates exceptional plasticity among other vdW inorganic semiconductors.This behavior can be attributed to the low resistance of InSe with respect to the interlayer gliding.The origin of such unusual behavior has been discussed in Ref. 10 in terms of the cohesion characteristics of vdW inorganic semiconductors, and anomalously low * a.rudenko@science.ru.nlYoung's modulus of individual layers [9].As suggested in Ref. 10, soft intralayer bonding between In and Se atoms may be responsible for rather low energy barrier to slipping between the {001} layers.At the same time, the microscopic mechanism of such an exceptional behavior is not yet understood, and a quantitative theory of the high plasticity observed in InSe is still missing.
It should be noted that a similar mechanism of plasticity is realized in MAX-phases [13], in which slipping of close-packed planes is governed by the motion of dislocations in the basal plane {001}, facilitated by weak interlayer bonding.Although MAX-phases are characterized by the localization of deformations along with the formation of kink bands, the dislocation motion plays a key role in this process.At the same time, the criteria determining the regime of plasticity typical to true bulk materials remain applicable even for layered MAX-phases [15].
In this study, we employ a multiscale theoretical approach to analyze the structure and mobility of dislocations as well as brittle vs. ductile behavior of InSe.We start from an ab initio description of the interatomic interactions to retrieve the generalized stacking fault (GSF) energy surface, which is further used to calculate the dislocation properties within the generalized Peierls-Nabarro (PN) model.
We show that the GSF energy surface of InSe is characterized by a very low energy barrier and existence of the stacking fault with almost zero energy.This leads to a large splitting of the dislocation core and exponential suppression of the Peierls stress.This is sufficient to provide a ductile behavior of a single-crystalline layered InSe.

II. THEORETICAL APPROACH AND COMPUTATIONAL DETAILS
Peierls-Nabarro (PN) model provides a basis of our understanding of dislocation properties which determines strength and plasticity of materials [16].A combination of ab initio calculations and the PN model [17,18]  from accurate electronic structure calculations of specific materials.The key point is to calculate the energetics of GSF from the first principles.This approach has been successfully used to analyze of the structure and mobility of dislocations in metals and alloys in a large number of realistic situations [18][19][20][21][22][23][24][25].In the framework of the generalized PN model, the total energy of the crystal with an infinite straight dislocation can be written as [16] where E el is the elastic energy term, and E mis is the nonlinear atomistic misfit energy term, both depending on relative displacements u(x) of atom rows below and above the cut plane, with x being the distance from the dislocation axis in the slip plane The misfit energy is defined by a periodic energy profile Φ(u) where n is the number of atomic rows, h is the distance between them, and l is the position of the dislocation center in the lattice.The function Φ(u) is usually approximated by the GSF energy or the γ-surface [26], which is associated with shearing half of the crystal by a vector u in the slip plane.Within the PN model, the displacement field u(x) is determined from the balance condition for the elastic stress and the atomistic forces originating from the interaction between two half-spaces of the crystal (see Ref. 16 for details).Within the PN model, this balance condition is determined by the minimum of the functional Eq. ( 1), δE tot /δu(x) = 0.In order to estimate the resistance of InSe with respect to shear deformation, we follow the GSF scheme and construct a series of different stacking configurations of bilayer InSe corresponding to a shearing of one layer relative to the other.We then calculate the energy profiles Φ[u(x)] along two sliding pathways: [100] (zigzag) and [120] (armchair), as it is shown in Fig. 1.The energies of various structural configurations were calculated using density functional theory.The calculations were carried out using the projected augmented wave formalism [27] as implemented in the Vienna ab initio simulation package (vasp) [28,29].The exchange-correlation effects were taken into account by using the dispersioncorrected nonlocal functionals, vdW-DF [30] and opt-B88-vdW [31].For comparative purposes, we also used the generalized gradient approximation (GGA) [32] and Hartree-Fock approximation (HF) in the form of the exact exchange energy [33] to calculate the interlayer energies.We used hard pseudopotentials, which include 4s, 4p, 4d, 5s, 5p valence electrons for In, and 2s, 2p, 3d, 4s, 4p electrons for Se.A 900 eV energy cutoff for the planewaves and a convergence threshold of 10 −8 eV were used.The Brillouin zone was sampled by a (16×16) k-point mesh.In order to avoid interactions between the supercell images in the nonperiodic direction, a 30 Å thick vacuum slab was added in the direction normal to the InSe sheet.In the calculations of slipping energies, the atomic positions were relaxed in the vertical direction until the residual forces were less than 10 −3 eV/ Å.In the calculations with variable interlayer separation, one of the atoms in each layer was fixed.The lateral lattice constant was fixed to the equilibrium value a = 4.0 Å.

A. Interlayer interaction
In this study, we consider InSe crystal, which is composed of vertically stacked weakly interacting InSe layers.Each layer adopts a crystal structure with two vertically displaced 2D buckled honeycombs (Fig. 1).Symmetry of each layer is characterized by the point group D 3h .In accordance with Ref. 10, we assume that InSe crystallizes in the β-InSe phase with the space group P6 3 /mmc (No. 193), in which In and Se atoms in the adjacent layers are vertically aligned.
The calculated separation energy E sep as a function of the distance between the two crystal layers cleaved along the plane {001} is presented in Fig. 2. The results in Fig. 2 are shown for two typical dispersion-corrected exchange-correlation functionals, vdW-DF and optB88-vdW, in comparison with the standard generalizedgradient approximation (GGA) and Hartree-Fock (HF) calculations.Both GGA and HF binding energy curves demonstrate a shallow minimum, indicating the absence of any sizeable binding between the layers.This demonstrates that the interlayer binding in InSe is mainly of the vdW origin.On the contrary, nonlocal vdWfunctionals predict a considerable interlayer interaction, typical to layered materials.The resulting cleavage energy G c = −E sep (0) amounts to 0.26 and 0.29 J/m 2 for vdW-DF and optB88-vdW functionals, respectively.Both functionals result in similar E sep (d) behaviour.In the discussion that follows, we consider the results obtained using the vdW-DF functional only.

B. Generalized stacking fault energy
The energy profiles Φ(u) calculated along two sliding pathways, corresponding to the [120] (armchair) and [010] (zigzag) cross-sections of the GSF surface, are presented in Fig. 3.One can see that the sliding along [120] between the AA and AB1 configurations is the most en-ergetically favorable.The corresponding energy barrier is around ∼14 mJ/m 2 .Interestingly, the energies of the initial (AA) and final (AB1) configurations are very close to each other.This means that the AA-AB1 sliding results in the formation of a stacking fault with nearly zero energy.This is due to the fact that the main contribution to the energy upon sliding of InSe layers is given by the adjacent Se atoms, which have equivalent positions in both AA and AB1 stacking configurations.As a consequence, the most plausible mechanism of plasticity in InSe is the motion of partial dislocations with the Burgers vector b par =AA-AB1 in the {001} plane.
It is convenient to perform a Fourier expansion of the GSF energy over the reciprocal vectors of a 2D lattice in the {001} plane, which can be written as The second and third terms in Eq. ( 3) describe contributions from two triangular lattices shifted by the vector u 0 = (1, 0, 0) relative to each other, which ensures a hexagonal symmetry of the function Φ(u) in the {001} plane.The coefficients A k have been determined by fitting to the results of ab initio calculations presented in Fig. 3.The expansion length has been set to k max = 3, which is sufficient to obtain a reliable approximation of the calculated points.
The resulting approximation to the entire GFS surface is presented in Fig. 4, with the individual crosssections along the [010] (zigzag) and [120] (armchair) directions shown by the solid lines in Fig. 3. From Fig. 4 one can see that a direct translation by vector c 2 = (0, √ 3, 0) along [010], which was considered in Ref. 10, is not the most energetically favorable slipping path.Instead, the minimum energy path corresponds to a decomposition of c 2 into two partial translations as (0, /2, 0).Since the partial translation (1/2, √ 3/2, 0) leads to the formation of a stacking fault with nearly zero energy (see Fig. 3), we expect that the motion of partial dislocations with the Burgers vector b par = (1/2, √ 3/2, 0) in the {001} plane is the main mechanism of plastic deformations in InSe.

C. Dislocation structure and mobility
We now analyze the structure and mobility of partial dislocations with the Burgers vectors b par = (1/2, √ 3/2, 0) in the {001} plane.To determine the dislocation core structure, we employ a modified PN model as it is discussed in Refs.17 and 18.We perform a minimization of the total energy functional, Eq. ( 1), with a discrete representation of the misfit energy, Eq. ( 2), using trial functions, u(x), defined from the Laurent expansion [18,34] of their derivatives where µ is the shear modulus, D is a parameter equal to 1 for a screw dislocation, and to 1/(1-ν) for an edge dislocation with ν being the Poisson ratio, and ζ determines the width of the dislocation core.It is sufficient to use p = 2 for the case under consideration.Using Eq. ( 4), the expressions for displacements u(x) and restoring forces ∂E mis [u(x)]/∂u can be represented in the conventional form as where θ = 2 arctan(x/ζ), ζ = αω 0 , ω 0 = dD/2 is the width of the dislocation core defined within the original PN model [16], and d is interplane distance.The parameters α and β have been determined from fitting Eq. ( 5) to the calculated values of E mis (u) in the interval AA-AB1 (Fig. 3).The resulting distribution of displacements u(x) and ρ(x) for screw dislocation are shown in Fig. 5.One can see that the dislocation under consideration has a rather wide core, facilitating gliding processes.
Let us estimate the Peierls stress, σ P , which determines the lattice resistance to the dislocation motion.Within the generalized PN model considered in this work, σ P can be represented as [17] where S 0 is the area per atom in the slip plane.The integral Eq. ( 6) can be approximated using the method of steepest descent, yielding Using the calculated GSF energy, we obtain σ P /µ ≈ 4.7×10 −5 .Therefore, the dislocations under consideration are characterized by considerably low Peierls stress and, as a consequence, high mobility typical to ductile metals (see Ref. 19).

D. Cleavage fracture condition
Rice and Thomson [35] suggested an approach to characterise whether a material is ductile or brittle by considering the balance of the crack opening and dislocation emission from the crack tip, processes that characterize the microscopic properties of materials.The emission of a dislocation can cause crack blunting, lowering the strain energy release rate and thus leading to ductile behavior.According to the Rice-Thompson criteria, brittle failure is realized when µb/γ > 7.5 − 10, where γ = G c /2 is the surface energy.The applicability of this criterion for the ranking of materials with different crystal structures has been widely discussed (see Ref.  From the results presented in Fig. 2, we have γ = 0.13 J/m 2 , which is considerably smaller compared to typical hexagonal close-packed metals with γ = 1-2 J/m 2 [36].However, because of the low shear modulus in the {001} plane µ = 1.75 GPa, we obtain µb/γ ≈ 3.2.Such a small ratio µb/γ indicates that brittle fracture along the {001} plane is highly unlikely. Another important parameter of fracture is the ratio γ/γ us [37], where γ us is the so-called unstable stacking fault energy, which characterizes a barrier for sliding.In the case under consideration γ/γ us has its maximum along the AA -AB1 path in Fig. 3(b).An estimation gives rather high ratio γ/γ us ≈ 8.3, which is similar to the values typical for ductile face-centered cubic metals Cu and Ni [37].It should be noted that the plasticity of layered crystals was discussed in Ref. 10 in terms of the relation of slipping and cleavage energies, which is similar to the γ/γ us ratio.
Strictly speaking, the criteria discussed above should be applied with caution to layered crystals such as InSe or MAX phases (see discussion in Ref. 15).For a correct description of plasticity in such materials, the dislocation mobility must be taken into account.Since we have found that the intrinsic dislocation mobility in InSe is rather high, one can expect that the Rice-Thomson criterion correctly predicts the ductile behavior in this compound.

IV. DISCUSSION AND CONCLUSION
To elucidate the origin of ductile behavior of the 2D layered compound InSe an ab initio based analysis of the dislocation structure, mobility, and cleavage decohesion was performed.We found that the AA-AB1 sliding in the 120 (armchair) direction is characterized by a very low energy barrier and results in the formation of stacking faults with nearly zero energy (see Fig. 3).As a result, partial dislocations with the Burgers vector b par along 120 can easily glide in the basis plane.We note that a similar behavior of the GSF energy surface was reported earlier for graphite [38] where high mobility of dislocations in the basis plane is also to be expected [39].
Low energy barrier for the AA-AB1 sliding results in rather small stacking fault energy γ us , resulting in the large ratio γ/γ us , indicating a ductile character of InSe.This allows us to conclude about the intrinsic plasticity of InSe in the situation when the deformation is applied along the basal plane; a similar conclusion can be drawn considering the ratio µb/γ that appears in the Rice-Thomson criterion of the fracture condition.Thus, monocrystalline InSe is microscopically ductile material.This behavior does not appear surprising for layered materials provided there is a preferred gliding direction (see discussion in Ref. 15).At the same time, this does not guarantee macroscopic plasticity of polycrystalline samples, since only one slipping plane cannot provide an arbitrary shape change of material.According to the von Mises criterion, five independent slip systems are required to provide an arbitrary deformation [16] which is highly improbable for layered materials with geometrically determined special directions.Therefore our analysis is sufficient only to explain ductility of single-crystalline materials under special kind of deformation, exactly what

FIG. 1 .
FIG. 1. Schematic representation of the InSe crystal structure (AA-stacking) shown in two projections.The arrows in the (001) plane denote [100] (zigzag) and [120] (armchair) nonequivalent directions in which the slipping of one layer is considered relative to the other.

FIG. 2 .
FIG. 2. The energy of bilayer InSe calculated as a function of the interlayer separation d using different exchangecorrelation functionals.d = 0 corresponds to the equilibrium distance.Gc = −Esep(0) denotes the cleavage energy.

FIG. 3 .
FIG. 3. Cross-sections of the GSF energy along (a) [100] (zigzag) and (b) [120] (armchair) directions in the (001) plane.Points correspond to the results of ab initio calculations within the nonlocal vdW-DF functional.Solid lines show the approximation obtained by fitting with Eq. (3).Arrows denote symmetric configurations corresponding to AA-and AB-types of stacking of the two layers. 15).

FIG. 4 .
FIG. 4. The GSF surface obtained from the Fourier expansion Eq. (3) with parameters fitted to the calculated GSF crosssections, presented in Fig. 3.The black arrows correspond to [010] (zigzag) and [120] (armchair) directions.The red lines correspond to the Burgers vector of partial dislocations discussed in the text.