Mechanical properties and stabilities of g-ZnS monolayers

Abstract

We investigate the mechanical properties and stabilities of the planar graphene-like zinc sulfide (g-ZnS) monolayers under various large strains using electronic structure calculations. The g-ZnS has a low in-plane stiffness, about 1/8 of that of graphene. The potential profiles and the stress–strain curves indicate that the free standing g-ZnS monolayers can sustain large tensile strains, up to 0.16, 0.22, and 0.19 for armchair, zigzag, and biaxial deformations, respectively. However, both the strength and flexibility are reduced compared to graphene-like zinc oxide monolayers. The third, fourth, and fifth order elastic constants are indispensable for accurate modeling of the mechanical properties under strains larger than 0.02, 0.04, and 0.08 respectively. The second order elastic constants, including in-plane stiffness, are predicted to monotonically increase with pressure, while the trend in the Poisson ratio is reversed. Our results imply that g-ZnS monolayers are mechanically stable under various large strains. The elastic limits provide a safe-guide for strain-engineering the g-ZnS based electronics.

---

1 Introduction

As one of the first discovered semiconductors, bulk zinc sulfide (ZnS) exists naturally in two polymorphs, zinc blende and wurtzite, with a band gap about 3.7 eV.^1,2 ZnS has important applications in ultraviolet light-emitting diodes and injection lasers,³ flat panel displays,⁴ sensors,⁵ and infrared optical windows,⁶ in addition to photocatalysis⁷ and triboluminescence.⁸ There are extensive interests in nanostructures of ZnS including nanotubes, nanowires, nanobelts, and nanosheets.^1,9 Experimentally, the width-to-thickness ratio dependence of the photoplastic effect of the ZnS nanobelts was reported.¹⁰ Nanosheets with a thickness of about 11 Å of wurtzite ZnS were obtained.¹¹ However, to the authors' best knowledge, the fabrication of mono-atomic-layer graphene-like hexagonal ZnS (g-ZnS) has not yet been reported, as contrasted to a few theoretical investigations.

Theoretical studies show that the ZnS monolayers have planar and buckled structures.¹² Once the out-of plane displacements of atoms are smaller than 0.03 Å, the planar structure (Fig. 1a) is reserved during the energy minimization. However, when the out-of-plane displacements of atoms are large, it will re-constructed to various non-planar structures with lower energies, about 0.08 eV/ZnS.¹² The planar g-ZnS monolayers have a direct band gap over 2.0 eV (ref. 12) and have a transition to indirect band gap when a biaxial strain of 0.028 is applied.¹³ The 2D bulk modulus of g-ZnS monolayers is predicted to be 23.94 N m⁻¹ using density functional theory calculations.¹³

Fig. 1 Geometry and structure. (a) Top-view and (b) side-view of the g-ZnS monolayer plane; (c) the unit cell with simulation box and atoms; (d) the orientations of the system.

Mechanical forces are much more tangible, reliable, and widely applicable than other stimuli to materials. Macroscopic mechanical stimuli are used for control of molecular systems and molecular machines, in addition to micro-device fabrication.¹⁴ The knowledge of mechanical properties of a material^15,16 is the base for development of mechanically responsive nanomaterials which will make a significant contribution to improvements in our lifestyles and our society.¹⁴ It is desirable to know the mechanical properties of g-ZnS monolayers for three reasons. First, it is critical in designing parts or structures regarding their practical applications of atomic 2D materials.¹⁷ For instance, applications in high-end bendable electronics require the integration of g-ZnS monolayers with stretchable polymer substrates. Second, strain engineering is a common and important approach to tailor the functional and structural properties including band gaps¹⁸ and radiation hardness¹⁹ of nanomaterials.^20,21 Finally, nanomaterials is vulnerable to strain with or without intent because of its monatomic thickness.^22,23 For instance, there are strains because of the mismatch of lattice constants. In addition, elastic limits and tolerances are desirable for functional manipulation.²⁴

The goal of this paper is to study the mechanical stabilities and properties of g-ZnS monolayers. We focus on the planar structure of g-ZnS in this study because it is more interested in applications as buckled structure can be suppressed under constraints such as substrate, sandwiching, or heterogeneous multilayers,^25–27 or applied tensile strains in their applications. We use density functional theory (DFT) calculations to model their responses under various mechanical loadings. The strain potential profile, the stress–strain relationships, the high order elastic constants, and the pressure dependent properties are studied. Our results are compared to that of graphene and graphene-like ZnO (g-ZnO). The results of graphene²⁸ and g-ZnO²⁹ were reported previously. They serve as references to g-ZnS monolayers because they have similar structures. In addition, graphene/ZnS and ZnS/ZnO composites were synthesized for superior photoelectric applications²⁶ and unusual photoluminescence emissions,²⁷ respectively. Therefore, a comparative study could be helpful in the material design of these composites.

