Giant piezoelectricity of monolayer group IV monochalcogenides: SnSe, SnS, GeSe, and GeS

Piezoelectric materials, which convert mechanical energy to electrical energy, have the advantages of large power densities and ease of application in sensors and energy harvesting.^1,2 For example, a widely used piezoelectric material is lead zirconate titanate Pb[Zr_xTi_1−x]O₃, a piezoceramic known as PZT.^3–5 However, the piezoceramic's brittle nature causes limitations in the sustainable strain.⁶ Meanwhile, non-centrosymmetric wurtzite-structured semiconductors, such as ZnO, GaN, and InN, are wildly used in the piezotronic and piezo-phototronic devices.^7–9 In particular, their nanowires or nanobelts^10–12 are expected to be useful for electromechanical coupled sensors, nanoscale energy conversion.^10–13 However, the much smaller piezoelectric coefficients of wurtzite semiconductors limit the mechanical-electrical energy conversion efficiency.^7,8

Recently, two-dimensional (2D) materials have sparked interest for the piezoelectric applications because of their high crystallinity and ability to withstand enormous strain. For those hexagonal structures with a D_6h point group, such as boron nitride (h-BN) and many transition-metal dichalcogenides (TMDCs), as well as layered orthorhombic structure with a D_4h point group, such as group-III monochalcogenides, their symmetry is reduced to the D_3h group when thinned down to monolayer. This breaks the inversion symmetry, as shown in Figs. 1(a) and 1(b), giving rise to piezoelectricity. They were predicted to be intrinsically piezoelectric,^14–21 and this idea has been demonstrated by experiments on the MoS₂ monolayer.^22–24 Unfortunately, the piezoelectric effect is rather small, e.g., the measured piezoelectric coefficient e₁₁ of monolayer MoS₂ is only around 2.9 × 10⁻¹⁰ C/m,²³ and the mechanical-electrical energy conversion rate is limited to be about 5%.²² 

Therefore, finding flexible, stable, and efficient 2D piezoelectric materials is crucial. This motivates us to study another family of 2D semiconductors, group IV monochalcogenides (MX, M=Sn or Ge, M=Se or S), i.e., SnSe, SnS, GeSe, and GeS. Their atomic structure is presented in Figs. 1(c) and 1(d), which exhibit a C_2v point group. We expect an enhanced piezoelectricity due to the following reasons: (1) As shown in Fig. 1(c), their stable monolayer structures are non-centrosymmetric, allowing them to be piezoelectric. (2) Their puckered C_2v symmetries are much more flexible (softer) along the armchair direction. This can further enhance the piezoelectric. (3) Significant advances in fabrication techniques have been achieved. For example, few-layer SnSe has been fabricated recently.²⁵ 

In this letter, we employ the first-principles density functional theory (DFT) simulations to calculate the piezoelectric effects of monolayer group-IV monochalcogenides. The piezoelectric effect of these monolayer materials is dramatically enhanced and anisotropic, and the most important piezoelectric coefficient d₁₁ is about two orders of magnitude larger than that of 2D and bulk materials,^7,14,16 which have been widely used in the industry. These anisotropic, giant piezoelectric materials represent a new class of nanomaterials that will allow for the next generation of ultra-sensitive mechanical detectors, energy conversion devices, and consumer-touch sensors.

The DFT calculations with the Perdew-Burke-Ernzerh (PBE) functional²⁶ have been carried out by using the Vienna Ab initio Simulation Package (VASP) with a plane wave basis set^27,28 and the projector-augmented wave method.²⁹ The Brillouin zone integration is obtained by a 14 × 14 × 1 k-point grid. The convergence criteria for electronic and ionic relaxations are 10⁻⁶eV and 10⁻³eV/Å, respectively. We use the “Berry-phase” theory of polarization to directly compute the electric polarization.^30–32 The change of polarization (ΔP) occurs upon making an adiabatic change in the Kohn-Sham Hamiltonian of the crystal.

The DFT-optimized monolayer and bulk structure parameters, i.e., the in-plane lattice constants a and b, are listed in Table I. The corresponding experimental or the previous DFT results of the bulk phase are listed as well.^33–42 We observe a similar trend as that in Ref. 42, in which the lattice constant a increases and the constant b decreases with increasing the number of layers for most group IV monochalcogenides, except for GeS. These monolayers are metastable. This is evidenced by recent experimental fabrications²⁵ and theoretical phonon calculations.⁴³ 

We have calculated the electronic structure of group IV monochalcogenides in the supplementary material.⁴⁴ All these materials exhibit an indirect band gap at the DFT level. We list the values of band gaps in Table I. These DFT gap values are for reference purposes only, as excited-state calculations are needed to get the reliable band gap of MXs. According to our experience,^45,46 the quasiparticle band gaps of MX range from 1.2 eV to 2.7 eV, which are within a very useful range for electronic applications. Moreover, huge excitonic effects are expected, which can substantially lower the optical absorption edge, promising for solar energy applications.^25,47

