Magnetic properties and critical behavior of magnetically intercalated WSe2: a theoretical study

Transition metal dichalcogenides, intercalated with transition metals, are studied for their potential applications as dilute magnetic semiconductors. We investigate the magnetic properties of WSe2 doped with third-row transition metals (Co, Cr, Fe, Mn, Ti and V). Using density functional theory in combination with Monte Carlo simulations, we obtain an estimate of the Curie or Néel temperature. We find that the magnetic ordering is highly dependent on the dopant type. While Ti and Cr-doped WSe2 have a ferromagnetic ground state, V, Mn, Fe and Co-doped WSe2 are antiferromagnetic in their ground state. For Fe doped WSe2, we find a high Curie-temperature of 327 K. In the case of V-doped WSe2, we find that there are two distinct magnetic phase transitions, originating from a frustrated in-plane antiferromagnetic exchange interaction and a ferromagnetic out-of-plane interaction. We calculate the formation energy and reveal that, in contrast to earlier reports, the formation energy is positive for the intercalated systems studied here. We also show that in the presence of W-vacancies, it becomes favorable for Ti, Fe, and Co to intercalate in WSe2.


Introduction
Transition metal dichalcogenides (TMDs) are layered materials that have a wide variety of interesting physical properties, and are being studied for applications in various technological fields [1]. TMDs exhibit a large range of electronic and magnetic properties. The presence of an electronic band gap, combined with the possibility for strong spin-orbit coupling, makes doped TMDs potential candidates for spintronic applications [2][3][4][5]. The fabrication of fieldeffect devices using TMDs [6][7][8][9], has encouraged research on TMD-based spintronic devices.
The intercalation of elements or even molecules between the layers of TMDs has been studied extensively and has yielded pathways to modify the properties of the material, using, for example, organic molecules [10]. It is possible to use intercalation of organic molecules to modulate the strength of charge density waves in for example 2H-NbSe 2 to increase the stability of the superconducting phase [11]. Other properties that can be tuned using intercalation are, for example, the electrochemical tuning of 2H MoS 2 using Li + ions [12,13]. Dopants can also affect the structure of the TMD itself, and can be used in phase engineering [14].
Recently, there have been reports of roomtemperature ferromagnetism in monolayer TMDs, for example in metallic monolayer VSe 2 [15] and in ultrathin VS 2 flakes [16]. Doping with magnetic atoms, for example third-row transition metals such as Fe or Mn, enhances the magnetic properties of the host TMD. Theoretical studies predict more TMD ferromagnets, such as doped WSe 2 [17], and doped MoS 2 [3,4]. Recent experimental reports discuss the effect of increased doping concentration on the strength and nature of the magnetic ordering in magnetically doped TMDs [18,19].
Antiferromagnetic materials are being studied in the context of spintronics as well. The absence of a net magnetic moment in antiferromagnets provides protection from external magnetic fields, and results in small parasitic fields, with little 'cross-talk' between antiferromagnetic memory elements [20,21]. Because of the insensitivity to external magnetic fields, antiferromagnetic memory elements need to be manipulated using electrical current, which must be able to switch the state of the material and perform a non-destructive read on the memory element. Recently, Nair et al demonstrated electrical switching in Fe-intercalated Fe 1/3 NbS 2 [22].
In this work, we investigate the magnetic properties of WSe 2 intercalated with Ti, V, Cr, Mn, Fe and Co. We calculate their Curie/Néel temperatures using density functional theory (DFT) and Monte Carlo simulations. To gauge the viability of these doped structures, we compute the formation energies of the TMDs intercalated with Ti, V, Cr, Mn, Fe and Co, which determines their thermodynamic stability.

Computational methods
We use a combination of Density Functional Theory (DFT) calculations and Monte Carlo simulations of the spin dynamics to obtain the magnetic transition temperatures of the materials under study. For all DFT calculations, we use the Vienna Ab-initio Simulation Package (VASP) [23][24][25][26]. We employ the projector-augmented wave (PAW) method, [27], in combination with the generalized gradient approximation as proposed by Perdew, Burke and Ernzerhof (PBE) [28] for the exchange and correlation functionals. We use the DFT-D3 method of Grimme et al [29], to account for the van der Waals interactions between the layers and the intercalants. We employ a plane-wave cut-off of 500 eV for all our DFT calculations, to ensure the results are accurate. We use a Gamma-centered 6 × 6 × 4 k-grid for the relaxation of the 2 × 2 × 1 supercell with intercalant atoms, and a 3 × 6 × 4 k-grid for the 4 × 2 × 1 supercell used in the magnetic calculations. We use an energy cutoff of 10 −5 eV for all DFT calculations. During the structural optimization of our systems, the lattice parameters and ion positions are changed until the forces ions are all lower than 0.005 eV Å −1 .

