Piezoelectricity in two-dimensional materials: a comparative study between lattice dynamics and ab-initio calculations

TL;DR

This study compares lattice dynamical theory and ab-initio calculations to analyze piezoelectric properties of 2D materials like h-BN and MoS2, revealing microscopic factors influencing their piezoelectric behavior.

Contribution

It provides a detailed comparison of lattice dynamical and ab-initio methods for calculating piezoelectric constants in 2D materials, clarifying the role of inner strains and microscopic quantities.

Findings

Inner strain contributions depend on ionic displacements

Lattice dynamics expresses inner strains via microscopic quantities

Differences in piezoelectric behavior between h-BN and MoS2 explained

Abstract

Elastic constant C_{11} and piezoelectric stress constant e_{1,11} of two-dimensional (2D) dielectric materials comprising h-BN, 2H MoS2 and other transition metal dichalcogenides (TMDCs) and -dioxides (TMDOs) are calculated using lattice dynamical theory. The results are compared with corresponding quantities obtained by ab-initio calculations. We identify the difference between clamped-ion and relaxed-ion contributions with the dependence on inner strains which are due to the relative displacements of the ions in the unit cell. Lattice dynamics allows to express the inner strains contributions in terms of microscopic quantities such as effective ionic charges and optoacoustical couplings, which allows us to clarify differences in the piezoelectric behavior between h- BN versus MoS2. Trends in the different microscopic quantities as functions of atomic composition are discussed.

Tables

Table 1. Table 1: c 11 subscript c 11 \mathrm{c}_{11} in units of N/m, e 1 , 11 subscript 𝑒 1 11 e_{1,11} in units of 10 -10 C/m. Inner strain values ( i ​ s 𝑖 𝑠 is ) follow from differences of relaxed-ion ( r ​ i 𝑟 𝑖 ri ) and clamped-ion ( c ​ i 𝑐 𝑖 ci ) values, superscripts (1) and (2) refer to input values taken from Refs. 12, and 16, , respectively. ω E ′ subscript 𝜔 superscript E ′ \omega_{\mathrm{E}^{\prime}} in units of cm -1 .

Compound	${c_{11}\mid }_{ri}$	${c_{11}\mid }_{ci}$	${c_{11}\mid }_{is}$	${e_{1,11}\mid }_{ri}$	${e_{1,11}\mid }_{ci}$	${e_{1,11}\mid }_{is}$	${\omega }_{E^{\prime }}$
2D h-BN	291	300	-9	1.38	3.71	-2.33	1310
2H-MoS2	130	153	-23	3.64	3.06	0.58(1)	390
	133	157	-24	4.93	3.20	1.73(2)
2H-MoSe2	108	131	-23	3.92	2.80	1.12(1)	278
	107	133	-26	5.26	2.92	2.34(2)
2H-MoTe2	80	101	-21	5.43	2.98	2.45(1)	231
	84	106	-22	6.55	2.75	3.80(2)
2H-WS2	144	170	-26	2.47	2.20	0.27(1)	348
	146	175	-29	3.76	2.33	1.43(2)
2H-WSe2	119	147	-28	2.71	1.93	0.78(1)	240
	120	147	-27	3.98	2.05	1.93(2)
2H-WTe2	89	116	-27	3.40	1.60	1.80(1)	192
	89	115	-26	5.02	1.75	3.27(2)

Table 2. Table 2: Surface mass density ρ 2 ​ D subscript 𝜌 2 D \rho_{\mathrm{2D}} (units 10 -6 kg/m 2 ); optical mode frequency ω E ′ subscript 𝜔 superscript E ′ \omega_{\mathrm{E}^{\prime}} (units cm -1 ); Opto-acoustic coupling D 11 , 1 ⋅ ⁣ ⋅ ( 1 ) subscript superscript D ⋅ ⋅ absent 1 11 1 \mathrm{D}^{\mathrm{\cdot\cdot}(1)}_{11,1} (units 10 19 cm s -2 ); effective charges e M ∗ subscript superscript 𝑒 M e^{*}_{\mathrm{M}} of metal ions (units elementary charge e 𝑒 e = 1.602 × \times 10 -19 C), (a) present theory, superscripts (1) and (2) refer to e 1 , 11 ∣ i ​ s ( 1 ) evaluated-at subscript 𝑒 1 11 𝑖 𝑠 1 e_{1,11}\mid_{is}^{(1)} and e 1 , 11 ∣ i ​ s ( 2 ) evaluated-at subscript 𝑒 1 11 𝑖 𝑠 2 e_{1,11}\mid_{is}^{(2)} from Table 1 , respectively, (b) effective charge calculated with DFPT.

Compound	${\rho }_{2D}$	${\omega }_{E^{\prime }}$	$D_{11,1}^{\cdot \cdot \left(1\right)}$	$e_M^{\ast \left(a\right)}$	$e_M^{\ast \left(b\right)}$
2D h-BN	0.76	1310	-8.44	+2.78	2.71
2H-MoS2	3.14	390	-1.99	-0.41(1)	-1.06
				-1.21(2)
2H-MoSe2	4.42	278	-0.91	-0.70(1)	-1.84
				-1.45(2)
2H-MoTe2	5.34	231	-0.63	-1.59(1)	-3.28
				-2.46(2)
2H-WS2	4.70	348	-1.14	-0.17(1)	-0.53
				-0.92(2)
