Unveiling the multi-level structure of midgap states in Sb-doped MoX₂ (X = S, Se, Te) monolayers

We begin by recalling some properties of pristine 2H-MoX₂, as they will be central to our discussion of the electronic properties of the doped structures. In Table 1, we present the equilibrium properties of the three materials as obtained from our calculations. The calculated lattice constants, bond lengths and thicknesses of the monolayers are in excellent agreement with previous reports [32, 61, 62, 63, 35]. In Fig. 1, we present the band structures in the absence of spin-orbit coupling (SOC), together with the projected and total density of states (PDOS and DOS). As we can see, all gaps are direct, with the valence band maximum (VBM) and conduction band minimum (CBM) lying at the K point. Their magnitudes are reported in Table 1 and agree with previous reports [35]. The PDOS calculations show that the low-energy states (close to the Fermi energy) are composed mostly of $d$-orbitals from the Mo atoms and $p$-orbitals from the X atoms. Additionally, the contributions from $d_{xy}$ and $d_{x^{2}-y^{2}}$ orbitals of Mo are identical, as are those from the orbitals $d_{yz}$ and $d_{xz}$ of Mo and $p_{x}$ and $p_{y}$ of X. This is a consequence of the crystal field of the material, as the 2H-MoX₂ structure has the point group symmetry $D_{3h}$ and this group has only one and two-dimensional irreducible representations. The orbitals with identical contributions correspond to basis functions of the same two-dimensional representation and, as such, transform as linear combinations of the same functions when subject to all symmetry operations of the group. Therefore, as the PDOS is calculated through a sum over all k-points in the BZ, all symmetry-equivalent points in a band are included and the sum of the weights of each basis function of the representation over all equivalent points is identical. Moreover, in the case of the X atoms, the orbitals from the top and bottom layers are arranged in symmetrical and anti-symmetrical combinations which are either even or odd with respect to the mirror plane symmetry in the central Mo layer. This symmetry is present at every k-point of the BZ, despite the general symmetry being reduced at arbitrary points.

The inclusion of SO coupling introduces important changes to the electronic structure, as can be appreciated in Fig. 2. The most important effects are a reduction of the band gap and the splitting of the valence band, particularly at the K point. These values are also reported in Table 1 and, again, an excellent agreement is found with previous reports at similar levels of calculation [64, 65, 66, 67, 68] . In particular, the expectation value of the $z$-component of the spin in each of the split bands is opposite at the $K$ and $K^{\prime}$ points of the BZ, which is a signature of the spin-valley coupling. A similar effect is also observed at the conduction band, but the splittings are much smaller. Our calculated values range from $3$ meV in MoS₂ to $35$ meV in MoTe₂.

In order to assess the stability of Sb doping, we have evaluated the formation energies of Sb(Mo) and Sb(S) substitutional impurities in MoX₂. They are defined as:

We begin by discussing the energetics of defects in MoS₂, as it is the most widely studied material. From our fully relaxed structures, the formation energies of the Sb(Mo) and Sb(S) defects under Mo rich conditions are $4.68$ and $0.80$ eV, respectively, which suggests that the Sb(S) is the most favorable substitution in that scenario. In contrast, under S rich conditions, the formation energies are $0.72$ and $2.69$ eV, thus reversing the trend. The latter scenario is compatible with a recent experimental observation in which the Sb(Mo) substitution is favored in homogeneously Sb-doped MoS₂ nanosheets grown by CVD on a SiO₂/Si substrate [53]. Other growth conditions could also influence this result, such as temperature and the substrate itself. In particular, the formation of single vacancies could also aid the formation of the Sb(Mo) defect. For reference, we have calculated the formation energies of single vacancies in MoS₂ at the same level of theory. Under Mo rich conditions, the values we have obtained are $6.95$ eV for a Mo vacancy and $1.27$ eV for a S vacancy. Under S rich conditions, the formation energies become quite close, with values of $2.99$ and $3.26$ eV, respectively. Therefore, both vacancies are likely to form in this case and could be stabilized by being filled with Sb atoms, with a preference for Mo vacancy sites.

