Substitutional transition metal doping in MoS2: a first-principles study

Single-layer MoS2 is a direct-gap semiconductor whose band edges character is dominated by the d-orbitals of the Mo atoms. It follows that substitutional doping of the Mo atoms has a significant impact on the material’s electronic properties, namely the size of the band gap and the position of the Fermi level. Here, density functional theory is used along with the G0W0 method to examine the effects of substituting Mo with four different transition metal dopants: Nb, Tc, Ta, and Re. Nb and Ta possess one less valence electron than Mo does and are therefore p-type dopants, while Re and Tc are n-type dopants, having one more valence electron than Mo has. Four types of substitutional structures are considered for each dopant species: isolated atoms, lines, three-atom clusters centered on a S atom (c3s), and three-atom clusters centered on a hole (c3h). The c3h structure is found to be the most stable configuration for all dopant species. However, electronic structure calculations reveal that isolated dopants are preferable for efficient n- or p-type performance. Lastly, it is shown that photoluminescence measurements can provide valuable insight into the atomic structure of the doped material. Understanding these properties of substitutionally-doped MoS2 can allow for its successful implementation into cutting-edge solid state devices.


Introduction
The insertion of charge carriers into a semiconducting material by extrinsic doping is necessary for its implementation into functional solid-state devices. In particular, substitutional doping of a semiconductor with a species of higher or lower valence than the host can introduce charge carriers in the form of electrons in the conduction band (n-type doping) or holes in the valence band (p-type doping). Joining two regions of opposite doping type forms a p-n junction, the principal element of electronic and photovoltaic technology. Meeting modern demands for both higher performance and broader application of solid-state devices requires fundamental changes to the incorporated materials. One prominent driver of this is Moore's law, which predicts that the packing density of transistors on a chip should double every two years [1]. This expectation has propelled rapid miniaturization of the p-doped channels that comprise the transistor gates, and recently has brought transistor technology to the point at which conventional semiconductors suffer from short channel effects at the cutting-edge length scales [2]. On another technological front, the demand for alternative energy sources, including solar energy, continues to grow. Expanded use of solar power would be vastly facilitated by appliances that are flexible and transparent, traits possessed by none of the currently prevalent semiconductors. Features such as these motivate the search for new semiconductors that can handle these evolving demands. As a result, the materials science community is looking to two-dimensional (2D) semiconductors to fill this void.
Perhaps the most aggressively studied materials to this end are the group-VI transition metal dichalcogenides (TMDs), having the chemical formula MX 2 , where M denotes a transition metal (e.g., Mo or W) and X represents a chalcogen (e.g., S or Se). These TMDs are layered semiconducting materials that have direct band gaps in their single-layer form [3,4]. Monolayer TMDs boast impressive mechanical strength and flexibility, [5] strong exciton binding energies, [6] and theoretical carrier mobilities comparable to that of Si [7,8]. They also exhibit a diverse set of tunable properties which depend on the transition metal and chalcogen species. For example, the band gaps of monolayer TMDs span a broad range of energies, meaning that heterostructures can be used to absorb light from a large portion of the visible spectrum [9]. Meanwhile, both vertical and in-plane combinations of monolayer TMDs of differing carrier type and concentration can produce atomically sharp p-n junctions [10]. TMDs are also attractive for photovoltaic device engineering, as their direct band gaps allow for fast radiative recombination and light absorption. They can absorb 5%-10% of incident sunlight, an absorption rate one order of magnitude higher than those of the currently used photovoltaic materials with subnanometer thickness, namely GaAs and Si [11]. These features have naturally led to the swift realization of TMD-based ultrasensitive photodetectors [12] and solar cell devices [13][14][15]. Simultaneously, extensive research has been devoted to implementing TMDs into field effect transistors (FETs). For instance, Radisavljevic et al have fabricated MoS 2 -based FETs possessing on/off ratios of over 10 8 [16]. Meanwhile, Nourbakhsh et al reported producing MoS 2 -based sub-10nm FETs, opening new possibilities for miniaturization [17].
To date, most of the aforementioned devices rely on the intrinsic doping behaviors of TMDs, which are believed to spawn from defect formation during the fabrication processes [18][19][20]. Unfortunately, such defects are difficult to control, and it has been reported that, within a single sample, the Fermi level can vary by up to 1eV between regions separated by only tens of nanometers [18]. Thus, recent efforts have been made to develop methods by which TMDs can be doped to reliably induce predictable doping behavior. One promising approach is substitutional doping of the transition metal atom. For a given TMD, both the conduction band minimum (CBM) and valence band maximum (VBM) belong to d-orbitals on the transition metals [21]. Therefore, one can anticipate that substituting a transition metal atom with a species possessing one more(fewer) valance electron than the host can induce n(p)-type doping. If the dopant's atomic radius is similar to that of the host, then substitutional doping places only minor strain to the crystal structure, while the host material retains a clean, planar surface with no dangling bonds. This allows for the construction of electrical contacts that are both atomically sharp and extremely stable [22].
While there have been numerous studies on transition metal doping in TMDs, substitutional transition metal doping has only recently been achieved. Gao et al have successfully doped WS 2 with Nb using chemical vapor deposition (CVD) with a concentration of about 6.7% by simultaneously introducing NbCl 5 and S into the CVD furnace. In the same study, the authors also demonstrated Re doping MoS 2 with a concentration of about 1% using an ReO 3 precursor [22]. Density functional theory (DFT) calculations show that Nb produces an acceptor state just above the VBM, which is indicative of successful p-type doping, while Re produces a donor state just below the CBM, demonstrating n-type doping. Other substitutional transition metal doping combinations have been achieved experimentally. At 1% concentration, transistors comprised of Ta-doped WS 2 and Nb-doped MoSe 2 show p-type performance and have been incorporated into transistors with on/off ratios greater than 10 6 [23,24]. Enhanced carrier concentrations have also been demonstrated with both Nbdoped and Re-doped WSe 2 [25]. Recently, pulsed laser deposition has been used to produce p-type Nb-doped WS 2 at concentrations of 0.5%-1.1%, which could allow for spatial control of doping concentrations [26]. These results suggest that transition metal doping is a viable method for production of p-n junctions, and thus, for device engineering.
Here, first-principles DFT and G 0 W 0 calculations are performed to examine the effect of transition metal doping in monolayer MoS 2 . The stabilities of several doping configurations are investigated to determine what sort of doping structures one can expect. For each configuration, the effect of doping concentration on the electronic properties, including band gap, electron mobility, and photoluminescence (PL) spectra are examined. The aim is to shed light on how various doping structures and concentrations can accommodate particular functions. While first-principles studies on substitutional transition metal doping in MoS 2 have been carried out before, [27][28][29][30][31][32][33] no other study has systematically investigated the effects of doping concentration and configuration on p-and n-type doping performance. A careful understanding of the atomic-scale doping structure and the corresponding properties can pave the way for controlled synthesis of doped TMD monolayers for the purpose of targeted functionality. Furthermore, one can expect that these results can be naturally extended to other semiconducting group-IV TMDs, whose pristine electronic structures are very similar to that of MoS 2 [9].