2H-WSe2	5.60	240	-0.68	-0.49(1)	-1.22
				-1.21(2)
2H-WTe2	6.60	192	-0.51	-1.11(1)	-2.6
				-2.01(2)

Table 3. Table 3: Symbols have the same meaning and unit as in Table 2 .

Compound	${\rho }_{2D}$	${\omega }_{E^{\prime }}$	$D_{11,1}^{\cdot \cdot \left(1\right)}$	$e_M^{\ast \left(a\right)}$	$e_M^{\ast \left(b\right)}$
CrS2	2.41	415	-1.44	-1.79	-2.44
CrSe2	3.10	317	-0.82	-2.33	-3.06
CrTe2	4.87	262	-0.58	-3.19	-4.13
CrO2	2.34	591	-2.01	-0.189	1.334
MoO2	3.06	522	-1.91	0.215	2.875
WO2	5.15	491	-1.44	0.382	3.326

Table 4. Table 4: Trends in effective ionic charges e M ∗ subscript superscript 𝑒 M e^{*}_{\mathrm{M}} (units elementary charge e 𝑒 e ) for 2H-MX 2 ; integers 16, 24 ,… are charge numbers Z.

	Cr24	Mo42	W74
S16	-2.44	-1.06	-0.53
Se34	-3.05	-1.84	-1.22
Te52	-4.19	-3.28	-2.60

Table 5. Table 5: Trends in opto-acoustic couplings D 11 , 1 ⋅ ⁣ ⋅ ( 1 ) subscript superscript D ⋅ ⋅ absent 1 11 1 \mathrm{D}^{\mathrm{\cdot\cdot}(1)}_{11,1} in units 10 19 cm s -2 .

	Cr24	Mo42	W74
S16	-1.44	-1.43	-1.04
Se34	-0.82	-0.91	-0.68
Te52	-0.58	-0.63	-0.51

Table 6. Table 6: Opto-acoustic coupling D 11 , 1 ⋅ ⁣ ⋅ ( 1 ) subscript superscript D ⋅ ⋅ absent 1 11 1 \mathrm{D}^{\mathrm{\cdot\cdot}(1)}_{11,1} in units 10 19 cm s -2 ; ω E ′ subscript 𝜔 superscript E ′ \omega_{\mathrm{E}^{\prime}} in units of cm -1 ; e 1 , 11 ∣ i ​ s evaluated-at subscript 𝑒 1 11 𝑖 𝑠 e_{1,11}\mid_{is} (units 10 -10 C/m); surface mass density ρ 2 ​ D subscript 𝜌 2 𝐷 \rho_{2D} (units of 10 -6 kg/m 2 ); effective charges e M ∗ subscript superscript 𝑒 M e^{*}_{\mathrm{M}} of metal ions (units e 𝑒 e ), (a) present theory, and (b) ab-initio calculations using density functional perturbation theory.

Compound	$D_{11,1}^{\cdot \cdot \left(1\right)}$	${\omega }_{E^{\prime }}$	$e_{1,11}$${\mid }_{is}$	${\rho }_{2D}$	$e_M^{\ast \left(a\right)}$	$e_M^{\ast \left(b\right)}$
CrS2	-1.440	415	1.30	2.408	-0.971	-2.443
CrSe2	-0.817	317	1.80	3.100	-1.329	-3.056
CrTe2	-0.582	262	2.68	4.872	-1.950	-4.130
MoS2	-1.429	377	0.56	3.029	-0.376	-1.063
MoSe2	-0.906	278	1.09	4.421	-0.677	-1.840
MoTe2	-0.626	231	2.16	5.344	-1.396	-3.277
WS2	-1.135	348	0.28	4.695	-0.179	-0.526
WSe2	-0.681	240	0.73	5.600	-0.460	-1.224
WTe2	-0.508	192	1.74	6.671	-1.071	-2.600
CrO2	-2.010	591	-0.38	2.335	0.300	1.334
MoO2	-1.910	522	-0.98	3.063	0.658	2.875
WO2	-1.440	491	-1.15	5.153	0.745	3.326

Equations

$${\mathrm{D}}_{ij}^{\kappa {\kappa }^{\prime }}(\overrightarrow{q})={\mathrm{D}}_{ij}^{\kappa {\kappa }^{\prime }(0)}+i\sum _k{\mathrm{D}}_{ij,k}^{\kappa {\kappa }^{\prime }(1)}{q}_k+\frac{1}{2}\sum _{kl}{\mathrm{D}}_{ij,kl}^{\kappa {\kappa }^{\prime }(2)}{q}_k{q}_l+....$$

$${\hat{\mathrm{D}}}_{ij}(\overrightarrow{q})=\frac{1}{{\rho }_{2D}}\sum _{kl}\{[ij,kl]+(ik,jl)\}{q}_k{q}_l.$$

$$[ij,kl]=\frac{1}{2A_{2D}}\sum _{\kappa {\kappa }^{\prime }}(m_{\kappa }{m}_{{\kappa }^{\prime }}{)}^{1/2}{\mathrm{D}}_{ij,kl}^{\kappa {\kappa }^{\prime }(2)},$$