Fortunately, piezoelectric properties are ground-state properties and the DFT calculation is a suitable tool shown to reliably predict the values. For example, the DFT-calculated piezoelectric coefficients are in excellent agreement with the experimental values for the bulk GaN.⁴⁸ Very recently, experiments measured the piezoelectric coefficient e₁₁ = 2.9 × 10⁻¹⁰ C/m for monolayer MoS₂, which is close to the DFT results (3.6 × 10⁻¹⁰ C/m).^14,23 Therefore, we employ the same theoretical approach in this work.

We first obtained the planar elastic stiffness coefficients C₁₁, C₂₂, and C₁₂ of the MX monolayer by fitting the DFT-calculated unit-cell energy U to a series of 2D strain states (ε₁₁,ε₂₂), based on the formula

where A₀ is the unit-cell area at the zero strain. Due to the existence of mirror symmetry along the zigzag direction (y direction), at the small strain limit, we can write

where Δu(ε₁₁, ε₂₂) = [U(ε₁₁, ε₂₂) − U(ε₁₁ = 0, ε₂₂ = 0)]/A₀ is the change of unit-cell energy per area. We carry out the strain energy calculation on an 11 × 11 grid with ε₁₁ and ε₂₂ ranging from −0.005 to 0.005. The atomic positions in the strained unit cell are allowed to be fully relaxed. Following the definitions of previous works,²³ the coefficients C₁₁, C₂₂, and C₁₂, which are calculated using a fully relaxed atomic configuration, are called relaxed-ion stiffness coefficients, which are experimentally relevant. In contrast, if the atomic positions are held when applying strain, the so-called clamped-ion coefficients, which represent the piezoelectric effect from the electronic contribution,⁴⁹ can be calculated as well.

Table II summarizes the clamped and relaxed-ion stiffness coefficients for the four types of C_2v symmetry MX monolayers. Additionally, we have also listed the elastic stiffness of another two typical D_3h symmetry piezoelectric materials, MoS₂¹⁴ and GaSe.¹⁶ According to the structures shown in Figs. 1(c) and 1(d), group IV monochalcogenides are soft along the armchair (x) direction. This is consistent with our DFT results in Table II. In particular, for both clamped and relax-ion cases, the elastic stiffnesses (C₁₁) of group-IV monochalcogenide are about 4–6 times smaller than that of MoS₂ and GaSe. This will significantly enhance the piezoelectric effects. An unexpected result from Table II is that the elastic stiffness (C₂₂) along the zigzag (y) direction is also substantially smaller (around 2–3 times) than that of MoS₂ and GaSe. This may be attributed to the intrinsic electronic properties of group IV monochalcogenides, whose covalence bonds are weaker than those of hexagonal TMDCs and group III monochalcogenides. This is also reflected in their longer bond lengths (2.50–2.89 Å), compared with those of GaSe (2.47 Å) and MoS₂ (1.84 Å).^14,16

Recently, puckered 2D structures, such as few-layer black phosphorus (phosphorene) and similar isoelectronic materials,⁵⁰ have attracted significant research interests. Due to their novel structure, phosphorene exhibits an unusually negative Poisson ratio.⁵¹ Here, we have calculated the Poisson ratio ν⊥ obtained directly from relaxed ion coordinates by evaluating the change of layer thickness in response to in-plane hydrostatic strain Δh/h = −ν_⊥(ε₁₁ + ε₂₂). The Poisson ratio ν⊥ is investigated by averaging the results of the armchair direction and zigzag direction for small stress (−0.8%–0.8%). Interestingly, our calculated value is positive and similar to those of the TMDCs and group III monochalcogenides. This differs from the result of phosphorene, in which the Poisson ratio is evaluated by the value only from the armchair direction within a much larger stress range (−5%–5%).⁵¹ 

Next, we calculated the linear piezoelectric coefficients of the group IV MX monolayers by evaluating the change of unit-cell polarization after imposing uniaxial strain. The linear piezoelectric coefficients e_ijk and d_ijk are third-rank tensors as they relate polarization vector P_i, to strain ε_jk and stress σ_jk, respectively,

Because of the mirror symmetry along the zigzag (y) direction, the independent piezoelectric coefficients are {e₁₁₁, e₁₂₂, e₂₁₂ = e₂₂₁} and {d₁₁₁, d₁₂₂, d₂₁₂ = d₂₂₁}. Indices 1 and 2 correspond to the x and y directions, respectively. The reason that e₂₁₂ = e₂₂₁ and d₂₁₂ = d₂₂₁ is because the strain tensor is usually defined to be symmetric, namely, ε_jk = ε_kj. The piezoelectric coefficients e₂₁₂ and d₂₁₂ describe the response of polarization to shear strain ε₁₂. In the following, we particularly focused on {e₁₁₁, e₁₂₂} and {d₁₁₁, d₁₂₂}, as well as the relationship between the e_ijk and d_ijk.

