Limited Ferromagnetic Interactions in Monolayers of MPS3 (M = Mn and Ni).

We present a systematic study of the electronic and magnetic properties of two-dimensional ordered alloys, consisting of two representative hosts (MnPS3 and NiPS3) of transition metal phosphorus trichalcogenides doped with 3d elements. For both hosts, our DFT + U calculations are able to qualitatively reproduce the ratios and signs of all experimentally observed magnetic couplings. The relative strength of all antiferromagnetic exchange couplings, both in MnPS3 and in NiPS3, can successfully be explained using an effective direct exchange model: it reveals that the third-neighbor exchange dominates in NiPS3 due to the filling of the t2g subshell, whereas for MnPS3, the first-neighbor exchange prevails, owing to the presence of the t2g magnetism. On the other hand, the nearest neighbor ferromagnetic coupling in NiPS3 can only be explained using a more complex superexchange model and is (also) largely triggered by the absence of the t2g magnetism. For the doped systems, the DFT + U calculations revealed that magnetic impurities do not affect the magnetic ordering observed in the pure phases, and thus, in general in these systems, ferromagnetism may not be easily induced by such a kind of elemental doping. However, unlike for the hosts, the first and second (dopant-host) exchange couplings are of similar order of magnitude. This leads to frustration in the case of antiferromagnetic coupling and may be one of the reasons of the observed lower magnetic ordering temperature of the doped systems.

Introduction

Nonmagnetic van der Waals layered materials such as transition metal dichalcogenides have been extensively studied over the last several years.[1] Just recently, the intrinsic ferromagnetism in the true 2D limit has been reported,[2] initiating increasing excitement in the spintronics and 2D material communities. In particular, the long-range magnetic order has been observed for insulating CrI3,[3] semiconducting Cr2Ge2Te6,[4] and metallic Fe3GeTe2[5] compounds. In addition, topological spin structures have been predicted for 2D materials.[6]

Currently, the attention is on transition metal phosphorus trichalcogenides (MPX3, where M stands for a transition atom and X = S and Se), which could be easily exfoliated down to monolayers[7,8] and are semiconducting materials with a wide range of band gaps.[9] The MPX3 structures exhibit various antiferromagnetic (AFM) arrangements within the magnetic ions, which are theoretically expected to be measurable using different light polarizations.[10] Interestingly, the metal to insulator transition and superconductivity phase have been observed in this compound family.[11−13] In particular, the AFM insulator phase in bulk FePS3 can be melted and transformed into the superconducting phase under high pressure, providing similarity to the high-Tc cuprate phase diagram.[14] In addition, theoretical predictions point to the existence of a large binding energy of excitons in MnPS3, whereas the experimental reports have observed excitons in few layers of NiPS3 strongly related to magnetic order.[13,15,16] Recent reports have demonstrated an all-optical control of the magnetic anisotropy in NiPS3 by tuning the photon energy in resonance with an orbital transition between crystal field-split levels.[17] The aforementioned demonstrates that this family of compounds is an ideal platform to study correlation effects in the true 2D limit.

In contrast to ferromagnetic (FM) materials, the antiferromagnets exhibit limited applications, mostly in the terahertz regime as ultra-fast components or specialized embedded memory–logic devices.[18−20] Most of the current applications of magnetic crystals are based on the FM semiconductors, for which the band gap and FM order are crucial factors. The magnetic-phase transitions for the MPX3 materials can be accomplished by applying stress,[21] changing the carrier concentration or applying voltage.[22] In addition, the “M” atoms in MPX3 crystals might be substituted with other transition metal atoms inducing the FM order, as recently reported for a particular concentration of the non-magnetic dopants in the CrPSe3 host, resulting in a half-metallic FM state.[23] Moreover, a series of mixed systems in a bulk form[24−35] has been experimentally realized, thus entering a new playground of magnetic phases.

In addition, the magnetic and electronic properties of MPX3 materials and critical Néel temperature (from TN = 78 K for MnPS3[36−38] up to TN = 155 K for NiPS3[36,37]) strongly depend on the type of magnetic ion in the host. Thus, an elemental substitution could be an efficient way to tune the magnetic properties of atomically thin layers, by changing the lattice parameters and magnetic moments. These quantities result in manipulation of the exchange interactions that could be an effective way to engineer highly functional materials, similar to magnetic heterostructures.[39] In particular, the main reason behind the idea that adding ions with the partially filled d-shells into the system might lead to the enhancement of the FM interactions is related to the so-called double-exchange mechanism:[40] the hopping between two correlated ions with different valences (and thus comparable Coulomb interactions) and relatively strong Hund’s exchange is energetically favored provided that the spins on the neighboring ions are aligned in a parallel fashion. Such a mechanism might for instance be at play along the Cr–Mn bond: here the hopping from the manganese eg orbitals (eg2 configuration) to the chromium eg orbitals (eg1 configuration) lowers the total energy of the system provided that the spins on Cr and Mn are aligned ferromagnetically (for the antiparallel configuration, the strong Hund’s exchange makes such a hopping energetically unfavorable). Note that such a double-exchange mechanism is similar to the one observed in the doped manganites which are FM—albeit here it is a bit more complex, for it involves two different magnetic ions (nevertheless, the latter should not be that important, for Cr and Mn have comparable values of the Coulomb interactions, cf. Table II of ref (41)).