Hubbard U correction
To account for the correlation effects of electrons in the 3d-shell of the intercalants, we adopt the Hubbard U model within the DFT+U framework [30,31].
We estimate the value of U with the linear response method of Cococcioni and de Gironcoli [32]. We calculate the different linear response values of U for each different dopant/host combination. The obtained value of U is used in all subsequent spin-polarized DFT+U calculations for that system.

Energy calculations of the magnetic states
The magnetic behavior of a material depends on the relative stability of different ferromagnetic, ferrimagnetic and antiferromagnetic states. In order to get a good estimate of the Curie temperature for ferromagnets, or the Néel temperature for antiferromagnets, we take different magnetic states into consideration for each material. We take the ferromagnetic state, along with different ferrimagnetic and antiferromagnetic spin configurations. To find the different magnetic configurations, we first generate all permutations of spin up or down on each intercalant site. We obtain a list of all possible magnetic configurations. For numerical efficiency, we determine which magnetic structures are symmetrically equivalent, simulating them only once. We perform a collinear spin-polarized DFT+U calculation for each of the magnetic structures. In the subsequent Monte Carlo calculations, we take into account the multiplicity of the symmetry-equivalent states. We apply an initial magnetic moment to each intercalant, corresponding to the number of unpaired electrons in the dshell of the respective intercalant element. Due to the formation of bonds, the final spin states of the intercalant exhibit a reduced magnetic moment compared to their atomic configuration, as shown in table 1.

Exchange parameters
We model the magnetic phase change by calculating the energy of the system using a local Heisenberg Hamiltonian, where the indices i and j indicate the magnetic intercalant sites and S i is the magnetic moment of the intercalant at site i. Since our magnetic calculations are collinear, S i and S j are scalars. Using the method developed by Tiwari et al [33], we obtain the magnetic exchange parameters J ij from the DFT+U calculations. The Heisenberg Hamiltonian of equation (1), with the parametrized exchange interactions J ij , serves as the basis for the Monte Carlo simulations. For each intercalated WSe 2 material, we report the in-plane and out-of-plane nearest-neighbor exchange parameter.

Monte Carlo
For the Monte Carlo calculations, we employ a supercell with 9 × 9 × 8 intercalant sites. For each temperature, we perform 2000 equilibration steps and 2000 subsequent steps. For better statistics, we average observables (e.g., magnetization, total energy) over six independent Monte Carlo runs, each with a different initial state. The magnetic phase-transition temperature, i.e. Curie temperature T C for ferromagnets and Néel temperature T N for antiferromagnets, is defined as the peak of the specific heat, T C/N = argmax(c(T)). Following [34], the magnetic susceptibility and specific heat are calculated from the variance of the magnetization and energy, respectively: where k B is the Boltzmann constant, T is the temperature, N is the number of magnetic dopants in the supercell and ⟨x⟩ is the expectation value of the quantity x.
We compare the results from our more refined method to the Curie (Néel) temperatures for our (anti-) ferromagnetic materials found using the crude mean-field estimate [35] where ∆E is the energy difference between the ferromagnetic state and the most stable antiferromagnetic state in for each host-dopant system.

Formation energy
The formation energy per unit cell is calculated using where E system is the total energy of the intercalated system as obtained by DFT. E pure is the total energy of the pure TMD material, without any intercalants. E intercal is the total energy per atom of the pure intercalant, in its metallic form. M is the number of unit cell repetitions in the supercell, and N is the number of intercalant atoms in the supercell. Figure 1 shows the supercells we use in our calculations. We start from a 2 × 2 × 1 supercell of 2H WSe 2 with AB stacking, and we put two intercalant atoms per supercell, one between each layer, as illustrated in figure 1. After relaxation of the atomic positions, we find that the 2H AB stacking depicted in figures 1(a) and (b) is the most stable structure for Fe and Co intercalation. However, for Ti, V, Cr and Mn intercalation, we find that the distorted Bernal structure shown in figures 1(c) and (d), is more stable.