It should be noted that selecting a suitable reference structure for the evaluation of the chemical potential of S under S rich conditions could be a complicated task due to the existence of various S allotropes depending on temperature. For instance, the $\alpha$-S allotrope, which is a solid made of S₈ rings loosely bound, can be found below 100°C. The liquid-phase can be found beyond 120°C and the gas-phase is found by crossing 720°C, comprising of numerous S molecules. Nevertheless, S exists as a diatomic gas (S₂) above 880°C [69]. Therefore, considering that the CVD method is usually employed for growing quality materials in which the temperature range varies between 650 and 970°C, S can be found in gaseous form at typical CVD temperatures [24, 70, 71, 72, 73, 36]. In that scenario, our results for the energetics of vacancies do not fully agree with those from Ref. 35, in which the S vacancy is favored in both Mo and S rich conditions. This difference can be attributed to a different choice of reference structure for S in that calculation, namely, the bulk crystal. By using as reference the bulk structure of S instead of the S₂ gas phase in our calculations, we find formation energies of $2.39$, $1.35$, $4.66$ and $2.42$ eV for Sb(Mo), Sb(S), Mo vacancy and S vacancy, respectively, under S rich conditions. The values for the vacancies agree with those from Ref. 35, thus highlighting the importance of the choice. Similar outcomes for vacancies in MoS₂ are reported in Ref. 36, which considered both reference structures as the chemical potential of S is varied. However, for experimental growth conditions such as those of Ref. 53, where temperatures as high as $955$ K are achieved, the gas phase of S₂ should be predominant. This could explain the preference for the Sb(Mo) substitution in the experiment, as supported by our calculations.

Naturally, similar trends should also be observed in defective MoSe₂ and MoTe₂ layers with appropriate choices of reference structures for the Se and Te atoms. For MoSe₂, the formation energy of the Sb(Mo) substitution is $3.98$ eV under Mo rich conditions. By using the gas phase of Se₂ as reference, the formation energy of this defect reduces to $0.47$ eV under Se rich conditions, indicating that it is likely to be formed. If we use the most stable bulk phase of Se instead, the formation energy is $2.01$ eV. For MoTe₂, the formation of the Sb(Mo) defect is $3.15$ eV under Mo rich conditions. Now, considering the boiling temperature of Te is higher than the temperatures commonly found in (CVD) growth conditions, a gas phase should be absent and we only use its most stable bulk phase as a reference. In that case, the formation energy under Te rich conditions is $2.27$ eV, thus making the Sb(Mo) defect less likely to occur in MoTe₂. Finally, we have also cross-checked our calculations with the aid of the VASP code and a good overall agreement is observed for the formation energies. A complete table of the values can be found in the supplemental file (see Table S1). Cohesive energies are also included for comparison.

We now focus our discussion on the Sb(Mo) substitution in MoX₂ structures, which was observed experimentally for X = S [53]. We begin with the structural properties of the doped systems. Our variable-cell relaxations indicate that the introduction of the Sb induces important structural distortions in its environment. In Table 2, we report a few parameters that describe these modifications. First, notice that the average lattice constant, defined as the lattice constant of the supercell divided by its dimension (5), is slightly increased in all cases in comparison with the lattice constants of the pristine materials. This highlights the importance of a full relaxation, including the supercell vectors. Next, notice that the neighboring atoms move away from the impurity, resulting in an increased Sb-X bond length in comparison with the value for the Mo-X bond in the pristine system. The distance between the first-neighboring X atoms is also increased in comparison with the original value, which is the lattice constant of the pristine system. In contrast, the local thickness of the monolayer, which is the distance between first-neighboring X atoms in opposing layers, is slightly reduced. This thickness increases with the distance from the impurity up to the value found in the pristine system. In addition, we observe that the trigonal prismatic coordination is preserved in the vicinity of the impurity. This is expected since the impurity does not remove the mirror symmetry found in the original structure, as we discuss in more detail below. Therefore, a structural transition from the 2H phase to the 1T or 1$T^{\prime}$ phases upon Sb doping is not expected. However, the picture could be different in less symmetrical scenarios, such as different impurity configurations in multilayers and heterostructures. As an example, a recent report indicates that Cu/Co intercalation in multilayer SnS₂ induces a stacking-type modification in the structure [74]. A similar outcome may result from Sb intercalation in MoX₂ multilayers.