Altogether, in this work, we study the magnetic properties of the MPX3 monolayers and try to understand (i) how FM interactions can be stabilized in these compounds and (ii) whether one can easily modify these compounds by elemental doping so that the FM interactions can be strengthened. To this end, we examine two representative hosts MnPS3 and NiPS3 in the undoped case and at a particular doping concentration of magnetic ions. Using first-principles calculations, we give qualitative and quantitative explanations of the origin of the exchange coupling strengths up to the third nearest neighbors and their respective ratios for MnPS3 and NiPS3 structures. Considering the model Hamiltonians with ab initio parametrization, we discuss the competition between the direct exchange and the superexchange mechanisms for the host structures. Next, we study various substantial sites of dopant atoms with mixed spins and mixed nearest-neighbor magnetic interactions. Here, again using first-principles calculations, we examine in detail the mixed exchange coupling parameters between the metal host and dopant atoms.

In the next section, we describe the computational details. In the third section, we will present our Results and Discussion while the last section is devoted to the Conclusions.

Methods

DFT Computational Details

The first-principles calculations are performed in the framework of spin-polarized density functional theory (DFT) as is implemented in the VASP code.[42,43] The electron–ion interaction is modeled by using PAW pseudopotentials[44,45] with 3s and 3p states for P and S atoms and 3d and 4s states for Mn, Ni, and Cr being treated as valence states. The Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional is employed.[46] The kinetic energy cutoff for the plane-wave expansion of the wave functions is set to 400 eV. A k-mesh of 10 × 6 × 2 is taken to sample an irreducible first Brillouin zone of the rectangular planar cell (see Figure S1 in Supporting Information) containing 20 atoms including four transition metal atoms. The lattice parameters have been fully optimized within the PBE + U approach for the magnetic ground state of the monolayers, assuming the rectangular supercell. In particular, the magnetic ground state (Hubbard U parameter) for MnPS3 and NiPS3 is AFM-N (U = 5 eV) and AFM-z (U = 6 eV), respectively. In the case of (2 × 2) supercell, which consists of four primitive hexagonal unit cells, the 5 × 5 × 1 k-mesh is chosen to obtain the optimized position of the atoms. Considering density of states (DOS) calculations, the denser k-mesh equal to 10 × 10 × 2 k-points is taken into account. The convergence criteria for the energy and force are set to 10–5 eV and 10–3 eV/Å, respectively. In order to properly model a monolayer system, 20 Å of a vacuum is added to neglect the spurious interaction between the image cells. Note that the standard exchange–correlation functionals are insufficient to account for a non-local nature of dispersive forces, which are crucial for layered materials and adsorption molecules on the surfaces.[47−50] Thus, the semi-empirical Grimme method[51] with a D3 parametrization is applied.[52] For the 2D materials, we can use the HSE for a better description of the gap;[53] however, in the case of magnetism also, GGA + U provides the same effect. We employ the GGA + U formalism proposed by Dudarev[54] to properly account for on-site Coulomb repulsion between 3d electrons of transition metal ions, by using effective Hubbard U parameters.

Note that the proper choice of the U values is not straightforward due to the lack of accurate experimental information on electronic properties. Also, the common choice to compare the band gaps obtained in DFT with experiments to judge the U values also need caution due to the fact that one-particle Kohn–Sham DOS cannot be directly compared to the measured data. Thus, we decided to compute the Hubbard U using the linear response method proposed by Cococcioni[55] for the monolayers of MnPS3 and NiPS3, and we obtained 5.6 eV for the Ni and 5.3 eV for the Mn. Therefore, we have used UCr = 4 eV, UMn = 5 eV, and UNi = 6 eV similar to the linear response results, which are typical U values reported for MPX3 materials in previous reports.[21] Moreover, our Coulombic repulsions are close to the typical values of U in semiconductor compounds. Indeed, we find U = 6.4 for Ni2+ and U = 3.5 eV for Cr in oxides,[56] while the typical value of U for Mn2+ is U = 5 eV.[57,58]

To calculate the AFM direct exchange for MnPS3 and NiPS3, we extracted the real-space tight-binding Hamiltonian with atom-centered Wannier functions with d-like orbital projections on the transition metals using the Wannier90 code.[59,60] We calculated separately the hopping parameters for the orbitals, symmetric and antisymmetric with respect to the basal plane. The different symmetry and the separation in energy help to disentangle the two subsectors of the d manifold.[61,62] The calculation of the hopping parameters was carried out in the non-magnetic case to get rid of the magnetic effects and evaluate just the bare-band structure hopping parameters, and then, the hopping parameters will be used for the model Hamiltonian part.[63] In order to have parameters to use for the model Hamiltonian, we do not perform the maximum localization so as to have the Wannier function basis of our tight-biding model as close as possible to the atomic orbitals.

Results and Discussion

The results are presented as following: first, we present the results for the magnetic ground state of the hosts (pure MnPS3 and NiPS3 systems). In particular, we consider the exchange couplings within the DFT + U approach. Next, the AFM exchange mechanism is discussed within the minimal many-body model, and the AFM exchange coupling strengths are evaluated numerically using the Wannier basis with ab initio parametrization. Next, the qualitative explanation of the FM superexchange is presented. Finally, we present comprehensive studies of electronic and magnetic properties of benchmark alloys with a fixed concentration of dopants. The elemental substitution is employed at various atomic sites of the honeycomb lattice.

Undoped Hosts MnPS3 and NiPS3

Magnetic Couplings within the DFT + U Approach