Our results for the continuum formulation could also be useful in finite element modeling of the multi-scale calculations for mechanical properties of g-ZnS monolayers at the continuum level. The organization of this paper is as follows. Section 2 presents the computational details of DFT calculations. The results and analysis are in Section 3, followed by conclusions in Section 4.

2 Density functional theory calculations

We consider a conventional unit cell containing six atoms (3 zinc atoms and 3 sulfur atoms) with periodic boundary conditions (Fig. 1c) to model g-ZnS monolayers. This 6-atom conventional unit cell is chosen to capture the “soft mode”, which is a particular normal mode exhibiting an anomalous reduction in its characteristic frequency and leading to mechanical instability. This soft mode is a key factor in limiting the strength of monolayer materials and can only be captured in unit cells with hexagonal rings.³⁰

Following convention, we define the armchair direction as the nearest neighbor direction and the zigzag direction is perpendicular to it within the atomic plane formed by nearest neighbors (Fig. 1d). Periodic boundary conditions along the plane are applied. The total energies of the system, forces on each atom, stresses, and stress–strain relationships of g-ZnS monolayers under the desired deformation configurations are characterized via DFT. The calculations were carried out with the Vienna Ab initio Simulation Package (VASP)^31,32 which is based on the Kohn–Sham Density Functional Theory (KS-DFT)³³ with the generalized gradient approximations as parameterized by Perdew, Burke, and Ernzerhof (PBE) for exchange–correlation functions.³⁴ The electrons explicitly included in the calculations are the 3d¹⁰4s² electrons for zinc atoms and 3s²3p⁶ for sulfur atoms. The core electrons are replaced by the projector augmented wave (PAW) and pseudo-potential approach.³⁵ A plane-wave cutoff of 500 eV is used in all the calculations. The convergence of the total energy and forces is 10⁻⁵ eV and 10⁻³ eV Å⁻¹, respectively. The calculations are performed at zero temperature.

The atomic structures of all the deformed and undeformed configurations were obtained by fully relaxing a 6-atom-unit cell where all atoms were placed in one plane. The simulation invokes periodic boundary conditions for the two in-plane directions. There is a 15 Å thick vacuum region to reduce the inter-layer interaction to model the single layer system. The irreducible Brillouin Zone was sampled with a gamma-centered 24 × 24 × 1 k-mesh. Such a large k-mesh was used to reduce the numerical errors caused by the strain of the systems. To eliminate the artificial effect of the out-of-plane thickness of the simulation box on the stress, we use the second Piola–Kirchhoff stress³⁶ to express the 2D forces per length with units of N m⁻¹.

For a general deformation state, the number of independent components of the second, third, fourth, and fifth order elastic tensors are 21, 56, 126, and 252 respectively. However, there are only fourteen independent elastic constants that need to be explicitly considered due to the symmetries of the atomic lattice point group D_6h which consists of a sixfold rotational axis and six mirror planes.³⁶ The fourteen independent elastic constants of g-ZnS monolayers are determined by least-squares fit to the stress–strain results from DFT calculations in two steps, detailed in our previous work,³⁶ which have been well used to explore the mechanical properties of 2D materials.^{28,29,37–46} A brief introduction is that, in the first step, we use a least-squares fit of five stress–strain responses. Five relationships between stress and strain are necessary because there are five independent fifth-order elastic constants (FFOEC). We obtain the stress–strain relationships by simulating the following deformation states: uniaxial strain in the zigzag direction (zigzag); uniaxial strain in the armchair direction (armchair); and equibiaxial strain (biaxial). From the first step, the components of the second-order elastic constants (SOEC), the third-order elastic constants (TOEC), and the fourth-order elastic constants (FOEC) are over-determined (i.e., the number of linearly independent variables are greater than the number of constrains), and the fifth-order elastic constants are well-determined (the number of linearly independent variables are equal to the number of constrains). Under such circumstances, the second step is needed: least-square solution to these over- and well-determined linear equations.

3 Results and analysis

3.1 Atomic structure