The resulting band-structures in the absence of spin-orbit coupling are shown in Fig. 3, together with the total DOS and the PDOS from the orbitals of the impurity atom and the sum of contributions from surrounding atoms. In the bands, we can clearly see five impurity levels inside the band gap (highlighted in orange), two of which are near-degenerate. Naturally, their signature is also present in the total DOS, where we see four peaks inside the gap and one of them is associated with the near-degenerate levels. We number these levels in order of increasing energy ($E_{1}$ to $E_{5}$), measured with respect to the valence band edge, and report the position of the corresponding DOS peaks in Table 3. The band gaps and the Fermi energy are also included ⁰⁰footnotetext: Here we report the Fermi energy as given by the calculations, which include the effects of smearing. The actual values, as predicted by electron counting, may differ by about $0.10$ eV, which is width of the smearing.. Note that the band gaps are slightly increased in comparison with the pristine values from Table 1. This result apparently contrasts with the calculations from Ref. 53, where a gap reduction is observed and attributed to a possible explanation of the redshift of the excitonic peaks induced by the impurity. However, we believe that both results can be reconciled with a proper assignment of the impurity levels, as we discuss below. Finally, our calculations yield Fermi energies about $0.32$ - $0.42$ eV above the valence band edge, suggesting a predominant $p$-type behavior in a similar fashion to other dopants such as Nb [33, 41, 38].

The electronic structure of the impurity levels can further be investigated through the analysis of their wavefunctions. We carry this analysis through PDOS calculations, shown in Fig. 3, and a calculation of the local density of states (LDOS) integrated over small energy ranges that span each level or set of levels in the case of near-degeneracy. The LDOS plots for Sb-doped MoS₂ are shown in Fig. 4, while those of MoSe₂ and MoTe₂ are shown in the supplemental file and display similar features (see Figs. S1 and S2). From these plots, we can see that all levels are localized around the impurity and display quite unique features. First, note that only the $E_{1}$ and $E_{4}$ levels contain contributions from the $5s$ orbital of Sb itself. On the other hand, all levels contain contributions from $p$ orbitals of neighboring X atoms, but these are particularly stronger in the $E_{4}$ level. Finally, all levels also contain contributions of comparable magnitudes from $4d$ orbitals of Mo atoms that are second neighbors to the impurity and smaller contributions from more distant neighbors. Interestingly, the $5p$ orbitals from Sb are not involved in any impurity level and, although we have not included the $4d$ electrons of Sb as semicore electrons in our calculations, earlier reports indicate that their contribution is very small in the energy range of interest [53]. In that sense, the Sb(Mo) substitution effectively reduces the number of available valence orbitals involved in the hybridization, suggesting that some of these levels could display similar features to that of impurity levels induced by a single Mo vacancy. To investigate this possibility, we have calculated the electronic properties of the Mo vacancy in MoS₂ at the same level of theory as our Sb(Mo) calculations. The band structure and PDOS/DOS plots are shown in Fig. 5 and the LDOS plots for the impurity levels are shown in Fig. 6. As in that case, we find five impurity levels inside the band gap, but now we find an additional near-degenerate pair right above the valence band edge. The positions of the corresponding DOS peaks are $0.32$, $0.64$ and $0.96$ eV, as measured from the valence band edge, and agree quite well with previous calculations [32, 35]. The band gap is $1.70$ eV, which is also slightly larger than the pristine value from Table 1. If we think of this system as a starting point to understand the electronic properties of Sb(Mo) doping, we may argue that the defect states of the vacancy system rehybridize with the $5s$ orbital of Sb and give rise to the defect states observed in Sb(Mo). By comparing the PDOS and LDOS profiles, it appears that the structure of the second near-degenerate pair in the vacancy system (levels $E_{3}$ and $E_{4}$) is preserved by the introduction of Sb, resulting in the near-degenerate pair observed in Sb(Mo) (levels $E_{2}$ and $E_{3}$). However, the contribution from the S $p$-orbitals is reduced. The first degenerate pair (levels $E_{1}$ and $E_{2}$) and the non-degenerate level in the vacancy system ($E_{5}$) appear to strongly interact with the $5s$ orbital of Sb, resulting in the remaining levels found in Sb(Mo). In particular, the near-degeneracy of the pair is lifted, resulting in only one level remaining inside the gap, while the other is mixed within the occupied extended states, as we can identify through the PDOS peak at about $-1.0$ eV in Fig. 3 (top). Notice that this peak moves towards the valence band edge as we change X from S to Te. Simultaneously, the contributions of Sb $5s$ orbitals to levels 2, 3 and 4 are reduced, which suggests that the strength of the interaction diminishes. This could be related to the reduction of the band gap and to an increase of the structural deformation induced by the impurity, as we can see from the increase of the $d_{\rm{X-X}}/a$ ratio in Table 2.

