First-principles Simulation of Piezoresistivity of Transition Metal Dichalcogenide Monolayers

Longitudinal gauge factors of transition metal dichalcogenide (TMDC) monolayer models have been evaluated by first-principles calculation. TMDC monolayers have a multivalley/ multipeak electronic structure, and uniaxial tensile strain causes a significant change in electrical conductivity through carrier redistribution in conduction-band valleys or valenceband peaks for some species of doped TMDC monolayers. A high piezoresistivity was observed selectively depending on the combination of transition metal and chalcogen atoms. In this simulation, a p-doped MoS2 monolayer model gave high longitudinal gauge factors of 94.2 and 87.4 at 20 °C in the carrier concentration range of 1017–1019 cm−3.


Introduction
The exploration of novel high-sensitivity piezoresistive materials is indispensable for expanding the usability and operability of force-sensing micro/nano-electromechanical system (MEMS/NEMS) devices. Transition metal dichalcogenide (TMDC) layers have attracted much attention as new candidates of piezoresistive materials composing MEMS mechanical sensors. In particular, TMDC monolayers containing Group VI transition metals, such as molybdenum (Mo) and tungsten (W), have been the subject of a great deal of interest as a post-graphene twodimensional nanomaterial for application to MEMS devices in various fields of electronics, photonics, and spintronics. (1)(2)(3)(4)(5)(6)(7)(8)(9)(10) In contrast to graphene, the TMDC is a typical semiconducting material with a clear band gap, and the effective masses of carrier electrons and holes are nonzero. (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24) From the viewpoint of piezoresistivity, the measurement or evaluation of gauge factors (GFs) of TMDCs is a key issue for future application to mechanical sensors. Recently, the evaluation of GFs of Mo disulfide (MoS 2 ) with 1-3 layers has been reported by the finite-element method with AFM measurement data. (25) On the other hand, the author has already developed the theory and simulation method to evaluate the piezoresistivity of any semiconducting materials on the basis of first-principles calculation. (26)(27)(28) The electronic structures of TMDC monolayers have been investigated by many researchers, and the author has also presented the details of band structures of TMDC monolayers for the simulation of thermoelectric properties. (11) It is well known that the electronic structures of TMDC monolayers containing Group VI transition metals are affected significantly by spin-orbit coupling. (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24) In the analysis of TMDC properties regarding stress or strain such as piezoresistivity, the influence of spin-orbit coupling on stress or strain should be sufficiently incorporated.
In this study, longitudinal GFs in doped TMDC monolayer models MX 2 (M = Mo and W; X = S, Se, and Te) have been simulated in terms of first-principles band structures, and the potential of TMDC monolayers for piezoresistive materials shall be discussed.

First-principles calculation
The electronic band structures of the TMDC monolayer models have been calculated using the VASP program package (29,30) based on the density functional theory (DFT). (31) For DFT exchange-correlation interaction, the generalized-gradient approximation method was used with the Perdew-Wang (PW91) functional. (32,33) The projector augmented wave (PAW) method was adopted, (34) and the cutoff energy for wave functions of electrons was set at 40 Ry (544 eV). Spin-orbit coupling was included in all DFT calculations. The k-point sampling was set by the Monkhorst-Pack grids (35) of 6 × 6 × 1 for geometry optimization and of 48 × 48 × 1 for carrier integration in the Brillouin zone. Figure 1 shows unit cells of the TMDC monolayer models in this study. The configuration of Mo and W atoms is trigonal prismatic for MoX 2 and WX 2 . For the direction vertical to the two-dimensional monolayer, the vacuum space was introduced for establishing a threedimensional periodic boundary condition by setting the fixed vertical cell parameter to 20 Å. Table 1 shows the optimized geometrical parameters a and d shown in Fig. 1 for the TMDC monolayer models, (11) and the volume per unit cell was defined by , where h is the layer thickness and is conveniently regarded as 2d in this study. The effect of uniaxial tensile strain on the TMDC plane was represented by partial optimization with a fixed scale along the tensile direction. In this study, 1% tensile-strained models were introduced by applying the uniaxial tensile strain ε = 0.01 along the [112 __ 0] and [11 __ 00] directions. As reported by Yu et al., (36) Poisson's ratios were obtained as negative values, ranging from −0.02 to −0.10, for all TMDC monolayer models.

