A computational study of MoS2 for band gap engineering by substitutional doping of TMN (T = transition metal (Cu), M = metalloid (B) and N = non-metal (C))

Tunable electronic properties of two dimensional Molybdenum disulfide (MoS2) make it a potential material. In this study, we inspect electronic and structural properties of TMN-doped MoS2 (T = Transition metal (Cu-copper), M = Metalloid (B-boron) and N = Nonmetal (C-carbon)) by using first principles DFT (density functional theory) calculations. Cu is substituted by Mo with varying concentration, which ranges from 2.08 to 8.33%, whereas B and C are replaced by S atoms with varying concentration of 2.08 to 4.16%. The substitutions result into significant variations in electronic and structural properties of MoS2. Moreover, the importance of substitutional site has been elaborated. The substitution of these impurities, variation in concentration and the replaced sites of MoS2 cause to modify the structure and energy gaps. Resulting bandgap fluctuates remain between 0.16 eV to 0.48 eV relative to 1.95 eV of pristine MoS2. The PDOS calculations show good bonding relation among the host MoS2 and the foreign impurity TMN. Therefore, substitution of impurities gives the opportunity to vary the bandgap as required for its valuable applications as semiconducting materials.


Introduction
Being most attractive class of materials, three dimensional materials exhibit vast range of applications that involve super lubricants [1], solar converters [2], half metallic magnets [3] and catalyst in redox based reactions [4]. Besides three dimentional materials, unique and regular atomic arrangement of monolayered graphene exposes its outstanding properties. A new scope to the 2D material research has been introduced by the isolation of monolayer of graphene that is expected to facilitate the resembling transition metal dichalcogenides (TMDs) [5,6]. A representative role in nanoscience has been played by graphene because of its physiochemical properties. However scholars are taking a stance on graphene to overcome its defects including thermal stability [7], vacancies [8], weaker anchoring of shape of atoms and clusters [9]. These defects are necessary to be overthrown just to beat the restrictions that lemmatize its applications. Basically 2D materials have drawn a considerable interest as potential basis for next generation of electronics with exceptional properties. Nonmetallic materials generates 2D semiconducting and insulating material due to their intrinsic chemical properties. Ab-initio calculations have confirmed that the semi-conductive monolayer MoS 2 possess direct bandgap of nearly 1.95 eV, whereas in bulk form it has indirect bandgap of 1.2 eV [10,11].
Graphene because of its bounded properties, is replaced by MoS 2 . MoS 2 being transition metal dichalcogenide material, with properties of large surface area, high mechanical robustness has opened a new gateway for the 2D material usage in energy storage applications [12]. In addition, MoS 2 has attained great interest in the fields of sensors‚ organic light emitting diodes‚ memory and power nano devices [13]. It can also produce oscillating piezo electric voltage and current outputs because of odd number of layers [14].
As TMDs exhibits unique electrical and optical properties that has triggered the researchers towards 2D materials like graphene. Out of all TMDs, MoS 2 has attained a distinctive position because of its good mobility, high on/off ratio, large optical absorption and much more. Hence, the use of single 2D material is lemmatized because of its inability to meet up the functional performance requirements for 2D practical applications. Graphene has good conductive ability but the zero band gap is the obstacle in its use for switch control. MoS 2 possess inconsistent band gap that is commanded by several layers. Although it's matchless electron mobility as compared to graphene make it ineligible for transparent electrode usage. The absorption and electronic properties of MoS 2 are modified through bandgap tuning via quantum confinement. Thus for practical application of MoS 2 based nano materials, more struggle is required from integrating to modernizing the physiochemical properties.
To enhance the significant potential of 2D materials for use in optoelectronic, sensing and memory devices, energy storage and catalysis 'doping' is one of the important processes to modify the band gap of 2D materials. Different schemes have been considered for this purpose as for nano devices chemically doped/modified carbon nanotubes (CNTs) and graphene and n-p type doped graphene via chemical vapor deposition (CVD) [15,16]. Basically substitutional and surface charge doping are implemented in case of 2D semiconductors [17,18]. Dolui et al [19]. Theoretically proved that transition metal elements prefer to substitute at Mo site instead of S site. Substitution of some of transition metal like as Pd, Ag, Cd, Nb, Zr, and Y on MoS 2 categorize the Mo as n-type (Pd, Ag and Cd) or p-type (Nb, Zr and Y). One of electron is added or removed from Mo site by this substitution. The transition metals such as Fe, Ni, Co and Cu have been doped at Mo and S edge site by Wang et al [20]. Besides this, the doping of the transition metals on both Mo and S edge sites of MoS 2 has been observed by them. Ramasubramaniam et al [21]. Applied DFT to examine the electronic properties of Mn-doped monolayer MoS 2 . Lu et al [22]. Inspected the electronic characteristics of deformities of the doped MoS 2 . They discussed therein the band gap opening of monolayer MoS 2 due to adatoms of various number of impurities. Cheng et al [23]. Calculated the structural, electronic and magnetic properties of doped MoS 2 with numerous transition metals by first-principles based calculations. Wang et al [24]. Explored the magnetic and electronic properties of Co-doped monolayered MoS 2 by first-principles calculations. Andriotis et al [25]. Carried out the electronic and magnetic properties of transition metal doped and co-doped MoS 2 using ab-initio calculations. As a result a noticeable energy gap reduction has been observed by TM doping in MoS 2 , above all this in TM dopants magnetic features as ferromagnetic coupling has been observed. Hussain et al [26]. Studied the effect of dopant site and dopant concentration of Cu on the band gap of MoS 2 using DFT calculations. They revealed that electronic and structural properties of MoS 2 are modified due to the Cu doping. Joseph et al [27]. Used hydrothermal technique and analyzed the influence of Cu doping on thermoelectronic properties of MoS 2 . X-ray diffraction (XRD) and high resolution transmission electron microscope (HRTEM) investigations verified the interaction between Cu and MoS 2 . Xia et al [28] also used the hydrothermal process for synthesis of Cu-doped MoS 2 nanosheets and characterized the binding properties of Cu with Mo and S through energy dispersive X-ray (EDX) mapping and X-ray photoelectron spectra (XPS). Thus, MoS 2 being a nonzero bandgap material with good thermal stability, has dragged a noticeable interest. Moreover MoS 2 based devices exhibits great prospective to use in energy harvesting and semiconducting applications.
In this paper, we have provided a systematic study of electronic properties of TMN-doped MoS 2 . Cu, B and C are chosen as dopant which are experimentally [27,29] and computationally [30] proved to have a good chemistry with MoS 2 . Furthermore, TMN assists in fine bandgap tuning by changing its substitutional site and concentration variation which advocates for its good conducting role close to metallic character. The purpose of our study is to ascertain the outcomes of TMN concentration, effect of variation of substitutional site on the bandgap of MoS 2 and the chemistry of foreign element with the host that suggests it practicable as semiconductor. Further, we have reported the stability of the system based on cohesive energy. As an inference a significant variation in bandgap value has been observed which is advantageous to enhance the electronic properties for its use in the applications towards on/off and semiconducting gadgets.

