Semimetal transition in curved MoS2/MoSe2 Van der Waals heterojunction by dispersion-corrected density functional theory

We present a theoretical study for MoS2/MoSe2 Van der Waals heterojunction in the armchair direction, and periodicity in the y-direction, under the mechanical deformation process to explore electronic structure vs. curvature angle. Our findings reveal that the heterojunction maintains chemical stability, even under high deformation, and the bandgap of the heterojunction is inversely proportional to curvature angle; the shift from semiconductor—with a bandgap of 0.8 eV—to semimetal occurs at deformation angles as low as 5°, having a gapless material. The mentioned transition corresponds mainly to distortion of half-filled molybdenum d-orbitals and chalcogen–chalcogen p-orbitals overlapping near the Fermi level.


Introduction
Transition metal dichalcogenides (TMD) exhibit novel electrical and mechanical properties like tunable bandgap in addition to a transition from an indirect to a direct bandgap, [1] electron mobility up to 200 cm 3 V −1 s −1 , [2] experimental values of elastic modulus superior to 130 GPa, [3] and a remarkable chemical stability. [4]hus, layered TMD including MoS 2 , WS 2 , and WSe 2 can be used for atomic or monolayer deposits with potential use, especially in flexible electronics targeting the Internet of Things (IoT), and highthroughput applications, [5,6] furthermore, the ability to adjust the electronic properties of these materials allows new functionalities like tuned electronic structure as occurred at twisted angles between two-dimensional layers of MoS 2, which conforms Moiré patterns. [7] recent paper published by Pu et al. shows that electrical properties in double-layer MoS 2 transistor subjected to mechanical bending, still keeps its operating performance up to a bending radius of 0.75 mm [8] ; in agreement with theoretical insights as presented by Sharma et al. and others, [8,9] computational calculations, especially the density functional and molecular dynamics schemes, opens an opportunity to explore key aspects of charge migration and orbitals distribution under strain and compression conditions in the plastic regime.On the other hand, a state-of-the-art experiment presented by Casillas et al. with aid of Cs-corrected Transmission Electron Microscopy (TEM) and combined molecular dynamics simulations, determined the resilient character of MoS 2 nanosheets due to structural restitution performance even at a curvature radius of 0.5 nm at applied pressure in the GPa regimes. [10]Additionally, nowadays there is also an increased focus on the integration of TMD heterojunctions for the design and fabrication of highthroughput electronic devices, [11] targeting large-scale fabrication, and tunability with notable results and progress. [12]Li et al. reported MoSe 2 /MoS 2 heterojunctions as-fabricated using epitaxial growth techniques, indicating high-performance for hydrogen evolution catalytic reactions as mainly caused by band alignment and charge transport over its interface. [13]Hosseini et al. reported that electron mobility (μ s ) is proportional to an applied strain and multilayer thickness with values of μ s ~ 300 cm 2 V −1 s −1 at 3% tensile strain in multilayer MoS 2 the surface. [14]And Cui et al. [15] indicated a stable resistive switching behavior under bending conditions for Ni/TiO 2 , making a higher p-n heterojunction in comparison with ITO/Si heterojunction as reported by Yao et al. [16] for tunable band structure as a function of elastic strain.However, there are a small number of reports regarding the variation of the electronic properties of TMD heterojunctions -especially important MoS 2 and MoSe 2 -as a result of mechanical bending; thus, its understanding is crucial to take advantage of TMD heterojunctions in electronic applications.Here, we present a theoretical study to determine the electronic structure of Van der Waals MoS 2 /MoSe 2 heterojunction under mechanical deformation from 0° to 45° bending curvature.Our results indicate an overlap of chalcogen p-orbitals and distortion of half-filled metallic d-orbitals around the Fermi level caused by mechanical deformation; such distortion of the d-orbitals leads to a transition from semiconductor to semimetal character of the heterojunction at bending angles starting at ~ 5°.This theoretical study provides a new platform for designing bandgap-engineering devices employing TMDC and using new degrees of freedom.

