Flexoelectric and piezoelectric coupling in a bended MoS$_2$ monolayer

Low-dimensional (LD) transition metal dichalcogenides (TMDs) in the form of nanoflakes, which consist of one or several layers, are the subject of intensive fundamental and applied research. Due to the size-induced transition from a bulk to nanoscale, they can be both nonpolar, piezoelectric or even ferroelectric. Also, in terms of electronic properties, they can be direct-band semiconductors, semi-metals or even metals. The tuning of the electronic properties in the LD-TMDs are commonly related with applied strains and strain gradients, which can affect strongly their polar properties via the piezoelectric and flexoelectric couplings. Using the density functional theory (DFT) and phenomenological Landau approach, we studied the bended 2H-MoS$_2$ monolayer and analyzed its flexoelectric and piezoelectric properties. The dependences of the dipole moment, strain and strain gradient on the coordinate along the layer were calculated. From these dependences the components of the flexoelectric and piezoelectric tensors have been determined and analyzed. Obtained results are useful for applications of LD-TMDs in strain engineering and flexible electronics.


I. INTRODUCTION
Layered transition metal dichalcogenides (TMDs) in the form of bulk materials are typically non-polar centrosymmetric semiconductors with a relatively wide band gap ~(1.1 -2) eV and specific form of Fermi surface [1,2]. On transition from the bulk to the nanoscale additional long-range orderings and physical properties, such as piezoelectric, or even ferroelectric [3,4,5], semiconductive, semimetallic or metallic [6,7,8] were found in different structural phases (polymorphs) of TMD monolayers [9]. In particular, the polar and semiconducting properties of the low-dimensional (LD) TMDs with a chemical formula MX2 (M − metal Mo, V, W; X − chalcogen S, Se, Te) [10,11] and Janus-compounds (JC) with a chemical formula MXY (X, Ychalcogens) [12,13] are varying from non-polar to ferroelectric state, and from direct-band semiconductor to metallic conductivity.
The strain and strain gradient impact on LD-TMD polar and electronic properties can be principally important for the properties control, and therefore LD semiconductor materials, such as graphene, MX2 and MXY monolayers, are ideal candidates for the strain engineering [14] and recently introduced "straintronics" [15]. Their strain-induced conductive domain walls can act as mobile charged channels, similarly to the "domain wall nanoelectronics" in multiferroic thin films [16,17,18] and graphene-on-ferroelectric nanostructures [19,20]. There are remarkable possibilities for tuning the structural, polar and electronic properties of LD-TMDs by application of either homogeneous elastic strains [21,22] or inhomogeneous curvature-induced strain gradients [23,24]. In particular, Duerloo et.
al. [6] predicted a strain-induced phase transition from a semiconducting 2H to a metallic 1T' phase in various MX2. Later on, and Song et. al. [25] observed a room temperature semiconductor−metal transition in thin MoTe2 films induced by a homogeneous tensile strain of 0.2%.
A number of first-principles studies explored the surface-induced piezoelectricity [26,27] and ferroelectric polarization [5,28] in various MX2 and MXY. Later on, the possible mechanism of the ferroelectric state appearance in LD-TMDs was described by the Landau-Ginzburg-Devonshire (LGD) continuous approach [29]. Since only in-plane ferroelectricity can exist in a geometrically flat pure centrosymmetric MX2 layer, LGD analysis suggests that the switchable out-of-plane ferroelectric polarization emerges due to the rare-earth doping, and predicts that the domain walls in LD-TDMs should become conductive above a certain strain threshold.
The physical origin of the bending-induced changes of LD-TMDs in electronic and polar properties can be similar to the ones in a bended graphene [30,31] and boron nitride [32], at that the flexoelectric effect plays a very important role. Indeed, it was predicted theoretically, that the bending can induce an out-of-plane electric dipole moment of carbon nanoshells [33]. As anticipated, the bending-induced dipole moment is curvature-dependent, and strong curvatures can lead to a relatively high polarization of LD-TMDs, in fact induced by a flexoelectric effect [34,35]. Actually, the bendinginduced out-of-plane dipole moment with density ~(0.01 − 0.4)C/nm and flexoelectric polarization ~(0.1 − 2) C/cm 2 were calculated from the first principles for MoS2 [5], WS2 [35], and WTe2 [36,37] single-layers. In contrast to dielectrics and wide gap ferroelectrics, where the flexoelectricity is determined by the lattice deformation, in centrosymmetric narrow gap semiconductors and semimetals, such as MX2, contains a significant contribution from the deformed electronic density [32,38,39] that, in principle, can dominate over the lattice-mediated ionic contribution especially under photoexcitation [40,41].
The tunable out-of-plane piezoelectricity and enhanced conductivity, both induced by a flexoelectricity, were observed by Kang et al. [23] in semiconducting 2H-MoTe2 flakes by creating the surface corrugation. The experimental results were corroborated by their ab initio calculations [23], analytical calculations performed within LGD approach [42] and finite element modeling (FEM) [43].
Specifically, LGD approach explores the flexoelectric origin of the polarization induced by a spontaneous bending and by inversion symmetry breaking due to the interactions with substrate. FEM allows calculating the elastic and electric fields, flexoelectric polarization, and its correlation with free charge density for a TMD (or JC) nanoflake placed on a rough substrate with a sinusoidal profile of the corrugation [43].
Despite the progress, the complete information about the piezoelectric and flexoelectric coupling tensors in LD-TMDs is lacking, and the influence of the tensor symmetry and numerical values of its components on polar and electronic phenomena in LD-TMDs has not been studied. Since the couplings are This knowledge is required to control and predict the physical properties of LD-TMDs and JC for their novel applications in nanoelectronics and advanced memories. Using the density functional theory (DFT) and phenomenological Landau approach, in this work we consider the curved monolayer of 2H-MoS2 in order to calculate the nonzero components of its flexoelectric and piezoelectric coupling tensors and to analyze their possible dependence on the layer corrugation and surface-induced symmetry lowering.