$$(ik,jl)=-\frac{1}{2A_{2D}}\sum _{\kappa {\kappa }^{\prime }}\sum _{hp}{\Gamma }_{hp}^{\kappa {\kappa }^{\prime }}[\sum _{{\kappa }^{\prime \prime }}(m_{{\kappa }^{\prime \prime }}{)}^{1/2}{\mathrm{D}}_{hi,k}^{\kappa {\kappa }^{\prime \prime }(1)}][\sum _{{\kappa }^{\prime \prime \prime }}(m_{{\kappa }^{\prime \prime \prime }}{)}^{1/2}{\mathrm{D}}_{pj,k}^{{\kappa }^{\prime }{\kappa }^{\prime \prime \prime }(1)}].$$

$${\Gamma }_{hp}^{\kappa {\kappa }^{\prime }}=\sum _{\lambda }\frac{{\xi }^{(\lambda )}(\kappa h){\xi }^{(\lambda )}({\kappa }^{\prime }p)}{{\omega }_{\lambda }^2},$$

$${\mathrm{c}}_{11}={\mathrm{C}}_{11,11}=[11,11]+(11,11),$$

$${\mathrm{c}}_{12}={\mathrm{C}}_{11,22}=[11,11]-2[11,22]+(11,22),$$

$${\mathrm{c}}_{11}={\mathrm{c}}_{11}{\mid }_{ri}={\mathrm{c}}_{11}{\mid }_{ci}+{\mathrm{c}}_{11}{\mid }_{is},$$

$${\mathrm{c}}_{11}{\mid }_{ci}=[11,11]$$

$${\mathrm{c}}_{11}{\mid }_{is}=(11,11)$$

$$e_{i,jl}{\mid }_{is}=\frac{1}{A_{2D}}\sum _{\kappa {\kappa }^{\prime }{\kappa }^{\prime \prime }}\sum _{hk}\sqrt{m_{\kappa }}{\mathrm{D}}_{jh,l}^{\kappa {\kappa }^{\prime }(1)}{\Gamma }_{hk}^{{\kappa }^{\prime }{\kappa }^{\prime \prime }}\frac{Z_{ik}^{(0){\kappa }^{\prime \prime }}}{\sqrt{m_{{\kappa }^{\prime \prime }}}}.$$

$$e_{1,11}=e_{1,11}{\mid }_{ri}=e_{1,11}{\mid }_{ci}+e_{1,11}{\mid }_{is}.$$

$${\xi }_1^{{\mathrm{E}}^{\prime }}=(\sqrt{\frac{\mu }{m_{\mathrm{N}}}},-\sqrt{\frac{\mu }{m_{\mathrm{B}}}}),$$

$${\mathrm{c}}_{11}{\mid }_{is}\equiv (11,11{)}_{2D}=-{\rho }_{2D}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}({\mathrm{D}}_{11,1}^{\mathrm{NB}(1)}{)}^2,$$

$$e_{1,11}{\mid }_{is}\equiv [1,11]={\rho }_{2D}{\mathrm{D}}_{11,1}^{\mathrm{NB}(1)}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}\frac{e_{\mathrm{B}}^{\ast }}{\sqrt{m_{\mathrm{B}}{m}_{\mathrm{N}}}}.$$

$${\xi }_1^{{\mathrm{E}}^{\prime }}=(\sqrt{\frac{\mu }{4m_{\mathrm{S}}}},\sqrt{\frac{\mu }{4m_{\mathrm{S}}}},-\sqrt{\frac{\mu }{m_{\mathrm{Mo}}}}),$$

$${\mathrm{c}}_{11}{\mid }_{is}\equiv (11,11{)}_{2D}=-2{\rho }_{2D}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}({\mathrm{D}}_{11,1}^{\mathrm{S}\mathrm{Mo}(1)}{)}^2,$$

$$e_{1,11}{\mid }_{is}=[1,11]=\sqrt{2}{\rho }_{2D}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}{\mathrm{D}}_{11,1}^{\mathrm{S}\mathrm{Mo}(1)}\frac{e_{\mathrm{Mo}}^{\ast }}{\sqrt{2m_{\mathrm{S}}{m}_{\mathrm{Mo}}}}.$$

$${\mathrm{c}}_{11}{\mid }_{ri}-{\mathrm{c}}_{11}{\mid }_{ci}=-2{\rho }_{2D}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}({\mathrm{D}}_{11,1}^{\mathrm{XM}(1)}{)}^2,$$

$$e_{1,11}{\mid }_{ri}-e_{1,11}{\mid }_{ci}=\sqrt{2}{\rho }_{2D}\frac{1}{({\omega }_{{\mathrm{E}}^{\prime }}{)}^2}{\mathrm{D}}_{11,1}^{\mathrm{XM}(1)}\frac{e_{\mathrm{M}}^{\ast }}{\sqrt{2m_{\mathrm{X}}{m}_{\mathrm{M}}}}.$$

Full text

I Introduction

Piezoelectricity is the manifestation of electro-mechanical coupling that is present in non-centrosymmetric dielectric crystals1 . An electric polarization occurs in response to macroscopic strains, its converse is the change of shape of the crystal upon application of an electric field. Numerous technological applications are based on the use of bulk (three-dimensional) crystals and ceramics. Piezoelectricity in two-dimensional and in layered crystals is still at the stage of fundamental research.

The synthesis of hexagonal boron-nitride (h-BN) nanotubes2 ; 3 has stimulated theoretical work on polarization and piezoelectricity in BN nanotubes and its dependence on their topology4 ; 5 ; 6 ; nw1 . Even more the discovery of graphene and other 2D crystals like h-BN and MoS2 has opened the road for the study of mechanical and electric properties of a new class of materials with controlled composition and number of atomic layers (i.e. its thickness)7 ; 8 .

