Electronic Properties of Defective MoS$_{2}$ Monolayers Subject to Mechanical Deformations: A First-Principles Approach

Monolayers (ML) of Group-6 transition-metal dichalcogenides (TMDs) are semiconducting two-dimensional materials with direct bandgap, showing promising applications in various fields of science and technology, such as nanoelectronics and optoelectronics. These monolayers can undergo strong elastic deformations, up to about 10\%, without any bond breaking. Moreover, the electronic structure and transport properties, which define the performance of these TMDs monolayers in nanoelectronic devices, can be strongly affected by the presence of defects (point or extended), which are often present in the synthetic samples. Thus, it is important to understand both effects on the electronic properties of such monolayers. In this work, we have investigated the electronic structure and energetic properties of defective MoS$_{2}$ monolayers, as subject to various strains, using density functional theory simulations. Our results indicated that strain leads to strong modifications of the defect levels inside the bandgap and their orbital characteristics. Strain also splits the degenerate defect levels up to an amount of 450~meV, proposing novel applications.

cal treatment [7,23,24,25,26]. These defects have advantages and disadvantages based on the desired application. For example, depending on the concentration of these defects, the performance of MoS 2 -based FETs may differ by several orders of magnitude [27,19,28]. On the other hand, some studies imply that the vacancy creation can extend the application of MoS 2 nanosheets [29,30,31,32,26], e.g., as single-photon emitters, due to localized states of the isolated defects [31,26]. In a recent study, Mo vacancies were generated site-selectively to write optically active defect states in TMDs MLs [26]. The mid-gap localized levels have also observed to improve the photoresponsivity of MoS 2 MLs by trapping the photoexcited charge carriers, leading to a growth of the photocurrent in photodetectors [33,34,35]. Furthermore, intrinsic defects, especially molybdenum (Mo) vacancies, may improve the contact resistance and the carrier transport efficiency of devices based on Au-MoS 2 top contacts [36].
Semiconducting 2D materials exhibit high resilience towards mechanical deformations in comparison to the conventional three-dimensional (3D) semiconductors [37,38,39]. While, for example, silicon tends to crack under 1.5% tensile strain, MoS 2 MLs can withstand about 10% tensile strain [40,41]. This capability allows for a high degree of flexibility [42]. However, the electronic structure of the 2D materials, in particular TMDs ML, is very sensitive to the applied strain [43]. In the case of MoS 2 , about 1.5% uniaxial tensile strain results in the direct-to-indirect bandgap transition, while 10 − 15% of biaxial tenisle strain leads to the semiconductor-to-metal transition [41]. Thus, mechanical deformations offer rapid and reversible tuning of the bandgaps in 2D TMDs of Group-6. These strain-engineered properties lead to new potential applications for TMDs MLs, such as piezoelectricity in MoS 2 , broad-spectrum solar energy funnel, and flexible transparent phototransistors [9,8,10,44,45]. It was also observed that biaxial strain can tune the properties of photodetector devices based on MoS 2 MLs [45].
In this paper, we scrutinize the influence of compressive and tensile strains on the electronic and energetic properties of defective TMD MLs. The exploitation of these phenomena may allow building blocks for novel applications. In the present study, we show the effect of various strain situations on formation energies, orbital characteristics, and defect levels (DLs) of intrinsic point vacancies in MoS 2 using density functional theory (DFT). We find noticeable splittings of up to 450 meV for the intrinsic DLs under strain. Furthermore, mid-gap DLs are shifted up(down) as compressive (tensile) isotropic biaxial strains are applied to the MLs. We also investigate the sensitivity of the formation energies to mechanical deformations. For the MLs, uniaxial and isotropic biaxial strains modify their energetic properties more profoundly than shear T1 strain.
The paper is organized as follows: the computational methods are explained in section 2, the modulations of the DLs, energetic properties, and orbital characteristics with respect to applied strain in section 3, and, finally, our conclusions are given in section 4. Complementary plots and images are provided in the Supplementary Information (SI).