Computational methods
The performed density functional theory (DFT) calculations were completed using CASTEP© code [17] as part of the Materials Studio® package, using the revised Perdew-Burke-Ernzerhof (RPBE) as exchange-correlation functional as part of the Generalized Gradient Approximation (GGA) to describe particle's interactions in all our systems during single-point energy calculations.Band structure calculations were performed by setting a self-consistent field (SCF) energy value of 1 × 10 -5 eV atom −1 , using a k-point mesh with a minimum separation of 0.07 Å −1 .We set the cutoff energy at 435 eV and we employed the long-range dispersion correction, known as the DFT-D2, as described by Grimme et al., [18] setting the parameter values of s 6 and d to 0.5 and 10, respectively, following previous reports. [19]The partial density of states (pDOS) was calculated using a k-point mesh of 9 × 2 × 2 sampling in the reciprocal space and using the same exchange-correlation functional as before.
Previous studies show that there is no relation between bending stiffness and the armchair or zigzag edges in MoS 2 , [20] therefore, MoS 2 nanoribbon was built in the armchair direction with a length of 8-unit cells or approximately 2.2 nm and a thickness of one layer, approximately 0.06 nm.The MoS 2 nanoribbon was created using a MoS 2 unit cell (space group P6 3 /mmc) with lattice parameters of a = b = 0.31 nm, and c = 1.84 nm [Fig.1(a)].MoSe 2 nanoribbon consists also in 8-unit cells created from the armchair direction of a MoSe 2 unit cell (space group P6 3 /mmc) with lattice parameters of a = b = 0.33 nm, and c = 1.84 nm-both MoS 2 and MoSe 2 unit cells were previously optimized with the DFT parameters mentioned above.The thickness of the MoSe 2 nanoribbon is 0.06 nm corresponding to a one-layer thickness.Both MoS 2 and MoSe 2 nanoribbon consist of fifteen molybdenum atoms and thirty chalcogen atoms, either sulfur or selenium, and both have periodic conditions over y-direction [Fig.1(b)].All proposed models were placed in large crystallographic lattice with dimensions of 0.58 nm, 3.1 nm, and 2 nm in the x, y, and z-direction, respectively with sufficient space to prevent atomistic interaction with near neighbor atoms.The heterojunction was manually bent from 0° to 2°, 5°, 10°, 25°, 30°, and 45°.For comparison, all calculated density of states and band structures were compared with the pristine (unbend or flat) model.

Geometric optimization of bending structures
In our study, we started by computing the relative structural stability of MoS 2 /MoSe 2 Van Der Waals heterojunction ( E vdw ) using the following expression, where E ht is the computed total energy of MoS 2 /MoSe 2 heterojunction, n is the number of atoms in the structure, and E MoS and E MoSe are the computed total energy of the MoS 2 and MoSe 2 nanoribbon, respectively.Using the computed values of total energy and curvature (Table S1), the estimated value of E vdw is 0.06 eV atom −1 and this indicates that the heterojunction with no bending curvature is geometrically stable, in agreement with experimental evidence, [21] with a binding distance between MoS 2 and MoSe 2 of 0.32 nm.
By estimating the energy density (E A ) from the computational calculations, which provides valuable energetic information, we were able to compute structural stability for heterojunctions using Eq. 2 which reads where E bend corresponds to the total energy of the proposed molecular model (in eV), E flat corresponds to computed energy of pristine (or flat) model, and A is the cross-sectional area in nm 2 .As observed, the energy density of a proposed theoretical heterojunction is proportional to bending curvature, with (1) maximum value at a 45° bending angle (~ 0.04 nm −1 ) as shown in Fig. 1(c).It is observed that the variation of energy density does not have similar values in comparison with previous reports for bending characteristics of MoS 2 nanosheets. [20]n our work, energy density values remain with variations of 8.712 × 10 3 eV nm −2 at 5°, to 8.731 × 10 3 eV nm −2 at 45°, indicating a ΔE A of 0.019 × 10 3 eV nm −2 .Here, as curvature increases, variation of energy density also increases, attributed to deformation effects of molecular flat nature of TDMC to reallocation of energy orbitals, expanding and contraction of Mo-S and Mo-Mo bond length, mostly at the edges. [22]ectronic structure The calculations of the partial density of states (pDOS), which express several available states per energy unit (eV), reveal how orbital distribution varies as a function of bending curvature.This allows us to determine the contributions of d-orbitals near the Fermi energy level, with a semiconductor to metal transition, as shown in Fig. 2(a-g).For comparison, it is clear that in the pristine model all density of states shows typical contributions from p-and d-orbitals corresponding to sulfur-selenium and molybdenum at valence and conduction band near Fermi, in agreement with previous reports. [23]However, for bending scenarios, it was detected that contribution of p-and d-orbitals have an increased presence near the Fermi level.The bending effect induces a reduction of energy bandgap provoking a semiconductor to semimetal transition for maximum bending angle.The latter has been reported previously for twisting layers of 2H-MoS 2 , [7] or under compressive strain conditions in 2H-MoS 2 layers, [24] attributed to orbitals' interaction and indicating the high flexibility of TMD in terms of electronic modulation.Previous reports suggest the increase of orbitals contribution near the Fermi level is linked to an overlap of p x orbitals from sulfur-sulfur bonding, as caused by bending and structural deformation. [7]However, from pDOS plots we detect a strong metallic transition for large bending curvature values near energy Fermi level (E F ), from 0 to 1.5 eV, in comparison to s-and p-orbital contributions.Contrary, for all bending scenarios the d-, p-and s-orbitals have a similar distribution regardless of the bending curvature between 0 eV and − 1.0 eV near E F .The half-filled d-orbitals, as observed, presented the higher variation as a function of bending curvature and can be understood as easily perturbated energy levels, mainly between K and Γ points in the Brillouin zone for strained 2H-MoS 2 sheets, [25] which derived in the contraction of the bandgap.A higher orbital weight at the conduction band and valence band of these d-orbitals represent an easier distortional energy level; as shown in Fig. 2, we have a gapless material after 5° of bending, in agreement with the bandgap contraction occurring on strained MoS 2 . [24]The analysis of pDOS results for all density of states computations indicates an increase of orbital contribution near E F could be attributed to the interaction between d-orbitals from molybdenum atoms in our proposed model, [26] causing irreducible points values of k = ~ 0.7 nm −1 .
For 0° and 45° bending angles, the molybdenum d-character and sulfur-selenium p-character orbitals distribution near E F reveals metallic transition, caused by d-orbitals coming mainly from molybdenum atoms [Fig.2(h) and (i)] which is similar to metal-semiconductor interfaces as reported by Farmanbar et al. in similar models. [27]The p-orbitals from sulfur and selenium atoms present a partial overlapping near the Fermi level, represented by the similar distribution between 0 eV and 3 eV, however, the signal is low when compared to molybdenum d-orbitals.The latter serves as theoretical insights that indicate semiconductor to metallic transitions happening at low angles values of 5°, meaning its relevance and extraordinary behavior for mechanical switching engineering applications in devices such as micro-electrical mechanical systems. [13]