2D h-BN, the structurally most simple dielectric crystal, has been used by as a prototype in studies of piezo- and flexoelectricity by ab-initio calculations6 ; 9 . Lattice dynamical theory has been used to study piezoelectricity in single and multilayered crystals with application to h-BN as a specific example10 ; 11 .

First-principles calculations on 2H-MoS2 and other transition metal dichalcogenides (TMDCs) have revealed that piezoelectricity in these materials is even more pronounced than in 2D h-BN12 . Notice that h-BN and MoS2 have the same layer number $N$ dependent symmetry, D3h and D3d, for odd and even $N$, respectively13 . Piezoelectricity has been measured recently in MoS2 layered crystals14 ; 15 . In agreement with predictions originally made by analytical theory11 it was found that only crystals with an uneven number of layers are piezoelectric with strength of the piezoelectric stress coefficient decreasing as 1/$N$ with increasing $N$. Piezoelectric and elastic properties of a broad range of transition metal (group IV B, VI B) dichalcogenides and dioxides monolayers have been investigated by means of first-principles calculations. It was found that Ti-, Zn-, Sn-, and Cr-based TMDCs and TMDOs have much better piezoelectric properties as compared to Mo- and W-based materials16 .

The analytical microscopic theory of piezoelectricity goes back to Born and is based on the same lattice dynamical concepts as the theory of elastic properties in ionic crystals17 ; 19 . On the other hand, nowadays first-principles density functional investigations are the method of choice for studying piezoelectric and elastic properties of solidsabinitio ; abinitio2 ; abinitio3 . The aim of the present paper is to confront the results of analytical theory with first-principles calculations in 2D materials. Such an investigation will allow us to elucidate trends in the evolution of physical properties as a function of ionic composition in classes of similar materials with the same crystal structure such as TMDCs and TMDOs.

The content of the paper is as follows. In Sec. II we recall basic concepts of the analytical theory of elastic constants and piezoelectric moduli in 2D crystals. Next (Sect. III) we apply the analytical theory to 2D crystals with D3h point group symmetry, treating h-BN and 2H-MoS2 as specific examples. In Sect. IV we will give a discussion of numerical results obtained from $ab$-initio calculations and analytical theory for a broad range of 2D TMDCs and TMDOs materials.

II Basic theory

II.1 Elastic constants

The analytical theory of elasticity of crystals of composite structures is based on Born’s method of long waves 17 . In three-dimensional (3D) ionic crystals the formulation is complicated by the appearance of divergent results due to the long range Coulomb forces. In 2D crystals these divergences are absent 18 and we can restrict ourselves to the simple treatment. In particular, the corrections to the elastic constants due to piezoelectricity, well known in 3D crystalslandau , vanish in the 2D case 10 . The basic theoretical quantity in lattice dynamics is the dynamical matrix $\mathrm{D}(\overrightarrow{q})$, here $\overrightarrow{q}$=($q_{x}$,$q_{y}$) is a wave vector in the 2D Brillouin zone (BZ). One expands the dynamical matrix $\mathrm{D}(\overrightarrow{q})$ in powers of components of the wave vector $\overrightarrow{q}$. The dynamical matrix of order 3$s$ has the elements $\mathrm{D}_{ij}^{\kappa\kappa^{\prime}}(\overrightarrow{q})$, where $i$ ($j$) = $x,y,z$ labels the cartesian displacements and $\kappa$=1, 2,…$s$ refer to the particles in the unit cell. The long wave expansion for non-primitive crystals reads

[TABLE]

The matrices $\mathrm{D}^{(0)}$, $\mathrm{D}^{(1)}$ and $\mathrm{D}^{(2)}$ represent the moments of the atomic force constants 17 ; 19 and have the following meaning: $\mathrm{D}^{(0)}$ determines the optical phonon frequencies at the $\Gamma$ point of the BZ, $\mathrm{D}^{(1)}$ accounts for the coupling between optical and strains, $\mathrm{D}^{(2)}$ accounts for the coupling between strains. Notice that the elements $\mathrm{D}_{ij,k}^{\kappa\kappa^{\prime}(1)}(\overrightarrow{q})$ are only different from zero if the ions are not centers of symmetry, and that $\mathrm{D}_{ij,k}^{\kappa\kappa^{\prime}(1)}(\overrightarrow{q})$=-$\mathrm{D}_{ij,k}^{\kappa^{\prime}\kappa(1)}(\overrightarrow{q})$. By perturbation theory one eliminates the optical displacements and obtains the acoustic dynamical matrix

[TABLE]

Here $\rho_{2D}$ is the surface mass density, $[ij,kl]$ stands for

[TABLE]

and

[TABLE]

Here $A_{2D}$ is the area of the unit cell, $m_{\kappa}$ is the mass of particle $\kappa$ and $i,j,k,l$ are Cartesian indices. The quantity $\Gamma_{hp}^{\kappa\kappa^{\prime}}$ is given by

[TABLE]

where $\omega_{\lambda}$ and $\xi^{(\lambda)}$ are the optical eigenfrequencies and eigenvectors of the matrix $\mathrm{D}^{(0)}$, respectively. As is obvious from Eqs. (3) and (4), the quantities $[ij,kl]$ refer to the situation where centers of mass displacements of the unit cell contribute to homogenous crystal strains, while $(ik,jl)$ account for relative displacements of ions within the unit cells, also called internal strains 17 or inner displacements 19 .

