Modulating the magnetic properties of MoS2 monolayers by group VIII doping and vacancy engineering

In this work, density functional theory is adopted to study the electronic and magnetic properties of MoS2 monolayers combined with a single S vacancy defect and a group VIII (G8) atom dopant, in which the dopant is incorporated via Mo substitution. The calculated results show that the magnetic properties of monolayer MoS2 can be tuned by changing the distribution of the G8 atom and S vacancy. The S vacancy tends to decrease the net magnetic moment of the doped system when these two defects are in their closest configuration. By adjusting the distance between the dopant and the S vacancy, the doped MoS2 monolayer may show a variable net magnetic moment. In particular, all of the Ni-doped MoS2 monolayers show zero magnetic moment with or without an S vacancy. The mean-field approximation is used to estimate the Curie temperature (TC). Our results show that Fe, Co, Ru, Rh, Os and Ir-doped MoS2 monolayers are potential candidates for ferromagnetism above room temperature. The density of states calculations provide further explanations as to the magnetic behavior of these doped systems. These results provide a new route for the potential application of atomically thin dilute magnetic semiconductors in spintronic devices by employing monolayer MoS2.


Introduction
Dilute magnetic semiconductors (DMSs) have been the focus of extensive research over the last decade, driven by the prospect of achieving spintronic devices, which can exploit both the charge and spin freedom. [1][2][3][4][5] Meanwhile, twodimensional (2D) transition metal dichalcogenides (TMDCs) have also demonstrated great potential for the new generation spintronics devices due to their unique structural and electronic properties. A signicant amount of theoretical and experimental effort has been devoted in the eld of TMDCs to comprehending the role of magnetic impurities such as V, Cr, Mn, Fe, Co, Ni, and Cu as discussed in several reviews. [5][6][7][8][9][10] Despite the great advances in this eld, producing strong ferromagnetic interactions and the stability of such ferromagnetic order is still the most important and most difficult part.
Recently, many studies, theoretical and experimental, have focused on the magnetic properties of 1H-MoS 2 . [11][12][13] Magnetic doping through either adsorption or substitution is found to be an efficient method to introduce magnetism into MoS 2 . [5][6][7] Magnetic interactions can be tuned by carriers and strain. At the same time, sulfur vacancies have also been found to be related to the magnetic properties of MoS 2 . 14 Vacancies not only inuence the magnetic properties, but also change the carrier density of the semiconductor. Obviously, S vacancies can inuence the magnetism in the transition metal doping case. Lots of studies have been done on each factor. However, it is not clear how these two factors interplay and then what the inuence of this interplay is on the magnetic properties.
In this paper, theoretical methods are used to study the effect of substitutional group VIII (G8) atom doping in 1H-MoS 2  Table 1 Optimized parameters for the G8-doped pristine and defective 1H-MoS 2 systems Mo 1Àx A x S 2Ày with x ¼ 6.25% (y ¼ 3.13%). The rootmean-square deviation (RMSD) of the X-S bond length (d RMSD inÅ), the formation energy (E form in eV), the total magnetic moments of the PBE results (M tot in m B ), the Bader charge (Q in e) of the G8 atom, and the energy band gap (E g in eV)

Computational methods
First principles calculations were performed using the Vienna ab initio simulation package (VASP) 15 on the basis of density functional theory (DFT). The electron-ion interactions were described by the projected augmented wave (PAW) method, 16 and the electronic exchange-correlation potential employed the generalized gradient approximation with the Perdew-Burke-Ernzerhof functional (GGA-PBE). 17 In order to verify our results, the Hubbard-U correction method 18 was also used. Different U values were assigned to the G8 impurities, while the U parameterization was not used for the host materials since they have little impact on the magnetic ordering, as suggested by many authors in ref. [19][20][21]. A 4 Â 4 Â 1 MoS 2 supercell structure containing 32 S and 16 Mo atoms was constructed as the pristine model in this calculation. Moreover, a vacuum region of 15 A was added along the c plane to minimize the interaction between the adjacent periodic images. The cutoff energy for the plane-wave expansion is set at 500 eV aer extensive convergence analysis. The Brillouin zone (BZ) is sampled using a 3 Â 3 Â 1 gamma-centered Monkhorst-Pack grid. The valence electron congurations in this calculation are Mo 4p 6 4d 5 5s 1 , S 3s 2 3p 4 , Fe 3d 6 4s 2 , Co 3d 7 4s 2 , Ni 3d 8 4s 2 , Ru 4d 7 5s 1 , Rh 4d 8 5s 1 , Pd 4d 10 , Os 5d 6 6s 2 , Ir 5d 7 6s 2 , and Pt 5d 8 6s 2 , respectively. All of the structures are fully relaxed using the conjugate gradient method, and during the structural relaxation, the energy convergent criterion is 10 À6 eV per unit cell, and the forces on all relaxed atoms are less than 0.02 eVÅ À1 .