Methods
The local density approximation (LDA) Ceperley 1980 to DFT, as implemented in the Vianna ab initio Simulation Package (VASP), [34] was used to optimize the lattice vectors, ion positions, and electronic densities of all systems considered in this study. Structural relaxation calculations utilized the projector augmented wave (PAW) method, [35] with a basis set of 400eV, for which iterations persisted until all interatomic Hellmann-Feyman forces settled below 0.01 eV/Å. Simulation of doped systems with various doping concentrations adhered to the same convergence criteria. The Brillouin zone of the pristine 1L-MoS 2 primitive cell was sampled at k-points positioned on a 6×6×1 Γ-centered Monkhorst-Packgrid [36]. The k-point sampling for larger supercells was scaled down so that the number of subdivisions along a reciprocal vector scaled inverseproportionally to the scaling of the corresponding lattice vector. Dopant concentrations were mediated by the size of the supercell in which they were inserted. Supercell dimensions included all hexagonal Bravais lattices that can host an integer number of primitive cells.
A vacuum height of 18 Å served to minimize artificial interactions between periodic layers. Meanwhile, relaxation of in-plane lattice parameters and ionic positions gave the lattice constants listed in the DFT column of table 1. Considering that the LDA tends to underestimate lattice parameters by overestimating covalent bonding energies, [37], the calculated values are consistent with those of experiment [38,39].
The single shot G 0 W 0 was used to more accurately calculate the energies of the donor and acceptor states [40][41][42]. To this end, the response function is described by plane waves with energies up to 80 eV. Precise calculation of the frequency-dependent dielectric matrix requires the summation of a sufficiently large number of unoccupied bands. Therefore, the number of unoccupied bands always exceeded the number of occupied bands for all G 0 W 0 calculations. Lastly, PL spectra were obtained by calculating the matrix elements of the k-space derivative operator between each pair of Kohn-Sham orbitals.