Computational Details
We performed DFT calculations using numerical atomic orbitals (NAOs) basis sets to construct Kohn-Sham orbitals as implemented in the SIESTA code [46,47]. We used relativistic, normconserving pseudopotentials including the correction from core electrons, which were obtained by the Troullier-Martin method [48,49]. The exchange and correlation interactions were described by the Perdew-Burke-Ernzerhof (PBE) functional in the generalized gradient approximation (GGA) [50]. In all the geometric optimizations and electronic calculations, we employed a double-zeta basis set with one polarization function (DZP) and 4p diffusive orbitals. The Energy-Shift and the Split-Norm were set to 0.02 Ry and 0.16, respectively. We used the energy cut-off of 450 Ry to calculate hartree, exchang,e and correlation contribution to the total energy. A vacuum of 40Å along the out-of-plane (c) axis was used to make the structures effectively isolated as 2D layers. The standard conjugategradients (CG) technique was used to optimize the atomic positions and lattice vectors of equilibrium and strained configurations. The lattice constants along with the atomic positions of the unit cell were optimized until the Hellman-Feynman forces are below 10 meV/Å. Keeping the optimized lattice parameters, the same criterion was chosen to find the equilibrium atomic position of the pristine and defective supercells. Applying Monkhorst-Pack method, the Brillouin zone (BZ) was sampled with 25×25×1 and 5×5×1 k-points for the unit cell and supercell calculations, respectively. Spin-orbit coupling was not included in the present work, due to the computational resources. The total energies were considered converged when the difference between two consecutive self-consistent field steps was less than 10 −4 eV.  We used periodic boundary conditions to simulate the vacancies inside MLs, thus, in order to lower the defect-defect interactions in the com-Copyright line will be provided by the publisher putational model, we created supercells with sizes ranging from 6×6×1 up to 8×8×1. Accordingly, chalcogen vacancies were studied in a 7×7×1 supercell representation, while supercells of 8×8×1 were considered to scrutinize the transition metal vacancies and vacancy complexes.   [15,51,52,53]. S and Mo vacancies could also be produced under ion-irradiation [23,26]. Most of the defects kept the C 3v symmetry of the monolayer, except for the sulfur-pair vacancy in the top atomic layer, V 2S−par , as depicted in Figure 2(c). We illustrate the displacement map of the atoms surrounding these point defects in Figure S1. The displacement maps emphasized the significance of the local strain caused by creating intrinsic vacancies. Considering the fact that the electronic properties of 2D materials are very sensitive to strain, selecting the largest possible lattice size and performing atomic optimization are in fact crucial in further studying their characterization.
We applied four different in-plane strain variations to the defective MoS 2 MLs, as shown in Figure 3. Uniaxial strains in X-and Y-directions as well as isotropic biaxial strain in the XY-plane simulated the effects of simple mechanical deformations on the electronic structure, formation energies, and orbital characteristics of these MLs. We also considered an inhomogeneous shear type (T1), which only changed the angle between the in-plane lattice vectors, while keeping their magnitude constant. It has been reported that the electronic structure of pristine TMD MLs can be tuned in a controlled manner via strain [43,41,54]. We considered a wide range of strain values from 3.0% compression up to 5% tension, which were below the breaking point of the MoS 2 ML estimated from experiments [37]. In order to understand how the stability of these vacancies change under various types of deformations, we calculated their formation energies, E f , for all equilibrated and strained structures, as following: where E f,α , E d,α , and E p,α are the formation energy, the total energy of a defective structure and the corresponding pristine ML at strain α, respectively. If α = 0, this equation gives the formation energies for unstrained structures.
Here, n and m variables indicate the number of vacancy atoms. The chemical potential of S and Mo atoms are denoted with µ S and µ M o , respectively. In this study, we assumed that the chemical potentials of Mo and S are in a thermal equilibrium with MoS 2 , meaning: In our calculations, we considered the Mo-rich limit and the body-centered-cubic structure of metal atoms at 0 K temperature as a reference. While it has been reported that intrinsic defects in MoS 2 ML may have charge states, it was shown that the neutral defect states are the most stable over a wide range of Fermi-level positions [51]. Therefore, we focus on the neutral defects in this paper.

