Flexoelectric and Piezoelectric Coupling in a Bended MoS₂ Monolayer

Abstract

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. The tuning of the electronic properties of the LD-TMDs are commonly related with applied strains and strain gradients, which can strongly affect their polar properties via piezoelectric and flexoelectric couplings. Using the density functional theory and phenomenological Landau approach, we studied the bended 2H-MoS₂ 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. Our results revealed that the contribution of the flexoelectric effect dominates over the piezoelectric effect in both in-plane and out-of-plane directions of the monolayer. In accordance with our calculations, a realistic strain gradient of about 1 nm⁻¹ can induce an order of magnitude higher than the flexoelectric response in comparison with the piezoelectric reaction. The value of the dilatational flexoelectric coefficient is almost two times smaller than the shear component. It appeared that the components of effective flexoelectric and piezoelectric couplings can be described by parabolic dependences of the corrugation. Obtained results are useful for applications of LD-TMDs in strain engineering and flexible electronics.

1. 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 shape of the Fermi surface [1,2]. On transition from the bulk to the nanoscale, additional long-range orderings and physical properties, such as piezoelectric, or ferroelectric [3,4,5], semiconductive, semi-metallic 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 MX₂ (M—metal Mo, V, W; X—chalcogen S, Se, Te) [10,11] and Janus-compounds with a chemical formula MXY (X, Y—chalcogens) [12,13] vary 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, MX₂ 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 MX₂. Subsequently, Song et al. [25] observed a room temperature semiconductor−metal transition in thin MoTe₂ films induced by a homogeneous tensile strain of 0.2%.

A number of first-principle studies explored surface-induced piezoelectricity [26,27] and ferroelectric polarization [5,28] in various MX₂ and MXY. Subsequently, 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 MX₂ 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 electronic and polar properties can be similar to the ones in a bended graphene [30,31] and boron nitride [32]; moreover, the flexoelectric effect plays an important role [33,34]. It was predicted theoretically that the bending can induce an out-of-plane electric dipole moment of carbon nanoshells [35]. As anticipated, the bending-induced dipole moment is curvature-dependent, and strong curvatures can lead to a relatively high polarization of LD-TMDs, induced by a flexoelectric effect [36,37]. The bending-induced out-of-plane dipole moment with density $p~\left(0.01-0.4\right)$ C/nm and flexoelectric polarization $P~\left(0.1-2\right)\mu$C/cm² were calculated from the first principles for MoS₂ [5], WS₂ [35], and WTe₂ [38,39] single-layers. The flexoelectricity is determined by the lattice deformation in dielectrics and wide gap ferroelectrics. In contrast to this, centrosymmetric narrow gap semiconductors and semimetals, such as MX₂, contain a significant contribution from the deformed electronic density [32,40,41], and thus the electronic contribution to their flexoelectric response can dominate over the lattice-mediated ionic contribution, especially under photoexcitation [42,43].

The tunable out-of-plane piezoelectricity and enhanced conductivity, both induced by flexoelectricity, were observed by Kang et al. [23] in semiconducting 2H-MoTe₂ flakes by creating surface corrugation. The experimental results were corroborated by their ab initio calculations [23], analytical calculations performed within LGD approach [44], and finite element modeling [45]. Specifically, the 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. Finite element modeling allows calculating the elastic and electric fields, flexoelectric polarization, and its correlation with free charge density for a TMD nanoflake placed on a rough substrate with a sinusoidal profile of the corrugation [45].

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 important for LD objects, their knowledge is required to control and predict the physical properties of LD-TMDs 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-MoS₂ 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.

2. Theoretical Formalism

The static electric polarization of a LD-TMD, $P_i(\overrightarrow{x})$, has the form:

$$P_i(\overrightarrow{x})\cong f_{ijkl}\frac{\partial u_{kl}}{\partial x_j}+\frac{e_{ijk}^s}{t}{u}_{jk}-{\epsilon }_0{\chi }_{ij}\frac{\partial \phi }{\partial x_j}.$$

(1)

