Magnetic order and critical temperature of substitutionally doped transition metal dichalcogenide monolayers

Using first-principles calculations, we investigate the magnetic order in two-dimensional (2D) transition-metal-dichalcogenide (TMD) monolayers: MoS2, MoSe2, MoTe2, WSe2, and WS2 substitutionally doped with period four transition-metals (Ti, V, Cr, Mn, Fe, Co, Ni). We uncover five distinct magnetically ordered states among the 35 distinct TMD-dopant pairs: the non-magnetic (NM), the ferromagnetic with out-of-plane spin polarization (Z FM), the out-of-plane polarized clustered FMs (clustered Z FM), the in-plane polarized FMs (X–Y FM), and the anti-ferromagnetic (AFM) state. Ni and Ti dopants result in an NM state for all considered TMDs, while Cr dopants result in an anti-ferromagnetically ordered state for all the TMDs. Most remarkably, we find that Fe, Mn, Co, and V result in an FM ordered state for all the TMDs, except for MoTe2. Finally, we show that V-doped MoSe2 and WSe2, and Mn-doped MoS2, are the most suitable candidates for realizing a room-temperature FM at a 16–18% atomic substitution.


INTRODUCTION
The recent experimental realization of two-dimensional (2D) magnetic crystals like CrI 3 1-3 , CrGeTe 3 4 , and VSe 2 5 , has sparked great interest for possible applications like spintronics 6,7 , valleytronics 8 , and skyrmion 9 -based magnetic memories 10,11 . However, magnetic order in 2D magnetic crystals is hampered because of low magnetic anisotropy and weak exchange interaction strength 12 , resulting in low Curie temperatures (e.g., 45 K for CrI 3 and 42 K for CrGeTe 3 ), which limits their use in commercial applications.
Semiconducting 2D materials doped with impurity atoms, e.g., doped graphene 13 and metal-doped transition-metal dichalcogenides (TMDs) [14][15][16][17][18] , have emerged as promising candidates for high-temperature 2D magnetic order. Semiconductors doped with transition metal impurities couple the properties of semiconductors and magnets and are called dilute magnetic semiconductors (DMS). The ability to control magnetic order through charge transfer in a DMS [19][20][21] opens up the possibility for realizing magnetic devices because their magnetic state can be controlled using an external electric field 22 .
Among semiconducting 2D materials, which can be used as a base material for 2D DMSs, TMDs are of special interest. The heavy atomic mass of TMDs can lead to larger magnetic anisotropy, which is necessary for the existence of magnetic order in 2D 23 . The interest in TMDs is further fueled by recent experimental results that have demonstrated the existence of stable magnetic order in TMDs doped with a transition metal impurity [24][25][26] .
For the technological application of magnetically doped TMDs, it is necessary to find the optimal combination of a TMD and a dopant. However, the number of possible TMDs and dopant combinations is too large for a comprehensive experimental investigation, and theoretical guidance is desired. There have been previous theoretical attempts at modeling doped TMDs and calculating their critical temperature [14][15][16]27 . However, previous theoretical predictions of the Curie temperature of doped TMDs have predicted unrealistically high Curie temperatures in excess of 1000 K at low doping concentrations (≈5%) [14][15][16]27 , whereas experimental observations to date suggest a Curie temperature below 350 K at such doping levels 24,25,28,29 .
The reason for the discrepancy between the experimental and the theoretical work is that previous theoretical works have ignored the effect of magnetic anisotropy 30 and have used the collinear magnetic approximation 31 . Most of the previous theoretical works have either calculated the Curie temperature using the Ising model 15 , or the mean-field approximation 14 . Unfortunately, both methods (Ising and mean-field) result in an overestimation of the actual Curie temperature 32 . Moreover, previous theoretical works have taken into account only a few combinations of dopants and TMDs, and a comprehensive study for a vast set of TMDs doped with period four transition metals is still missing.
In this work, we calculate the critical temperature (Curie temperature or Kosterlitz-Thouless (K-T) transition temperature) for a set of 2H TMDs: MoS 2 , MoSe 2 , WSe 2 , WS 2 , and MoTe 2 , substitutionally doped with all the period four transition metals starting from Ti to Ni (Ti, V, Cr, Mn, Fe, Co, Ni). To model the magnetic structure of doped TMDs, we use our method, developed in ref. 12 . We model the magnetic exchange interaction (J) using a parametrized functional form (J(r)) and fit its parameters to first-principles calculations. Finally, we calculate the concentration-dependent critical temperature using the Monte-Carlo method. First, we apply our method to one of the TMDs, MoSe 2 , doped with all period four transition metals, and show all the possible magnetic ordered states originating from different dopants. Next, we present the magnetically ordered states and the critical temperature for all the TMDs doped with period four transition metals. We find that out of the thirty-five material combinations investigated, ten are non-magnetic (NM), nine are anti-ferromagnetic (AFM), and sixteen are ferromagnetic (FM). Out of the 16 FMs, 6 are FMs with an in-plane magnetic easy-axis, and the other 10 have an out-of-plane easy-axis. We find that the most promising FMs can be realized by doping MoSe 2 and WSe 2 with V along with doping MoS 2 with Mn, which have a Curie temperature of approximately 200 K at an atomic substitution of 15%. Figure 1 illustrates our computational model, which has two parts: DFT and Spin-Model. To reduce the computational cost and weed out the candidate ferromagnets, in the first part (DFT), we determine the magnetic ground state of the doped TMD by substituting the transition metal with two dopants in a supercell of size 3 × 3 × 1. If out of all possible magnetic configurations, the ferromagnetic state has the lowest energy, we make bigger supercells (4 × 4 × 1, 5 × 5 × 1, and 7 × 7 × 1) of the corresponding TMD and dopant. We substitute two transition-metal atoms (W/ Mo) in the TMD supercells with dopant atoms separated at distances ranging from the nearest neighbor to the fifth neighbor. We calculate the total energy for both the FM and the AFM magnetic orders with both the in-plane and the out-of-plane magnetic easy axis for the bigger supercells of TMDs 12 .