There are two independent elastic constants in the case of monolayer hexagonal crystals. They have the dimension of surface tension coefficient. Using Voigt’s notation ($xx$=1, $yy$=2), one has 10

[TABLE]

[TABLE]

with $\mathrm{c}_{11}$-$\mathrm{c}_{12}$=2$c_{66}$.

In ab-initio calculations of the elastic constants one distinguishes clamped ion and relaxed ion terms which we denote with subscripts $ci$ and $ri$, respectively. With the foregoing comments about the physical origin of the square brackets and the round bracket terms, we identify the square brackets terms on the right hand side of Eqs. (6) and (7) with the clamped ion contributions while the sum of square and round brackets terms corresponds to the relaxed ion terms, i.e. the experimentally measured elastic constants. Writing

[TABLE]

where the subscript $is$ stands for internal strains, we identify

[TABLE]

[TABLE]

II.2 Piezoelectric constants

Born’s long wave method 17 ; 19 allows one to calculate the piezoelectric moduli in ionic crystals that do not possess a center of inversion symmetry. Within lattice dynamical theory 19 we write the internal strain contribution (suffix $is$) to the piezoelectric stress tensor in 2D as

[TABLE]

Here $Z_{ik}^{(0)\kappa^{\prime\prime}}$ is the Fourier transform of the transverse effective charge tensor19 taken at the $\Gamma$ point of the 2D BZ, the other quantities have the same meaning as in the case of Eq. (4). Obviously, the right hand side of Eq. (11) takes only into account the inner displacements (strains) of the ions. It is therefore also called ionic contribution to the piezoelectric modulusabinitio , in contradistinction to the clamped ion or electronic contribution (see below). In the case of a 2D hexagonal crystal with $D_{3h}$ symmetry there exists only one independent nonzero piezoelectric stress constant $e_{1,11}$. The internal strain contribution $e_{1,11}\mid_{is}$=[1,11] to $e_{1,11}$ is evaluated within analytic theory by means of Eq. (11). In addition to the internal displacement term, the piezoelectric tensor is made up of a second contribution due to the redistribution of the electronic charge cloud upon application of a homogeneous macroscopic strainmartin ; baroni ; 19 . In analogy with Eq. (8) we then write

[TABLE]

In $ab$-initio calculations one then distinguishes again relaxed ion and clamped ion contributions. The former should be compared with experimentally measured quantity, the latter is obtained from calculations in absence of internal ionic displacements. While the inner displacements always lead to a reduction of the elastic constants in comparison with the clamped ion contribution, this is not necessarily so for the piezoelectric constants, as we will show in Sect. III.

In first-principles calculations the elastic stiffness tensor and piezoelectric tensor coefficients, $e_{ijk}$, are obtained by using density-functional perturbation theory (DFPT)baroni as implemented in the Vienna Ab initio Simulation Package (VASP) codevasp1 ; vasp2 ; vasp3 ; vasp4 . Here, a highly dense $k$-point mesh, 36$\times$36$\times$1, is used to accurately predict these tensor components. The clamped-ion elastic and piezoelectric coefficients are obtained from the purely electronic contribution and the relaxed-ion coefficients are obtained from the sum of ionic and electronic contributions. Within DFPT the VASP code gives electronic and ionic contribution to the piezoelectric tensor directly. A different approach, also implemented in the VASP code, is based on the Berry phase conceptabinitio2 ; berry2 . Here one calculates the polarization for a particular strain, the piezoelectric tensor then follows by calculating the change in polarization due to a strain change. We have used the Berry phase approach with applied uniform strain. At this point, in order to apply strain in a desired direction, the hexagonal primitive cell structure of each material is transformed to a tetragonal one composed of two hexagonal primitive cells12 ; 16 . A 24$\times$24$\times$1 $k$-point mesh is used to calculate the change in polarization. For all the calculations, the exchange-correlation interactions are treated using the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) formulationcem3 . The single electron wave functions are expanded in plane waves with a kinetic energy cutoff of 600 eV. For the structure optimizations, the Brillouin-zone integrations are performed using a $\Gamma$-centered regular 26$\times$26$\times$1 $k$-point mesh within the Monkhorst-Pack schemecem4 . The convergence criterion for electronic and ionic relaxations are set as 10*-7* and 10*-3* eV/Å , respectively. In order to minimize the periodic interaction along the $z$-direction the vacuum space between the layers is taken at least 15 Å.

III Composite 2D materials

Two-dimensional hexagonal boron nitride (Fig. 1(a)), the structurally most simple non centrosymmetric dielectric crystal with point group D3h and two atoms per unit cell, has been used as a prototype for studies of piezoelectricity by ab-initio methods6 ; 9 . Lattice dynamical theory has been applied to study elastic and piezoelectric effects in single and multilayered crystals10 ; 11 . Although h-BN was considered as a specific example, the analytical results are general and are readily applied to materials with the same crystal symmetry such as MoS2 and other TMDCs and TMDOs.

Since there are only two atoms per unit cell, in case of B and N, the x$\equiv$1 component of the eigenvector of the optical mode of symmetry E*′* (Fig. 2) is given by

[TABLE]