Results and Discussion
For the fully optimized unit cell of MoS 2 ML, we obtained the in-plane (xy-plane) lattice constant of 3.176Å. The Mo-S bond-lengths and the S-Mo-S angles are equal to 2.427Å and 81.9 • , respectively. Our calculated direct bandgap for this geometry is 1.73 eV at the K point, which is in accordance with other GGA-based results [16]. It should be emphasized, that the true quasi-particle bandgap is about 2.4 to 2.9 eV [16,55,56]. Our results differ due to the well-known underestimation of the bandgap in GGA. The computed lattice parameters are also in good agreement with experimental measurements and other theoretical results [17,18,16]. Previously, a detailed description of the orbital characteristics of the bands has been obtained by means of a 11-band tightbinding model Hamiltonian [57]. All occupied and unoccupied bands around the Fermi level consist of hybridized : Mo-4d and S-3p orbitals, as illustrated in Figure 4. Here, we display the band structure along with the partial density of states (PDOS) of a unit cell of the MoS 2 ML. The main contributions to the band edges come from the d z 2 , d x 2 −y 2 , d xy , p x , and p y orbitals which agree very well with previous studies [58,57].

Influence of strain on the formation energy
The calculated formation energies of defective MoS 2 MLs are shown in Figure 5. In line with previous studies, we find that the most probable defect in MoS 2 MLs is a single S vacancy,V S [51,52,58,59]. The formation energies of V 2S−top and V 2S−par are very close to each other. The case of the vacancy complex V M o+3S is more likely to happen than a single Mo vacancy, V M o , due to the coordination of the metal atoms and the fact that they are sandwiched between two S-atom layers. Hence, when creating Mo vacancies for single-photon emitters at selective sites [26], care has to be taken not to generate vacancy complexe. The energy required to form a molybdenum vacancy with its six neighboring sulfur atoms, V M o+6S , is the highest among all the suggested types of defect. Figure 6 shows the formation energies of the six studied vacancies as function of four different mechanical deformations applied to the defective MoS 2 MLs. In earlier reports, the case of V S vacancy in MoS 2 ML under uniaxial and biaxial strain was considered [11,60], corresponding to pluses, crosses, and triangles in Figure 6(a). While our results agree very well with these works, we study several other sulfur based vacancy complexes and strain situations in Figure 6  comparison with the other defects with 4d orbitals of Mo neighbors. For the case of V 2S−par , the strain in X-and Ydirections do not cause identical energy shifts, due to the broken C 3v symmetry. The geometry modifications also lead to a similar behavior for a Mo vacancy under ±1.5% of uniaxial strains. Applying shear T1 tensile and compressive deformations, the formation energies of all the defects are reduced.