Computational model
In the second part (Spin-model), we model the magnetic structure of doped TMDs using a classical Heisenberg Hamiltonian, which features parameterized exchange tensor J i,j , and onsite anisotropy D determined from DFT 12 . Specifically, we approximate the elements of J i,j tensor as a continuous function of distance J(r i − r j ) (see Eq. (2) in the Methodology section). We go beyond the nearest-neighbor interaction because long-range interaction plays a decisive role in determining the magnetically ordered state of doped materials 33,34 . We take into account the exchange interactions up to the 5th neighbor (N = 5), beyond the 5th nearest neighbor the exchange interaction (J(r)) is numerically truncated.
To obtain the parameters for J(r), we fit the parameterized Heisenberg Hamiltonian to the total energy obtained in the DFT step. The details of our fitting procedure are outlined in ref. 12 . We study the phase change of the parameterized Heisenberg Hamiltonian for large (40 × 40) supercells with an atomic substitution ranging from 6% to 18%, using the Monte-Carlo (MC) algorithm. We obtain the median critical temperature (Curie/ K-T) from the peak of the average susceptibility for each percent atomic substitution, obtained from the MC simulations. For each percent substitution, we average over 20 different substitutional configurations to account for configurational entropy. We provide further computational details and parameters at the end of this article.