We first optimized the geometry of the 6-atom-unit cell. The relaxed structure shows that the six atoms are still coplanar, as shown in Fig. 1b. The Zn–S bond length is 2.248 Å (Fig. 1c), agreeing well with previous DFT studies.^12,13 The S–Zn–S angle is 120°. The lattice constant a = 3.89 Å, which is larger than the experimental measurement (3.75 Å) of wurtzite ZnS,⁴⁷ where zinc and sulfur atoms are in different planes. Our planar structure agrees with previous DFT calculations.^12,13 This strain-free structure is then set as the reference when mechanical strain is loaded.

3.2 Strain energy profile

When strain is applied, the system will be disturbed away from the equilibrium state. Since the configuration energy of the strain-free configuration is the minima of the potential well, any strain will increase the system's energy. By applying different amounts of strain along different directions, the potential well can be explored. Once the strain is applied, all the atoms of the system are allowed full freedom of motion. A quasi-Newton algorithm is used to relax all atoms into equilibrium positions within the deformed unit cell that yields the minimum total energy for the imposed strain state of the super cell.

Both compression and tension are considered with Lagrangian strains ranging from −0.1 to 0.3 with an increment of 0.01 in each step for all three deformation modes. It is important to include the compressive strains since they are believed to be the cause of the rippling of the free standing atomic sheet.⁴⁸ It was observed that a graphene sheet experiences biaxial compression after thermal annealing,⁴⁹ which could also happen with g-ZnS monolayers. Such an asymmetrical range was chosen due to the non-symmetric mechanical responses of material, as well as its mechanical instability,⁵⁰ to the compressive and the tensile strains.

We define the strain energy per atom E_s = (E_tot − E₀)/n, where E_tot is the total energy of the strained system, E₀ is the total energy of the strain-free system, and n = 6 is the number of atoms in the unit cell. This size-independent quantity is used for comparison between different systems. The E_s of the g-ZnS monolayer as a function of strain in uniaxial armchair, uniaxial zigzag, and biaxial deformation are plotted in Fig. 2a. The result of the potential profile of g-ZnS is compared with that of g-ZnO (Fig. 2b). E_s is anisotropic with strain direction. E_s is non-symmetrical for compression (η < 0) and tension (η > 0) for all three modes. This non-symmetry indicates the anharmonicity of the g-ZnS monolayer structures. The harmonic region where the E_s is a quadratic function of applied strain can be taken between −0.02 < η < 0.02. The stresses, derivatives of the strain energies, are linearly increasing with the increase of the applied strains in the harmonic region.

Fig. 2 Energy profile. The strain energy per atom under uniaxial strain in armchair and zigzag directions, and equibiaxial strains of g-ZnS (a), compared with g-ZnO (b).

The anharmonic region is the range of strain where the linear stress–strain relationship is invalid and higher order terms are not negligible. With even larger loading of strains, the systems will undergo irreversible structural changes, and the systems are in a plastic region where they may fail. The width of the stable region is the opening width of the potential energy well (Fig. 2). The opening width and depth of the potential energy well η_s could serve as a scale to quantify the flexibility and strength of a nanostructure, respectively. The average width of the stable regions of the three deformation modes (i.e., the opening width of the potential energy wells) is a reasonably good scale for the mechanical stabilities of the nano structures. As a result, from the point view of potential energy, we conclude that g-ZnS is mechanically stable. However, it is less stable than g-ZnO, because both the depth and width of the potential energy well are smaller than that of g-ZnO.

The ultimate strains are determined as the corresponding strain of the ultimate stress, which is the maxima of the stress–strain curve, as discussed in the following section. It is worth noting that in general the compressive strains will cause rippling of the free-standing thin films, membranes, plates, and nanosheets.⁴⁸ The critical compressive strain for rippling instability is much less than the critical tensile strain for fracture, for example, 0.0001% versus 2% in graphene sheets.⁵⁰ However, the rippling can be suppressed by applying constraints, such as embedding (0.7%),⁵¹ substrates (0.4% before heating),⁴⁹ thermal cycling on SiO₂ (0.05%)⁵² and BN (0.6%),⁵³ and sandwiching.⁵⁴ Our study of compressive strains is important in understanding the mechanics of these non-rippling applications. The rippling phenomena are interesting and important, however, they are outside the scope of this study.

3.3 Stress–strain relationships

The second P-K stress versus Lagrangian strain relationships of g-ZnS sheets for uniaxial strains along the armchair and zigzag directions, as well as biaxial strains, are shown in Fig. 3a, compared with that of g-ZnO (Fig. 3b), since these two nano materials were extensively studied and well known. The results show that the g-ZnS sheets can sustain large strains. The ultimate tensile strength is the maximum stress that a material can withstand while being stretched, and the corresponding strain is the ultimate strain. Under ideal conditions, the critical strain is larger than the ultimate strain. The systems of perfect g-ZnS sheets under strains beyond the ultimate strains are in a meta-stable state. The ultimate tensile strain, which reflects the intrinsic bonding strengths and acts as a lower limit of the critical strain, should be considered when exploring the potential applications.