Another interesting feature of the PDOS calculations in both Sb(Mo) and Mo vacancy systems is that the paired orbital contributions observed in the pristine calculations are preserved. This is a consequence of the fact that the $D_{3h}$ point group symmetry is preserved with respect to the defect site. With the supercell construction, the periodic structure also preserves that symmetry. This property is even observed in the case of Sb(S) doping and the S vacancy, where the mirror symmetry and related operations are lost. In these cases, the symmetry with respect to the defect site is reduced to $C_{3v}$, but the irreducible representations of this group have a similar paired structure of basis functions. The electronic structure of these two cases is briefly discussed in the supplemental file (see Figs. S3 and S4). In particular, the results for the S vacancy agree quite well with a previous report from Ref. 35. Meanwhile, the Sb(S) substitution behaves as an acceptor with a near-degenerate ground state. Naturally, these symmetries will further be reduced in less idealized doping scenarios with disordered configurations, which may lead to an even richer electronic structure. Since such scenarios cannot be address through ab-initio calculations, future investigations with semi-empirical models such as the tight-binding method can contribute to increase our understanding of the electronic properties of these defects in MoX₂ and other semiconducting TDMCs.

Another important point to address is the effect of the supercell size in our calculations, which allows us to identify possible residual interactions between the impurity and its periodic images. To that end, we have performed a calculation in a larger $10\times 10$ supercell for two cases of interest, namely, the Sb(Mo) impurity in MoS₂ and the Mo vacancy in MoTe₂. These situations corresponds to atomic defect concentrations of 1% with respect to the total number of Mo atoms. A comparison of the band structures for the $5\times 5$ and $10\times 10$ supercells can be found in Figs. S5 and S6 of the supplemental material. We also show LDOS profiles for some impurity levels of the larger cell in Fig. S5. We find that the structure of the impurity levels is preserved, including their energies with respect to the valence band edge and the corresponding wavefunctions. In particular, for the Sb(Mo) defect, the weak dispersion of the $E_{1}$, $E_{2}$ and $E_{3}$ levels observed in Fig. 3 is eliminated in the larger cell, while that of the $E_{4}$ is reduced. A similar outcome is observed for the Mo vacancy. These results provide further evidence for the identification of these levels as localized states induced by the impurity. Additionally, they confirm that the main features of these levels, including their energetics and wavefunctions are well described by the smaller $5\times 5$ cell. We have also performed cross-check calculations with the VASP code for all of the structures considered in this work in the same $5\times 5$ cell (see Figs. S7 - S10 of the supplemental file) and we find an excellent agreement with the band structures from QE, thus confirming the robustness of our results.

We now discuss the effects of the spin-orbit interaction on the electronic structure of the Sb(Mo) doped systems. The band structures in the presence of SOC are shown in Fig. 7. As we can see, the main effect of the interaction is the introduction of small splittings in the impurity levels, which are highlighted in orange and ordered following the notation introduced in Fig. 3. In particular, the levels $E_{2}$ and $E_{3}$ remain near-degerate and split into two pairs by about $0.05$ eV, which is the largest splitting observed. On the other hand, no splitting is observed in the valence band and a very small splitting is observed in the conduction band, in contrast with the features observed in the pristine materials. Moreover, note that the nature of the band gap is also modified by the presence of the impurity. In Sb-doped MoS₂, the VBM is shifted to the $\Gamma$ point of the supercell BZ. In MoSe₂, the valence band energies at $\Gamma$ and $K$ are nearly identical and in MoTe₂ the VBM lies at the $K$ point. In all cases, the CBM lies at a point at the $\Gamma-K$ line, but note that the conduction band is particularly flat throughout a portion of the BZ in MoS₂ and MoSe₂.