where $\mu$=m_{B}$$m_{N}/M is the reduced mass, with M=$m_{B}$+$m_{N}$ the total mass per unit cell. The corresponding eigenfrequency is the doubly degenerate mode \omega_{E^{\prime}}$$\equiv$$\omega_{TO}=$\omega_{LO}$ at $\Gamma$.

The internal strain contribution to $c_{11}$ follows from Eq. (4) and is given by10

[TABLE]

where $\rho_{2D}$=M/$A_{2D}$ is the surface mass density. The unit cell area is $A_{2D}$=a^{2}$$\sqrt{3}/2, where $a$=$\mid\overrightarrow{a_{1}}\mid$=$\mid\overrightarrow{a_{2}}\mid$ is the length of the in-plane basic vectors of the hexagonal lattice. Likewise the internal strain contribution to the piezoelectric stress constant $e_{1,11}$ follows from Eq. (11) and is given by

[TABLE]

Here the quantity $e^{\ast}_{\mathrm{B}}$, to be called effective charge of boron ion, stands for the effective charge tensor component $Z_{11}^{(0)\mathrm{B}}$. In the following we will use the notation $e^{\ast}_{\kappa}$ for $Z_{11}^{(0)\kappa}$.

We now extend these results to transition metal dichalcogenides and dioxides MX2 where M=Mo, W, Cr and X=O, S, Se, Te. These materials have the same point group symmetry as 2D h-BN.

As a specific example we consider MoS2, see Fig. 1(b). We recall that the nomenclature 2H refers to the most abundant polytype with trigonal prismatic coordination between an Mo center surrounded by six sulfide ligands. Each sulfur center is pyramidal and connected to three Mo centersSchonfeld , Fig. 1(b) The internal displacements that contribute to $c_{11}$ and $e_{1,11}$ are of E*′* symmetry where the two S ions on top of each other in each unit cell move in unison and Mo in opposite directionmolina (Fig. 2). We then can treat them as an effective particle with mass 2$m_{\mathrm{S}}$ and charge 2$e^{\ast}_{\mathrm{S}}$=-$e^{\ast}_{\mathrm{Mo}}$. The $x\equiv$1 component of the optical eigenvectors of E*′* symmetry reads

[TABLE]

where $\mu$=2m_{\mathrm{S}}$$m_{\mathrm{Mo}}/(2$m_{\mathrm{S}}$+$m_{\mathrm{Mo}}$) is the reduced mass. It can be shown that the E” mode, where the two S ions move in opposite direction while Mo stays at restmolina , does not contribute to $c_{11}\mid_{is}$ or to $e_{1,11}\mid_{is}$ . Writing $\mathrm{D}^{\mathrm{Mo}\mathrm{S}(1)}_{11,1}$ for the matrix element $\mathrm{D}^{\kappa\kappa^{\prime}(1)}_{11,1}$ where $\kappa$=Mo and $\kappa^{\prime}$=S, we obtain by means of Eqs. (4), (5) and (16)

[TABLE]

where $\rho_{2D}$=M/$A_{2D}$ is the surface mass density, M=2m_{S}+$$m_{Mo} the total mass per unit cell, and $A_{2D}$ the basis area of the hexagonal unit cell. Similarly we find by means of Eqs. (5), (11) and (16)

[TABLE]

The translation of Eqs. (17) and (18) to other MX2 compounds is straightforward. Using Eqs. (8) and (17) we can write

[TABLE]

and similarly, combining Eqs. (12) and (18) we get

[TABLE]

Here $\rho_{2D}$ and $\omega_{\mathrm{E}^{\prime}}^{(0)}$ are the surface mass density and E’ the optical mode frequency specific for the given MX2 compound. In the next section we will use Eqs. (19) and (20) for a numerical evaluation of the internal strain coupling $\mathrm{D}^{\mathrm{X}\mathrm{M}(1)}_{11,1}$ and of the effective charge $e^{\ast}_{\mathrm{M}}$.

IV Numerical results

In the previous section we have derived from lattice dynamical theory analytical expressions for the strain contributions $c_{11}\mid_{is}$ and $e_{1,11}\mid_{is}$ to the elastic and piezoelectric constants in 2D hexagonal crystals. On the other hand numerical values of the internal strain quantities can be obtained by taking the differences of the corresponding relaxed-ion and clamped-ion quantities that are calculated by $ab$-$initio$ methods. Whenever the optical phonon frequency $\omega_{\mathrm{E}^{\prime}}$ is known from $ab$-$initio$ or from experiment, Eqs. (19) and (20) allow us to determine the value of the opto-acoustic coupling $\mathrm{D}^{\mathrm{X}\mathrm{M}(1)}_{11,1}$ and of the effective charge $e^{\ast}_{\mathrm{M}}$ respectively for a given material. Such an analysis will also allow us to understand the magnitude of macroscopic quantities such as $c_{11}$ and $e_{1,11}$ on the basis of atomistic concepts that are specific for various materials. The results for 2D h-BN and a series of TMDCs and TMDOs will be presented in Tables 1, 2 and 3.