Fig. 3 Stress–strain responses. The stress–strain relationships of the g-ZnS sheets under the armchair, zigzag, and biaxial strains. Σ₁(Σ₂) denotes the x(y) component of stress. “Cont” stands for the fitting of DFT calculations (“DFT”) to continuum elastic theory. The compressive domain is η < 0 (cyan) and the tensile domain is η > 0 (green). The harmonic region is η ≤ η_h and the anharmonic region is η_h < η ≤ η_m. The mechanically unstable region is η > η_m (yellow region), and the mechanically stable region is η ≤ η_m.

The ultimate tensile strains of g-ZnS monolayers are 0.16, 0.22, and 0.19 along the armchair, zigzag, and biaxial, respectively (Table 1). The ultimate tensile stresses are 5.0, 5.2, and 5.7 N m⁻¹ for the armchair, zigzag, and biaxial strains, respectively. The g-ZnS monolayer exhibits small ultimate tensile stresses compared to g-ZnO. They are about one sixth that of graphene. The biaxial ultimate tensile stress is 17.9 GPa if we take the monolayer height of 3.19 Å, which is the separation of layers in a wurtzite ZnS.⁴⁷ The positive ultimate tensile stresses could serve as evidence of its stability and bond strengths. Opposed to the ultimate stresses, the ultimate strains of g-ZnS are smaller, but comparable to that of g-ZnO and graphene. Both the strength and flexibility of g-ZnS are reduced compared to that of g-ZnO monolayers. This could be understood from the viewpoint of bond length. The Zn–S bond length is 2.248 Å, 18% larger than that of Zn–O bonds (1.900 Å). The Zn–S bonds can be viewed as being stretched in prior by the introduction of sulfur atoms, in reference to Zn–O bonds. These stretched bonds are weaker than those un-stretched, resulting in a reduction of the mechanical strength.

Table 1 Elastic limits. Ultimate strengths (Σ^a_m, Σ^z_m, Σ^b_m) in units of N m⁻¹ and ultimate strains (η^a_m, η^z_m, η^b_m) of g-ZnS monolayers under uniaxial (armchair and zigzag) and biaxial strain from DFT calculations, compared with that of g-ZnS and graphene

	g-ZnS	g-ZnO^a	Graphene^b
a Ref. 29.b Ref. 28.
Σ^am	5.0	8.4	28.6
η^am	0.16	0.17	0.19
Σ^zm	5.2	8.7	30.4
η^zm	0.22	0.24	0.23
Σ^bm	5.7	9.3	32.1
η^bm	0.19	0.20	0.23

---

The g-ZnS sheet behaves in an asymmetric manner with respect to compressive and tensile strains. With increasing strains, the Zn–S bonds are stretched and eventually rupture. The positive slope of the stress–strain curves and the positive ultimate tensile stresses indicate that this structure is mechanical stable. Our results show that the g-ZnS monolayers are stable under various strains.

Note that the softening of the g-ZnS monolayers under strains beyond the ultimate strains only occurs for ideal conditions. The systems under this circumstance are in a meta-stable state, which can be easily destroyed by long wavelength perturbations and vacancy defects, as well as high temperature effects, and enter a plastic state.⁵⁵ Thus only the data within the ultimate strain has physical meaning and was used in determining the high order elastic constants in the following subsection.

3.4 Elastic constants

The elastic constants are critical parameters in finite element analysis models for mechanical properties of materials. Our results of these elastic constants provide an accurate continuum description of the elastic properties of g-ZnS monolayers from ab initio density functional theory calculations. They are suitable for incorporation into numerical methods such as the finite element technique.

The second order elastic constants model the linear elastic response. The higher (>2) order elastic constants are important to characterize the nonlinear elastic response of g-ZnS monolayers using a continuum description. These can be obtained using least-squares fit of the DFT data and are reported in Table 2. Corresponding values for graphene are also shown.

Table 2 Elastic constants. Nonzero independent components for the SOEC, TOEC, FOEC, and FFOEC tensor components, Poisson's ratio ν, and in-plane stiffness Y_s of g-ZnS monolayers from DFT calculations, compared with g-ZnO and graphene

		g-ZnS	g-ZnO^a	Graphene^b