Structural optimization
For each intercalant, we use the most stable relaxed structure in the calculations that follow.
To assess the validity of our relaxed structures, we perform an additional relaxation step where we include spin polarization. In the spin-polarized relaxations, we apply initial magnetic moments to the dopants corresponding to the most stable magnetic configuration of each material. Comparing the lattice parameters with and without spin polarization, our results remain qualitatively the same and the quantitative change is very low (<1%). Only the lattice parameter c shows appreciable change for Mnintercalated WSe 2 (6.61%) and Fe-intercalated WSe 2 (1.95%), which are antiferromagnetic. These results indicate that structural changes due to magnetism can be neglected in our calculations. Section 3 of the supplementary information (available online at stacks.iop.org/2DM/8/025009/mmedia) discusses the spin-polarized relaxations in detail. Table 1 provides an overview of our simulation results. Monte Carlo simulations yield the transition temperatures between an ordered magnetic state and a paramagnetic state for each of our materials. For our Monte Carlo simulations, we use a Heisenberg Hamiltonian with, for each material, an in-plane (J ∥ ) and an out-of-plane (J ⊥ ) nearest-neighbor magnetic exchange parameter between magnetic dipoles located on the dopants. We obtain these parameters from Table 1. We report the calculated Hubbard U value, the calculated magnetic moment on the intercalant atoms, the exchange parameters used in the Monte Carlo simulations, and the transition temperatures of the magnetic states. J ⊥ is the out-of-plane nearest-neighbor exchange parameter, while J ∥ is the in-plane exchange parameter. For each host-dopant combination, we include the transition temperature as calculated using our Monte Carlo method and the mean-field approximation. We calculate the exchange interaction in the in-plane and the out-of-plane direction, using the method described in section 2.3. For Ti and Cr intercalation, we find ferromagnetic (J > 0) exchange interactions in both the in-plane and out-of-plane direction. We find indeed that DFT+U predicts ferromagnetic ground states for both Ti and Cr intercalation. In the case of Mn intercalation, both the in-plane and the out-of-plane exchange interactions are antiferromagnetic (J < 0). V, Fe and Co intercalants feature a ferromagnetic interaction out-of-plane (J ⊥ > 0) while the in-plane interaction is antiferromagnetic (J ∥ < 0). V-intercalated WSe 2 has in-and out-ofplane exchange interactions that are similar in magnitude but opposite in sign, which leads to complex phase transition behavior. Figure 2 shows the magnetic susceptibilities of each of the investigated systems. The magnetic susceptibility of each material shows a peak where the magnetization changes abruptly in the system. At temperatures below the peak, the system's magnetic moments are ordered, which means the system is in a ferromagnetic, antiferromagnetic or ferrimagnetic state. At temperatures above the peak, the magnetic moments are randomly ordered, and the system is in a paramagnetic state. Figure 3 shows the specific heat curves of each material. The specific heat of the material goes through a peak when the internal energy changes abruptly, marking a phase transition. In our case, the phase transition arises from the Heisenberg Hamiltonian and is therefore a transition from an ordered magnetic state to a paramagnetic state.

Magnetic ordering and transition temperature
We extract the transition temperatures from the specific heat peaks. For Ti and Cr intercalation, we find Curie temperatures of 34 and 225 K, respectively. For Mn, Fe and Co intercalation, we find Néel temperatures of 88, 327 and 48 K, respectively. In the case of V intercalation, the specific heat exhibits a curious shape: it has a strong peak at 148 K and a cusp at 65 K.
Fe and Co have the largest out-of-plane exchange parameters of all materials. At first glance, we would expect the large values of the out-of-plane exchange parameter to lead to high magnetic transition temperatures. For Fe intercalation, it does. However, for Co, the transition temperature is low despite the large out-of-plane exchange constant, because the low magnetic moment on the Co atoms bring down the transition temperature. Due to the significantly larger magnetic moment of the Fe atoms, together with the large out-of-plane exchange interaction, the Feintercalated WSe 2 has the highest transition temperature of all materials.
To verify the nature of the magnetic phases, figure 4 shows the magnetization against temperature to further discuss the nature of the phase transitions in the different intercalated materials. The magnetization is extracted at each temperature step and subsequently averaged over the different runs. We then normalize the magnetization curves to the saturation magnetization, which is the maximum achievable magnetization in the cell used in the Monte Carlo runs. The magnetization of Ti and Cr has a high value below T C , which drops, essentially to zero, at high temperatures. This is in line with the ferromagnetic ground state transitioning into a paramagnetic state at high T. Mn, as the only material which has antiferromagnetic interactions in the in-plane and outof-plane directions, shows a vanishing magnetization below T N , indicative of an anti ferromagnetic ground state.
For V, Fe and Co, all of which have a ferromagnetic in-plane and an antiferromagnetic out-of-plane interaction, we see evidence of frustration in the system, as the magnetization retains a nonzero value . The specific heat of magnetically intercalated WSe2 with respect to temperature. We normalize the specific heat curves, by dividing them by their respective maximum values. From these curves, we extract the transition temperature. below T N . Frustration appears in magnetic materials, when competing interactions prevent the system from reaching a ground state where all interactions are fulfilled. In our case, the triangular lattice formed by intercalants features antiferromagnetic interactions between the spin sites. However, each triangle can never have three antiferromagnetically ordered spins. In our study, we see that for V-, Co-and Fedoped WSe 2 , where the ferromagnetic interaction is dominant, the low temperature state is frustrated with a non-zero magnetization, transitioning to a paramagnetic state at high temperatures. The remaining magnetization at 0 K of V-, Co-and Fe-doped WSe 2 is due to the finite size of our simulation. In frustrated systems, the antiferromagnetic interactions cannot be completely satisfied, leading to an imbalance in the magnetic moments at short length-scales. In infinite systems, the magnetic moments cancel out, leading to a material with a net zero magnetization at 0 K.
In special cases, frustration leads to unexpected phenomena, such as the appearance of multiple phase transitions and partially ordered phases [36,37]. For V-doped WSe 2 the frustration yields a clear double-phase transition visible as a peak (at 148 K) and cusp (at 65 K) in the specific heat, shown in figure 3. This double phase transition is also apparent in as a complicated signature in the magnetization of V-doped WSe 2 . As discussed, below 65 K, the system is a frustrated antiferromagnet, in-plane, while being  ferromagnetically ordered out-of-plane, yielding a non-zero magnetization. For temperatures above 65 K but below 148 K, the system finds itself in partially ordered state, with an in-plane frustrated anti ferromagnetic state, while the out-of-plane ferromagnetic order is not yet fully established.

