Stability of Strained Stanene Compared to That of Graphene

Stanene, composed of tin atoms, is a member of 2D-Xenes, two-dimensional single element materials. The properties of the stanene can be changed and improved by applying deformation, and it is important to know the range of in-plane deformation that the stanene can withstand. Using the Tersoff interatomic potential for calculation of phonon frequencies, the range of stability of planar stanene under uniform in-plane deformation is analyzed and compared with the known data for graphene. Unlike atomically flat graphene, stanene has a certain thickness (buckling height). It is shown that as the tensile strain increases, the thickness of the buckled stanene decreases, and when a certain tensile strain is reached, the stanene becomes absolutely flat, like graphene. Postcritical behaviour of stanene depends on the type of applied strain: critical tensile strain leads to breaking of interatomic bonds and critical in-plane compressive strain leads to rippling of stanene. It is demonstrated that application of shear strain reduces the range of stability of stanene. The existence of two energetically equivalent states of stanene is shown, and consequently, the possibility of the formation of domains separated by domain walls in the stanene is predicted.


Introduction
In addition to graphene [1], other single-element two-dimensional materials with a similar hexagonal structure have been obtained and studied. These are group IV elements, such as silicene, germanene, and stanene (2D-Xenes). In the past few years, they have attracted the attention of researchers due to a combination of electronic and optical properties that can be used in nanotechnology [2][3][4][5].
Unlike atomically flat graphene, all two-dimensional group IV materials (silicene, germanene, and stanene) are buckled and have a certain thickness (buckling height [6,7]). The free-standing stanene has a thickness of h = 0.85 ÷ 0.90 Å according to density functional theory (DFT) calculations [8,9]. Stanene grown on a substrate may be thinner or thicker than a free-standing stanene.
The first successful attempt to synthesize stanene was made on the Bi 2 Te 3 (111) surface in 2015 [10] and later on many other substrates using the epitaxial methods of growth [4,11]. The thickness of stanene on Bi 2 Te 3 (111) is reported to be h = 1.2 Å [10], which is larger than that of free-standing stanene. From the DFT calculations, the structure of the stanene with dumbbell units with a thickness of 3.41 Å was predicted [12]. Large area stanene can be epitaxially grown on Ag(111) with h = 0.12 Å [13,14]. A nearly flat stanene-like honeycomb structure with three Sn atoms per unit cell was epitaxially grown on an Au(111) surface [15], and ultraflat stanene (h = 0) was epitaxially grown on a Cu(111) surface [16]. Strain-free Elastic properties and crack propagation in single and bilayer stanene were analyzed based on molecular dynamics simulations using Tersoff bond order potential [41]. At the strain rate of 10 8 s −1 , failure of stanene was observed at about 38% strain, and 15 GPs stress for tension along the armchair direction and lower critical values were found for tension along the zigzag direction (33% strain and 14 GPa stress). A comparative analysis of the mechanical properties of 2D-Xenes suggests that stanene is the weakest and softest member of this family [42].
Formulas for the shear modulus of hexagonal nanostructures were derived based on a physically motivated model [43]. In the works [44][45][46], concrete examples of materials with unusual properties are considered, and the origin of these properties at the discrete level is described. Electronic properties of bilayer stanene can be improved via modulating the stacking order and angle of the bilayers [47].
In order to take full advantage of the possibilities of elastic strain engineering, it is necessary to know the stability region of the stanene in the space of the strain tensor components ε xx , ε yy , and ε xy . For graphene, this problem was solved in [48] using the potentials developed by Savin [49]. This problem is difficult to solve experimentally because it is not easy to apply high enough stress to break a defect-free stanene. In such circumstances, theoretical studies make it possible to shed light on many properties of twodimensional structures, the study of which by experimental methods is difficult. However, theoretical studies also encounter certain difficulties. For example, ab initio methods, due to the high demands on computer resources, consider a small number of atoms in short periods of time [6,[50][51][52], and the molecular dynamics method gives results that depend on the interatomic potentials [53][54][55][56]. Despite this, modeling allows at least a qualitative characterization of many important properties of nanomaterials.
In this work, we aim to calculate the range of stability of plane stanene modelled with the Tersoff potential in the space of the strain tensor components ε xx , ε yy , and ε xy .

Materials and Methods
Molecular dynamics is used to simulate the stability of stanene.

Structure of Stanene
A free-standing, infinitely large, defect-free, single-sheet stanene is considered. According to theoretical and experimental data [57][58][59][60][61], undeformed, free-standing stanene has a hexagonal buckled structure. A primitive translational cell of stanene defined by the translation vectors a 1 and a 2 , is shown in Figure 1a projected onto the XY plane. The X and Y axes of the Cartesian coordinate system are oriented in the zigzag and armchair directions, respectively. Atoms 1 and 2 represent two triangular sublattices colored black and red; generally speaking, they have different Z-coordinates. Figure 1b shows three valence bonds, 1, 2 and 3, and six valence angles, ϕ i , and τ i , i = 1, 2, 3. These structural parameters will be analyzed for homogeneously deformed stanene. Figure 1c shows the side views of buckled stanene having thickness h projected onto the XZ and YZ planes. In Figure 1d, similar side views are shown for the absolutely flat stanene.

Structure Generation and Homogeneous Deformation Application
The structure of stanene is generated by translation of a pair of atoms by the vectors a 1 = a(1, 0, 0) and a 2 = (a/2)(1, √ 3, 0), where a is the lattice parameter (see Figure 1a). The length of the translation vectors |a 1 | = |a 2 | = a can be obtained by projecting the equilibrium bond length ρ 0 onto the XY plane. The radius-vectors of Sn atoms in unstrained structure are R 0 m,n,k = ma 1 + na 2 + s k , where the integers m, n define the primitive cell number, and k = 1, 2 is the number of atoms in a primitive translational cell. The vectors s k , k = 1, 2, in unstrained stanene are taken as s 1 = (0, 0, 0), s 2 = (a/2, a/2 √ 3, h), where h is the structure thickness depicted in Figure 1c.
Stanene is considered under homogeneous plain stress conditions achieved by application of in-plane strain with the components ε xx , ε yy , ε xy . The equilibrium positions of atoms in homogeneously strained stanene are R s m,n,k = mp 1 + np 2 + q k , where vectors p 1 = a 1 + a 1 T, p 2 = a 2 + a 2 T are defined through the matrix Vectors q k describe the position of atoms in a primitive cell after deformation. By the rigid shift of the structure in space, one can always take q 1 = (0, 0, 0). On the other hand, the equilibrium position of atom 2 in the primitive cell, q 2 , must be obtained by relaxation of the structure. As a result of relaxation, the equilibrium thickness of the deformed structure will also be determined.
The relaxation is performed by the steepest descent method using self-made code written in the C++ algorithmic language. Relaxation stops when the value of the maximum force acting on the atoms becomes less than 10 −8 eV/Å.

Tersoff Potential
The many-body Tersoff potential [62] is used to describe interactions between tin atoms in stanene. In Tersoff's models, the potential energy of the system is taken as where r ij = |R i − R j | is the distance between atoms i, and j (R i is the radius-vector of atom i). The function f c (r ij ) limits the range of interaction to the nearest neighbors. This function is taken in the form where R and D are the constants describing the cut-off radius of the potential. The potential energy (2) includes the embeddings of repulsive and attractive forces, where A, B, λ 1 , and λ 2 are the potential parameters describing the dimer strength and Pauling constants [63,64]. The component a ij for the repulsive force is taken equal to 1 [65]. The b ij term in (2) describes the weakening of the bond between atoms i and j due to the presence of another bond between atoms i and k. This term describes the energy of valence angles, Parameters β, N, λ 3 , c, d, and δ are adjustable, while M = 3.0 and γ = 1.0 are constants.
In this work, a more symmetrical form of the potential is used in which b ij = b ij +b ji 2 . Parameters of the Tersoff potential are listed in Table 1. The Hamiltonian of the system (total energy of the computational cell is) where U s is the potential energy of the system defined by Equation (2), and the second term on the right-hand side of the Hamiltonian gives the kinetic energy (M is the mass of the stanene atom). The equations of motion are derived from the Hamiltonian Equation (8) using the Hamilton's principle.
Let the computational cell include N x × N y primitive translational cells, each of which contains two atoms. Atoms are conveniently numbered with three indices m = 1, ..., N x , n = 1, ..., N y , and k = 1, 2, as described in Section 2.2. Table 1. Parameters of the Tersoff potential [66]. 638.6396

Check for Stability of Strained Stanene
For strained stanene, we check its stability by calculating the frequency of phonon modes. To do so we write down the equations of motion for the two atoms of the (m, n)th translational cell Then we linearize these equations and look for their solutions in the form In Equations (9) and (10), u m,n,k is the displacement vector of (m, n, k)th atom from equilibrium position, and it is assumed that |u m,n,k | ρ 0 ; M is the mass of tin atom; P m,n is the potential energy of (m, n)th primitive cell. Six components of the displacement vectors u m,n,1 and u m,n,2 define positions of two atoms in a primitive translational cell; e is the normalized six component eigenvector (e, e) = 1; i is an imaginary unit; 0 ≤ q x , q y < 2π are the components of the wavevector; ω is the phonon frequency.
Substitution of Equation (10) into linearized equations of motion results in the eigenproblem whose solution gives eigenvalues ω 2 j and corresponding eigenvectors e j , j = 1, ..., 6. A flat stanene is stable if all eigenvalues ω 2 j in the entire first Brillouin zone in the phonon spectrum are positive. The first Brillouin zone is scanned with ∆q x and ∆q y steps equal to π/100. The eigenvalue problem mentioned above is solved at each point, and if all ω 2 j , j = 1, ..., 6 are all positive, then the flat structure is stable. If the spectrum has negative eigenvalues, then the planar state is unstable.
The calculation of dispersion curves and stability analysis for stanene are performed using a home-made C++ code that implements the solution of the eigenvalue problem based on the iterative Jacobi algorithm for calculating the eigenvalues and eigenvectors of a real symmetric matrix. The iterations stop when the maximum off-diagonal element of the stiffness matrix becomes less than 10 −5 eV/Å 2 .

Results and Discussion
Let us present the numerical results obtained for the equilibrium structure of stanene and its stability.

Equilibrium State of Stanene
Before applying deformation, it is necessary to find the equilibrium state of the structure relative to which deformation will be applied. To do this, we change the lattice parameter a and independently the height of the structure h in order to find the state of minimum energy. The process of finding this state is illustrated in Figure 2. The equilibrium state of stanene is characterized by the equilibrium bond length ρ 0 ≈ 2.88 Å, structure height h = 0.99 Å, potential energy per atom U ≈ −3.31 eV, and the lattice parameter a = 2.70 Å. Recall that according to the DFT calculations [8,9], the free-standing stanene has somewhat smaller thickness of h = 0.85 ÷ 0.90 Å.
Let us estimate the Poisson's ratio of the undeformed mill. To do this, we calculate the energy of elastic deformation per atom for three deformation states: Let us estimate the Poisson's ratio of the undeformed stanene. To do this, we calculate the elastic strain energy per atom for the three strain states, state 1: ε xx = ε, ε yy = 0; state 2: ε xx = 0, ε yy = ε; and state 3: ε xx = ε yy = ε, with ε xy = 0 in all cases. The elastic energy density for the plain stress conditions can be expressed through the elastic constants and strain components as follows where E and ν are the Young modulus and Poisson's ratio. Substituting the values of strain for the three considered strain states and eliminating E, one finds The result for ε = 10 −3 is U 1 = 3.87 × 10 −6 , U 2 = 3.86 × 10 −6 , and U 3 = 9.61 × 10 −6 eV. The almost equal values of U 1 and U 2 confirm that the stanene is elastically isotropic at small strains. Substituting the values of energy into Equation (12), one finds ν ≈ 0.242 ≈ 0.245. DFT results give larger values of Poisson's ratio. In the work [67], ν zigzag = 0.42, ν armchair = 0.36, and in [68], ν = 0.395.
It should be noted that in molecular dynamics calculations, Poisson's ratio is often found with a larger uncertainty than Young's modulus. For example, in the review [69], dozens of papers on the elastic properties of graphene are analyzed, and the following ranges of calculated Young's modulus E = 312 ÷ 384 N/m and Poisson's ratio ν = 0.16 ÷ 0.46 are indicated.

Phonon Dispersion Curves
Our stability analysis is based on the calculation of the frequencies of small-amplitude oscillations of stanene. It is useful to calculate the phonon dispersion curves of undeformed stanene and compare them with the available results. Our results for the Tersoff potential are shown in Figure 3 and can be compared with the DFT-based results reported in [32], see Figure 1(c) of their work. There are three acoustic and three optic branches of the dispersion curves. We note a reasonably good agreement of the results. The maximum phonon frequency is 6.2 THz for the Tersoff potential and about 6.0 THz from ab initio calculations [32]. The ZA and ZO modes with out-of-plane displacements of atoms have the lowest acoustic and optical frequencies, respectively. Longitudinal phonons (LA and LO) have frequencies higher than transverse phonons (TA and TO, respectively) because the shear modulus is less than the tensile modulus.  Figure 4 shows the stability region of a planar stanene subjected to biaxial normal strain ε xx and ε yy at ε xy = 0. A flat stanene is stable inside the region shown by the solid line ABCDEFA and unstable outside this region. The postcritical behavior depends on which boundary of the stability region is crossed when leaving the stability region. When the CD line is crossed, atomic bonds of type 1 and 3 are broken, shown in Figure 1a. If the DE line is crossed, then the type 2 bonds break as they are most loaded in this case. If line ABC is crossed, then the stanene sheet experiences planar compression along the Y axis; hence ripples are formed that are elongated along the X axis as the planar configuration becomes unstable. Similarly, when crossing the AFE line, ripples oriented along the Y axis appear in the stanene compressed along the X axis.

Stability Region of Planar Stanene under In-Plane Tensile and Compressive Strain
Interestingly, within the region of stability, stanene exists in two different states. For relatively small tensile strain, stanene is buckled (h > 0), and it becomes absolutely flat (h = 0) under sufficiently large tensile strain. These two regions are separated by the dashed line ε yy = −1.03ε xx + 0.272 in Figure 4.
The change in the height of stanene in the range of stability of the planar stanene is shown in Figure 5a. It can be seen that h decreases with increasing tensile strain ε xx and ε yy and is equal to zero in the BCDEFB region.
Possible values of stanene height along the line ε xx = ε yy are shown in Figure 5b. The negative value of h means that the red and black sublattices shown in Figure 1 change the sign of the Z coordinate, Z → −Z. Suppose we have obtained an absolutely flat stanene with a sufficiently large tensile strain. Then a decrease in the tensile strain will lead to an approach to the bifurcation point in which the system has two energetically equivalent paths of development. In the region of nonzero h, atoms belonging to the red and black sublattices can have coordinates Z > 0 and Z < 0, respectively, or vice versa. In the first case, h > 0, and in the second case, h < 0. The red dashed line in Figure 5b shows an absolutely flat stanene, which is unstable after crossing the bifurcation point as the tensile strain decreases.
The fact of the existence of two energetically equivalent states implies the possibility of the existence in stanene of domains separated by domain walls.

Comparison of the Stability Regions of Stanene and Graphene
In the literature, the stability regions of graphene were described [48,70] using the Savin potential [49] and the modified Brenner potential [65,71], and here we compare our results for stanene with existing results for graphene. Note that the accuracy of various interatomic potentials for graphene have been assessed in the work [54]. It was proved that the Savin potential reproduces many properties of carbon structures well [72][73][74][75].
A comparison of the stability regions of stanene and graphene is shown in Figure 6. At first glance, the stability regions of stanene and graphene are similar, but there are some important differences between them. The most striking difference is that stanene exists in a buckled and absolutely flat configuration, while graphene is flat throughout its stability range. Another important difference is that buckled stanene allows very little negative strain, unlike graphene, which remains flat in a range of negative values of ε xx at large positive values of ε yy and in a range of negative values of ε yy for large positive ε xx . Note that absolutely flat stanene, like graphene, has regions of stability under negative strains. Based on this, it can be concluded that absolutely flat stanene has similar properties to graphene, while buckled stanene behaves differently. This can be explained by the difference in Poisson's ratio of buckled stanene and graphene.
The fact that the stanene is stable under greater biaxial tension with ε xx = ε yy can be explained by stanene stretching due to flattening.

Region of Stability of Planar Stanene and Graphene in the Presence of Shear Strain
The effect of shear strain on the stability regions of graphene and stanene is presented in Figure 7a and Figure 7b, respectively. The borders of stability regions of planar structures are shown as black, red, and blue lines for ε xy = 0.0, 0.15 and 0.3, respectively.
Both graphene and stanene stretched along the armchair direction (i.e., along the Y axis) retain greater stability when shear strain is applied. However, as ε xy increases, the stability region of graphene narrows along the strain component ε xx and shrinks towards the center, while the stability region of stanene reduces from the right and from the bottom.
It is also interesting to look at the structure of strained stanene; for this, two points, A and B, shown in Figure 7b were chosen. Point A corresponds to biaxial strain ε xx = ε yy = 0.118 and point B to uniaxial strain ε xx = 0, ε yy = 0.3. The effect of addition of shear strain ε xy on the stanene structure will be analyzed.
In Figure 8a, it can be observed that the structure at ε xy = 0 retains the hexagonal symmetry under biaxial tension. This symmetry is lost after application of shear stain, see (c) and (e). Figure 8 also shows that the height of the structure decreases with increasing shear strain, cf. (b) with (d) and (f). Parameters of the structures at points A and B marked in Figure 7b are presented in Table 2 for different values of shear strain. Valence bonds and valence angles are defined in Figure 1b. It can be seen that at point A for ε xy = 0 all valence angles are equal and all valence bonds have the same length. The equality disappears when a shear strain is applied. At point B at ε xy = 0 all three bond angles are different, and all three bonds have different lengths.
The penultimate column of Table 2 shows that the potential energy per one atom is greater for stanene at point B compared to point A. The last column of Table 2 shows the thickness of the stanene for all considered strain states. In point B, stanene is absolutely flat for all three considered strain states.

Conclusions and Future Work
By finding the vibration frequencies of tin atoms in equilibrium stanene, the stability region of a planar structure was found in the space of strain components ε xx , ε yy , and ε xy , see Figures 4 and 7b. The results were compared to the stability region of homogeneously strained graphene reported in [48,70], see Figures 6 and 7a. Stanene was simulated with the use of the Tersoff potential [62] and graphene with the Savin [49] and Brenner [65,71] potentials. The accuracy of the potentials was discussed in [54,[72][73][74][75]. The accuracy of the potential at not very large atomic displacements can be estimated by comparing the phonon dispersion curves shown in Figure 3 with those calculated on the basis of the DFT theory [32]. As already mentioned, the agreement between the dispersion curves calculated using the Tersoff potential and those calculated ab initio is quite good. The accuracy of the Tersoff potential at large values of strain may not be as good because the dispersion curves are calculated for the linearized equations of motion.
It was demonstrated that the stanene is buckled, has a nonzero thickness h under a small stretch, and becomes absolutely flat under a sufficiently large stretch, namely, for ε yy > −1.03ε xx + 0.272, see Figure 4. Stanene height as the function of ε xx and ε yy is presented in Figure 5. Actually, Figure 5b shows the bifurcation diagram where negative h means that the red and black sublattices shown in Figure 1 change the sign of the Z coordinate, Z → −Z. The fact of the existence of two energetically equivalent states (with positive and negative h) implies the possibility of the existence in stanene of domains separated by domain walls.
In a future work, the combined effect of deformation and temperature on the stability of stanene will be analyzed (in this work, the temperature effect was neglected). The structure and energy of domain walls, the existence of which was predicted in this work, will be investigated. Postcritical ripples in stanene under in-plane compression will be studied in accordance with the works [36,70,76]. Of interest are also the mechanical properties of stanene nanotubes and nanotube bundles similar to their carbon counterparts [77][78][79]. The propagation of a shock wave in stanene may have interesting features due to the nonzero thickness [80].

Data Availability Statement:
The data and the home-made codes of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.