Magnetic order in MoSe 2
We first apply our computational method to MoSe 2 doped with all period four transition metals. We present all possible magnetically ordered states in doped MoSe 2 : the non-magnetic state (NM), the ferromagnetic state with out-of-plane spin polarization (Z FM), the out-of-plane polarized clustered FMs (clustered Z FM), the inplane polarized FMs (X-Y FM), and the anti-ferromagnetic state (AFM). We then apply our method to all the TMDs doped with period four transition metals and calculate their critical temperature and uncover their magnetically ordered phase at the atomic substitution of 15%.
We perform DFT calculations on a 3 × 3 × 1 cell of MoSe 2 doped with two dopant atoms (Ti, V, Mn, Cr, Fe, Co, Ni), which amounts to 22.22% atomic substitution. We define the atomic substitution using, where N dopant is the number of dopant atoms and N transition-metal is the total number of transition metal atoms in the TMD supercell. We find that Ni and Ti as dopants result in a negligible magnetic moment of 0.09 μ B per dopant. Hence, Ni and Ti as dopants result in a weakly-magnetic/non-magnetic (NM) state in MoSe 2 . On the other hand, V, Cr, Mn, Fe, and Co have a magnetic moment of 1.7 μ B , 3.0 μ B , 3.7 μ B , 3.1 μ B , and 1.6 μ B per dopant, respectively. We find that the magnetic easy-axis is inplane for Fe and Mn, whereas it is out-of-plane for V and Codoped MoSe 2 . Figure 2a shows the concentration-dependent critical temperature (Curie temperature for out-of-plane FM and Kosterlitz-Thouless transition temperature for in-plane FM) of MoSe 2 doped with V, Mn, Fe, and Co. We observe that V-doped MoSe 2 exhibits room-temperature out-of-plane FM at an atomic substitution of about 16.5%, and Fe-doped MoSe 2 exhibits room-temperature inplane FM at an atomic substitution of 16%. We also observe that the variance in the obtained critical temperatures is very low across different substitutional configurations (<30 K) except for Co doping. Low variance implies that the critical temperature is robust to the random position of dopants in V, Fe, and Mn-doped MoSe 2 . Figure 2b shows the saturation magnetization (M = √(S x 2 + S y 2 + S z 2 )) per dopant atom in V, Mn, Fe, and Co-doped MoSe 2 , obtained from the MC simulations at an atomic substitution ranging from 6% to 18%. The dotted lines show the starting magnetization of each magnetic dopant at the start of the MC simulation (also shown in Table 1). We observe that the saturation magnetization in V-doped MoSe 2 remains almost flat, and decreases slightly at a percent substitution below 7%. For Fe, the saturation magnetization increases with increasing atomic substitution and saturates to its maximum magnetization at around 18% atomic substitution. However, for both Mn and Co-doped MoSe 2 , the saturation  Figure 3a shows the magnetically ordered state in a sample of V, Cr, and Co-doped MoSe 2 at an atomic substitution of 15%, at a temperature of 5 K. We plot the magnetization in the out-of-plane direction ðŜ z ¼ S z =jSjÞ. The magnetically ordered state in V and Co-doped MoSe 2 is FM with an out-of-plane easy axis, whereas, for Cr substitution, the magnetically ordered state is AFM. We observe that the magnetically ordered state of V-doped MoSe 2 saturates to a perfect FM state with an out-of-plane easy-axis. For Co dopants, we find that there are clusters of FM-oriented Co ions, but the long-range order is missing. In the case of Cr, we observe that the magnetically ordered state has a randomized magnetic order with spins orienting randomly. The reason for such a randomized magnetically ordered state is that the AFM order leads to a magnetic frustration for some clusters, which leads to a randomized orientation for the magnetic moments. Figure 3b shows the magnetically ordered state in a sample of Fe and Mn-doped MoSe 2 at an atomic substitution of 15% at a temperature of 5 K. Because the magnetic easy-axis is in-plane, we plot the in-plane angle (ϕ) of each dopant atom: ϕ ¼ cos À1 ðS x =jS k jÞ, where, S ∥ is the in-plane magnetization. We observe that for both Mn and Fe, the orientation of the magnetic moment remains randomized with some short-range (<8 Å) order. For both Fe and Mn, we observe two effects. First, we observe domain formation with FM clusters. Second, the FM clusters are not perfectly ferromagnetic because their magnetic order breaks at slightly longer distances. We will discuss domain formation separately later in this section. The slight breaking of magnetic order at a longer distance appears due to the Kosterlitz-Thouless transition behavior 35 , where, even at temperatures below the K-T transition temperature, long-range magnetic order does not exist. Interestingly, we see some level of quasi vortex formation both for Fe and Mn-doped MoSe 2 . The observation is in line with the K-T physics for magnets with in-plane anisotropy 35 . However, due to broken translational symmetry because of the random placement of defects, as well as finite in-plane anisotropy, vortex formation is imperfect. Nevertheless, observation of quasi vortices in doped MoSe 2 implies that the topological vortices in the K-T transition are quite robust to lattice imperfections. However, an in-depth analysis of this phenomenon is beyond the scope of this article. Figure 4 shows the exchange function J(r) for doped MoSe 2 . The solid line shows the average between the exchange between electrons with moments along the x and z-direction (J(r) = (J x (r) + J z (r))/2), and the dots show the calculated discrete J parameters obtained from the DFT calculations using the J 1 − J 2 model 31 . We define a distance >7 Å as long-range, and <7 Å as short-range (up to the third-nearest neighbor).
We observe that J(r) is the strongest for V, Fe, and Co-doped MoSe 2 in the short range. For Mn, J(r) is weaker in the short range. Whereas, for Cr, J(r < 7 Å) is negative, signifying an AFM interaction. Looking at the long-range interaction beyond 8 Å (inset of Fig. 4), we find that the long-range interaction is strongest in V-doped MoSe 2 and significantly weaker in Co and Mn-doped MoSe 2 . In Fe-doped MoSe 2 , the long-range interaction is the second highest.
Analyzing the ordered magnetic states for each dopant shown in Fig. 3a and b, their saturation magnetization in Fig. 2b, and their J(r), we find that a strong short-range interaction results in FM cluster formation, e.g., in Co, Fe, and Mn-doped MoSe 2 , and the clusters then orient randomly due to weak long-range inter-cluster interaction. At higher percent atomic substitution in Fe-doped MoSe 2 , we find that bigger clusters start forming as seen by the increase in saturation magnetization for Fe-doped MoSe 2 at higher atomic substitution, as shown in Fig. 2b. However, the X-Y nature of magnetism prohibits the long-range order in Fe-doped MoSe 2 . It should also be noted that the appearance of a saturation magnetization at these concentrations in the X-Y magnets is due to the finite size of the lattice, the random position of the magnetic ions, and the long-range behavior of J(r).
To summarize this section, we find that five magnetically ordered states are possible, depending on the J(r) in a doped TMD including,  Critical temperature of MoS 2 , WS 2 , MoSe 2 , WSe 2 , and MoTe 2 Figure 5 shows the magnetically ordered state, as well as the critical temperature at an atomic substitution of 15% for the selected combination of the TMD and the dopant. Material combinations indicated with "Z" are ferromagnetic (FM) with their magnetic easy-axis in the out-of-plane direction, while "X-Y" indicates the X-Y magnets, with an in-plane easy axis. The TMD dopant combinations, which are shaded, are clustered Z FM. Materials indicated as AFM have an anti-ferromagnetic (AFM) ordered state, while (NM) represents a non-magnetic (paramagnetic) state. As discussed above, FMs with an out-of-plane magnetic easy axis have a Curie temperature, while, FMs with an in-plane magnetic axis have a Kosterlitz-Thouless phase transition, which results in a quasi-ordered magnetic state, but the long-range order remains randomized. We provide the full atomic substitution-dependent (6% to 18%) critical temperature for all the Z and the X-Y ferromagnets in the supplementary documentation (Supplementary Table. 2).  Fig. 2).  Some general trends can be extracted from Fig. 5. We observe that Cr dopants result in an AFM-ordered state for all the TMDs. Ti and Ni as dopants result in a non-magnetic state. Interestingly, V, Fe, and Mn always result in an FM ordered state for all the TMDs except for MoTe 2 , for which only Fe and Co result in an FM state. Co dopants result in an AFM ordered state for disulfides (MoS 2 , WS 2 ), but they result in an FM ordered state for diselenides (MoSe 2 , WSe 2 ), and MoTe 2 . V dopants always yield an out-of-plane FM. Whereas, Mn and Fe substitution result in an X-Y magnetic order for WS 2 , MoSe 2 , and WSe 2 .
The main highlights of Fig. 5 are the material combinations that result in an out-of-plane FM with strong long-range interaction. We find five combinations with V as a dopant for all the TMDs except MoTe 2 and Mn-doped MoS 2 . Mn substitution in MoS 2 results in an FM-ordered state, with high median Curie temperature of 190 K at 15% atomic substitution. Also, V as a dopant in MoSe 2 and WSe 2 results in an out-of-plane FM with a high Curie temperature measuring ≈200 K.
Finally, we would like to mention that the electronic origin of magnetism in doped TMDs is a result of the super-exchange interaction in the short-range 14,16 and carrier-mediated interaction in the long range. For example, the electronic origin of magnetism and the doping stability in MoSe 2 are briefly discussed in the supplementary document ( Supplementary Fig. 4).