The magnetic ground states of MnPS3 and NiPS3 exhibit AFM Néel (AFM-N)- and AFM zigzag (AFM-z)-type ordering (see Figure ), respectively. The neutron diffraction data predicted that the Mn2+ (3d5, S = 5/2) spins are slightly tilted (around 8°) from the perpendicular direction of the honeycomb lattice,[68] whereas the Ni2+ (3d8, S = 1) spins are aligned within the honeycomb plane. In addition, due to the different filling of the 3d orbitals for various metals ions (Fe, Mn, Ni etc.), the size of the magnetic exchange coupling (J) also changes. The existence of magnetic ordering at finite temperature in 2D limit requires the magnetic anisotropy in accordance to the Mermin–Wagner theorem. Currently, only FePS3, which possesses a strong out-of-plane easy axis, has been experimentally reported to exhibit AFM order in the monolayer up to 118 K.[8] In order to explain the anisotropic order, the dipolar interactions between the magnetic moments or a single-ion anisotropy resulting from a nonzero spin–orbit coupling should be considered. The recent experimental reports, such as magnon band measurements,[69] support the claim that the dipolar interactions should be the leading term, whereas the electron spin resonance[70] and critical behavior measurements[71] indicate that the single-ion anisotropy might come into play. In addition, experimental observations demonstrated that the AFM ordering persists down to bilayer samples and is suppressed in the monolayer.[72] Notably, the recent theoretical report[73] has questioned the Raman criterion used for the monolayer studies therein, suggesting that NiPS3 magnetic ordering could be presented in monolayer samples, as also indicated by strong two-magnon continuum existing in thin samples of NiPS3.[72] In the case of MnPS3, the magnetic order has been presented down to the bilayer and was reported to be absent in the monolayer.[74] The suppression of the Néel temperature in thin samples can be associated with reducing the interplanar coupling in atomically thin samples.[3,72] Aforementioned results demonstrate that the magnetic ordering in monolayers of MPX3 is still under a hot debate and many experiments are being carried out to verify the theoretical predictions. Similar effects of strong spin fluctuation and absence of interlayer exchange coupling that weaken the long-range spin order in the 2D limit have been reported in other layered magnets such as Cr2Ge2Te6 and CrI3.[3]

Figure 1

Spin arrangements of the metal atoms in (a) AFM Néel (AFM-N), (b) AFM stripy (AFM-s), and (c) AFM zigzag (AFM-z) states. (d) Top and (e) side views of the crystal structure of the MPX3 system with the metal atoms denoted in blue and surrounded by the sulfur atoms (yellow balls). The gray, green, and violet lines indicate the NN, 2NN, and 3NN distance between the metal atoms, respectively. The dotted and solid lines indicate that the sulfur atoms are below and above the metal layer, respectively.

Here, we focus on the rationally invariant Heisenberg contribution, whereas the dipolar and single-ion contributions are out of the scope of the present work. The latter has been estimated to 0.3 and 0.009 meV for NiPS3 and MnPS3, respectively, and discussed theoretically in ref (66). Here, we present the results of exchange couplings up to the third nearest neighbor (for clarification, see Figure S1 in Supporting Information).

The exchange interaction up to the third nearest neighbor (JM) between the metal atoms of the hosts has been widely studied in a series of previous studies[21,66,67,75−77] (see Table ). Note that the prediction of the magnetic ground state within the DFT calculations does not depend on Hubbard interaction, whereas it is well known that the exchange coupling strength is sensitive to both U and the lattice parameters. We set the Hubbard U parameter to U = 5 eV and U = 6 eV for the 3d orbitals of Mn and Ni atoms, respectively. These values are calculated from first principles using the Cococcioni approach.[55] Our predicted JM values are in good agreement with neutron diffraction experiments (see Table ). The dominant exchange coupling is J3Ni, which is much stronger than J1Ni. Note, that the experimental ratio of critical temperatures is reflected in the relation of dominant contributions J3Ni > J1Mn. In both cases, J2 is much smaller than the other two exchange couplings. In particular, J2 and J3 couplings might be considered as superexchange interactions involving the atoms in the path M–S1···S4–M for J2 (see Figure d,e), where the S atoms are located in different sublayers, whereas the J3 interaction is mediated by S atoms located in the same sublayer through the bridge M–S1···S2–M. One could expect stronger hybridization of the S p states and M 3d states within the same sublayer of S atoms. In addition, the calculations reveal that J1 is AFM and FM for MnPS3 and NiPS3, respectively. Note that for both MnPS3 and NiPS3, the angle between the M–S–M atoms is close to 90° (83.4° for MnPS3 and 85.4° for NiPS3) for pointing to FM superexchange according to Goodenough–Kanamori–Anderson rules.[78,79] The direct M–M and indirect M–S–M (superexchange) mechanisms are of crucial importance to understand the differences between these two systems.

Table 1

Exchange Coupling Strengths JM of MnPS3 and NiPS3 Systems Using Model Calculations and Various Lattice Parameters and On-Site Coulomb Repulsion U in the DFT + U Calculationsa

structure	J1M	J2M	J3M	|J2M/J1M|	|J3M/J1M|
MnPS3 (model calculations, see caption above)	19.5	0.35	7.76	0.02	0.40
MnPS3 (PBE, a = 6.00 Å, AFM-N)[66]	1.21	0.18	0.54	0.15	0.45
MnPS3 (U = 5 eV, a = 6.11 Å, AFM-N) this work	1.04	0.05	0.53	0.05	0.51
MnPS3 (U = 5 eV, a = 5.88 Å, AFM-N)[67]	1.58	0.08	0.46	0.05	0.29
MnPS3 (U = 4 eV, a = 6.00 Å, AFM-N)[21]	0.4	0.03	0.15	0.08	0.38
MnPS3 (U = 3 eV, a = 6.00 Å, AFM-N)[66]	1.42	0.08	0.52	0.06	0.37
MnPS3 (experiment, AFM-N)[64]	1.54	0.14	0.36	0.09	0.23
NiPS3 (model calculations, see caption above)	–4.9	0.06	14.8	0.01	3.0
NiPS3 (PBE, a = 5.78 Å, AFM-z)[21]	–11.3	–0.12	36	0.01	3.2
NiPS3 (U = 6 eV, a = 5.84 Å, AFM-z) this work	–3.34	–0.19	13.7	0.06	4.1
NiPS3 (U = 4 eV, a = 5.78 Å, AFM-z)[21]	–4.11	1.95	17.4	0.47	4.2
NiPS3 (U = 3 eV, a = 5.78 Å, AFM-z)[66]	–2.6	–0.32	14	0.12	5.4
NiPS3 (experiment, AFM-z)[65]	–3.8	0.2	13.8	0.05	3.6