a Ref. 29.b Ref. 28.
a (Å)		3.894	3.291	2.468
Ys (N m^−1)		43.6	47.8	340.8
ν		0.508	0.667	0.178
SOECs	C11 (N m^−1)	58.8	86.0	352.0
	C12 (N m^−1)	29.9	57.3	62.6
TOECs	C111 (N m^−1)	−457.7	−525.4	−3089.7
	C112 (N m^−1)	−240.5	−456.3	−453.8
	C222 (N m^−1)	−414.1	−452.7	−2928.1
FOECs	C1111 (N m^−1)	2893	2000	21 927
	C1112 (N m^−1)	2526	3504	2731
	C1122 (N m^−1)	816	2694	3888
	C2222 (N m^−1)	1456	−554	18 779
FFOECs	C11111 (N m^−1)	−10 661	−3296	−118 791
	C11112 (N m^−1)	−9343	−7685	−19 173
	C11122 (N m^−1)	−16 016	−20 648	−15 863
	C12222 (N m^−1)	−20 448	−17 926	−27 463
	C22222 (N m^−1)	−22 473	1607	−134 752

---

The in-plane Young's modulus Y_s (or stiffness) and Poison's ratio ν are obtained from the following relationships: Y_s = (C₁₁² − C₁₂²)/C₁₁ and ν = C₁₂/C₁₁. For g-ZnS, we have Y_s = 43.6 N m⁻¹, and ν = 0.508. The in-plane stiffness is smaller (about 90%) than that of g-ZnO and much smaller (about 1/8) than that of graphene (341 N m (ref. 28)). Taking the experimental value of 3.19 Å (ref. 47) as the height of g-ZnS monolayer, this stiffness is equivalent to 137 GPa, which is comparable to the bulk modulus (around 100 GPa) of a metal, for example Zirconium. Our result also agrees with a previous DFT study.¹³ The Poisson ratio of a g-ZnS monolayer is much larger than that of graphene, but much smaller than of g-ZnO, indicating less shear motion in g-ZnS than in g-ZnO monolayers under strains.

Besides the second order elastic constants, the higher order (>2) elastic constants are also important quantities.^56–58 Experimentally by measuring the changes of sound velocities under the application of hydrostatic and uniaxial stresses, these high order elastic constants can be determined.^59,60 The high order elastic constants are very important in studying the nonlinear elasticity,⁶¹ thermal expansion (through gruneisen parameter),⁶² temperature dependence of elastic constants,^62,63 harmonic generation,⁶⁰ phonon–phonon interactions,⁶⁴ photon–phonon interactions,⁶⁵ lattice defects,⁶⁶ phase transitions,⁶⁷ echo phenomena,⁶⁸ and strain softening,⁶⁹ and so on.²⁸ In addition, with the higher order elastic continuum description utilizing these elastic constants, one can model the stress and deformation state under uniaxial stress, rather than uniaxial strain. Explicitly, when pressure is applied, the pressure dependent second-order elastic moduli can be obtained from the high order elastic continuum description.³⁶ The third-order elastic constants are important in understanding the nonlinear elasticity of materials, such as changes in acoustic velocities due to finite strain. As a consequence, nano devices (such as nano surface acoustic wave sensors and nano waveguides) could be synthesized by introducing local strain,^17,70 as discussed in detail later in this paper.

Stress–strain curves in the previous section show that they will soften when the strain is larger than the ultimate strain. From the view of electron bonding, this is due to the bond weakening and breaking. This softening behavior is determined by the TOECs and FFOECs in the continuum aspect. The negative values of TOECs and FFOECs ensure the softening of g-ZnS monolayer under large strain.

A good way to check the importance of the high order elastic constants is to consider the case when they are missing. With the elastic constants, the stress–strain response can be predicted from elastic theory.³⁶ When we only consider the second order elasticity, the stress varies with strain linearly. Here we take the biaxial deformation of the g-ZnS sheet as an example. As illustrated in Fig. 4, the linear behaviors are only valid within a small strain range, about −0.02 ≤ η ≤ 0.02, as the same result obtained from the energy versus strain curves in Fig. 2. With the knowledge of the elastic constants up to the third order, the stress–strain curve can be accurately predicted within the range of −0.04 ≤ η ≤ 0.04. Using the elastic constants up to the fourth order, the mechanical behaviors can be well treated up to a strain as large as 0.08. For the strains beyond 0.08, the fifth order elastic are required for an accurate modeling. The analysis of the uniaxial deformations comes to a similar result.

Fig. 4 Order effects. The predicted stress–strain responses of biaxial deformation of ordered graphene oxide monolayer from different orders (second, third, fourth, and fifth order) of elastic constants of g-ZnS sheet, compared to that from the DFT calculations (circle line), using biaxial deformation as an example.