Results and discussion
A. Tendency of the formation energy The fully relaxed lattice constant of monolayer MoS 2 is 3.19Å and the Mo-S bond length is 2.42Å, which agrees with the experimental value (3.16Å). 22 Defects play an important role in the electronic and magnetic properties of materials. Doping atoms or vacancies in 1H-MoS 2 are introduced by replacing or removing a single host atom in the pristine system. In this   This corresponds to simulating a Mo 1Àx A x S 2Ày periodic system with x ¼ 6.25% (y ¼ 3.13%). Furthermore, to characterize the deviation of the impurity atom from the original Mo position, we calculated the root-mean-square deviation (RMSD) of the X-S bond length using the following formula: where N represents the number of X-S bonds around the dopants. The value of N is 6, except for the nearest neighboring structures which have ve X-S bonds. d 0 is the Mo-S length in pristine MoS 2 , while d i is the X-S bond length in the doped model. The results are listed in Table 1.
In order to inspect the stability and feasibility of all optimized geometrical structures, the formation energies were obtained utilizing the following formula: 23-25 where m doped and E pure represent the total energies of the G8doped monolayer MoS 2 and pure one. m Mo , m S and m G8 are the chemical potentials for the Mo host, S host and G8 dopant atoms, respectively. All of the chemical potentials of the G8 elements are obtained from faced-centered cubic (fcc) structures except for Fe, which uses a body-centered cubic (bcc) structure. The formation energy of MoS 2 itself, E form (MoS 2 ), can be calculated from the expression: where m MoS 2 is equal to E pure per MoS 2 formula unit, and m 0 Mo /m 0 S is the total energy per atom of Mo/S in its reference phase. For Mo, the reference phase is the bulk bcc metal. The reference phase for S is the S8 ring, which is the most stable state at room temperature. Thus, the value of E form (MoS 2 ) is À2.70 eV in DFT-GGA, consistent with the value of À2.84 eV reported in ref. 26.
The values of m Mo and m S in eqn (2) depend on the experimental growth conditions. For the Mo-rich case, the Mo chemical potential is equal to the bulk Mo value, m Mo-rich Mo ¼ m 0 Mo , and the S chemical potential can be obtained from m MoS 2 ¼ m Mo + 2m S on the basis of thermodynamic equilibrium. Hence, combined with eqn (3), the chemical potentials for the Mo-rich limit can then be written as: Likewise, under S-rich conditions, the values are: All of the formation energies are listed in Table 1.
Monolayer MoS 2 with a single S vacancy (V S ) or Mo vacancy (V Mo ) was fully relaxed. The optimized structure shows that the neighboring Mo and S atoms have slight displacements with respect to the vacancy site V S and V Mo , which is different from the obvious reconstruction in a graphene sheet with a single C vacancy. 27 For V S -MoS 2 , this is more likely to occur under S-rich conditions and the formation energy is 6.57 eV, which in agreement with the previous values of 5.89 (ref. 28) and 5.72 eV, 7 is much lower than the formation energy of the Mo vacancy V Mo (14.09 eV) in Mo-rich conditions. The previous studies have reported that the substitution of an Mo site is more stable than that of an S site. 7,29 Experimentally, S vacancies are more common than Mo vacancies. 30 Therefore, in this paper, the co-doped congurations mainly consist of an S vacancy and G8 impurity substitution of an Mo site. The C 3v symmetry of pristine monolayer MoS 2 is destroyed aer G8 element doping, and the distances between the impurities and the nearest S atom change with different distributions of the dopants and the vacancy. All the calculated formation energies are summarized in Table 1. The distances between the dopants atoms and the S vacancy are taken from their original positions in pristine 1H-MoS 2 . Fig. 2 shows our calculated formation energies as a function of the distance between the S vacancy and dopant atoms, which varies from 2.42 to 7.53Å. The rst data "0" represents the doped MoS 2 without an S vacancy. According to our theoretical results, Fe doping is most favorable energetically among all considered impurities. For all of the G8 elements, the formation energies from the doped defective 1H-MoS 2 congurations of the second-, third-, fourth-, and h-nearest neighbor are very close, and are 1.50 eV larger than the corresponding nearest neighboring case. The relevant data for the Co atom are consistent with ref. 31.
The positive formation energy indicates that the formation of the vacancy defect and the doping of the transition metal are endothermic processes. The equilibrium concentrations of the vacancies are usually very low because of their high formation energies. Nonetheless, new techniques have been developed to create these defects. The generation of nanomesh size vacancies in graphene has been reported. 32,33 Vacancy engineering of a doped MoS 2 monolayer can also be achieved.