Positive (negative) JM indicate AFM (FM) correlations, respectively. Note, that in refs (64) and (65), different conventions of the exchange couplings JM were used. The AFM-N and AFM-z indicate the magnetic ground state of the system. The exchange couplings in the model calculations are obtained within the direct exchange mechanism [for MnPS3 using eqs and 4 whereas for NiPS3 using eq ]—except for J1Ni for NiPS3, which contains also an important superexchange contribution and follows from eq ; see text for further details. More details on the calculation of the magnetic exchanges within the DFT + U approach are present in Supporting Information. All values are given in meV.

Effective Direct Exchange and AFM Couplings

In order to understand the origin of the exchange couplings and their relative strengths, we consider model Hamiltonians for the direct and superexchange interactions. Note, that the direct exchange discussed in this work is a (second order) kinetic exchange process, which involves only the hopping between the transition metal ions, with the ligands not explicitly involved, cf. ref (80), whereas a superexchange term is a fourth- (or higher-) order kinetic exchange process, which explicitly involves the hopping over the ligands, cf. ref (80). In addition, we are pointing out below the mechanisms which could impact the sign of J1.

To understand the origin of the magnetic couplings, we first write down a minimal many-body model which solely contains the valence electrons of the transition metal ion[81]—the multi-band Hubbard model. From this, using the second-order perturbation theory that is valid in the Mott insulating limit, we derive the (effective) direct exchange processes. By construction, all obtained spin couplings have to be AFM. Therefore, while surprisingly successful, the following simple analysis will not be able to explain the onset of the FM couplings in NiPS3 (more on this at the end of the section).

In the monolayers of MnPS3 and NiPS3, the metal atoms are surrounded by six sulfur atoms (MS6 octahedron) and exhibit D3 point group symmetry in trigonal anti-prismatic environment of the ligands (sulfur atoms), see Figure a. Hence, due to this trigonal crystal field effect, the d manifold splits in two disentangled subsets of bands. The only coupling between these subsets is the spin–orbit coupling.[82] The bands lower in energy are even (d, d, and d) with respect to the basal plane, while the bands higher in energy are odd (d and d) with respect to the basal plane. The Mn ion is d5 and Mn d bands split into half-filled even and half-filled odd bands; therefore, the magnetic coupling acquires contributions both from the even and the odd orbitals. Instead, because the Ni ion has a d8 configuration, the even orbitals are fully occupied, and therefore, there is no magnetic contribution from these orbitals. Altogether, the minimal model is the five-band Hubbard model with a simplified structure of the Coulomb interactions (no spin on-site spin exchange and pair-hopping terms, cf. ref (83))In this model, c† creates an electron with spin σ = ↑, ↓ in a Wannier orbital |m⟩ = |x2 – y2⟩, |xy⟩, |xz⟩, |yz⟩, or |3z2 – r2⟩ at site i, and n = c†c. ↑ (↓) indicates the spin up (down). The parameter t is the hopping integral from orbital m at site i to orbital m′ at site i′. The on-site terms t = ε give the crystal field splitting. U and JH are the direct and (Hund) exchange terms of the screened on-site Coulomb interaction.

Figure 2

(a) Trigonal anti-prismatic and (b) octahedral geometries. On the right side, the d-orbital splitting of the metal atom in the corresponding crystal fields is shown. M and S denote the metal atom and sulfur atom, respectively.

We applied second-order perturbation theory in t/U and for the commensurate electron filling of Mn and Ni ions. The direct exchange constant for the valence electrons occupying the Mn odd orbitals iswhere a and b are the dxz and d orbitals, respectively, while i and i′ are the Mn lattice sites. In this formula, we take in consideration both the cases with i = 1 ≠ i′ = 2 and i = i′ = 1. UMn and JHMn are the Coulomb repulsion and Hund coupling, respectively, in the case of the Mn atoms. The Hund’s rule interaction between odd and even electrons yields a magnetic coupling between these electrons; therefore, the denominator depends on the occupancy of the even orbitals. By symmetry, the on-site energies are ϵMn1 = ϵMn2 = ϵMn1 = ϵMn2.

Similarly, we obtain the direct exchange constant for the valence electrons occupying the odd orbitals of the Ni atomswhere a and b are the d and d orbitals, respectively, while j, j′ are the Ni lattice sites. In this formula, we take in consideration both the cases with j = 1 ≠ j′ = 2 and j = j′ = 1. UNi and JHNi are the Coulomb repulsion and Hund coupling, respectively, in the case of the Ni atoms. By symmetry the on-site energies are ϵNi1 = ϵNi2 = ϵNi1 = ϵNi2.

Finally, the direct exchange constant for the even Mn orbitals is as followswhere c, d and e are the orbitals d, d, d respectively, while i and i′ are the metal lattice sites. In this formula, we take in consideration both the cases with i = 1 ≠ i′ = 2 and i = i′ = 1. UMn and JHMn are the Coulomb repulsion and Hund coupling, respectively, in the case of the Mn atoms. By symmetry the on-site energies are ϵMn1 = ϵMn2 and ϵMn1 = ϵMn2 = ϵMn1 = ϵMn2.

The obtained (total) direct exchange, J = Jodd + Jeven, is positive for both Mn and Ni ion and for any distance j – j′. Hence, as already mentioned, it is always AFM by construction.

Numerical Evaluation of the AFM Couplings