Our results illustrate that the monatomic layer structures possess different mechanical behaviors in contrast to the bulk or multi-layered structures, where the second order elastic constants are sufficient in most cases. The second order elastic constants are relatively easier to be calculated from the strain energy curves,^55,71 however, they are not sufficient for monatomic layer structures. The high order elastic constants are required for an accurate description of the mechanical behaviors of monatomic layer structures since they are vulnerable to strain due to the geometry confinements.

Our results of mechanical properties of g-ZnS monolayers are limited to zero temperature due to current DFT calculations. Once finite temperatures are considered, the thermal expansions and dynamics will in general reduce the interactions between atoms. As a result, the longitudinal mode elastic constants will decrease with respect to the temperature of the system. The variation of shear mode elastic constants should be more complex in response to the temperature. A thorough study will be interesting, which is, however, beyond the scope of this study.

3.5 Pressure effect on the elastic moduli

With third-order elastic moduli, one can study the effect of the second-order elastic moduli on the pressure p acting in the plane of g-ZnS monolayers. Explicitly, when pressure is applied, the pressure dependent second-order elastic moduli (C˜₁₁, C˜₁₂, C˜₂₂) can be obtained from C₁₁, C₁₂, C₂₂, C₁₁₁, C₁₁₂, C₂₂₂, Y_s, and ν as:^39–41

The general trend is that the second-order elastic moduli increase linearly with the applied pressure (Fig. 5). C˜₁₁ is asymmetrical to C˜₂₂ unlike the zero pressure case. C˜₁₁ = C˜₂₂ = C₁₁ only occurs when the pressure is zero. This anisotropy could be the outcome of anharmonicity. For Poisson ratio, Fig. 5 depicts a trend that decreases monotonically with the increase of pressure, opposite to that of the second-order elastic moduli. Our results show that the g-ZnS monolayers are mechanically stable under the various pressures.

Fig. 5 Pressure dependent elastic moduli. Second-order elastic moduli and Poisson ratio as function of the pressure for the g-ZnS (a) and g-ZnO (b) monolayers from DFT predictions.

Compared to g-ZnO monolayers, g-ZnS monolayers exhibit smaller pressure dependent second-order elastic moduli, as expected. The shear component C˜₁₂ is more insensitive to pressure. The Poisson's ratio is also smaller in the examining range of −4.0 N m⁻¹ to 4.0 N m⁻¹. However, the change in the Poisson's ratio is more sensitive to the variation of pressure, indicating a more severe shear movement in g-ZnS monolayers in responds to the change of pressures.

4 Conclusions

In summary, by applying various mechanical strains, we studied the mechanical properties and stabilities of the g-ZnS monolayers. The potential profiles, the stress–strain relationships, the in-plane stiffness, Poisson's ratio, the second, third, fourth, and fifth order elastic constants, the ultimate stresses, ultimate strains, critical strains, and the pressure effect on the elastic moduli are studied. g-ZnS monolayers have a low in-plane stiffness, about 1/8 that of graphene, and 90% of g-ZnO monolayers. The potential profiles and the stress–strain curves indicate that the free standing g-ZnS monolayer can sustain large tensile strains, up to 0.16, 0.22, and 0.19 for armchair, zigzag, and biaxial deformations, respectively. However, both the strength and flexibility are reduced compared to graphene-like zinc oxide (g-ZnO) monolayers. According to the results of the positive ultimate strengths and strains, second order elastic constants, and the in-plane Young's modulus, we conclude that g-ZnS structures are mechanically stable under various strains and pressures.

We obtained an accurate continuum description of the elastic properties of this structure by explicitly determining the fourteen independent components of high order (up to fifth order) elastic constants from the fitting of stress–strain curves obtained from DFT calculations. This data is useful to develop a continuum description which is suitable for incorporation into a finite element analysis model for its applications at large scale. We also determined the valid range using these elastic constants in different orders. The harmonic elastic constants are only valid with a small range of −0.02 ≤ η ≤ 0.02. With the knowledge of the elastic constants up to the third order, the stress–strain curve can be accurately predicted within the range of −0.04 ≤ η ≤ 0.04. Using the elastic constants up to the fourth order, the mechanical behaviors can be accurately predicted up to a strain as large as 0.08. For the strains beyond 0.08, the fifth order elastic constants are required for accurate modeling. The high order elastic constants reflect the high order nonlinear bond strength under large strains. We predicted that both the second order elastic constants and the in-plane stiffness monotonically increase with elevating pressure, while the trend of Poisson ratio is reversed. Our results could serve as a safe-guide for the strain-engineering of g-ZnS monolayers for their promising applications and advanced function modification.