DISCUSSION
We have presented the magnetic order in TMDs doped with period four transition metals. We have determined the nature of the magnetically ordered states, as well as their critical temperature as a function of percent atomic substitution. We showed that there are five possible magnetically ordered states for doped TMDs, depending on the nature of their exchange interaction J(r), the magnetic anisotropy, and the atomic substitution. The possible magnetically ordered states are non-magnetic (NM), perfectly ferromagnetic (FM Z), clustered ferromagnetic (clustered FM Z), X-Y ferromagnetic (X-Y FM), and randomized anti-ferromagnetic (AFM).
We have shown that Ti and Ni dopants always result in a nonmagnetic state. Moreover, Cr dopants result in an AFM configuration for all the TMDs. Both Mn and Fe dopants result in an X-Y magnet for MoSe 2 , WS 2 , and WSe 2 . From this study, we conclude that the best chance of realizing a 2D DMS using TMDs with roomtemperature Curie temperature is found in Mn-doped MoS 2 and Vdoped MoSe 2 and WSe 2 at an atomic substitution in excess of 16.5%.
We have provided a generalized method of modeling the magnetic interaction in doped 2D materials. For further usability of our method, the parameters of the functional form for all the TMD and dopant combinations are provided in a supplementary document (Supplementary Table 1).
There have been recent experimental reports regarding magnetic order in TMDs [24][25][26] , and FM clusters have been detected in V-doped WSe 2 , using magnetic scanning tunnel microscopy (MTM) 26 . The magnetically ordered states presented in this work for transition-metal doped TMDs, and their critical temperature can be verified experimentally using a similar procedure as used in 26 . Moreover, recent experimental reports have shown that the transition-metal substitution is often accompanied by vacancies 36,37 , and a possible future extension of our work will be to include the impact of structural defects on the magnetic order of TMDs.