Influence of strain on defect levels Sulfur Vacancies
The electronic band structure of the MoS 2 ML containing sulfur vacancies is shown in Figure 7, revealing the position of several localized DLs (green lines). As shown in Figure 7(a), sulfur vacancy of the V S type introduces a localized band near the valence band maximum (VBM) and a double degenerate empty mid-gap state. A sulfur divacancy, V 2S−top , introduces the same localized levels together with another doubledegenerate state around the conduction band minimum (CBM), as depicted in Figure 7(b). Figure 7(c) shows that in the case of a S-pair vacancy on the top atomic layer, V 2S−par , five non-degenerate localized levels occur, one occupied and four unoccupied. The 4d orbitals of Mo neighbors close to the vacancy form the DLs which are in very good agreement with previous reports for the unstrained ML [58,20,25].
To study the effect of strain on the electronic properties of defective MoS 2 MLs, we have investigated the change in the positions of the DLs and their orbital composition for all studied strain situations. Figure 8 and 9 show the results obtained for the V 2S−top and V 2S−par defects, respectively. In addition, the corresponding results obtained for a single sulfur vacancy, V S , are shown in Figure S2 and S3 in the supporting information. Green dashed-lines and black lines indicate the position of the Fermi energy and band edges, respectively. The occupied DLs in all the defective structures are shifted significantly less than the unoccupied DLs under any of the applied strain situations.
In the relaxed MoS 2 ML with V 2S−top , based on analysis of the orbital characteristics, DL2, DL3, DL4, and DL5  levels are mostly composed of d x 2 −y 2 , d xy , d xz , and d yz orbitals of the neighboring Mo atoms, respectively. For strain in X(Y)-direction, the DL2(DL3) level is strongly tuned, as illustrated in Figure 8(a) and (b). We speculate that this opposite shift of degenerate DLs is the consequence of the relative direction of each uniaxial strain to the nodal planes of the orbitals involved in the bands. Shown in Figure 8(c), since biaxial isotropic strain does not break the C 3v symmetry, the degeneracy remains intact, but the bands are shifted up(down) by compressive (tensile) strains. The hexagonal symmetry is removed via uniaxial and shear T1 strains which leads to breaking of the degeneracy of both the deep levels and the states around the CBM. Although the orbital composition can be uniquely identified for the relaxed MLs, shear T1 strain mixes the orbital contributions from Mo neighbors into the DLs resulting in a shift in the opposite directions, as shown in Figure 8(d). While a composition of d x 2 −y 2 and d xy makes up the states DL2 and DL3, localized level DL4 and DL5 are a mixture of d xz and d yz orbitals. In Figure S4, the change in the orbital components of the mid-gap DLs as function of applied strains is demonstrated. The degeneracy splitting of DL2 and DL3, ∆E, is shown as function of compressive and tensile strains in Figure 8(e). The splitting reaches to almost 200 meV and 450 meV for 5% of tensile uniaxial and shear T1 strain, respectively. As illustrated Figure S2, the DLs for a V S inside ML MoS 2 demonstrate similar behavior under strain. However, compare to Figure 8(e), the splitting of the deep degenerate levels of V S is about half for the same amount of strains.
Due to the absence of the C 3v symmetry in the MoS 2 ML with V 2S−par , no degenerate levels are present in the band structure of the unstrained defective monolayer, as shown in Figure 7(c). These bands are labeled DL1 to DL5 in Figure 9(a)-(d). In Figure S5, we show the principal orbitals constituting these DLs of relaxed and strained defective ML. Two defect states closer to the Fermi level, DL2 and DL3, mostly consist of d x 2 −y 2 and d xy orbitals of neighboring Mo, respectively. Moreover, main orbitals in DL4 (DL5) are d xy and d z 2 (d x 2 −y 2 and d xz ). The uniaxial strain in X-direction (Y-direction) tunes DL3 (DL2) much more than DL2 (DL3) in such a way that bands cross each other at +1% (−1.5%) strain. This tendency is related to the directional influence of uniaxial strains on the nodal planes of the orbitals corresponding to these localized states. As it shows in Figure 9(c), the isotropic biaxial compressive and tensile strains move the deep DLs relative to each other as the overlap between their orbitals varies. Shear T1 strain combines the orbital components of DLs in a way that a mixture of d xz and d yz (d x 2 −y 2 and d xz ) orbitals is added to d x 2 −y 2 (d z 2 ) orbital to make DL2 (DL3) state, as shown in Figure 9(d). The orbital characteristics of other two bands, DL4 and DL5, are also mixed. Thus, they move away for all strain values. The absolute splitting of DL2 and DL3 under four types of strain is shown in Figure 9(e). For the case of shear T1 strain, even though the band edges are not modified as much as the other deformations, the separation of the localized bands is up to 290 meV for 5% of tensile strain.
Since these defects are optically active, the degeneracy splitting could be a way to nondestructively identify the type of defects as well as to measure the applied strain.  This should also shift and broaden the optical spectra of defective MLs at low temperature. Furthermore, due to high resilience of MoS 2 MLs to the mechanical deformations, it is possible to use such splittings as a switch in desired devices. Besides, as mid-gap states trap the charge carriers, the shift in the DLs under deformations, in particular isotropic biaxial strains, results in changes of the photoresponsivity and other characteristics of the photoconductor devices based on MoS 2 MLs.