The Wannier basis provides us with ab initio values of the hopping integrals and crystal field splitting. We calculate the hopping parameters and the on-site energies using the interpolated band structure of the Wannier functions of the d-subsector. The on-site energies for Mn and Ni are: ϵMn1 = −1197.7 meV, ϵMn1 = −2082.6 meV, ϵMn1 = −2179.3 meV, and ϵNi1 = −1768.5 meV. We have three first nearest neighbors (1NN), six second nearest neighbors (2NN), and three third nearest neighbors (3NN). In the case of odd orbitals, the 3NN couplings are greater than the 1NN couplings, even by an order of magnitude; therefore, it is very important to consider these hopping amplitudes in the calculations. On the other hand, for even orbitals, the 1NN couplings are greater than the 3NN couplings. The 2NN couplings are always smaller with respect to the 1NN couplings and 3NN couplings. In the case of even orbitals, we neglect the difference between the on-site energies assuming the following approximation ϵMn – ϵMn, ϵMn – ϵMn ≪ (UMn + 4JHMn).

Now, we will numerically evaluate the second- and third-neighbor direct exchange as a function of the first-neighbor direct exchange. For the odd subsector of MnPS3, we obtain J2Mn,odd = 0.037J1Mn,odd and J3Mn,odd = 2.026J1Mn,odd. For the even subsector of the MnPS3, we obtain J2Mn,even = 0.015J1Mn,even and J3Mn,even = 0.016J1Mn,even. Considering that JNi,even = 0 due to fully occupied even orbitals, for the odd subsector of NiPS3, we obtain J2Ni,odd = 0.047J1Ni,odd and J3Ni,odd = 11.09J1Ni,odd. If we consider the total direct exchange value as the sum of the odd and even sector, we have J1Mn = (152.4 eV·meV)/(UMn + 4JHMn) for the first-neighbor coupling in MnPS3. For the second and third neighbors, we obtain J2Mn = 0.018J1Mn and J3Mn = 0.397J1Mn, and therefore, for MnPS3, the dominant direct exchange comes from the first-neighbor coupling. Using a Coulomb repulsion of 5 eV for Mn and a Hund coupling of 0.7 eV, we obtain the numerical values equal to J1Mn = 19.5 meV, J2Mn = 0.35 meV, and J3Mn = 7.76 meV, (see Table ). Note that the values are overestimated in comparison to experimental values (see Table ); however, the signs and the dominant contributions are the same.

When we numerically evaluate the direct exchange of NiPS3, we obtain J1Ni = (8.957 eV·meV)/(UNi + JHNi) for the first-neighbor coupling, J2Ni = 0.047J1Ni and J3Ni = 11.09J1Ni for the second and third neighbors, respectively. Remarkably, the 3NN magnetic direct exchange is larger than the 1NN exchange in the odd case. Using a Coulomb repulsion of 6 eV for Ni and a Hund coupling of 0.7 eV, we obtain the numerical values equal to J1Ni = 1.3 meV, J2Ni = 0.06 meV, and J3Ni = 14.8 meV, where last two are reported in Table . Note, that J2Ni and J3Ni are in good agreement with experimental values (see Table ), while J1Ni is significantly different and has a wrong sign. Therefore, in the case of J1Ni, a more complex superexchange model will be presented further in the text.

Altogether, we obtain that the simple direct exchange scheme gives J3Ni ≫ J1Ni ≫ J2Ni and J1Mn > J3Mn ≫ J2Mn. In the Ni case, the leading term is J3Ni, while in the Mn case, the leading term is J1Mn. The reason for this different behavior comes from the different filling that produces JNi,even = 0. The calculated direct exchange couplings are qualitatively in agreement with the magnetic couplings obtained experimentally or using DFT—except that J1Ni and J2Ni are FM due to more complex magnetic exchange not taken into account by the (simple) direct exchange scheme (see below). Even though the latter coupling is relatively small, note that considering in more detail, such discrepancy is relevant for an accurate description of the magnetic coupling and, hence, of the magnetic critical temperature. The different leading exchange terms in Mn and Ni compounds open the way to manipulate the magnetism by tuning the concentration of Mn and Ni compounds or by adding new magnetic materials as dopants.

Superexchange and FM Coupling in NiPS3

So far, we considered an effective direct exchange model, which solely contained the exchange processes due to the electrons hopping between the transition metal (Mn or Ni) ions. Note that such a model should be considered an effective one, for in reality, the hopping between the neighboring Mn or Ni ions is predominantly mediated by the sulfur ions. Hence, a natural extension of the direct exchange model should explicitly contain the exchange processes on the sulfur ions—in fact, it is such a nearest neighbor superexchange model[84] that is studied below to explain the onset of the nearest neighbor FM exchange in NiPS3. Note that for consistency, we comment at the end of this subsubsection why the more complex superexchange model is not needed to understand the other, the first- and third-neighbor,a magnetic couplings in MPS3—that is, the effective direct exchange model is enough.

We introduce the nearest neighbor superexchange model[84] for NiPS3 by considering two Ni ions connected via two sulfur atoms over two 90° bonds, see Figure . Note that, because Ni2+ is in a d8 configuration and sulfur ions are fully occupied, it is easier to consider the hole language as discussed below. Moreover, in what follows, we neglect the small trigonal distortions and we assume an octahedral crystal field with the division in t2g and eg orbitals [with the coordinate system defined in such a way that the xy plane coincides with the plane formed by the nearest-neighbor transition metal ion and sulfur, see Figures b and 3]. Hence, we begin by considering a fully atomic limit without hopping (zeroth order of perturbation theory in the small kinetic energy) in which there are two holes localized in two distinct Ni eg orbitals (d, d3; in what follows assumed to be energetically degenerate, see also above) and two (lying higher by the charge transfer energy Δ) empty p orbitals (p, py) on sulfur. As before, due to the strong Hund’s rule JH, the two Ni2+ holes form a high spin S = 1 state. Now, let us perform a perturbation theory in the kinetic energy (over Coulomb repulsion U and charge transfer energy Δ) and consider the possible exchange processes—which are of two kinds:

Figure 3

Schematic comparison between the AFM direct exchange (top panels) and the FM superexchange (bottom panels) processes in NiPS3. Finite direct-AFM exchange processes due to (a1) nonzero hopping elements ∝ ddπ between the nearest-neighbor d orbitals on nickel and (a2) nonzero hopping elements ∝ ddσ between the nearest-neighbor d3 orbitals on nickel. Lack of direct-AFM exchange processes due to vanishing hopping elements between the nearest-neighbor d and d3 on nickel. Finite FM superexchange processes due to (b1) nonzero hopping elements ∝ pdσ between the nearest-neighbor d orbitals on nickel and the nearest-neighbor p/p orbital on sulfur, (b2) nonzero hopping elements ∝ pdσ between the d3 orbitals on nickel and the nearest-neighbor p/p orbital on sulfur, and (b3) nonzero hopping elements ∝ pdσ between the d orbital on nickel and the nearest-neighbor p/p orbital on sulfur and the d3 orbital on nickel and the nearest-neighbor p/p orbital on sulfur. See text for more details.

First, there are direct exchange processes between the nickel ions, see Figure a. By definition, these concern virtual occupancies of one of the nickel ions by three holes [with a relative energy cost of U + JH according to the simplified structure of the Coulomb interactions, see eq ] and are possible once the hole can directly hop back and forth between the nickel orbitals under consideration. According to the Slater–Koster scheme,[85] which is qualitatively confirmed by our DFT calculations, the latter is allowed between the pair of d orbital (ddπ hopping element) and between the pair of d3 orbitals (ddσ/4 hopping element; the small ddδ element can be neglected), cf.Figure (a1,a2); note that the hopping between the d and the d3 orbital vanishes in this geometry, cf.Figure (a3). Altogether we obtain the nearest-neighbor direct exchange contributionwhere we assumed a typical relation between the Slater–Koster hopping integrals ddσ ≈ 2ddπ. It is important to state here that the above direct exchange process is different than the effective direct exchanges process defined in the previous subsections—for the latter ones may include all indirect hoppings (e.g., via sulfur) between the two nearest neighbor nickel ions.

Second, there are superexchange processes between the nickel ions, see Figure b. By definition these concern virtual occupancies of one of the sulfur ions by two holes (with an energy cost associated with the charge transfer 2Δ for antiparallel spins or 2Δ – JH for parallel spinsb) and are possible once the two holes hop back and forth between the sulfur and nickel orbitals under consideration. Again, using the Slater–Koster scheme,[85] and considering the three distinct possible hopping processes shown in Figure (b1–b3), we obtain the nearest neighbor superexchange contributionwhich is negative (FM) due to the lowering of the energy for a FM (virtual) occupancy of sulfur by two holes in the superexchange process. Note that for a given Ni–Ni bond, there are always two superexchange processes: one over the top-left and one over the bottom-right sulfur (hence the factor of two within the square brackets in the above formula). Moreover, the process depicted in Figure (b3) has to be multiplied by a factor of two, for one can interchange the position of the d and the d3 orbital and in this way double the amplitude of this process. Finally note that, due to the fact that we have spins S = 1 on nickel, overall the above superexchange process is reduced by a factor 1/2 with respect to to an analogous one for the S = 1/2 on copper (the superexchange processes have to be projected on the high-spin S = 1 states on both nickel ions—hence a factor of reduction).

Let us now comment why the above FM superexchange mechanism is not important for the nearest neighbor exchange in MnPS3. The reason for this is that in the case of manganese ions, we are a bit closer to the situation discussed in ref (84), which, for instance, shows the AFM exchange coupling in the case of half-filled t2g subshells of Cr3+ in LiCrS2. More precisely, the situation for the Mn2+ ions is as follows. On one hand, the AFM direct exchange is much stronger for Mn2+ because one of the t2g electrons (the d) can hop over the ddσ bond. On the other hand, the superexchange also contains an additional strong[84] AFM contribution due to the superexchange processes over one p sulfur orbital—which strongly hybridizes with the two nearest t2g orbitals. Altogether, as confirmed by the effective direct exchange studies in the previous subsection, these two mechanisms originating in the t2g exchange processes easily overcome the (above-described) FM processes for the eg orbitals.

Finally, we mention that the (surprisingly) strong third-neighbor coupling in MPS3 can easily be explained using the superexchange model. First, for the t2g sector, the strong AFM third-neighbor coupling is already discussed in ref (84)—see Figure 7 of the reference. Second, one can easily imagine that a similarly strong AFM coupling can also be realized for the eg electrons: in order to understand it, one just needs to replace the two third-neighbor d orbitals with the d orbitals and rotate the sulfur p orbitals by 90° in the process shown in Figure 7 of ref (84). In fact, the latter one should have a very high amplitude and hence a really strong third-neighbor exchange in NiPS3.

Numerical Evaluation of the FM Coupling in NiPS3

Having derived the direct exchange and superexchange processes between the nearest neighbor nickel ions, we are now ready to estimate the contributions of these both (competing) spin interaction terms. To this end, we assume that Coulomb repulsion U = 6 eV (the chosen value of our DFT + U approach, see Table ) and take a typical (for 3d transition metal compounds) value of Hund’s exchange JH = 0.7 eV. Next, based on the DFT calculations, we estimate that (i) the charge transfer energy Δ = 3 eV, (ii) the hopping pdσ = 0.9 eV, and (iii) the hopping ddσ = 0.05 eV. From this and using eqs and 6, we can easily calculate the AFM and FM contributions to the spin exchangeAs the second exchange is larger (by absolute value), we conclude that it is a relatively strong superexchange along the 90° nickel-sulfur-nickel bonds, which triggers the FM exchange along the nearest neighbor nickel ions. The sum of the aforementioned contributions leads to FM exchange coupling (J1Ni = −4.9), which is in good agreement with experimental studies (see Table ).