Magnetic structure and the exchange interactions
We model the magnetic structure of doped TMDs using a parameterized Heisenberg Hamiltonian assuming a localized nature of the magnetic interaction 38 , The first term is the exchange term between the i th and the j th magnetic atom (dopant) with S = S x x + S y y + S z z, as the magnetic moment vector. J i,j is the strength of the exchange interaction between the i th and the j th magnetic atoms and is a tensor as described in ref. 12 . Because anisotropy plays an important role in determining the magnetic ground state of doped magnetic systems 39 , we take into account the J i,j tensor instead of an effective isotropic exchange. The second term is the onsite-anisotropy with strength D. We only use the diagonal elements of J i,j which are, J xx , J yy , and J zz for the magnetic axis in the x, y, and z direction, respectively. Because of the in-plane isotropy in TMDs, we modify the J i,j tensor by choosing, J xx = J yy = J ∥ and J zz = J ⊥ , with J ∥ being the in-plane exchange interaction, and J ⊥ being the out-of-plane exchange interaction. We approximate J i,j as a function of distance J(r) because we go beyond the nearest-neighbor interaction. The functional form is, c i B i ðrÞ hðr c À rÞ þ c ?=k expðÀr=λÞ hðr À r c Þ: Here, h(r) is the Heaviside step function. r c is the cut-off radius within which we approximate the J parameters using B-splines B i (r) 40 with order 3, and outside r c we approximate them using an exponential decay c ?=k expðÀr=λÞ 33,34 . Parameters A ⊥/∥ and c i are the free parameters. We choose r c to be within the third nearest-neighbor range, which is 7 Å for all the TMDs. Because of the continuity at the boundary r = r c , the parameter c ⊥/∥ and λ have an analytical form in terms of the spline functions, Note that, traditionally, going beyond the nearest-neighbor interaction increases the number of parameters as 2N, where N is the interaction range. For example, for next-neighbor interaction, N = 2, and the total number of J i,j parameters required to model the magnetic structure is 4. However, in our method, we take into account the exchange interactions up to the 5th neighbor (N = 5). Thanks to the functional form (J(r)) we use, the number of free parameters remains fixed to five. The parameters of the functional form for all the TMDs and the dopant combinations are provided in a separate supplementary document (Supplementary Table 1).
It should be noted that the generalized functional form of the Eq. (2) is useful for materials with in-plane rotational invariance. Materials with broken rotational symmetry require an additional parameter to account for the angular dependence of J(r).
The magnitude of the magnetic moment for the Monte-Carlo simulations are obtained from DFT using Here, M l;j DFT is the average of the magnetic moment of the jth magnetic atom of lth magnetic configuration obtained from DFT. N c and N are the total number of the magnetic configurations simulated and the magnetic atoms, respectively 12 .