Results and discussion
For each dopant species, four types of substitutional structures are considered: isolated atoms, lines, three-atom clusters centered on a S atom (c3s), and three-atom clusters centered on a hole (c3h) (see figure 1). For the remainder of this text, the label D:MoS 2 refers to MoS 2 doped with a particular dopant species D, where D=Nb, Tc, Ta, or Re. To begin, the stability of each doped structure is examined for each dopant species by calculating its formation energy with respect to pristine monolayer MoS 2 . Here, it is assumed that all transition metal atoms not doped in the MoS 2 crystal belong to their monolayer dichalcogenide crystal (e.g., Nb in pristine NbS 2 ). This is because, in a given growth environment, the precursors containing Mo, D, and S are all simultaneously available, allowing for the possibility of both MoS 2 and the DS 2 growth. Therefore, doping of D in MoS 2 to form D:MoS 2 comes at the cost of forming DS 2 . Thus, the formation energy per dopant atom is therefore calculated as dwhere E doped is the energy of the doped system, E pris is that of pristine monolayer MoS 2 , μ Mo is that of an Mo atom in pristine MoS 2 , μ dop is that of the dopant atom in its monolayer dichalcogenide form, and n is the number of substituted Mo atoms. While many conventions can be used for the assignment of the transition metals' chemical potentials, the authors feel that referencing the transition metals' monolayer dichalcogenide form most faithfully acknowledges the competition between the growth of pristine and doped MoS 2 in the growth process. Figure 2 shows that the formation energy of any doped structure is fairly insensitive to concentration. Meanwhile, as the formation energies of isolated dopants are relatively large for all concentrations, dopant atoms  [39] appear to be more stable when grouped together than when isolated, which indicates some degree of bonding interaction between neighboring dopant atoms. That is, the inter-dopant binding energy per dopant atom in the grouped structures can be taken as the difference between the formations energies of the grouped and isolated structures for a given concentration. This is further evidenced by the fact that the c3h structure is the most stable  for all dopant species and concentrations. This may be because the hole at the center of the dopant cluster can easily accommodate the optimal atomic spacing for inter-dopant bonding. By changing the bond angles between the dopants and neighboring S atoms, the three-atom cluster can expand or contract, allowing for the metallic inter-dopant bonds to assume their natural lengths. Provided that changes in bond angles have a lower energy cost than bond length variations, this flexibility should vastly reduce the formation energy of c3h structures. The results suggest that at higher dopant concentrations, cluster and line doping will be thermodynamically favored over isolated doping. Such clusters and lines have been previously observed experimentally [43] in Nb:WS 2 at high concentrations. It is therefore crucial that the electronic and optical properties of grouped dopant structures are well-understood. We note that the highest concentration of 33% simulated in this work is much higher than what has been achieved experimentally. In fact, at these concentrations, the simulated systems are more alike to distinct ordered crystalline solids than to doped MoS 2 . However, we have chosen to include these systems in our study to show how the material properties very smoothly with concentration all the way to the degenerate doping limit. The investigation of how each dopant species and structure affect the electronic properties of MoS 2 begins by examining the band structures of doped MoS 2 . As a prototypical example, the effects of isolated dopants are inspected in figure 3. By comparing the pristine band structure in figure 3(a) with those of isolated Tc and Re dopants in figures 3(e) and (f), it is seen that the n-type dopants each introduce a partially occupied band, highlighted in orange, residing just below the bulk conduction band minimum (CBM). In these cases, each dopant donates an extra electron to the conduction band that is weakly bound to the dopant nucleus by a screened Coulomb interaction. The weak binding slightly downshifts the donor state from the CBM of the pristine material. Such a small downshift is ideal for n-type doping for device applications, as the donated electron can be thermally excited into the bulk conduction band at room temperature and contribute to electrical conductivity. Meanwhile, figures 3(c) and (d) show that isolated Nb and Ta dopants each provide a partially filled state just above the valence band. This signifies that an state is weakly repelled from the dopant nucleus, producing an acceptor band above the valance band maximum (VBM), indicative of successful p-type doping. In all cases, isolated dopants appear to leave the direct gap nature of pristine MoS 2 intact. Here, we define the band gap as the energy difference between the highest energy state on the highest empty band with lowest energy state on the lowest filled band, as this transition should correspond to the major photoluminescnece (PL) peak if the doping is nondegenerate.
As grouped dopant structures are significantly more stable than isolated dopants, the electronic properties of the four dopant structures considered are compared in figure 4. In contrast to what is observed for isolated dopants, grouped dopant structures can change the local chemistry of the material and induce a richer variety of effects. In n-doped materials, the magnitude of the downshifts of the donor bands are qualitatively proportional to the binding energies between the donor electrons and the dopant nuclei. In the case of the three-atom clusters, the donated electrons are more tightly bound than they are for isolated dopants. This is because the clusters more closely mimic the electronic structure of bulk ReS 2 or TcS 2 , in which the transitional metals are metallically bonded to one another. Comparing figures 4(n) and (r) with (m) and (q) reveals this downshift. In the same light, the separation between the VBM and acceptor bands is larger for clusters than it is for the isolated dopant atoms, as evidenced in the upshift of the acceptor bands in figures 4(f) and (j) compared to those of in figures 4(e) and (i). The electronic structures shown in figure 4 allow for a straightforward calculation of each system's band gap by subtracting the highest energy from the highest fully occupied band from the lowest energy of the lowest fully unoccupied band. The results are plotted in figure 5. It is immediately seen that gap tends to grow as concentration increases. This can be understood by considering the inverse relationship between localization in position and momentum spaces. As the separation between doped regions increases (i.e., as the concentration decreases), the dopant bands become flatter, indicating that the corresponding orbitals are more localized on the dopant nuclei. This makes sense, as the bound donor electrons or acceptor holes are less likely to spread to neighboring dopant nuclei if the nearest neighbors are further away. It follows that higher doping concentrations produce more curvature in the bands that are affected by the presence of dopants. Bands with larger oscillations occupy a broader energy range and are more likely to cross the Fermi level, in which case they become partially filled (see figure S2). In doing so, the energy gap between the fully occupied and unoccupied states increases, so that the band gaps of the doped systems tend to gradually rise with growing concentration. As discussed below, these trends could be confirmed experimentally with PL measurements. Generally, an expansion of the band gap that accompanies increasing concentration would blue shift the PL peak, while the increase in the number of partially occupied bands would reduce the peak's intensity.
It is also evident in figure 5 that p-type dopants tend to increase the band gap from the pristine value. This can be explained by assessing the charges of the dopant nuclei. Nb for example possesses one fewer electron than Mo does, reducing the screening to the valence states. This increases the binding energy on the p-doped material's valence band states, reducing their energies and pulling them further from the conduction band, increasing the gap. The effect is amplified when dopant nuclei are grouped, as evidenced by the fact that the band gap of Nb:MoS 2 is minimized for most concentrations when the Nb atom is isolated.
The fact that doped systems considered possess an odd number of electrons per basic unit cell gives rise to the possibility of spin polarization. This is investigated in figure 4 where the spin-up and spin-down eigenvalues are plotted separately. In a given plot, the area of each blue circle is proportional to the projection of a spin-up state onto the dopant atom(s), while those of red circles are proportional to the corresponding spin-down projections. To ensure that spin-up states do not obstruct the spin-down states, for a given projection value, the area of a red circle is 50% larger than that of a blue circle. Therefore, a blue circle centered on top of a slightly larger red circle indicates that the state has spin degeneracy. That this is usually the case indicates that most of the studied systems have no spin polarization. However, Tc clusters are an exception, as illustrated by the fact that the red and blue circles within the pristine gap do not overlap in figures 4(n) and (o). This means that Tc clusters produce spin polarized mid-gap states. The origin of the spin polarization is somewhat mysterious for a few reasons. First, the Tc atoms in pristine single layer TcS 2 are not spin polarized (see figure S4). Additionally, it is somewhat surprising that spin polarization exists in clusters of Tc but not Re, as the two atoms belong to the same group in the periodic table, and are therefore expected to induce similar electronic properties. Lastly, the possibility that antiferromagnetic coupling between Tc atoms in the triplet leads to geometric frustration [44,45] is considered. However, Bader analysis [46] reveals that the magnetic moments of the Tc atoms are parallel (see figure S3), invalidating this hypothesis. Nonetheless, the possibility of spin-polarized dopant states are exciting for potential spintronic applications. Being localized to the dopant structure, the spin polarization of Tc-clusters also hints at the prospect quantum emitters with spin-control.
For the purpose of electronic device applications, whether or not a donor or acceptor state can effectively induce n-or p-type behavior is contingent on how tightly the state is bound to the impurity. For n-type doping, weaker binding of a donor state to the dopant reduces the energy gap between the donor state and the CBM, facilitating room-temperature excitation of the donor electron into the conduction band, increasing the material's conductivity. In the same way, a weakly bound hole state will lie close the VBM and contribute to conductivity as well. The greater the energy gap between the donor (acceptor) state and the CBM (VBM), the lower probability of excitation from those states to a conduction channel. It is therefore crucial to understand the precise position of the dopant states within the band gap to predict the dopant's effect on the electronic performance. To this end, one should acknowledge that the position of the mid-gap dopant states in figure 4 can be misleading, as DFT's prediction of unoccupied state energy are usually unreliable. To remedy this, the G 0 W 0 method is employed to better approximate the energies of the unoccupied dopant states [40]. From this, density of states (DOS) is plotted for both the LDA and G 0 W 0 methods in figure 6 for all of the systems considered in figure 4. As is typically reported in the literature, in all cases, the G 0 W 0 predicts a much larger band gap than does the LDA [9]. This improvement is due to a more accurate description of the electron-self energy [47] along with the proper treatment of the discontinuity in the energy derivative with respect to particle number [48]. Meanwhile, the behavior of the donor states is heavily dependent on the structure from which they arise. Figures 6(e), (i), (m), and (q) show that the gaps between the donor (acceptor) states and the CBM (VBM) are small in both the LDA and G 0 W 0 methods. Again, this suggests that n-and p-type doping will be successful if the dopant atoms remain isolated in the crystal. On the other hand, the dopant states for the three grouped structures all move towards the center of the gap when the G 0 W 0 approximation is included, limiting their effectiveness as n-or p-type dopants. Furthermore, close inspection of the mid-gap states in clustered dopant structures show that they are either completely filled or empty in the G 0 W 0 approximation. In these cases, the G 0 W 0 method has amplified the spin splitting of these states, so that one band resides entirely below the Fermi level and the other well above it. This is indicates that p-or n-type doping is ineffective in these structures, since the conductivity arising from partially filled bands is completely absent in these cases. Thus, in doping MoS 2 for electronic device applications, it would be beneficial to avoid grouping of dopant atoms.
The last feature investigated here is the presence of splitting in the midgap states of c3s Re:MoS 2 in figure 4(s). This arises because bulk ReS 2 is most stable in the distorted 1T phase, which has much less symmetry than the host MoS 2 [49,50]. Figure S5 in the SI shows that the bond angles subtended by the Re atom are distorted in the simulated structures containing c3s doped structures. The bending of bonds lifts the degeneracy in the donor states, causing a split in the donor bands. Interestingly, this distortion is only present at concentrations of 11% and 14%. This is likely an artifact of the particular configuration used to model these concentrations, as a precise reason as to why these concentrations produce a distortion while other concentrations did not cannot be pin pointed. However, while the circumstances in which Re clusters break the threefold symmetry of pristine MoS 2 remain unpredictable, it is clear that symmetry breaking is possible, and that one should therefore expect Re clusters to produce a relatively wide band of states inside of the pristine gap. It should be noted that bulk TcS 2 also adopts the distorted 1T phase in bulk form. However, the positioning of S atoms in Tc:MoS 2 seems to be dominated by the Mo-S interactions, which drives the system into an H-phase, so that the mid-gap states remain degenerate.
To connect these first-principles results with potential experimental observations, the emission spectrum for each dopant structure and species is calculated. The aim is to guide PL experiments to be able to establish the presence of various dopant structures within the material. After the excitation of an electron by photon absorption, it is assume that the excited electron relaxes non-radiatively to the minimum energy state of excited band before spontaneously emitting a photon via radiative relaxation to a lower band. The PL spectrum measures the frequency of this spontaneous emission (SE). The Einstein coefficient for the rate of SE from an initial state ñ i | to a final state ñ f | is given by [51] where w  fi is the energy difference between the states, c is the speed of light, α is the fine structure constant, and r is the position operator. As the emitted photons carry very little momentum compared to the dimensions of the the structures' Brillouin zones, it is assumed that the initial and final states ñ i | and ñ f | must possess the same crystal momentum. These emission rates are plotted as a function of the frequency of the emitted photon in figure 7. The rates of pristine MoS 2 shown in red are taken from the CBM to the VBM, while those of doped MoS 2 shown in blue are taken from all momentum-conserving transitions from the lowest energy states of any unfilled bands and the highest energy states of filled bands.
The first column of plots in figure 7 shows that in general, isolated dopant states will emit photons with energies very near the pristine PL peak. Meanwhile, the situation is more complicated for grouped structures, particularly the three-dopant clusters, whose states reside much deeper in the gap, increasing the number of possible radiative transitions. For example p-type dopant clusters, which have partially occupied midgap states to and from which an excited electron can relax. Relaxations involving these states entail smaller energy transitions, producing several red-shifted PL peaks, as shown in figures 7(f), (g), (j), and (k). This is consistent with previous experimental work [43], in which the PL peak of Nb:WS 2 was red shifted for a large doping concentration of about 6.7%. This high concentration allowed for the formation of grouped dopant structures, enabling emission of lower energy. In the same work, the PL peak of Re:MoS 2 at about 0.3% concentration was essentially unchanged from its pristine value. This is again consistent with the predictions here, as the low concentration of Re prohibited the grouping of dopant atoms, and isolated Re doping produces only a minor shift the PL peak in figure 7(q). This suggests that PL can be used to determine the structure of the substitutional doping and thus determine if the material is well-suited for electronic device applications.
One should bear in mind that the strong exciton binding energies of MoS 2 are not accounted for in these calculated emission spectra, and as a result, all experimental emission lines will be substantially red-shifted. Exciton binding energies can be treated accurately with the Bethe-Selpeter equation, [52,53] which has been shown to reduce calculated emission energies in MoS 2 by up to 1 eV [54]. Nonetheless, the results presented here should allow for a qualitative understanding of PL peak behavior in response to transition metal doping.

Conclusion
In summary, the dependence of electronic properties on the structure of the dopant atoms in transition metaldoped single-layer MoS 2 is demonstrated. It is found that c3h clusters are more stable than the other doping configurations examined at the same concentrations. Additionally, n-or p-type behavior is most effectively induced when the doping atoms remain isolated from one another. The combination of these two findings suggests that high doping concentrations should be avoided in order to prevent grouping from occurring at a significant rate. Meanwhile, emission calculations predict that isolated doping induces a blue shift in the PL peak, while grouped dopant atoms can prompt a red shift. Thus, it is proposed that PL can be a viable method for determining the merit of a transition metal-doped MoS 2 sample for electronic device fabrication. The authors believe that this information is valuable to those aiming to construct photovoltaic and electronic devices from doped MoS 2 .