Computational methodology
We performed DFT computations employing vienna ab-initio simulation package (VASP) [31], which implies projector augmented wave (PAW) potentials [32] using the scheme of exchange correlation functional Perdew-Burke-Ernzerhof (PBE) [33]. The cut off kinetic energy of plane wave basis set was set to be 400 eV. We used a 4×4 supercell of MoS 2 , comprising of 16 and 32 atoms of Mo and S respectively. To avoid interlayer interactions, the vacuum in z-direction is kept 15Å. A Monkhorst Pack mesh of 5x5x1 in brillion zone is used for geometry optimization [34], whereas a bigger grid of 15x15x1 is employed for density of states (DOS) calculations. Energy minimization computations are performed until the complete Hellmann Feynman forces were smaller than 0.02 eV Å −1 [35]. In all the calculations, before and after the substitutions, the ionic positions were allowed to relax in the supercell. Cohesive energy was determined as: Where E System is the all-out energy of the framework. E Mo , E S , E Cu, E B and E C are energies of confined Mo, S, Cu, B and C separately. N Mo , N S , N Cu , N B and N C are the quantity of Mo, S, Cu, B and C atoms of the framework separately and N is absolute number of atoms in supercell.  The formation energies reported in table 1 have been calculated employing the following formula: Where E(doped), E(pristine), E(dopant) and E(removed) are respectively, the calculated total energy of the doped sytem, pure system, dopant and removed host atom.