We first compare the elastic and piezoelectric properties of 2D h-BN with those of MoS2, where the last material is also representative for the other TMDCs. In Table 1 we have quoted values of $c_{11}$ and $e_{1,11}$ calculated by ab-initio calculations12 ; 16 under clamped-ion and relaxed-ion conditions. The piezoelectric coefficients $e_{1,11}\mid_{ri}$ and $e_{1,11}\mid_{ci}$ have been obtained12 ; 16 by the Berry phase method. The corresponding values of the internal strains contributions $c_{11}\mid_{is}$ and $e_{1,11}\mid_{is}$ are then obtained by means of Eqs. (8) and (12). Notice that the relative large difference between the values of $e_{1,11}\mid_{ri}$ obtained in Refs.12, and 16, respectively, leads to relatively larger difference in the values of $e_{1,11}\mid_{is}$. As mentioned in earlier work16 we attribute the difference in the values reported in Ref.12, to be likely due to the use of different pseudopotentials and other computational parameters.

As a general consequence of Eq. (4) the internal strains, irrespective of the material, yield a negative contribution to the elastic constants. In the present case, see expressions of Eq. (14) and (17) of $c_{11}\mid_{is}$ for 2D h-BN and 2H-MoS2 respectively and results in Table 1. On the other hand the situation is different for the contribution of internal strains to $e_{1,11}\mid_{is}$. As Eq. (11) suggests, the sign of the effective ionic charges determines the sign of $e_{1,11}\mid_{is}$. While the clamped-ion piezoelectric constant $e_{1,11}\mid_{ci}$ is larger for 2D h-BN than for 2H-MoS2, the relaxed-ion piezoelectric modulus $e_{1,11}\mid_{ri}$, which corresponds to the experimentally measured quantity 14 ; 15 , is larger for 2H-MoS2 than for 2D h-BN (see Table 1). In the latter material the internal strain component $e_{1,11}\mid_{is}$ leads to a reduction of $e_{1,11}$ and in the former to an increase. Notice that the large absolute value of $e_{1,11}\mid_{is}$ of 2D h-BN is in accordance with the large ionic contribution to the static dielectric response obtained from ab-initio finite electric field calculations in BN nanotubes nw1 . A positive value of $e_{1,11}\mid_{ci}$ has also been obtained by analytical calculations Droth using the Berry phase method. We observe that the opposite sign of the ionic contribution $e_{1,11}\mid_{is}$ in comparison with the electronic contribution $e_{1,11}\mid_{ci}$ as is here the case for 2D h-BN is not uncommon in other piezoelectric materials and has been found also in 3D III-IV semiconductors.abinitio

We now show that this different behavior of 2D h-BN and 2H-MoS2 is due to the opposite sign of the effective charges $e^{\ast}_{\mathrm{B}}$ and $e^{\ast}_{\mathrm{Mo}}$. We first calculate the effective charges by a semi-analytical method, inverting Eq. (11). Thereby we first determine the opto-acoustic coupling $\mathrm{D}^{(1)}$ from $ab$-initio values of the elastic constants. In case of 2D h-BN we start from Eq. (14), insert the numerical values for $\mathrm{c}_{11}\mid_{is}$ (Table 1), the optical mode frequency $\omega_{\mathrm{E}^{\prime}}$=1310 cm*-1*, the surface mass density $\rho_{2D}$=7.59$\times$10*-8* g/cm2, and obtain the value for $\mathrm{D}^{\mathrm{NB}(1)}_{11,1}$ quoted in Table 2. Here we have retained the negative value of $\mathrm{D}^{\mathrm{NB}(1)}_{11,1}$ which is consistent with force constant model calculations of the dynamical matrix10 . Notice that we have taken into account that the direction of the $x$-axis in the present paper is different from Refs.10 ; 11 . Likewise we proceed for 2H-MoS2 starting from Eq. (17), where we have used $c_{11}\mid_{is}$ from Table 1 and $\omega_{\mathrm{E}^{\prime}}$=390 cm*-1* and $\rho_{2D}$=3.14$\times$10*-7* g/cm2. See $\mathrm{D}^{\mathrm{..}(1)}_{11,1}$ in Table 2. The large difference in absolute value between $\mathrm{D}^{\mathrm{NB}(1)}_{11,1}$ and $\mathrm{D}^{\mathrm{SMo}(1)}_{11,1}$ is due to the fact that the interatomic forces in 2D h-BN are considerably stronger than in 2H-MoS2, in accordance with the different phonon spectra for h-BN10 ; serrano and MoS2molina ; waka . We next turn to Eqs. (15) and (18), insert the corresponding values for $e_{1,11}\mid_{is}$ from Table 1 as well as $\mathrm{D}^{\mathrm{\cdot\cdot}(1)}_{11,1}$ from Table 2, the corresponding frequencies $\omega_{\mathrm{E}^{\prime}}$, masses and densities, and solve with respect to $e^{\ast}_{\mathrm{M}}$. The results are shown for $e^{\ast}_{\mathrm{B}}$ and $e^{\ast}_{\mathrm{Mo}}$ in column 5 of Table 2. In column 6 we have quoted values of $e^{\ast}_{\mathrm{M}}$ obtained directly by DFPT calculations. Although the magnitude of $e^{*}_{\mathrm{Mo}}$ depends strongly on the values of the input quantity e_{1,11}$$\mid_{ri} (see Table 1), we obtain $e^{*}_{\mathrm{Mo}}$ $<$ 0 and $e^{*}_{\mathrm{B}}$ $>$ 0, in agreement with results from ab-initio calculations 25 . We attribute the positive value of $e^{*}_{\mathrm{B}}$ to the larger electronegativity 3.0 of N in comparison with 2.0 of B (see Ref.atkins, ), notwithstanding that there is an opposite electron transfer from N to B within the $\pi$ bond. In case of MoS2 we may assume that the 4$d^{5}$5$s^{1}$ valence electrons of Mo participate in the bonds with the six surrounding S atoms, each of which has valence configuration of 3$s^{2}$3$p^{4}$. Thereby the shielding effect at Mo is decreased and the excess electrons of S lead to an effective negative charge $e^{*}_{\mathrm{Mo}}$. Obviously the large difference between the values -0.41*(1)* and -1.21*(2)* of $e^{*(a)}_{\mathrm{Mo}}$ is due to the different values of $e_{1,11}\mid_{is}$, 0.58 and 1.73 obtained respectively from Refs.12, and 16, , see Table  1.