Here, $f_{ijkl}$ is the static flexoelectric tensor [46] determined by the microscopic properties of the material [47,48], $u_{ij}(\overrightarrow{x})$ is the elastic strain tensor, $e_{ijk}^s$ is the tensor of the surface-induced piezoelectric effect [49,50], and $t$ is the thickness of a nanoflake. The last term in Equation (1) is proportional to the gradient of electric potential $\phi (\overrightarrow{x})$, ${\epsilon }_0$ is a universal dielectric constant, and ${\chi }_{ij}$ is a real part of the TMD dielectric susceptibility. Einstein summation over repeated indexes is used hereinafter.

2.1. Ab Initio Calculations of the Flat and Corrugated MoS₂ Nanolayer

We performed the calculations of the atomic position in the MoS₂ nanolayer for the case of a flat layer and for 1–10% corrugated layers with a 1% step in corrugation within the DFT [51] in the generalized gradient approximation, implemented in the Quantum Espresso code [52]. We used ultrasoft Perdew-Burke-Ernzerhof pseudopotentials [53], which include 14 valence electrons for molybdenum and six valence electrons for sulfur. An integration of the Brillouin zone was performed using 1 × 14 × 1 k-points mesh centered on Γ in the Brillouin zone, generated by Monkhorst-Pack scheme [54], and Methfessel-Paxton smearing [55] 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 MoS₂ 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 (Figure 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⁻⁴ atomic units, and at this point, all the necessary quantities (atomic coordinates and wave functions) were extracted. The difference between the initial corrugation and that obtained for the corresponding relaxed system does not exceed few percentages of their values. Therefore, the calculations revealed that the bending in the z-direction is sinusoidal along the x-axis (Figure 2).

In addition to the position of the atoms, the Bader charges [56] were computed. The distribution of the charges appeared to be uniform in a flat system and equal to $Q_{Mo}=1.083e$ for molybdenum atoms and is $Q_S=-0.5415 e$ for sulfur atoms, where $e=1.60217662·{10}^{-19}$C is the elementary charge. For strongly curved layers with a bending of 8–10%, a harmonic dependence of the charges on the x coordinate was observed (as shown in Figure 3). However, the difference between the maximal and minimal charges appeared to be small (see the scale in Figure 3) and close to the magnitude of expected numerical error. The effect of such a distribution on the electrophysical parameters of the layer is neglected hereinafter. This means that we neglect the terms proportional to the deformation potential in Equation (1), as the first approximation, and extract the terms proportional to the net flexoelectric coefficient, $f_{ijkl}$, from the DFT results. The charge variations can be found as the second approximation, which becomes important for large curvatures only [32].

To the best of our knowledge, Figure 3 calculated for MoS₂ has only qualitative analogue in the literature. Namely, the redistribution of Bader charges caused by the bending of phosphorene [33] and graphene [34] nanoribbons was calculated, and it appeared that the symmetries of Bader charges are broken under the applied strain gradient. This symmetry breaking is almost absent for the studied MoS₂ (compare Figure 3 in this work with Figure 2 from [33] and Figure 4 from [34]), but evidently present for the phosphorene and graphene nanoribbons. The possible reason is the small discrepancy in the net Bader charges of sulfur in a bended MoS₂ monolayer, which are smaller than the charge discrepancy in phosphorene [33] and graphene [34].

2.2. Determination of the Flexoelectric Coefficients from Ab Initio Calculations

The components of the piezoelectric and flexoelectric tensors for the 2H-MoS₂ that belong to the point symmetry group 6/mmm can be determined from the equations:

$$d_z(x)=e_{zzx}{u}_{zx}(x)+e_{zxx}{u}_{xx}(x)+f_{zxzx}{u}_{zx,x}(x)+f_{zxxx}{u}_{xx,x}(x),$$

(2)

$$d_x(x)=e_{xzx}{u}_{zx}(x)+e_{xxx}{u}_{xx}(x)+f_{xxxz}{u}_{xz,x}(x)+f_{xxxx}{u}_{xx,x}(x),$$

(3)

where $d_i$ is projection of the elementary dipole moment on the axis i; $u_{ij}$ and $u_{ij,k}$ are strain tensor and strain gradient, respectively; $e_{ijk}$ and $f_{ijkl}$ are the components of the piezoelectric and flexoelectric tensors in Cartesian coordinates.

All nonzero components of $d_i$, $u_{ij}$ and $u_{ij,k}$ were calculated in Mathematica 12.2 ^®Wolfram Research notebook [57] based on the DFT results. In particular, the effective elementary dipole moment $\overrightarrow{d}$ of each dipole (see Figure 4) was calculated from the following equation:

$${\overrightarrow{d}}_i={\overrightarrow{r}}_{{\mathrm{Mo}}_i}\cdot \left|Q_{\mathrm{Mo}}\right|-{\displaystyle \sum }_{k=1}^6{\overrightarrow{r}}_{{\mathrm{S}}_{ik}}\cdot \left|\frac{Q_{\mathrm{S}}}{3}\right|,$$

(4)

where ${\overrightarrow{r}}_{{\mathrm{Mo}}_i}={\{x_{\mathrm{Mo}},y_{\mathrm{Mo}},z_{\mathrm{Mo}}\}}_i$ and ${\overrightarrow{r}}_{S_{ik}}={\{x_S,y_S,z_S\}}_{ik}$ is a radius vector of the Mo and S atoms, respectively.

Further dependences were plotted on the x coordinate of the mass center of the dipole, and its radius vector was determined by the equation:

$${\overrightarrow{r}}^{(j,i)}=\frac{1}{m_{\mathrm{Mo}}+6m_{\mathrm{S}}}({\overrightarrow{r}}_{\mathrm{Mo}}^{(j,i)}{m}_{\mathrm{Mo}}+{\displaystyle \sum }_{ k=1}^6{\overrightarrow{r}}_{{\mathrm{S}}_k}^{(j,i)}{m}_{\mathrm{S}}),$$

(5)

where “$j$” is the percentage of corrugation, “$i$” is the dipole number, $m_{Mo}=95.9\frac{\mathrm{g}}{\mathrm{mol}}$ and $m_S=32 \frac{\mathrm{g}}{\mathrm{mol}}$ are atomic masses of molybdenum and sulfur, respectively. The points in Figure 5b and Figure 6b show the dependence of x-axis and y-axis projection of the dipole moment on x-coordinate of the dipole mass center.

Using the radius vector, the elementary displacement of the mass centers is calculated from the equation:

$${\overrightarrow{U}}^{(j,i)}={\overrightarrow{r}}^{(0,i)}-{\overrightarrow{r}}^{(j,i)},$$

(6)

where ${\overrightarrow{r}}^{(0,i)}$ is the radius vector of i-th dipole of the undeformed layer, ${\overrightarrow{r}}^{(j,i)}$ 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 ${\overrightarrow{U}}^{(j)}$ on x coordinate of the mass center (5) are shown in Figure 5a and Figure 6a, respectively.

Solid curves in Figure 5a,b and Figure 6a,b are harmonic (sinusoidal and co-sinusoidal) functions, which were selected in the form:

$$d_z=d_{zc}\cos (kx)+d_{zs}\sin (kx),$$

(7a)

$$d_x=d_{xs}\sin (2kx)+d_{xc}\cos (2kx)+c.$$

(7b)

Here, the coefficients $d_{zc}$, $d_{zs}$, $d_{xs}$, $d_{xc}$ and $c$ are fitting parameters. Theirs numerical values are selected to reach the best fit of the points in Figure 5b and Figure 6b.

The dependence of the displacement projections on x-axis was approximated by the functions:

$$U_z(x)\approx U_0\cos (kx),$$

(8a)

$$U_x(x)\approx V_{0}x+W_0\cos (2kx).$$

(8b)

Here, the coefficients $U_0$, $V_0$, and $W_0$ are fitting parameters. Theirs numerical values are selected to reach the best fit of the points in Figure 5a and Figure 6a.

The displacement causes the strain tensor $u_{jk}=\frac{1}{2}(\frac{\partial U_k}{\partial j}+\frac{\partial U_j}{\partial k})$, whose components (shown in Figure 5c and Figure 6c) are given by the following equations:

$$u_{zx}(x)=u_{xz}(x)=\frac{1}{2}\frac{\partial U_z}{\partial x}=-\frac{U_0}{2}k\sin (kx),$$

(9a)

$$u_{xx}(x)=\frac{\partial U_x}{\partial x}=V_0+2W_{0}k\cos (2kx).$$

(9b)

Their nonzero derivatives (shown in Figure 5d and Figure 6d) are equal to:

$$u_{zx,x}(x)=u_{xz,x}(x)=-\frac{U_0}{2}{k}^2\cos (kx),$$

(10a)

$$u_{xx,x}(x)=\frac{\partial u_{xx}}{\partial x}=-4W_0{k}^2\sin (2kx).$$

(10b)

The values $d_{zs}$, $d_{zc}$, $d_{xs}$, $d_{xc}$, $c$, $U_0$, $V_0$, $W_0$ and $k$ in Equations (6)–(8) are constants, which are different for each layer. The values of these coefficients are given in

Table 1. Fitting parameters.

Corrugation	1%	2%	3%	4%	5%	6%	7%	8%	9%	10%
$k$, 1/nm	1.44	1.441	1.443	1.446	1.449	1.453	1.457	1.463	1.468	1.475
$d_{xs}$, 10−3e·nm	0.006	0.028	0.059	0.107	0.191	0.301	0.432	0.553	0.678	0.826
$d_{xc}$, 10−3e·nm	0.0001	−0.006	−0.004	−0.059	−0.062	−0.048	−0.039	−0.078	−0.142	−0.166
$c$, 10−3e·nm	−0.997	−0.979	−0.979	−0.982	−0.978	−1.018	−1.058	−1.043	−1.002	−1.017
$d_{zs}$, 10−3e·nm	0.066	0.122	0.174	0.325	0.341	0.382	0.448	0.530	0.596	0.67
$d_{zc}$, 10−3e·nm	0.481	0.959	1.441	1.906	2.389	2.867	3.334	3.797	4.258	4.707
$V_0$, 10−3 a.u.	−0.238	−0.97	−2.193	−3.905	−6.1	−8.775	−11.93	−15.55	−19.62	−24.16
$W_0$, 10−3 nm	0.02	0.0258	0.0449	1.296	2.06	2.883	3.852	5.019	6.356	7.791
$U_0$, 10−3 nm	20.71	41.35	61.94	82.76	103.3	123.7	144	164.1	184.2	204

, and the dependences given by Equations (9) and (10) are shown in Figure 5c,d and Figure 6c,d.

Analyzing the dependences shown in Figure 5a, we can conclude that the projection of the atomic displacements on the z-axis at the top and bottom of corrugation profile are maximal and minimal and equal to each other by the absolute value. This is due to the created curvature of the MoS₂ layer. Moreover, the bending is caused not only by a shift of unit cells relative to the initial state, but also by their deformation. Thus, the S atoms, which are inside the bend line, shift closer to each other under the action of mechanical forces, and accordingly they are closer to Mo atoms. Therefore, the S atoms outside the bend repulse and, as a result, move away from Mo atoms. Since the opposite electric charges are moving away from each other, the dipole moment in this direction also increases, which can be seen from Figure 5b.

The linear component of the projection of the offset on the x-axis (shown in Figure 6a) is explained by the conservation of the total length of the unit cell, and the projection of the length decreases in the x-direction. The presence of a harmonic component is explained by the fact that atoms located exactly on the top and bottom of the corrugation profile are exposed to electrostatic forces only in the z-direction and remain in place in the x-direction (excluding the linear shift). In addition, those atoms that are in places, where there is a change in the bend of corrugation, do not have a shift in either in the y-direction or in the x-direction (excluding the linear shift). All other atoms change their position under the action of electrostatic and elastic forces in order to achieve equilibrium. Based on the amplitude of the dipole moment (shown in Figure 6b), we can conclude that the deformation of the unit cell in x-direction is negligible, and the presence of the dipole moment in the absence of deformation is explained by the fact that the calculated dipole moment has the third order symmetry axis and, accordingly, is not symmetric in the x-direction in the Cartesian coordinate frame.

Substituting the functions from Equations (7)–(10) to Equations (2) and (3), we obtain:

$$\begin{array}{c}{d}_{zc}\cos (kx)+d_{zs}\sin (kx)=-e_{331}\frac{U_0}{2}k\sin (kx)+e_{311}\left[V_0+2W_{0}k\cos (2kx)\right]-\\ -f_{3131}\frac{U_0}{2}{k}^2\cos (kx)-f_{3111}4W_0{k}^2\sin (2kx),\end{array}$$

(11a)

$$\begin{array}{c}{d}_{xs}\sin (2kx)+d_{xc}\cos (2kx)+c=-e_{131}\frac{U_0}{2}k\sin (kx)+e_{111}\left[V_0+2W_{0}k\cos (2kx)\right]-\\ -f_{1113}\frac{U_0}{2}{k}^2\cos (kx)-f_{1111}4W_0{k}^2\sin (2kx),\end{array}$$

(11b)

The components of flexoelectric and piezoelectric tensors can be obtained from the above expressions in the following form:

$$f_{3131}=-\frac{2d_{zc}}{U_0{k}^2},\hspace{1em}{e}_{331}=-\frac{2d_{zs}}{U_{0}k},\hspace{1em}{e}_{311}=0,\hspace{1em}{f}_{3111}=0.$$

(12a)

$$f_{1111}=-\frac{d_{xs}}{4W_0{k}^2},\hspace{1em}{e}_{111}=\frac{d_{xc}}{2W_{0}k},\hspace{1em}{e}_{131}=0,\hspace{1em}{f}_{1113}=0.$$

(12b)

The values of $f_{3131}$, $f_{1111}$, $e_{111}$ and $e_{331}$, calculated according to Equation (12), considering the values of Table 1, are given in

Table 2. The components of flexoelectric and piezoelectric coefficients.

Curvature, %	${\mathit{f}}_{3131}, \mathit{e}/\mathbf{n}{\mathbf{m}}^2$	${\mathit{e}}_{331}, \mathit{e}/\mathbf{n}\mathbf{m}$	${\mathit{f}}_{1111}, \mathit{e}/\mathbf{n}{\mathbf{m}}^2$	${\mathit{e}}_{111}, \mathit{e}/\mathbf{n}\mathbf{m}$
1	−0.0224	−0.00443	−0.0338	0.00163
2	−0.0223	−0.00409	−0.131	−0.0766
3	−0.0223	−0.00389	−0.156	−0.0292
4	−0.022	−0.00544	−0.00983	−0.0159
5	−0.022	−0.00456	−0.011	−0.0103
6	−0.022	−0.00425	−0.0124	−0.00573
7	−0.0218	−0.00427	−0.0132	−0.00352
8	−0.0216	−0.00442	−0.0129	−0.00529
9	−0.0214	−0.00441	−0.0124	−0.00759
10	−0.0212	−0.00445	−0.0122	−0.00721

.

Dependences of the flexoelectric and piezoelectric tensors components on the corrugation of the MoS₂ monolayer, which are determined from the DFT calculations, are shown in Figure 7. Let us underline the dominant contributions of the flexoelectric effect over the piezoelectric effect in both x and z directions, namely $\frac{\left|f_{ijkl}\right|}{\left|e_{ijkl}\right|}\cong (5-30)\frac{1}{nm}$. This means that a realistic strain gradient of about $\frac{1}{nm}$ can induce an order of magnitude higher flexoelectric response in comparison with a piezoelectric reaction. In addition, the dilatational flexoelectric coefficient $\left|f_{1111}\right|$ is almost two times smaller than the shear component $\left|f_{1313}\right|$ for all corrugations of more than 1%; moreover, both values, $f_{1111}$ and $f_{1313}$, are negative (Figure 7a). The inequality $\left|e_{331}\right|>\left|e_{111}\right|$ is valid for all corrugations of more than 1% except for 7% (see Figure 7b). All components of the flexoelectric and piezoelectric tensors can be described by an empirical parabolic dependence, $A{\mu }^2+B\mu +C$ of the corrugation “$\mu$” (see Figure A1 in Appendix A). Components $f_{3131}$ and $e_{331}$ weakly depend on corrugation magnitude and slightly and linearly change with corrugation around a constant value. Components $f_{1111}$ and $e_{111}$ have an extremum at the corrugation 7% and reach high absolute values at small corrugations $\mu <4\%$. These unreasonably high values originated from the fitting error of the dipole moment projection $d_x(x)$ (see Figure A2 in Appendix A); thus, these points should be disregarded.

To the best of our knowledge, experimental results, which can reliably determine the flexoelectric and piezoelectric properties of a monolayer 2H-MoS₂, are lacking. Kang et.al. [23] studied MoTe₂ nanoflakes on a corrugated gold substrate and have shown that the effective piezoelectric response linearly depends on the tip bias (Figure 3 in ref. [23]). They also analyzed different contributions to the dependence of the effective piezoresponse on the flake thickness and concluded that the flexoelectric contribution dominates (Figure 4 in ref. [23]). These results are in a complete qualitative agreement with our calculations for a monolayer MoS₂.

In addition, there are few DFT calculations of the piezoelectric coefficients in a flat centrosymmetric MX₂. In particular, $e_{31}^s$ was estimated as several −2.4 pC/m for a reconstructed SL-MoS₂ in the 1T’ phase [5]. The value is much smaller than in-plane components, $e_{1jk}^s$~(300–500) pC/m, estimated for H-MoS₂ and MoSTe (see Table 1 in Ref. [28]). The values of $e_{331}^s$, calculated in this work, are around −0.7 pC/m, which is several times smaller than −2.4 pC/m [5]. The values of $e_{111}^s$, calculated in this work, vary from −1 to −12 pC/m for different corrugations, which is one to two orders smaller than the values listed in Table 1 from Ref. [28]. The discrepancy can be related with different symmetry (monolayer 2H-MoS₂ via H-MoS₂) and, more likely, with different configurations of the layer (flat configuration via bended layer). The components of the flexoelectric (in mC/m²) and piezoelectric (in mC/m²) coefficients, recalculated from Table 2, are listed in

Table A1. Components of flexoelectric (in mC/m²) and piezoelectric (in pC/m) coefficients, recalculated from Table 2.

Curvature, %	${\mathit{f}}_{3131}, \mathbf{m}\mathbf{C}/{\mathbf{m}}^2$	${\mathit{e}}_{331}, \mathbf{p}\mathbf{C}/\mathbf{m}$	${\mathit{f}}_{1111}, \mathbf{m}\mathbf{C}/{\mathbf{m}}^2$	${\mathit{e}}_{111}, \mathbf{p}\mathbf{C}/\mathbf{m}$
1	−3.59	−0.71	−5.42	0.261
2	−3.58	−0.655	−21	−12.3
3	−3.58	−0.624	−25.1	−4.67
4	−3.53	−0.871	−1.58	−2.55
5	−3.53	−0.731	−1.77	−1.65
6	−3.52	−0.681	−1.98	−0.918
7	−3.49	−0.685	−2.11	−0.563
8	−3.47	−0.708	−2.06	−0.848
9	−3.44	−0.706	−1.98	−1.22
10	−3.4	−0.713	−1.95	−1.15

in the end of Appendix A.

3. Conclusions

Using the phenomenological Landau approach combined with DFT ab initio calculations, we consider the atomic displacements and charge state of the curved 2H-MoS₂ monolayer and determine the effective components of its flexoelectric and piezoelectric coupling tensors, $f_{ijkl}$ and $e_{ijk}$, as a function of the layer corrugation varying from 0% to 10%.

It appears that the components of effective flexoelectric and piezoelectric couplings can be described by parabolic dependences of the corrugation. In particular, the components $f_{3131}$ and $e_{331}$ weakly depend on the corrugation magnitude and slightly change with the corrugation around an almost constant value. The components $f_{1111}$ and $e_{111}$ have an extremum at the corrugation 7% and reach unreasonably high absolute values at small corrugations $<4\%$, which are related with the fitting error of the dipole moment projection.

Our calculations reveal the dominant contributions of the flexoelectric effect over the piezoelectric effect in both in-plane and out-of-plane directions of the monolayer. This means that a realistic strain gradient of about 1 nm⁻¹ can induce an order of magnitude higher flexoelectric response in comparison with a piezoelectric reaction. The dilatational flexoelectric coefficient $f_{1111}$ is almost two times smaller than the shear component $f_{1313}$, while $\left|e_{331}\right|>\left|e_{111}\right|$ for the most of corrugations.

Obtained quantitative results can be useful for elaboration of nanoscale flexible electronic devices based on the bended MX₂ layers. In particular, the bended monolayers are promising candidates for the ultra-small diodes and bipolar transistor on MX₂, for which the principal schemes are presented in Ref. [45]. 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.