Results and discussions
We incorporated Mo atom(s) in MoS 2 with Cu and, S atom(s) with B and C and observed that the change in bandgaps depends upon the impurity substitution, variation in impurity concentrations and substitutional sites.  [26,36]. The Cu is displaced from the original position of Mo by 0.042 Å, whereas B and C by 0.98 Å and 0.88 Å respectively after geometry optimization process. The C-Cu-B, C-Cu-S and B-Cu-S bond angles are determined as 116°, 107°and 89°r espectively relative to their 82°values in pristine MoS 2 [26].

Substitution by single impurities
In order to find the energy gap of the model given in figure 1(a), we computed the band structure computations. The results obtained are depicted in figure 1(b). At Fermi level, the energy gap is observed to be 0.17 eV. A similar energy gap with Mg doping in graphene has been accounted in our recent research [37]. The presence of impurity lines divide the wide bandgap into three parts having widths 0.38, 0.31 and 0.39 eV. The impurity lines near to Fermi level appear due to Cu atom. Our outcomes are in agreement with pervious detailed examinations [26,[36][37][38]. To investigate the role of dopants, we performed the PDOS calculations. The consequences are represented in figure 1(c). The participation to DOS around the Fermi level is due to the orbitals of Cu, S and Mo. The hybridization appears just below the Fermi level due to p orbitals of B and d orbitals of Cu. A little participation comes to the PDOS from the p x of C, p z of S, and d of Cu at 0.5 eV. A sharp peak can be noted in conduction band, which mainly comes due the d-orbitals of Cu, p z of S, p orbitals of B, and p x of C. An intense peak in the conduction band shows that a strong overlapping among the B, Cu, S and Mo is observed around 1.2 eV energy. These results are consistent with our recent reported results [26].

Effect of concentration and substitutional site
In this subsection we systematically changed the sites of B and C atoms by keeping the impurity concentration of Cu (4.16%) same, and have observed a variation in bandgap values. Furthermore, bonding mechanism is also considered.

Two Cu substituted by Mo and, B and C by S
We substituted two Mo atoms with Cu, which made the Cu concentration 4.16%. Two S atoms were additionally replaced, one with B and other with C as shown in figure 2(a). The bond lengths of d Cu1-S , d Cu1-B , d Cu2-S , d Cu2-C , d Mo-B and d Mo-C are observed to be 2.6 Å, 2.1 Å, 2.7 Å, 1.9 Å 2.1 Å and 2 Å respectively. All d Mo-S bond lengths are found to be the same as the pristine MoS 2 [26]. The Cu is displaced from the original position of Mo by 0.16 Å, whereas B and C by 0.39 Å and 0.68 Å respectively after geometry optimization process. These findings are well agreed with our recent reported results [26,30].
The electronic band structure calculations showed the direct and indirect band gaps having values 0.26 and 0.46 eV respectively as depicted in figure 2   After electronic band structure computations, we noticed that the highest band gap is observed 0.38 eV which is below the Fermi level as indicated in figure 4(b). The observed band gap (0.38 eV) has smaller value as compared to the above reported band gaps. Emergence of impurity lines splits the bandgap into five parts having widths 0.38, 0.30, 0.22, 0.20 and 0.19 eV (see figure 4(b)). Thus, by changing the site of dopants a significant variation in band gap occurs.
The PDOS calculations results indicate that the main involvement to DOS originates from the orbitals of Cu, Mo, B, C and S in the energy area −3 to −0.8 eV as shown in figure 4(c). A peak is noticed due to hybridization of the orbitals of C, S, Cu and Mo at −0.45 eV. The hybridization is observed in the locality of Fermi level due to the the orbitals of S and Cu. A small overlapping is noted at Fermi level due to the orbitals of C and Mo, while a large overlapping is found due to p z of S and d of Cu, Mo and C at 0.55 eV. We also determined a strong hybridization due to the orbitals of B, Mo and Cu at 0.85 eV, which indicates a sigma bonding. In conduction band, an overlapping can be seen due to the orbitals of C, S, Cu and Mo at 1 eV. The main overlapping in the conduction band originates from the orbitals of B, S and Mo as depicted in figure 4(c).

Three Cu substituted by Mo and, B and C by S
In the present case, the concentration of B and C are kept same as discussed above and only the concentration of Cu has been increased to 6.25% as depicted in figure 5(a). After geometry optimization, the bond lengths d  The electronic band structure computations were performed for the present configuration. These results indicate that the direct and indirect bandgaps appear above and below the Fermi level. The highest value of indirect bandgap is found as 0.46 eV above the Fermi level. The same value of bandgap (0.46 eV) was observed in earlier reported section 3.2.1 below the Fermi level. The direct bandgap of 0.30 eV is determined below the Fermi level, while 0.26 eV direct band gap was determined above the Fermi level reported in previous section 3.2.1. Thus, by increasing the impurity (Cu) concentration the change of location of direct and indirect band gaps occur. We also observed a small band gap just below the Fermi level as shown in figure 5(b). We observed many impurity lines which cause the reduction in big band gap. Our results are well agreed with previous studies [26,28].
The PDOS calculations indicate that an overlapping is observed due to p orbitals of B, S and C, and d orbitals of Mo and Cu at −0.7 eV. We also observed hybridization among p y of C, p z of S and d of Cu and Mo around 0.25 eV. A little peak is can be seen at 0.7 eV due to overlapping of concerned p and d orbitals of the elements involved. An overlapping is is noted due to p orbitals of B, C and S, and d orbitals of Cu and Mo in energy area 0.8 to 1.2 eV as shown in figure 5(c).

Three Cu substituted by Mo and, one B and two C by S
Keeping the concentration of Cu and B same as reported in section 3.2.4, only variation in C concentration has been made to observe its effect on bond length, which comes to be 2.  arises from the p z of S, p of C and d of Mo and Cu, while in conduction band, a strong hybridization is observed due to these atoms in 0.6 to 1.1 eV range. A small peak is also determined by the overlapping of C, Cu, and S and Mo atoms at 1.3 eV.
3.2.6. Three Cu substituted by Mo and, one C and two B by S In this configuration, the concentration of Cu and C is kept constant as discussed in sections 3.2.4 and 3.2.5, only the variation in B concentration has been made to observe its effect on geometry parameters and other properties. The highest bond length value between Cu 1 and C is observed to be 2.18 Å. The bond lengths d Cu2-C and d Cu3-C are measured to be same being 2.06 Å. The atomic distance between Cu and B is determined as 2.01 Å. The bond lengths between impurity atoms and S are observed to be 2.24 Å, while bond length of 2.1 Å between Mo and B atoms is noticed. The bond length values between Mo and S are determined in ranging from 2.33 to 2.5 Å as shown in figure 7(a). The highest bond angle (Cu-B-Mo) is observed to be 121°. Thus, a significant variation in geometry parameters are observed on varying the concentration of B atoms.

Four Cu substituted by Mo and, two C by S
In this system, we substituted four Mo atoms with Cu, which made the Cu concentration 8.33%. Additionally, two S atoms were replaced with two C atoms and creation of a single vacancy of S atom as depicted in figure 8(a).   [39].
A large overlapping is determined among the d of Mo and Cu, and p of C and S in ranging from −3 to −0.8 eV as shown in figure 8(c). The hybridization is observed due to the p of S and C, and d of Mo and Cu in energy area −0.6 to −0.15 eV. At the Fermi level, a small peak is observed which is from p of C and S, and d of Cu and Mo. Many peaks are also determined by overlapping of the orbitals of Cu, Mo, S and C in conduction band. In conduction band, orbitals of B and Cu also overlapping in energy area 1 to 2 eV.

Four Cu substituted by Mo and, B and C by S
The concentration of Cu is kept same as reported in section 3.2.7, only the one C from two is replaced by one B including the site change and after this the change in geometrical structure has been noticed. The bond length between Cu 1 and B is observed to be 2.5 Å, while all bond lengths between impurity atoms are determined as 2.2 Å. The atomic distance between S and Cu is found to be 2.3 Å.
The highest value of indirect bandgap is observed to be 0.37 eV above the Fermi level. We also observed the smallest value of indirect bandgap 0.18 eV just below the Fermi level as shown in figure 9(b).
The PDOS calculations indicate that a large overlapping is determined among the d of Mo and Cu, and p orbitals of C, B and S in valance band as shown in figure 9(c). An intense peak is found at −0.45 eV which comes from the orbitals of B, C, S, Mo and Cu. This overlapping forms a bond. At the Fermi level, an overlapping is observed which arises from the orbitals of S, Cu and Mo, while a small overlapping of the orbitals of B also played a role. At 0.4 eV, the p y of C is overlapping with d of Mo and Cu at 0.45 eV. A sharp overlapping is determined between p x of C and B at 1.1 eV. The large overlapping is observed due to the orbitals of B, S and Mo in the conduction band ranging from 1.3 to 2.5 eV (see figure 9(c)).

Four Cu substituted hexagonally by Mo and, B and C by S
The concentration of Cu, B and C is kept the same as reported in section 3.2.8, only the site changing of dopants has been made to notice its effect on bond lengths and other properties. The results indicate the bond lengths between B and its nearest atoms (Cu and Mo) are 2.1 Å, while C and its nearest atoms (Cu and Mo) have 1.9 Å as depicted in figure 10(a). All other atomic distances are determined to be 2.4 Å. Hence, a small variation in bond lengths is produced by changing the dopant sites. A similar hexagonal doping of Cu in MoS 2 , Be in h-BN and Be in graphene have been studied in earlier studies [26,36,40].
The band structure computations results are depicted in figure 10(b). These results show that the bandgap value is observed to be 0.27 eV above the Fermi level. Fermi level moves 0.70 eV below the bandgap that indicates a p type doping. By keeping the same value of concentrations (Cu, C, B) as discussed in the previous section only the sites of dopant are changed to notice the effect on bandgap, a significant change is observed as outcome.
The PDOS calculations indicate that a large overlapping is determined among the orbitals of Mo, Cu and S in energy area −3 to 0.7 eV as shown in figure 10(c). In vicinity of Fermi level, the contributions from C, Cu, S and Mo atoms are populated in PDOS. A hybridization is observed due to p of B, C and S, and d of Mo and Cu at 1.1 eV.

Four Cu substituted hexagonally by Mo and, one B and two C by S
The hexagonal geometry of Cu is maintained by keeping the concentration and sites of Cu and B same as discussed in section 3.2.9. The variation in bond length value has been noted down by increasing C concentration. A 2.2 Å bond value is observed between Cu and B atoms. The bond length values between Cu and C atoms is determined as 1.97 Å. The bond length values between Mo and S are observed in ranging from 2.36 to 2.61 Å. Thus, a prominent variation in bond lengths can be observed by playing with the concentration of C atoms. A modification in band gap has been noticed by altering concentration as shown in figure 11(b). Firstly the band gap is calculated by keeping the concentration in same ratio as narrated in section 3.2.9. Then the concentration of C is changed irrespective of the Cu and B. Resultantly, a remarkable decrease in band gap from 0.27 to 0.18 eV has been perceived.
The PDOS calculations indicate that a large overlapping is determined among the d of Mo and Cu, and p orbitals of S and C in energy ranging from −3 to 1 eV, while a small hybridization of the orbitals of B is also observed as shown in figure 11(c). An overlapping is determined due to the p of S and C, and d of Mo and Cu in energy area 0.6 to 1 eV. A peak is observed due to the orbitals of C, B, S, Cu and Mo at 1.2 eV.

3.2.11
. Four Cu substitute hexagonally by Mo and, two B and two C by S In this model, the concentration and sites of Cu and C are kept same as reported in the last case. The variation in bond lengths have been observed by increasing the B concentration as shown in figure 12(a). After geometry optimization, the bond lengths of d Cu1-C2 , d Cu2-B2 , d Cu3-C1 , d Mo-B1 , d Cu4-S and d Mo-B2 are noted to be 2.04 Å, 2.21 Å, 2.04 Å, 2.11 Å, 2.5 Å and 2.09 Å respectively. All bond lengths of d Mo-S are found to same as in the pristine MoS 2 . A similar hexagonal doping of Cu in MoS 2 , Be in h-BN and Be in graphene have been reported in earlier studies [26,36,40].
Band structure results are displayed in figure 11(b). A variation in band gap has been seen by playing with concentration. Firstly, the bandgap is measured by keeping the concentration in same ratio as narrated in section 3.2.10. Then the concentration of B is changed irrespective of the Cu and C concentrations. The surprising results have been noted by this strategy, as no variation in bandgap value has been identified.
The consequences of PDOS computations are shown in figure 12(c) that points out large contribution to PDOS arises from d of Mo and Cu, and p of B, C and S atoms in energy ranging from −3 to 0.8 eV. In the conduction band, hybridization occurs between B and C atoms in energy ranging from 1 to 2.5 eV.
3.2.12. Four Cu substituted rectangularly by Mo and, one B and two C by S In this case, the concentration of dopants is same as discussed in section 3.2.10 only the dopants sites are changed. Cu-C bond distance is found 2.01 Å. We observed a bond separation between Cu and S of 2.5 Å. The bond distance between Mo and B is determined to 2 Å as shown in figure 13(a).   Electronic band structure outcomes are exhibited in figure 14(b). These outcomes indicate that the highest value of bandgap is found as 0.28 eV above the Fermi level. A variation in band gap is noticed by keeping the same concentration of dopants discussed in section 3.2.12 but the locations of dopants are altered. Hence, the significant variations in band gaps have been determined after changing the sites and concentration of dopants.
The PDOS calculations indicate that a very large overlapping is found among the d of Mo and Cu, and p of S, B and C in energy ranging from −3 to −1 eV as shown in figure 14(c). An overlapping is determined due to the orbitals of S, C, Mo and Cu in energy ranging from −0.6 to 0.8 eV. In the conduction band, the population of PDOS arises from the p of B and S, p y of C and d of Mo and Cu in energy ranging from 1.1 to 1.3 eV. The p of C, B and S, and d of Mo show main contribution in conduction band.
Inspired by the investigations carried out by Ouma et al [41]. where they have reported the lanthanide substitutions in host MoS 2 and resulting magnetic moments indicating dilute magnetic semiconductor behavior; we also repeated some of our configurations employing spin polarized DFT calculations. The configurations reported in figures 1, 2, 7, 9, 10 and 14 have been recalculated as these models contains all the concentrations of Cu substituted (i.e. Cu 1-4 atoms) in addition to hexagonal and rectangular patterns. It was interesting to note that only two systems described in figures 2 and 9 have in respective order, the total magnetic moments of 1.08 and 0.54 μB. However, the magnetic moments on Mo were negligibly small. Additionally, the band gaps being main focus of this study are not altered. However, this interesting behavior could be observed by devoting a comprehensive study on MoS 2 by varying substituents, concentration and positions to ascertain the reported alike characteristics in MoS 2 which we intend to report in future.

Conclusions
Via first principles DFT calculations, we studied the geometrical and electronic properties of TMN doped MoS 2 . We substituted Mo by Cu and, B and C by S. We determined the effect of variation of concentration of impurities in addition to replaced position. We do note the changes in geometrical parameters of the host material after optimization procedure. However, from PDOS calculations we noticed that the impurity atoms make bonding with host MoS 2 advocating good interaction of the foreign elements with the pristine. The maximum change in Cu-C bond distance is determined 0.51 Å comparative to pure MoS 2 . We also observed a variations in bond angles and buckling heights of dopants. However, despite these variations, overall geometry of the system remains intact. However, the impurities cause to decrease the E coh as low as − 4.56 eV relative to − 5.2 eV of MoS 2 . Due to doping the impurity lines appear in the big energy gap of MoS 2 causing to divide it into small energy gaps. For TMN-doped structures, the change in concentration and site of impurities in doped structure modifies energy gaps. These efforts exhibit a variety of band gaps in the range 0.16 to 0.48 eV. Thus having achieved the band gap of desired value, MoS 2 become a valuable material for its use as semiconducting materials.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).