B. Magnetic properties
It is known that MoS 2 has D 3h symmetry and the d-orbitals have a schematic band structure with a d 2 conguration of the metal atom, as shown in Fig. 3. The two valence electrons of the Mo Paper ion occupy the lowest d z 2 orbital, which is the reason for the lack of magnetism for MoS 2 . In our case, the sulfur atom vacancy introduces two defective bands in the band gap region, which mainly consists of the "unsaturated" d orbitals of Mo atoms. The relatively large energy gap of MoS 2 can almost accommodate the ve d states of TM magnetic impurities. Therefore, about seven defective bands may appear in the forbidden gap region of the pristine bands when the doped G8 atoms are incorporated with an S vacancy. The magnetic properties mainly depend on the correlation between these defective bands.
The Fe case, which corresponds to the d 4 conguration, is likely to have a net magnetic moment of 2m B , except for the A case. In the nearest conguration, the interaction between the Fe and the S vacancy (the unsaturated Mo atoms) leads to strong crystal eld splitting between the d z 2, d xy , and d x 2 Ày 2 orbitals, which causes the system to become a small band gap semiconductor. Among all of the ve pristine dopant-vacancy congurations, the B and E cases have more symmetry and the vertical mirror plane is preserved. It is important for the system to keep the energy levels of the hybridized states below the Fermi level in the close energy range. This can explain the net magnetic moment of the B and E cases of Fe, and Co. On the other side, the C and D cases are more likely to have zero magnetic moment. Our results also indicate that local stress appears to be a crucial factor in the development of magnetism as Jahn-Teller distortions that destroy the C 3v lattice symmetry lead to the disappearance of magnetism. For the Fe and Co cases, the C and D congurations may also have a net magnetic moment because of the coupling between the defective d orbital of the unsaturated Mo atoms and the lower d xy , d x 2 Ày 2 , and d z 2 orbitals of the impurity atom.
For 4d or 5d dopants such as the Ru-doped cases, only the structure of mono-doped Mo 15 RuS 31 generates a 2m B magnetic moment. The origin of the magnetism is attributed to the near degeneracy of the Ru d x 2 Ày 2 and d z 2 orbitals. The Os-doping has similar characteristics to the Ru doped cases. For the Pddoping, only the Mo 15 PdS 31 -B structure generates a 2m B magnetic moment. For the Rh-and Ir-doping, all of the structures generate a 1m B magnetic moment. For the Pt-doping structures, all are nonmagnetic. Taken as a whole, the net magnetic moments from the 3d elements are larger than those from the 4d and 5d elements, because strongly localized 3d orbitals are more likely to induce net magnetic moments due to strong Hund coupling, which competes with the ligand eld energy splitting.
Furthermore, the ligand eld energy splitting is connected with the distortion of the local structure around the dopants. The root-mean-square deviations of the X-S bonds have been obtained to quantify the local deformation caused by the dopants and are listed in Table 1. The results clearly show that the RMSDs of the models with smaller magnetic moments are much larger, which means a larger deviation from C 3v symmetry. This rule holds for the Fe, Co and Os cases.
To verify our results, the DFT+U method was also used. It was found that the magnetic properties with the GGA+U method are consistent with the results without the Hubbard-U parameter. The results with the corresponding U values are shown in ESI  Table S1. † There is no relevant inuence on our conclusions and therefore we discuss in the following the results without onsite interaction.
The system is stabilized by charge transfer. The Bader charge has also been obtained to analyze the charge transfer between the G8 atoms and MoS 2 , as shown in Table 1. In a perfect 1H-MoS 2 , the formal valences of Mo and S are +4 and À2, respectively. Since a Mo has six S neighbors, it contributes 2/3 electrons to each Mo-S covalent bond. Therefore, the charge transfer of the S and Mo atoms is approximately +0.67e and À1.33e, respectively. Taking Fe-doped 1H-MoS 2 as an example, in the monodoping case, the charge number of Fe will be 8 À 2 Â 2/3 ¼ 6.67 (8 is the valence electron number of Fe treated by the PBE pseudopotential). It will be 8 À 5/3 Â 2/3 ¼ 6.89 when one S atom is removed. The theoretical charge numbers of Ru and Os are equal to that of Fe due to their similar valence electron numbers. For the other G8 elements, the two values for mono-and co-doping are 7.67 and 7.89 for Co, Rh and Ir, and 8.67 and 8.89 for Ni, Pd and Pt. The more the Bader charge deviates from the ideal charge number, the more electron transfer there is from the host ions to the dopant ion (Fig. 4).
In order to understand the characteristics of the magnetic moment induced by impurity states in the G8-doped 1H-MoS 2 systems, the spin density distribution is plotted to visualize the distribution of the magnetic moments of the doped 1H-MoS 2 systems in Fig. 5. Taking Fe-doped 1H-MoS 2 as an example, Fig. 5(a)-(f) show that the spatial extensions of the spin polarizations have reached the rst-nearest S atoms and the secondnearest Mo atoms. The structure Mo 15 FeS 32 displays weak ferromagnetic and antiferromagnetic coupling between Fe and three neighboring S atoms, and ferromagnetic coupling among the six second-nearest Mo atoms. Fig. 5(b) shows the nearest conguration of the co-doping system, which is typical for the no magnetism case. The antiferromagnetic alignment of d electrons is observed.