DFT calculations
All the ab-initio DFT calculations reported in this work were performed using the Vienna ab-initio simulation package (VASP) 41,42 . The ground state self-consistent field (SCF) calculations were performed using a projectoraugmented wave (PAW) potential 41 with a generalized-gradient approximation as proposed by Perdew-Burke-Ernzerhof (PBE) 43 . We have used a kinetic energy cut-off of 450 eV for our DFT calculations. The Brillouin zones were sampled using a Γ-centered k-point mesh of size 5 × 5 × 1 points for 4 × 4 × 1 supercells, and 3 × 3 × 1 points for 5 × 5 × 1 and 7 × 7 × 1 supercells. The TMD supercells doped with transition metals were relaxed until the force on each of the ions was below 10 meV/Å. The energy convergence criterion for the subsequent SCF calculations was set to 10 −4 eV.
We have used the Hubbard U model within DFT + U 44 to take into account the electron-electron interaction in the d orbital of the magnetic S. Tiwari et al.
transition-metal atoms. We use the linear response method 45 to determine the Hubbard U parameter for the d orbitals of the dopant atoms. The U values we obtain from the linear response calculation are in the range 4 − 6 eV for Ti, V, Cr, Mn, Fe, Co, and Ni-doped TMDs. For the transitionmetal and the chalcogen atoms of the TMDs, we use a U = 0 for all their orbitals. We have verified our results by applying a U on the d-orbital of the base transition metal atoms for the TMDs, and our result does not change qualitatively and quantitative changes were small. For example, for Crdoped MoS 2 the near-neighbor anti-ferromagnetic interaction increased by a mere 4% when we applied a U = 4 eV on the Mo atoms.
To obtain Eq. (2) parameters, we place two dopant atoms in a supercell of the base TMD at various positions (near, next, and next-next neighbor) and calculate the total energy of various magnetic configurations. We use supercells of sizes 4 × 4 × 1, 5 × 5 × 1, and 7 × 7 × 1 to calculate the total energies of various magnetic configurations.

Monte-Carlo simulations
We simulate the phase change of the Heisenberg Hamiltonian using Monte-Carlo (MC) simulations. We obtain the critical temperature from the peak of the susceptibility obtained from the MC simulations. The critical temperature was calculated only for the samples with M sat ≥ 0.33 × M start , where M sat is the saturation magnetization per dopant atom obtained from MC, and M start is the starting magnetization per dopant atom obtained from DFT. To ensure that we capture the effect of configurational entropy, we run 20 separate MC calculations for each material combination and percent atomic substitution, each starting from randomly doped configurations of a 40 × 40 supercell of a TMD. We use a pseudorandom number generator to generate random positions in a pristine TMD lattice, and place the dopant atoms at those positions. We investigate TMDs with an atomic substitution ranging from 6% to 18%. For each randomly doped configuration, we use 1000 equilibration steps, and 1000 MC steps for averaging observables, at each temperature step. For each equilibration and MC step, we perform N atom spin-flip steps, where N atom is the number of dopant atoms in the unit cell.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.