Evaluation of gauge factor
In the two-dimensional transport model of doped semiconductors with a multivalley and/or multipeak electronic structure, the electrical conductivity tensor G or electrical resistivity tensor ρ can be represented as a 2 × 2 tensor in terms of carrier densities and effective mass tensors by where n j,α is the jth conduction-band (CB) carrier electron density for the CB valley α, p j,β is the jth valence-band (VB) hole density for the VB peak β, m * j,α , and m * j,β are band effective mass tensors, τ j,α and τ j,β are relaxation time tensors, and e 2 is the square of the absolute value of the elementary electric charge. (37) An excess or lack of electrons as large as one particle per unit cell, in order to carry out regular first-principles calculation in the n-or p-doped semiconductor state, should lead to an enormous overestimation of carrier concentration. In actual doped TMDC monolayers, the total number of carrier electrons or holes per unit cell, δ, must be less by a few orders than 1. Under the condition that a small amount of carrier occupation causes no significant changes in band structure, δ can be given by an appropriate shift in the Fermi energy E F as Table 1 Optimized geometrical parameters for TMDC monolayer models.
in the n-or p-type carrier occupation with the temperature T, respectively, where E j,k is the intrinsic-semiconductor-state band energy of the jth subband at the k point, w k is the k-point weight for k, and k B is the Boltzmann constant. Practically, δ = NV was set to an appropriate constant that corresponds to the carrier concentration N ranging from 1 × 10 17 to 1 × 10 19 cm −3 , and then E F was solved according to Eq. (2) or (3). The values of n j,α and p j,β were respectively determined by regional partitioning as where R α is the region of the CB valley α and R β is the region of the VB peak β. (26) The band effective mass tensors in each region can be defined as 2 × 2 tensors, ( ) at the minimum band energy E j in R α or the maximum E j in R β . ħ is equal to Planck's constant divided by 2π and R is the planar rotation matrix around the k z || [0001] axis. On the righthand side in Eq. (6), a positive sign is adopted for carrier electrons (γ = α) and a negative sign for holes (γ = β). In this paper, the x-and y-directions can be respectively defined as the [112 __ 0] and [11 __ 00] directions, and R becomes a unit matrix in this definition. For the relaxation times, we have adopted the approximation that all of the band relaxation times are equal and constant regardless of strain, (26)(27)(28) because the variation rate of carrier conductivity can be easily and adequately represented by cancelling most of the band relaxation times.
Under the uniaxial tensile strain ε, the longitudinal GF K l is defined as where ρ l 0 is the longitudinal diagonal element of the resistivity tensor ρ without strain and Δρ l is the variation in resistivity in the direction toward the uniaxial tensile strain.

Carrier redistribution by uniaxial strain
The details of the multivalley/multipeak band structure for the TMDC monolayer models have been reported in the author's previous study. (11) The band diagrams in the vicinities of the CB bottom and VB top for each TMDC monolayer model are shown in Fig. 2. For the CB bottom of monolayers, there are deep valleys not only at the K points but at the ΓK points along the Γ-K paths. Valleys at the ΓM points along the Γ-M paths are actually very shallow for the direction perpendicular to the Γ-M paths toward the ΓK valleys, so the occupation of carrier electrons on the ΓM valleys can be ignored. For the VB top, there are peaks at the K and Γ points. Figure 3 shows the regional partitioning to define R α (R K , R ΓK1 , and R ΓK2 ) and R β (R K and R Γ ) for computing n j,α and p j,β for each valley or peak for Eqs. (4) and (5). Six K valleys and six K peaks on the corners of the Brillouin zone respectively remain equivalent with the application of any strain, (27) and they can be treated as one valley and one peak. On the other hand, the uniaxial strain gives a small difference in band energy level among six ΓK valleys.  The temperature dependences of the carrier electron occupation ratios for the strain-free and [112 __ 0] tensile-strained TMDC monolayer models are shown in Fig. 4. For n-doped MoSe 2 and WSe 2 , carrier electrons are distributed on both the ΓK and K valleys, because the ΓK valleys in MoSe 2 lie on almost the same band energy level of the K valleys, and the ΓK valleys clearly lie lower than the K valleys in WSe 2 at the calculation level in this paper. By applying the uniaxial strain, a remarkable redistribution of carrier electrons can be easily triggered from the ΓK valleys to the K valleys in both n-doped MoSe 2 and WSe 2 monolayer models. On the other hand, the redistribution of carrier electrons can hardly be observed in n-doped MoS 2 , MoTe 2 , WS 2 , and WSe 2 . The changes in carrier electron occupation ratio by the tensile strain for MoSe 2 and WSe 2 models can be supported by the band energy diagram for the corresponding tensilestrained models shown in Fig. 5. The relative band energy levels at the ΓK 1 and ΓK 2 valleys to the K valley are raised by 0.05 and 0.03 eV for MoSe 2 and WSe 2 , respectively, owing to the 1% tensile strains.
For the VB of monolayer models, all holes are located at the K point in p-doped TMDCs except for MoS 2 , and the redistribution of holes cannot occur in the presence of any strain. Only in the MoS 2 model, the band energy of the Γ peak is close to that of the K peak, and holes are distributed on both the K and Γ peaks for p-doped MoS 2 . Figure 6 shows the temperature dependences of the hole occupation ratios for the strain-free and [112 __ 0] tensilestrained monolayer models. The uniaxial strain can give rise to the redistribution of holes from the K point to the Γ point. The change in hole occupation ratio by the tensile strain for the MoS 2 model can also be supported by the band energy diagram for the corresponding tensilestrained models shown in Fig. 7. The difference in band energy level between the Γ and K