The magnitude of the gaps are $1.77$, $1.58$ and $1.17$ eV for X = S, Se and Te, respectively. They are larger than the corresponding pristine gaps in the presence of SOC, in agreement with our calculations without SOC. This behavior contrasts with a previous calculation from Ref. 53, in which a slightly different band structure and a gap reduction were reported for X = S. We believe these differences are related to two factors. First, a full cell relaxation may result in a slightly different band structure than one obtained in a calculation in which only internal degrees of freedom are relaxed. In Fig. 8, we display a comparison of these two situations both in the absence and the presence of SOC for Sb-doped MoS₂ (calculated with the VASP code). As we can see, the calculation without a full-relaxation, but with SOC agrees quite well with the band structure reported in Ref. 53. Second, and more importantly, Ref. 53 prescribes a different assignment for the valence band and the electronic structure of the impurity levels, despite the similarities to our results. In that work, the first bundle of six weakly-dispersive levels (shown in blue in Fig. 8) is assumed to consist entirely of extended states, such that the topmost band is assigned as the valence band. However, as we see in our analysis of the wavefunctions of these levels, shown in Figs. 3 and 4 in the absence of SOC and Fig. S11 (supplemental material) in the presence of SOC, these states, which correspond to $E_{1}$, $E_{2}$ and $E_{3}$, show a clear localization around the impurity so we assign them as impurity states. Additionally, as we have mentioned before, these results are consistent with those observed in a larger $10\times 10$ supercell (see Fig. S5 of the supplemental material). Therefore, we consider that the valence band actually lies further below, corresponding to the closest red line in Fig. 8.

In light of this discussion, if we use the valence band assignment of Ref. 53, we find band gaps between 1.4 and 1.5 eV for the calculations in Fig. 8. These are in good agreement with the theoretical value of 1.39 eV reported in that reference and would be consistent with a reduction when compared to the pristine values (see Table 1). However, if we follow our assignment as described above, the actual gaps are around 1.8-1.9 eV for the same calculations and are consistent with an increase. In that scenario, the redshifts of the A and B excitonic peaks observed in experimental photoluminescence measurements, which were previously assigned to a reduction of the band gap, could be related to other mechanisms. These include excitonic effects such as a different quasiparticle correction, a stronger electron-hole interaction in the presence of the impurity or a mixture of the original excitons with new transitions involving the impurity-induced levels. All these effects are not included in a ground state DFT calculation, so only a calculation that fully addresses them, such as a GW/BSE calculation, would be able to fully elucidate the optical properties of these materials.

Finally, we now briefly discuss other features that are not included in our calculations which could be relevant to the stability and the electronic properties of the substitutional impurities. First, note that we have only considered neutral charge configurations for all defect configurations. Previous studies of single vacancies with charge in MoS₂ [36] indicate that, for S vacancies, the neutral state is favored for Fermi energy values (charge doping) throughout most of the gap, followed by a transition to a -1 charge state in the vicinity of the CBM, which is accompanied by the occupation of the associated defect levels. In contrast, for the Sb(S) substitutional defect, we see in our calculations that the defect levels are already (partially) occupied in the neutral state (Figs. S3 and S8 in the supplemental file), so a +1 charge state could be favored for Fermi energy values in the vicinity of the VBM. For Mo vacancies, Ref. 36 indicates the stability of different charge states throughout the gap, which is a consequence of the presence of multiple defect levels inside the gap. Considering that we see a similar defect level structure for the Sb(Mo) substitutional defect in our calculations, we expect a similar outcome for the stability of charged states. We stress however, that additional calculations would be required for a definitive conclusion, which are beyond the scope of the present work.

Second, since the electronic and optical properties of 2D materials can strongly be influenced by their environment, the substrate can play an important role in the description of the defect-induced levels inside the gap. For instance, Refs. 75 and 76 have found that the quasiparticle band gaps of semiconducting TMDCs can decrease by up to $\sim 100$ meV in the presence of capping graphene layers. As such, we expect that the defect levels may experience renormalizations of the same order. On the other hand, the renormalizations due to quasiparticle and excitonic effects are found to partially cancel out, resulting in small changes to the optical properties. In Ref. [53], the homogeneously Sb-doped MoS₂ nanosheet is prepared on top of SiO₂/Si substrate. Considering that SiO₂ is less polarizable than graphene, we expect that these renormalizations will be weaker in this case. However, the scenario may be very different with other 2D semiconductors, such as phosphorene, where much larger renormalizations are seen in the presence of a quartz substrate and capping layers of h-BN [77].