Acknowledgements

The authors would like to acknowledge the generous financial support from the Defense Threat Reduction Agency (DTRA) Grant # BRBAA08-C-2-0130 and # HDTRA1-13-1-0025, and a NSF Career award under the Award number of DMR 1151028. Support by Department of Defense (Grant W911NF-12-1-0083) is gratefully acknowledged by Z.C. Computational resources were provided by Rensselaer Polytechnic Institute through AMOS, the IBM Blue Gene/Q system at the Center for Computational Innovations.

References

1. X. Fang, T. Zhai, U. K. Gautam, L. Li, L. Wu, Y. Bando and D. Golberg, Prog. Mater. Sci., 2011, 56, 175–287.
2. H. Xu, Y. Li, A. L. Rosa, T. Frauenheim and R. Q. Zhang, J. Phys. Chem. C, 2008, 112, 20291–20294.
3. T. Yamamoto, S. Kishimoto and S. Iida, Phys. B, 2001, 308310, 916–919.
4. M. Bredol and J. Merikhi, J. Mater. Sci., 1998, 33, 471–476.
5. C. N. Xu, T. Watanabe, M. Akiyama and X. G. Zheng, Appl. Phys. Lett., 1999, 74, 1236–1238.
6. P. Calandra, M. Goffredi and V. Liveri, Colloids Surf., A, 1999, 160, 9–13.
7. T. A. Kennedy, E. R. Glaser, P. B. Klein and R. N. Bhargava, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 52, R14356–R14359.
8. O. Agyeman, C. N. Xu, M. Suzuki and X. G. Zheng, J. Mater. Res., 2002, 17, 959–963.
9. X. Zhang, M. Zhao, S. Yan, T. He, W. Li, X. Lin, Z. Xi, Z. Wang, X. Liu and Y. Xia, Nanotechnology, 2008, 19, 305708.
10. Y. Wei, L. H. Xu, Y. W. Tao, X. J. Zheng, J. H. Yang, D. F. Zou and S. X. Mao, Appl. Phys. Lett., 2012, 101, 091904.
11. S. Yu and M. Yoshimura, Adv. Mater., 2002, 14, 296–300.
12. N. Krainara, J. Limtrakul, F. Illas and S. T. Bromley, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 233305.
13. H. Behera and G. Mukhopadhyay, J. Phys. D: Appl. Phys., 2014, 47, 075302.
14. K. Ariga, T. Mori and J. P. Hill, Adv. Mater., 2012, 24, 158–176.
15. Mechanical Properties of Engineered Materials.
16. J. F. Nye, Physical Properties of Crystals, Oxford Science Publications, Oxford, 1995.
17. Q. Peng, A. R. Zamiri, W. Ji and S. De, Acta Mech., 2012, 223, 2591–2596.
18. Q. Peng and S. De, Phys. E, 2012, 44, 1662–1666.
19. Q. Peng, W. Ji and S. De, Nanoscale, 2013, 5, 695–703.
20. F. Guinea, M. I. Katsnelson and A. K. Geim, Nat. Phys., 2010, 6, 30–33.
21. M. Topsakal, S. Cahangirov, E. Bekaroglu and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 235119.
22. Q. Peng, J. Crean, A. K. Dearden, X. Wen, C. Huang, S. P. A. Bordas and S. De, Mod. Phys. Lett. B, 2013, 27, 1330017.
23. Q. Peng, A. K. Dearden, J. Crean, L. Han, S. Liu, X. Wen and S. De, Nanotechnol., Sci. Appl., 2014, 7, 1–29.
24. Q. Peng and S. De, Nanoscale, 2014, 6, 12071–12079.
25. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499, 419–425.
26. Y. Lei, F. Chen, R. Li and J. Xu, Appl. Surf. Sci., 2014, 308, 206–210.
27. A. K. Kole, P. Kumbhakar and T. Ganguly, J. Appl. Phys., 2014, 115, 224306.
28. Q. Peng, C. Liang, W. Ji and S. De, Phys. Chem. Chem. Phys., 2013, 15, 2003–2011.
29. Q. Peng, C. Liang, W. Ji and S. De, Comput. Mater. Sci., 2013, 68, 320–324.
30. C. A. Marianetti and H. G. Yevick, Phys. Rev. Lett., 2010, 105, 245502.
31. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
32. G. Kresse and J. Furthuller, Comput. Mater. Sci., 1996, 6, 15.
33. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
34. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
35. R. O. Jones and O. Gunnarsson, Rev. Mod. Phys., 1989, 61, 689–746.
36. Q. Peng, W. Ji and S. De, Comput. Mater. Sci., 2012, 56, 11–17.
37. Q. Peng, C. Liang, W. Ji and S. De, Model. Numer. Simul. Mater. Sci., 2012, 2, 76–84.
38. Q. Peng, X.-J. Chen, S. Liu and S. De, RSC Adv., 2013, 3, 7083–7092.
39. Q. Peng, C. Liang, W. Ji and S. De, Mech. Mater., 2013, 64, 135–141.
40. Q. Peng, X.-J. Chen, W. Ji and S. De, Adv. Eng. Mater., 2013, 15, 718–727.
41. Q. Peng, C. Liang, W. Ji and S. De, Appl. Phys. A, 2013, 13, 483–490.
42. Q. Peng, Z. Chen and S. De, Mech. Adv. Mater. Struct., 2013 DOI:10.1080/15376494.2013.839067.
43. Q. Peng, X. Wen and S. De, RSC Adv., 2013, 3, 13772–13781.
44. Q. Peng and S. De, Phys. Chem. Chem. Phys., 2013, 15, 19427–19437.
45. Q. Peng and S. De, RSC Adv., 2013, 3, 24337–24344.
46. Q. Peng, L. Han, X. Wen, S. Liu, Z. Chen, J. Lian and S. De, Phys. Chem. Chem. Phys., 2015, 17, 2160–2168.
47. A. K. Kole and P. Kumbhakar, Results Phys., 2012, 2, 150–155.
48. E. Cerda and L. Mahadevan, Phys. Rev. Lett., 2003, 90, 074302.
49. W. Bao, F. Miao, Z. Chen, H. Zhang, W. Jang, C. Dames and C. N. Lau, Nat. Nanotechnol., 2009, 4, 562–566.
50. Y. Zhang and F. Liu, Appl. Phys. Lett., 2011, 99, 241908.
51. G. Tsoukleri, J. Parthenios, K. Papagelis, R. Jalil, A. C. Ferrari, A. K. Geim, K. S. Novoselov and C. Galiotis, Small, 2009, 5, 2397–2402.
52. D. Yoon, Y.-W. Son and H. Cheong, Nano Lett., 2011, 11, 3227–3231.
53. W. Pan, J. Xiao, J. Zhu, C. Yu, G. Zhang, Z. Ni, K. Watanabe, T. Taniguchi, Y. Shi and X. Wang, Sci. Rep., 2012, 2, 893.
54. R. Quhe, J. Zheng, G. Luo, Q. Liu, R. Qin, J. Zhou, D. Yu, S. Nagase, W.-N. Mei, Z. Gao and J. Lu, NPG Asia Mater., 2012, 4, E6.
55. M. Topsakal, S. Cahangirov and S. Ciraci, Appl. Phys. Lett., 2010, 96, 091912.
56. R. N. Thurston and K. Brugger, Phys. Rev. A, 1964, 133, 1604.
57. K. Brugger, Phys. Rev. A, 1964, 133, 1611.
58. Y. Hiki, Annu. Rev. Mater. Sci., 1981, 11, 51.
59. K. Brugger, J. Appl. Phys., 1965, 36, 768.
60. W. B. Gauster and M. A. Breazeale, Phys. Rev., 1968, 168, 655–661.
61. H. Kobayashi and Y. Hiki, Phys. Rev. B: Solid State, 1973, 7, 594–601.
62. Y. Hiki, J. F. Thomas and A. V. Granato, Phys. Rev., 1967, 153, 764–771.
63. J. A. Garber and A. V. Granato, Phys. Rev. B: Solid State, 1975, 11, 3990–3997.
64. W. P. Mason and T. B. Bateman, J. Acoust. Soc. Am., 1966, 40, 852.
65. J. Melngailis, A. A. Maradudin and A. Seeger, Phys. Rev., 1963, 131, 1972–1975.
66. A. D. Brailsfo, J. Appl. Phys., 1972, 43, 1380.
67. M. Wuttig and T. Suzuki, Acta Metall., 1979, 27, 755–761.
68. K. Fossheim, K. Kajimura, T. G. Kazyaka, R. L. Melcher and N. S. Shiren, Phys. Rev. B: Condens. Matter Mater. Phys., 1978, 17, 964–998.
69. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385–388.
70. Q. Peng, W. Ji and S. De, Phys. Chem. Chem. Phys., 2012, 14, 13385–13391.
71. M. Topsakal, E. Aktürk and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 115442.

---