II. THEORETICAL FORMALISM
The static electric polarization of a LD-TMD, ( ⃗), has the form: Here ( ⃗) is the elastic strain tensor, is the static flexoelectric tensor [44] determined by the microscopic properties of the material [45,46], is the tensor of the surface-induced piezoelectric effect [47,48], and is the layer or nanoflake thickness. The last term in Eq.(1) is proportional to the electric field, = − , where ( ⃗) is an electric potential, 0 is a universal dielectric constant, is a real part of the TMD dielectric susceptibility. Einstein summation over repeated indexes is used hereinafter.

A. Ab initio calculations of the flat and corrugated MoS2 nanolayer
We performed the calculations of the atomic position in the MoS2 nanolayer for the case of a flat layer and for (1 -10)% corrugated layers with a 1% step in corrugation within the DFT [49] in the generalized gradient approximation (GGA), implemented in the Quantum Espresso code [50]. We have used ultrasoft Perdew-Burke-Ernzerhof (PBE) pseudopotentials [51], which include 14 valence electrons for molybdenum and 6 valence electrons for sulfur. An integration of the Brillouin zone has been performed using 14×1×1 k-points mesh centered on Γ in the Brillouin zone, generated by Monkhorst-Pack scheme [52], and Methfessel-Paxton smearing [53] with a parameter of 0.005 Ry. We applied 50 Ry cutoff for smooth part of the wave function and 350 Ry for the augmented charge density to ensure a sufficient convergence of the results.
The MoS2 monolayer was modeled by a supercell approach with 25 Å of vacuum layer added to avoid a coulomb interaction between the periodic images. Initial (flat) layer was built with experimental value of the lattice constant of 3.161 Å and contained 48 atoms (Fig. 1). For the corrugated layers, we reduced the lattice constant of the supercell in x-direction and modulated the atomic z (out-of-plane) coordinates by suitable sinusoidal distribution to match an integer value of corrugation. After that, all the systems were relaxed through all the internal coordinates until the Hellmann-Feynman forces became less than 10 −4 atomic units, and at this point all the necessary quantities (atomic coordinates and wave functions) were extracted. It is worth to note, that the difference between the initial corrugation and that obtained for the corresponding relaxed system does not exceed few percent of their values. Therefore, the calculations revealed that the bending in z-direction is sinusoidal along the x-axis (Fig. 2). In addition to the position of the atoms, the Bader charges [54] were computed. The distribution of the charges appeared to be uniform in a flat system and equal to = 1,083 for molybdenum atoms, and is = -0,5415 for sulfur atoms (in the units of elementary charge e). For strongly curved layers with a bending (8-10)%, a harmonic dependence of the charges on the x coordinate was observed (as shown in Fig. 3). However, the difference between the maximal and minimal charges appeared to be very small (see the scale in Fig. 3) and close to the magnitude of expected numerical error. The effect of such distribution on the electrophysical parameters of the layer is considered neglected hereinafter.
This, in fact, means that we neglect the terms proportional to the deformation potential in Eq.(1), as the first approximation, and extract the terms proportional to the net flexoelectric coefficient, , from the DFT results. The charge variations can be found as the second approximation, which only becomes important for large curvatures [32].  The components of the electrostriction and flexoelectric tensors for the 2H-MoS2, that belongs to the point symmetry group 6/mmm, can be determined from the equations: where is projection of the elementary dipole moment on the axis i; and , are strain tensor and strain gradient caused by the layer bending; and are the components of the piezoelectric and flexoelectric tensors.
All nonzero components, , and , , were calculated based on the DFT results. In particular, the effective elementary dipole moment ⃗ of each dipole (see Fig. 4) was calculated from the following equation: where ⃗ = { , , } is a radius-vector of the corresponding X atoms.
where ⃗ (0, ) is the radius vector of i-th dipole of the undeformed layer, ⃗ ( , ) is the radius vector of i-th dipole of the layer with corrugation j%. The dependences of the x-axis and y-axis projections of the displacement vector ⃗ ⃗⃗ ( ) on x coordinate of the mass center (5) are shown in Fig. 5а and Fig. 6 a, respectively.
Solid curves in Fig. 5a-b and Fig. 6a-b are harmonic (sinusoidal and co-sinusoidal) functions, which were selected in the form: Here the coefficients , , and are fitting parameters, which numerical values are selected to reach the best fit of the points in Fig. 5b and Fig. 6b.
The dependence of the displacement projections on x-axis was approximated by the functions: ( ) ≈ 0 + 0 cos(2 ).
Here the coefficients 0 , 0 , and 0 are fitting parameters, which numerical values are selected to reach the best fit of the points in Fig. 5a and Fig. 6a.
The displacement causes the strain tensor = 1 2 ( + ), whose components (shown in Fig. 5c and Fig. 6c) are given by the following equations: Their nonzero derivatives (shown in Fig. 5d and Fig. 6d) are equal to: The values , , , , , 0 , 0 , 0 and in Eqs.(6)-(8) are constants, which are different for each layer. The values of these coefficients are given in Table I, and the dependences given by Eqs.(9)-(10) are shown in Fig. 5c-d and Fig. 6c- Table. II.     All dependences correspond to the effective elementary dipole mass center. The parameters are given in Table I.
Dependences of the flexoelectric and piezoelectric tensors components on the corrugation of MoS2 monolayer, which are determined from the DFT calculations, are shown in Fig. 7. All components of flexoelectric and piezoelectric tensors can be described by an empirical parabolic dependence, 2 + + of the corrugation " " (see Fig. A1 in Appendix A). Components 3131 and 331 weakly depend on corrugation magnitude, and slightly linearly changes with corrugation around a constant value.
Components 1111 and 111 have an extremum at the corrugation 7%, and reach high absolute values at small corrugations < 4%. These unreasonably high values originated from the fitting error of the dipole moment projection ( ) (see Fig. A2 in Appendix A), so these points should be disregarded. Obtained quantitative results can be useful for elaboration of nanoscale flexible electronic devices based on the bended MX2 layers. In particular, the bended monolayers are promising candidates for the ultra-small diodes and bipolar transistor on MX2, which principal schemes are presented in Ref. [43]. Here the bending profile of the layers controls the sharpness of p-n junctions between the regions with n-type (electron) and p-type (hole) conductivity.