In Table 3 we present results of similar calculations for Cr based dichalcogenides and some TMDOs. Although these materials have not been synthesized so far, their mechanical and dynamical stability have already been shown by first principles calculationsmx2-deniz .

From Table 2 and 3 if follows that for all TMDCs, notwithstanding quantitative differences, the effective charges $e^{*(a)}_{\mathrm{M}}$ and $e^{*(b)}_{\mathrm{M}}$ obtained by the present analytical calculations or by DFPT respectively, are negative. Notice also that the absolute values increase in the order X=S, Se, Te, in accordance with the charge number Z of these elements, see Table 4. For the same chalcogene, $e^{*}_{\mathrm{M}}$ decreases in absolute value in the order M= Cr, Mo, W, and we conclude that the decrease of shielding mentioned before for Mo is less (more) efficient for heavier (lighter) metals. Our results are in quantitative agreement with recent Born effective charge tensor calculations for MoS2, MoSe2, WS2 and WSe2 nw2 .

For TMDOs the $e^{*}_{\mathrm{M}}$ values obtained by analytical theory and direct ab-initio calculations differ by an order of magnitude, also the values of $e^{*}_{\mathrm{M}}$ obtained by $ab$-$initio$ are all positive. Obviously, the large electronegativity 3.4 atkins of O leads to a negative internal strain contribution to $e_{1,11}$ as is the case for 2D h-BN. The large quantitative difference between analytical and $ab$-$initio$ calculations suggests that the concept of ”effective” ionic chargew used in the analytical theory gives fair results for TMDCs but breaks down for TMDOs. In the latter case the deformability of the electron cloud19 of the metal ions upon inner displacements should not be neglected while in Eqs. (15), (18) and (20) we have assumed a rigid ion model.

We observe that the opto-acoustic coupling $\mathrm{D}^{\mathrm{XM(1)}}_{11,1}$ given in column 4 of Tables 2 and 3 are a measure of the X-M bond strength. We see that for a given metal M this quantity decreases in absolute value with increasing size (i.e. charge number Z) of the X atom, likewise for a given X, $\mathrm{D}^{\mathrm{XM(1)}}_{11,1}$ decreases with increasing size of the metal atom, see Table 5.

We have also carried out calculations of $e_{1,11}\mid_{ir}$ and $e_{1,11}\mid_{ci}$ by using density functional perturbation theory instead of the Berry phase method and then determined $e^{*}_{\mathrm{M}}$. The results are quoted in Table  6. Although the overall agreement with directly calculated values $e^{*}_{\mathrm{M}}$ is less satisfactory than in the case of Tables 2 and 3, the material dependent trends are similar.

V Conclusion

We have presented a synthesis of lattice dynamical theory and $ab$-$initio$ calculations results for the description of elastic and piezoelectric properties of two dimensional ionic crystals with hexagonal lattice structure. As specific examples we investigated 2D h-BN as well as 2H-TMDCs and 2H-TMDOs of composition MX2 where M is a transition metal ion and X a chalcogen or oxygen ion. Such a study allowed us to separate quantitatively electronic and ionic contributions to the elastic and piezoelectric constants. We have investigated the validity of microscopic concepts such as rigid ion model and effective ionic charges for various MX2 compounds.

Further we have been able to discern trends in the values of the opto-acoustic coupling $\mathrm{D}^{\mathrm{XM}(1)}_{11,1}$ and of the optical mode frequency $\omega_{\mathrm{E}^{\prime}}$ as function of the atomic composition. For a given metal M, $\mathrm{D}^{\mathrm{XM}(1)}_{11,1}$ and $\omega_{\mathrm{E}^{\prime}}$ decrease in absolute value with increasing atomic number Z of the chalcogen X ion, and similarly for a given X, these quantities decrease with increasing Z of M. These properties reflect the strength of the chemical bonds of the M ion to the 6 surrounding X ions. We find that the effective charge $e^{*}_{\mathrm{M}}$ of the metal ion is negative for all TMDCs while $e^{*}_{\mathrm{B}}$ in 2D h-BN is positive. This difference in sign entails that the inner strain contribution e_{1,11}$$\mid_{is} leads to a reduction of $e_{1,11}$ in case of 2D h-BN and to an increase in case of 2H-MoS2 and other TMDCs.

Acknowledgments. The authors acknowledge useful discussions with L. Wirtz, and A. Molina-Sanchez. This work was supported by the Methusalem program, and the Flemish Science Foundation (FWO-Vl). Computational resources were provided by HPC infrastructure of the University of Antwerp (CalcUA) a division of the Flemish Supercomputer Center (VSC), which is funded by the Hercules foundation.