Molybdenum Vacancy and Vacancy Complexes
The position of localized DLs (green lines), originating from Mo vacancy and vacancy complexes inside the MoS 2 ML, is demonstrated in Figure 10. There are a nondegenerate and two double-degenerate levels in the band structure of the MoS 2 ML with a single Mo vacancy. For the vacancy complexes such as V M o+3S and V M o+6S , there are an occupied localized band, two double-and a triple-degenerate, and a single unoccupied state in the band structure, as shown in Figure 10(b) and (c). The localized states are mainly originated from 3p orbitals of neighbor-ing sulfurs with a small contribution from Mo 4d orbitals in accordance with results of Refs. [58,20].
The change in the localized defect states of the MoS 2 ML with V M o as function of various types of compressive and tensile strain are shown in Figure 11(a)-(d). The defect states inside the bandgap are labeled as DL1-DL6. The DL1 state is consistuted of a mixture of p x and p y , while DL2 level is mostly made of p x orbitals of the six neighboring sulfurs. For the case of DL3 and DL4, d z 2 orbital of Mo is mixed with sulfur p y and p x , respectively, to construct the states. Both DLs also containing a small contribution from d x 2 −y 2 orbital of the neghiboring molybdenums. The non-degenerate DL5 state is mostly composed of d xy , d x 2 −y 2 , p x , and p y orbitals of atoms surrounding the vacancy. The orbital characteristics of the levels in defective structures, as well as geometry modifications under strains are presented in Figure S6 and S7 of supporting information. As shown in Figure 11 ometry modify the hybridization of the orbitals which, in turn, leads to an abrupt shift in the localized states of DL1-DL5. This will be discussed in detail below. Moreover, except for DL3 and DL4 bands under tensile strain in Xdirection, the degeneracy of DLs is removed by applying uniaxial and shear T1 compressions and tensions. Shown in Figure 11(c), isotropic biaxial strains shift the degenerate bands, but not separating them. The amount of degeneracy splitting of occupied (unoccupied) degenerate DLs are displayed with solid(dashed)-lines as function of four types of strain in Figure 11(e). Of significance, tensile shear T1 breaks the degeneracy of the occupied levels the most and up to 330 meV. The effect of various strains on the DLs of vacancy complexes, V M o+3S , and V M o+6S , are depicted in Figure S8 and S9. As localized states are shifted due to applied strains, we expect a drastic change in the position and height of the peaks in the optical spectra of defective MLs. Accordingly, mechanical deformations can impact the performance of flexible optoelectronic devices based on MoS 2 MLs, e.g. single-photon emitters.
The change in the charge density of the MoS 2 ML with V M o under strain in Y-direction is shown in Figure 12.
Here, the strain cases of −2.0% and +2.5% are presented as an example in Figure 12(a) and (c), respectively. At zero strain, the hexagonal symmetry is visible in the density profile of the defective monolayer, displayed in Figure 12(b). As compressive or tensile strains applied, the symmetry is broken resulting in two sets of neighboring sulfurs reducing their distance by up to 27.39% (29.34%) for +2.5% (−2.0%) strain and forming a charge density overlap. The remaining set is simultaneously pushed away from its equilibrium position. As a consequence, stepwise shifts of localized bands are observed in the electronic structure of the MoS 2 ML with V M o , as shown in Figure 11.

Conclusion
In this paper, applying first-principles calculations, we have scrutinized the influence of four dif-ferent compressive and tensile strains on the electronic and energetic properties of the MoS 2 ML with point defects: We have shown that applying strain is a simple yet powerful tool to tune defect properties in MoS 2 ML, e.g. changing significantly the formation energies. As an example, strain reduced the energy of formation for V M o and V M o+3S vacancies. In addition, shear T1 strains lowered the formation of all the point defects. Breaking the symmetry of the monolayers lead to considerable degeneracy splitting of the DLs, ranging from a few meV to more than 400 meV, depending on the vacancy and type of strain. These could be used as a noninvasive method to identify the type of defect. Since MoS 2 MLs are robust to mechanical deformations, such splittings could be used as switches in devices. It also allows for a measurement of strain via optical means. We observed stepwise shifts in the localized energy levels of the MoS 2 ML with Mo vacancies under strain. These shifts are shown to originate from the transition of the charge overlaps between neighboring atoms. The tunability of the photodetector devices properties via strain could stem from the shift in the localized DLs, acting as trapping sites for photoexcited charge carriers, under the applied deformation. Therefore, for flexible optoelectronic devices, the effect of strain on the localized DLs position needs to be considered. Due to the analogy of the properties and geometries of various compounds in the TMD family, we expect a similar response to strain from the intrinsic defects inside their MLs. These findings should stimulate further experimental investigations on strain and defect engineering of TMDs monolayers and their potential application in selfpowered nanosystems, electromechanical sensors, photovoltaic and flexible devices.