Formation energy and stability
In figure 5, we plot the formation energy, E form , of intercalated WSe 2 for every dopant element used.
Moving from Ti to Co along the 3d atomic series in the periodic table, E form peaks with Cr intercalation. Ti intercalated WSe 2 is the most stable material.
Our results disagree with previously reported results [38] that claim negative formation energies. However, upon closer inspection, Kumar et al. [38] calculated the formation energy by comparing the total energy of the combined system with the energy of pure WSe 2 and an atomic intercalant, i.e., a single metallic atom in vacuum. When using the atomic intercalant energy, the calculation yields the bonding energy, not the formation energy, of an intercalant in the lattice. However, thermodynamic stability is given by the formation energy, in reference to the stable metallic state of the intercalant. While the bonding energies in [38] are negative, mistakenly taken to indicate stability, our calculations reveal the intercalated structures are not stable, with all systems having positive formation energies.
To investigate a possible scenario of how intercalation could practically be realized, we investigate the formation energy in WSe 2 with a W vacancy. Details about this study are given in the supplementary information. We perform a structural relaxation for WSe 2 with W vacancy intercalated with Ti, V, Cr, Mn, Fe or Co, using the same method used in section 2. Figure 6 shows the calculated formation energies of intercalated WSe 2 with vacancy and lower doping concentration. Ti-intercalated WSe 2 now has a negative formation energy measuring −0.64 eV. The resulting formation energies for Co-and Feintercalated WSe 2 are 0.05 and 0.01 eV, respectively. Furthermore, intercalation increases entropy, strongly favoring intercalation at room temperature and above.

Conclusions
We have used density functional theory to investigate the magnetic properties of WSe 2 intercalated with Ti, V, Cr, Mn, Fe and Co. The magnetic properties of intercalated WSe 2 are highly dependent on the dopant type. We found that the spins align ferromagnetically for Ti and Cr intercalation with predicted Curie temperatures of 34 and 225 K, respectively. For Mn we found antiferromagnetic behavior that is predicted to persist up to 88 K. The high predicted Néel temperature for Fe-intercalated WSe 2 is appealing because of the potential applications for spintronics and other novel magnetic devices. For V, Fe and Co, the anti ferromagnetic interaction in-plane competes with the ferromagnetic interaction out-ofplane, resulting in a complex ground state that is a frustrated anti ferromagnet in-plane with ferromagnetic ordering out-of-plane. The frustration is due to anti ferromagnetic interactions in the triangular lattice formed by the intercalants in the plane. V-doped WSe 2 shows an interesting double phase transition, visible in the specific heat and the magnetization. At 65 K, the system undergoes a first phase transition to a partially disordered state, and at 148 K, the second phase transition occurs, and the magnetic ordering disappears, creating a paramagnetic state.
We have calculated the formation energies of our materials and found that, in contrast to [38], the formation energies are positive for all intercalants, indicating that the doped systems are not stable. Ti is the most stable intercalant, while Cr is the least stable. To address the question of stability, we have calculated the formation energies of our materials in the presence of W vacancies. We found that the formation energy of Fe-, Co-and Ti-intercalated WSe 2 reduces to 0.01, 0.05 and −0.64 eV, respectively. Taking the larger entropy of an intercalated structure into account, intercalation will take place at room temperature and above despite slightly positive formation energies for intercalation in pristine WSe 2 .
The magnetic properties of intercalated WSe 2 are favorable for future spintronic devices. We predict a Néel temperature for Fe-intercalated WSe 2 of 327 K, which offers the possibility for room-temperature applications. Additionally, we have shown that the unfavorable formation energies of our materials can be lowered by introducing W vacancies and lowering the doping concentration.
Finally, we would like to mention that addition of spin-orbit coupling and non-collinearity in modeling similar intercalated systems could be an interesting avenue for future exploration.