At this point, we would like to comment on the fact that for the edge-sharing copper chains, typically one entirely neglects the direct exchange cf. ref (86). The reason for this is that, in the copper chains, the direct copper hopping ddσ = 0.08 eV,[87] whereas usually one assumes that for the cuprates, pdσ ≈ 1.5 eV. With such hopping parameters (and assuming the typical cuprate values of U = 8 eV, Δ = 3 eV, and oxygen JH = 0.7 eV), one can immediately see that for the copper chains with 90° geometry, JCu,direct ≈ 3 meV ≪ |JCu,SE| ≈ 50 meV—which justifies why the direct copper exchange is typically neglected in such studies. In reality, the FM exchange is actually really small in such copper chains due to the angle along the copper–oxygen–copper bond not being strictly equal to 90°.

A strong suppression of the direct exchange is also observed in CrI3,[88] which is an experimentally confirmed 2D ferromagnet. CrI3 is isostructurally identical to MPS3 compounds when P atoms are removed and has partially-filled t2g shells, just as MnPS3. This might suggest that CrI3 should predominantly have AFM interactions, just as MnPS3. However, this is not the case because the relative strength of the direct and superexchange mechanisms is different in CrI3 and MnPS3. In CrI3, the Cr–I–Cr bond angle is very close to 90°, leading to a strong FM superexchange, while the Cr–Cr distance is relatively large (3.95 Å), giving rise to a very weak AFM direct exchange. Interestingly, due to the larger distances between the magnetic atoms in CrI3, the third-neighbor exchange is also less relevant.

Doped Systems (M3/4,X1/4)PS3

Our strategy is to try to induce long-range FM order via chemical doping, keeping the two-dimensional structure of the mother compounds intact. In principle, the chemical doping could influence the magnetic exchanges, especially the long-range ones. To this end, we consider doping of the host systems with distinct possible 3d elements. In particular, at low doping concentration, the main exchange couplings are the magnetic exchanges between the host magnetic atoms (JM for i = 1, 2, 3) and between the host impurity atoms (JMX for i = 1, 2, 3). Although for most of the 3d doping, the Mn–Mn and the Ni–Ni magnetic couplings remain AFM, in the case of Cr doping, the AFM exchanges between the magnetic atoms of the host turn into FM ones.[89] Thus, we performed DFT + U calculations for 25% concentration of the magnetic dopants (X = Cr, Mn, and Ni) in the hosts of NiPS3 and MnPS3. We have employed various structural arrangements of the atoms (see Figure ) and collinear spin configurations such as AFM Neel (AFM-N), zigzag (AFM-z), stripy (AFM-s), and FM case. For each of the configurations, the atomic positions have been fully optimized, keeping the lattice constants equal to the pure optimized monolayer structures (see Methods).

Figure 4

Crystal structure of 2 × 2 supercells consisting of fourfold primitive unit cells. The M and X are the transition metals of the host and dopant atoms, respectively. 1NN, 2NN, and 3NN denote the nearest neighbors, second neighbors, and third neighbors of the dopants, respectively.

Let us first discuss the impact of the dopant on the energetic and structural properties of the hosts. The energy difference between a particular configuration and the magnetic ground state for each considered alloy is presented in Figure . Note that each of the employed alloys preserves the magnetic ground state of the host, independent of the structural arrangement of the impurity and the type of the dopants (see Figure ). In addition, the favorite Mn and Ni dopant position is at 1NN (see Figure c,d), whereas the Cr dopants prefer to lie further apart, in particular, at 2NN (see Figure a) and 3NN (see Figure b) structural positions for NiPS3 and MnPS3 hosts, respectively. This reveals that the Mn and Ni dopants have tendency to cluster, while the Cr ions prefer to spread over the host. From now on, we only discuss the magnetic ground-state configurations of the alloys (the one which exhibits the lowest energy in Figure for particular alloy). The magnetization and the band gaps of these systems are collected in Table . All considered alloys exhibit a semiconducting behavior (see Table ). However, only (Ni3/4,Cr1/4)PS3 and (Ni3/4,Mn1/4)PS3 alloys have a nonzero net magnetization. Thus, we focus our discussion on these two ferrimagnetic systems. Owing to the different spins of the host and dopant atoms and particular arrangement of the dopants in the NiPS3 host, the ferrimagnetic state appears in these systems. The band structure and the orbital projections of these two ferrimagnetic alloys are presented in Figure . Note that in the (Ni3/4,Cr1/4)PS3 case, bands close to the Fermi level are mainly composed of 3d states of Cr dopants, causing a sizeable reduction in the energy gap compared to pure NiPS3 (2.3 eV for U = 6 eV). In the case of Mn impurities, the valence band maximum is mainly composed of Mn 3d states, whereas the conduction band minimum consists of very flat bands of Ni atoms (see Figure b).

Figure 5

(a–d) Energy difference (ΔE) between a particular magnetic configuration (gray rhombus FM, green square AFM-N, blue triangle AFM-s, and red circle AFM-z) and the magnetic ground state for each of the alloys (M3/4,X1/4)PS3. Note, that in each of the plots, the most energetically preferable arrangement of the dopants is at the lowest energy. The energy is given in eV per magnetic atom for various structural arrangements of the dopant in the host (1NN, 2NN, and 3NN).

Figure 6