C. Estimation of the Curie temperature
The Curie temperature (T C ) values can be roughly estimated from the mean-eld expression 34 using the relation where k B is the Boltzmann constant and D is the difference in the supercell total energy between the antiparallel and parallel alignments. The results are listed in Table 2.
Our theoretical results show that Fe, Co, Ru, Rh, Os and Irdoped MoS 2 are candidates for room-temperature ferromagnetic materials in DMSs. It is should also be noticed that the mean eld approximation used in this paper is known to overestimate the Curie temperature. 35 Nevertheless, a high D will imply a high T C and the trends presented here should hold.

D. Density of states
To gain further insight into the emergent magnetic behavior, the total density of states (TDOS) and projected density of states (PDOS) of all of the models are shown, which also provides a specic description of the ground state electronic structure. The upper and lower panels denote the TDOS and PDOS, and positive and negative values represent spin-up and spin-down channels, respectively. The Fermi level is set at zero energy, which is indicated by black dotted lines, to easily identify the band gap and the relative position of the states from the impurity atoms. Fig. 6 shows the DOS and PDOS of pristine 1H-  Fig. 7. It can be found that the impurity states presented in the band gap region are mainly contributed by the Mo 4d, Fe 3d, and S 3p orbitals. All of these systems are ferromagnetic with net magnetic moments of 2m B , except for the structure Mo 15 FeS 31 -A, which is a nonmagnetic semiconductor with a 0.3 eV band gap. In conguration A, the largest deformation of the Fe atom can be found from visualizing the optimized structure, which results in the removal of the degeneracy of the 3d orbitals and leads to the formation of the localized hydride states near the valance bands. As for the ve magnetic systems, as expected, the Fe dopant is a main contributor to the total magnetic moment. This feature applies to all G8 impuritydoped congurations in this calculation. For the Co-doped cases, these systems have odd electrons. They are all ferromagnetic with a magnetic moment of 3 or 1m B . Unlike the Fe case, the correlation between the Co atom dopant and the S vacancy causes the impurity states to spread in the whole band gap region. Especially for the congurations Mo 15 CoS 31 -A and Mo 15 CoS 31 -D, the strong interaction between the Co-Mo-S pairs leads to the formation of a localized impurity state near the valance bands. The decrease in the magnetic moment can be attributed to the pairing of the electrons in this mixed orbital. In other words, the magnetic properties of the codoped system depend on the competition between the ligand eld splitting and the Hund coupling (Fig. 8).
For Ni-doped 1H-MoS 2 systems, there are 4 unpaired electrons remaining, which would lead to about a 4m B magnetic moment. But from our calculations, all of the ground states of the Ni-doped models are nonmagnetic semiconductors where the spin-up and spin-down channels are completely symmetrical and they possess tiny band gaps of 0.28, 0.43, 0.15, 0.30, 0.27, and 0.27 eV, respectively. The strong crystal eld of neighboring Mo and S atoms leads to large energy splitting of the defective orbital. Therefore, four outer electrons form two electron pairs without the existence of isolated electrons. During our calculations, we found that the Mo 15 NiS 32 and Mo 15 NiS 31 -B, C and E congurations would relax to transition states that have high magnetic order (net magnetic moment 4m B ). These transition states are local minima which are 199, 341, 355, and 208 meV above the global minima. In a recent experimental study, 9 the author reported that 4% Ni doped MoS 2 has a paramagnetic phase at room temperature, and the paramagnetic phase may dominate at low temperature. Our theoretical results of the transition states with an S vacancy provide an explanation for this phenomenon (Fig. 9).
For the Ru-doped cases, substitution of Mo atoms by Ru increases the degree of p-d hybridization leading to a shi of the majority spin below the Fermi level. All congurations with the existence of S vacancies have semiconductor character with a band gap shown in Table 1, while the structure Mo 15 RuS 32 is a magnetic semiconductor with a net magnetic moment of 2m B . The magnetism of Mo 15 RuS 32 can be attributed to the formation of degenerated bands mainly consisting of the Ru d x 2 Ày 2 and d z 2, Mo 4d and S 3p orbitals (Fig. 10).
Rh-doping provides 3 more valence electrons than the host Mo atom and all structures generate a 1m B magnetic moment. In addition, the mono-doped structure is a half-metal and the others are magnetic semiconductors (Fig. 11).
For Pd-doped 1H-MoS 2 , the structure Mo 15 PdS 31 -B is a magnetic semiconductor with a 2m B magnetic moment, although Pd is a nonmagnetic element. The other systems are The reason for the magnetism is the Hund coupling between the two hybridized states from the Mo 4d, S 3p and Pd 4d orbitals (Fig. 12).
Os is a 5d element with less localized d orbitals than the 3d and 4d elements discussed above. For the Os-doped 1H-MoS 2 , the impurity states are located near the conduction bands. The  (Fig. 13).
For Ir-doping, the DOS are very similar to the cases of Rhdoping, since the Ir atom belongs to the same column of the periodic table and has the same valence electron conguration as Rh. The difference is that the structure Mo 15 IrS 32 is unambiguously half-metallic and can be utilized as a spin lter. Similarly, all structures possess a 1m B net magnetic moment (Fig. 14). From the DOS gures, we can nd that the d bands of the Pt ion mainly appear above the Fermi level and the strong hybridization of the S 3p and Pt 4d orbital, and the electrons tend to occupy the delocalized defective bands that are caused by the S vacancy (Fig. 15).