By definition, the tensors are related by

where C_mnjk are elastic constants. In 2D structures, an index can be either 1 or 2. Therefore,

Using the Voigt notation, we simplify it as e₁₁ = e₁₁₁, e₁₂ = e₁₂₂, d₁₁ = d₁₁₁, d₁₂ = d₁₂₂, C₁₁ = C₁₁₁₁, C₁₂ = C₁₁₂₂ = C₂₂₁₁. Then, we calculated d₁₁ and d₁₂ by e₁₁ and e₁₂ as

We have directly calculated the polarization of the MX monolayers by applying the uniaxial strain ε₁₁ and ε₂₂ to the orthorhombic unit cell along the x and y directions, respectively. The change of polarization along the y direction is zero, because the mirror symmetry still remains under an uniaxial strain for the C_2v point group. The values of e₁₁ and e₁₂ are evaluated by a linear fit of 2D unit-cell polarization change along the x direction (ΔP₁) with respect to ε₁₁ and ε₂₂. In Figs. 2(a) and 2(b), we use ε₁₁ and ε₂₂ ranging from −0.005 to 0.005 in steps of 0.001 in the champed-ion case and −0.01–0.01 in steps of 0.002 in the relax-ion case. The dense steps of 0.001 are required for monolayer SnSe and SnS because their linear polarization changes occurring in the strain region are very small, less than ±0.004, as shown in Figs. 2(c) and 2(d). For e₁₁, the trend of polarization vs. strain is different for the cases of GeS and GeSe, and for e₁₂, the trend is different for SnSe and SnS. This is because of the competition between ionic polarization contribution and electron polarization for e₁₁ of GeS and GeSe, and e₁₂ of SnS and SnSe. The relaxed-ion (or clamped-ion) d₁₁ and d₁₂ coefficients are finally calculated by the corresponding e_11,e₁₂ coefficients and elastic stiffness coefficients C₁₁, C₂₂, and C₁₂ based on Eqs. (8) and (9).

We have summarized the calculated e₁₁, e₁₂, d₁₁, and d₁₂ coefficients in Table III. The most useful piezoelectric coefficients (relaxed-ion d₁₁ and d₁₂), which reflect how much polarization charge can be generated with a fixed force and thus decide the mechanic-electrical energy converting ratio, are about 75–250 pm/V. Compared with those of frequently used bulk piezoelectric materials, such as α-quartz, wurtzite AlN, and ZnO,^7,14,52–54 and recently emerging 2D piezoelectric materials, such as MoS₂ and GaSe,^14,16,23 these values are about two orders of magnitude larger, as shown in Fig. 3.

Finally, we find that the relaxed-ion d₁₁ and d₁₂ coefficients in the MX monolayers obey a periodic trend, as shown in the inset of Fig. 3. GeS possesses the smallest piezoelectric effect (d₁₁ = 75.43 pm/V and d₁₂ = −50.42 pm/V), and moving upward in groups 14 (crystallogens) and 16 (chalcogenide) enhances the magnitude of the effect until SnSe, which has the largest coefficient (d₁₁ = 250.58 pm/V, d₁₂ = −80.31 pm/V), is reached. Interestingly, this trend is similar to that discovered in the hexagonal TMDCs.¹⁴ More calculations on similar puckered C_2v symmetry materials may be necessary to conclude the interesting trend of piezoelectric effects revealed in those honeycomb structures.¹⁸ 

These group IV monochalcogenides have highly desirable properties useful for a broad range of applications. On the other hand, for realistic devices, many other factors, in addition to the piezoelectric coefficients, will need to be considered. For instance, substrate effects and carrier mobilities are important for deciding the converting ratio in energy capture devices and the mechanical fatigue of these flexible materials has not been tested yet. It is also known that the layer number may dramatically impact the piezoelectric effect.²² These are beyond the scope of this letter, but further research is expected.

By reliable simulations, we have demonstrated that monolayer group IV monochalcogenides MX, GeS, GeSe, SnS, and SnSe are strongly piezoelectric. Their piezoelectric coefficients are surprisingly one to two orders of magnitude larger than other frequently used piezoelectric materials. Encouraged by experimental achievements of monolayer samples, we expect that the huge piezoelectric properties of these materials to provide new platforms for electronic and piezotronic devices, and enable previously inaccessible avenues for sensing and control at the nanoscale.

We acknowledge the fruitful discussions with Vy Tran and Anders Carlsson. R. Fei and L.Y. were supported by the National Science Foundation (NSF) Grant No. DMR-1207141 and NSF CAREER Grant No. DMR-1455346. J.L. and W.B.L. acknowledge the support by NSF DMR-1410636 and DMR-1120901. The computational resources have been provided by the Stampede of Teragrid at the Texas Advanced Computing Center (TACC) through XSEDE.