Band structure and projected density of states (PDOS) for 2D ferrimagnetic alloys: (a) (Ni3/4,Cr1/4)PS3 and (b) (Ni3/4,Mn1/4)PS3. On the left side of (a,b), the corresponding atomic structural arrangements and local spin magnetic configurations are presented. Note that the spin of the dopant and the host atoms has different magnitudes, which are not visible here (spin arrows are not in scale). For better visibility, the spin channels and projections of the p states of S and P atoms are presented only in PDOS. The energy is given in eV, and the Fermi level is placed at zero.

Table 2

Results for the Most Energetically Favorable Position of the Dopants in the (M3/4,X1/4) Alloysa

alloy (M3/4,X1/4)PS3	magn. state (2 × 2) sc	M [μB] per sc	band gap [eV]
(Mn3/4,Cr1/4)PS3	3NN AFM-N	0	1.7 (D, Γ)
(Mn3/4,Ni1/4)PS3	1NN AFM-N	0	2.16 (D at K)
(Ni3/4,Cr1/4)PS3	2NN AFM-z	3.9	1.1 (↑, D) 2.11 (↓, ID)
(Ni3/4,Mn1/4)PS3	1NN AFM-z	5.7	2.11 (for both spins)

M is the total magnetization of the cell.

Next, we examine the mixed exchange couplings between the metal atom of the host and the impurity JMX within the classical Ising Hamiltonian in the honeycomb lattice for the smallest possible cell containing 25% dopant concentration (for details see Supporting Information and Figure S1 therein). Note that the FM exchange couplings are obtained for the Mn and Ni nearest neighbors and Cr–Cr ions at the second-nearest-neighbor distance (see J2X in Table ). The similar results for Cr have been recently reported in ref (90). In addition, the different values for J1MnNi (−1 meV) and J1NiMn (−0.7 meV) stems from the different lattice parameters of the hosts (see Table ). Although the mixed exchange coupling obtained for (Ni3/4,Mn1/4)PS3 correctly reflected the 1NN AFM-z ground-state configuration for this system, the difference between the 1NN AFM-N and 1NN AFM-z is just 1.5 meV per magnetic ion (see Figure c), thus resulting in magnetic frustration between the Mn and Ni atoms at the honeycomb sites due to the competition between the Néel and zigzag configurations, similar to the results reported in ref (35). Moreover, the critical temperature is directly related to the strength of the exchange couplings. Generally, the mixed exchange couplings are smaller than those in the corresponding metal atoms of the hosts. Thus, the critical temperature of the mixed structure is expected to be smaller than that for the corresponding pure system, which is in line with recent experimental reports on the Ni1–MnPS3 alloy[35] and with series of the mixed-system studies, where the suppression of TN temperature with dopant substitution is reported.[27,29,34] In addition, one should expect the further reduction of TN temperature of the employed mixed systems due to the possible disorder of the dopant atoms in the host, as reported in refs (30) and (34), which is not accounted in our studies.

Table 3

Exchange Coupling Strengths J Calculated for Various Alloy Systems, Implied by the Ising Modela

(M0.75,X0.25)PS3	J1XM	J2XM	J2X	J3XM
(Mn3/4,Cr1/4)PS3	2.5	0.1	–1.1	–0.4
(Mn3/4,Ni1/4)PS3	–0.7	–0.5	–0.6	2.0
(Ni3/4,Mn1/4)PS3	–1.0	0.3	0.15	4.9
(Ni3/4,Cr1/4)PS3	2.2	0.4	–0.3	2.3

Positive and negative J indicate AFM and FM ordering, respectively.

Conclusions

First, we examined the magnetic and electronic properties of the AFM-ordered systems MnPS3 and NiPS3 (without magnetic impurities). We presented a qualitative explanation for the relative ratio of the different nearest neighbor (first, second, and third neighbor) exchange couplings by studying an effective direct exchange for both MnPS3 and NiPS3. In particular, we demonstrated that the third-neighbor exchange dominates in NiPS3 due to the filling of the t2g subshell, whereas for MnPS3, the first-neighbor exchange is prevailed, owing to the presence of the t2g magnetism.

We showed in this work that the onset of the nearest-neighbor FM coupling in NiPS3 is due to the relatively strong (FM) superexchange, which is enabled by the complete filling of the t2g shell and by the (close to) 90° nickel–sulfur–nickel nearest neighbor bond. Nevertheless, even these relatively “fortunate” circumstances do not guarantee that the strong further neighbor exchange can also be FM. The reason for this lies in the further neighbor bonds, which are no longer of a “90° variety”—and the latter bonds are, apart from specific Jahn–Teller effects, the only way to obtain FM couplings in Mott insulators. From a more general perspective, this study confirms the paradigm that the AFM couplings are natural to the Mott insulating compounds.

Second, we examined the properties of the MnPS3 and NiPS3 compounds doped with 25% impurities (Ni in MnPS3, Mn in NiPS3, and Cr in both hosts). It turned out that all of the investigated alloys are Mott insulating, albeit with generally smaller band gaps than those of the corresponding host. Crucially, we demonstrated an extreme robustness of the AFM phases against impurity doping of the MnPS3 and NiPS3 compounds. Hence, ferromagnetism cannot be easily stabilized by impurity doping, in agreement with the above paradigm. Nevertheless, as the dopants have different spins than the pure phases, the alloys exhibit ferrimagnetic properties for particular arrangements of dopants.

The Mn and Ni impurities prefer to form dimers within the host, whereas the Cr dopants prefer to be further apart. Interestingly, unlike for the hosts, the first and second (dopant–host) exchange couplings are of similar order of magnitude. The latter leads to frustration in the case of AFM couplings. We suggest that this may be one of the reasons of the observed lower magnetic ordering temperature of the doped systems.[35]

Our work sheds light on the origin of magnetism in the AFM family of transition metal phosphorus trichalcogenides by pointing out the mechanisms which govern the benchmark compounds, thus extending the fundamental knowledge of 2D magnetism.