Conclusion
In summary, based on rst-principles calculations, we have investigated the effect of co-doping of G8 atoms and an S vacancy on the magnetic and electronic properties of monolayer MoS 2 . All of the mono G8-doped systems have magnetism, except for the Ni, Pd and Pt cases. According to the formation energy calculation, it is found that the nearest neighboring congurations for every doped system are preferred, and have lower formation energy (about 1.50 eV) than the remaining four congurations. Moreover, from the magnetic calculations, we nd that the magnetic moment from the co-doping of Fe and an S vacancy (Fe-V S ) is 2m B , except for the nearest neighboring conguration, which is a nonmagnetic semiconductor. For the Co-V S , the total magnetic moments are 1 or 3m B . Particularly for the Ni-V S case, all of the doped models are nonmagnetic semiconductors. The high magnetic order structures with a net magnetic moment of 4m B are transition states that are about 200 meV above the ground state. In the case of Ru-V S , only the structure without an S vacancy has a net magnetic moment. For Os-V S , the Mo 15 OsS 31 and Mo 15 OsS 31 -C congurations have a net magnetic moment of 2m B . For Rh-V S and Ir-V S , all of the congurations are magnetic and the total magnetic moments are 1m B . Moreover, Ir-doped pristine 1H-MoS 2 is a half-metal. In the case of Pd-V S , most congurations are nonmagnetic semiconductors except for the second-nearest neighboring conguration, which possesses a 2m B magnetic moment. For Pt-V S , all of the models are nonmagnetic semiconductors. Furthermore, the Curie temperatures (T C ) are calculated within mean-eld approximation. Our theoretical results show that Fe, Co, Ru, Rh, Os and Ir-doped MoS 2 monolayers with particular distributions exhibit room-temperature ferromagnetism.
The magnetic properties of the co-doped system depend on the competition between the Hund coupling and the ligand eld splitting. Our results suggest that the co-doping of G8 atoms and S vacancies is an efficient way to modulate the magnetic properties. The main obstacle is the control of the distribution. Our further work will focus on methods to decrease the formation energy of the defects and the maintenance of the magnetic moment.

Conflicts of interest
There are no conicts to declare. This journal is © The Royal Society of Chemistry 2018