Longitudinal GFs
The calculated longitudinal GFs of the doped TMDC monolayer models are shown in Table 2. The carrier concentration dependence of GFs was hardly observed in the carrier concentration range of 10 17 -10 19    models are hardly affected by uniaxial strain or stress, so the piezoresistivity arises mainly from the carrier redistribution. An effect of the intervalley scattering is small enough to neglect for the models with the effective carrier redistribution. For the models where uniaxial strain or stress does not cause the carrier redistribution, the GFs are very small. As discussed above, the effective carrier redistribution can be observed in n-doped MoSe 2 and WSe 2 and p-doped MoS 2 monolayer models, and then significant GFs are obtained in these models. In particular, the p-doped MoS 2 monolayer model gives high longitudinal GFs, namely, 94.2 and 87.4 for the [112 __ 0] and [11 __ 00] directions, respectively. However, the GFs in n-doped MoSe 2 and WSe 2 models are not very high, because the effective masses at the K valleys are not far from those at the ΓK valleys. In contrast, the diagonal elements of the reciprocal effective-mass tensor in Eq. (6) at the Γ peaks are 0.21-0.33 times as much as those at the K peaks for the MoS 2 monolayer model, and therefore, a drastic decrease in electric conductivity should be observed by the hole transfer from the K peaks to the Γ point.
For each of n-doped MoSe 2 and WSe 2 and p-doped MoS 2 monolayer models, the longitudinal GF for the [112 __ 0] direction is slightly larger than that for the [11 __ 00] directions. This crystal anisotropy in GF is due to a tiny difference in the amount of carrier redistribution in the firstprinciples calculation. As discussed above, the mechanisms of the longitudinal piezoresistance effect in n-doped MoSe 2 and WSe 2 and p-doped MoS 2 monolayer models are common for the [112 __ 0] and [11 __ 00] directions; therefore, it can be regarded that the crystal anisotropy in piezoresistance should be small. The temperature dependence of longitudinal GFs in specific models is shown in Fig. 8. Similarly to general semiconducting materials, the absolute values of longitudinal GFs decrease as temperature increases.

Conclusions
The longitudinal GFs of MX 2 TMDC monolayer models (M = Mo and W; X = S, Se, and Te) have been evaluated by first-principles calculation. TMDC monolayers have a multivalley/ multipeak electronic structure, and effective regional partitioning for carrier redistribution has been presented. For some species of doped TMDC monolayers, uniaxial tensile strain causes a significant change in electrical conductivity through carrier redistribution in CB valleys or VB peaks. A high piezoresistivity was observed selectively depending on the combination of M and X. In this simulation, a p-doped MoS 2 monolayer model gave high longitudinal GFs of 94.2 for the [112 __ 0] direction and 87.4 for the [11 __ 00] direction at 20 °C in the carrier concentration range of 10 17 -10 19 cm −3 , based on the hole transfer from the K peaks to the Γ point by the uniaxial strain or stress with the change in effective mass.