Electrostatic potential
The electrostatic potential is understood as charge distribution perpendicular to (001)-basal plane obtained from energy Kohn-Sham calculations by CASTEP©.The energy electrostatic potential presents severe fluctuations as bending curvatures occur from 0° to 45°, with higher visible electrostatic perturbation for bending curvature at 45° (Fig. 3).The bending values of 1° to 5° reduces by about 40%, decreasing steadily until the heterojunction reaches a bending angle value of 45°.This change in the electrostatic potential is directly related to the orbital distribution as previously mentioned, due to the transition from semiconductor to semimetal character; p-orbitals overlap while d-orbitals distort due to mechanical bending, decreasing a total contribution to the density of states near Fermi energy level, and hence, the charge distribution.
Furthermore, the position of the Fermi energy level, corresponding to zero values, decreases as heterojunction bending curvature occurs (Table I) and having a drastic change from positive to negative potential inducing a shift from semiconductor to semimetal.The Fermi energy level position remains stable with only significant changes from 0° to 1° and 35° to 45°, suggesting the pinning phenomenon in agreement with other theoretical studies. [28]y these results, we detect that when the heterojunction is flat its bandgap is ~ 0.8 eV (Fig. 3), which agrees with previous studies on the band structure of TMDC heterojunctions. [29]evertheless, the intersection of energy bands with the Fermi level as the MoS 2 /MoSe 2 heterojunction bends confirms the transition from a semiconductor to a semimetal character predicted with the pDOS happening between the M, K, and the Γ points of the Brillouin zone, and attributed to metallic d-orbitals from molybdenum and metal-metal interaction. [30]The computed variation of the bandgap as the heterojunction bends, also indicates that band alignment could be achieved even at small deformation angles, which can induce large areas of electrondonor acceptor ideal for catalytic properties or electronic modulation for supercapacitor applications as reported by Li et al. [13]

Conclusions
We present a theoretical study by the meaning of dispersioncorrected DFT approximation for Van der Waals MoS 2 /MoSe 2 heterojunction.Our results indicate a distortion of electronic structure caused by angular mechanical deformation in the range of 0° to 45°.The partial density of states near the Fermi energy level reveals a strong interaction of half-filled metallic d-orbitals with the transition from semiconductor to semimetal character, due to expanding and stretching of bonds on the sandwich-like Van der Waals solids at bending angles ~ 5°.The presented band structure calculations agree with several reports on TMD as in the literature.Our theoretical study serves as a framework for the experimental design of bandgapengineering devices employing transition metal dichalcogenide heterojunctions.
Figure 1.(a) MoS 2 and MoSe 2 unit cells (P6 3 /mmc space group), (b) model of the MoS 2 /MoSe 2 Van der Waals heterojunction side and top view, respectively.Dashed lines correspond to the periodic conditions on the as-proposed model.(c) Energy density values as a function of bending curvature.The red line is the polynomial fit to emphasize nonlinear behavior.Color code: cyan ball represents Mo, yellow and orange balls represent S and Se.

Figure 2 .
Figure 2. (a)-(g) The partial density of states calculated for as-proposed theoretical MoS 2 /MoSe 2 heterojunction as a function of bending angle.(h) The partial density of states of at 0° curvature angles.(i) The partial density of states of 45° curvature angles.Color code: cyan molybdenum, yellow and orange sulfur, and selenium atoms.The perturbation of the d-and p-orbitals is indicated by the red square box and the red arrows.

Figure 3 .
Figure 3. (a) Electrostatic potential profile along the z-direction from 0° to 45° of bending angle.The dashed line indicates the location of the interface.The distribution of the electrostatic potential degrades as the heterojunction bends, showing a clear alteration of the charge distribution along the z-direction.(b-d) Band structure of the MoS 2 /MoSe 2 heterojunction under a mechanical bending situation at 0°, 1°, and 45°; the dashed orange line indicates the position of the Fermi level just at the top of the valence band.As the bending increase, energy bands intercross with the Fermi level giving the transition of the heterojunction from a semiconductor to semimetal.

Table I .
Computed Fermi level position, work function, and bandgap at each bending angle of the MoS 2 /MoSe 2 heterojunction.