Introduction {#Sec1}
============

Due to their importance in the field of energy-saving technologies for the next generation of electrical devices, permanent magnets that do not contain rare-earths or platinum are a current focus of research^[@CR1]^. *α*″-Fe~16~N~2~, which is a martensite, has been studied over many years, and several groups have claimed to observe giant magnetic saturation moments in it^[@CR2]--[@CR4]^. It was first experimentally prepared by K. H. Jack^[@CR5]^ during low-temperature annealing of *α*′-FeN. It appeared as a metastable phase that survived up to approximately 500 K, above which it decomposed into *α*-Fe and Fe~4~N.

The existence of a giant saturation moment in *α*″-Fe~16~N~2~ was first reported by Kim *et al*.^[@CR2]^ and later confirmed by Sugita *et al*.^[@CR3]^ in single-crystal Fe~16~N~2~ that was grown epitaxially on InGaAs.  Since then, several investigations have explored the validation of the giant saturation moment in this material.  However, these investigations resulted in more controversy rather than bringing the issue to a firm conclusion as some of these studies confirmed the existence of a giant moment^[@CR6],[@CR7]^ while others were unable to reproduce such results^[@CR8]^.

Therefore, research on Fe~16~N~2~ has always faced two challenges: (1) stabilizing the material at higher temperature so that it does not decompose into *α*-Fe and Fe~4~N and (2) understanding its complex magnetic behavior, which made the material a "*40 year old mystery*"^[@CR9]^ among the magnetic materials.  In a recent work, we have demonstrated that^[@CR10]^, Fe~16~N~2~ can be stabilized via vanadium (V) doping.  Recently.  Ke *et al*.^[@CR11]^ studied the exchange interaction and ferromagnetic transition temperature for Fe~16~N~2~ and Co- and Ti-doped Fe~16~N~2~ using a first-principles-based method.  They obtained the Curie temperature (T~*C*~) via both mean field approximation (MFA) and random phase approximation (RPA).  Direct measurement of T~*C*~ of Fe~16~N~2~ is limited due to the decomposition of the material into Fe~4~ and Fe at temperatures above 500 C^[@CR2],[@CR5]^. The Curie temperatures that were obtained theoretically by Ke *et al*.  were almost three times higher than the temperature that was obtained by Sugita *et al*.  by extrapolating from experimental data that were obtained at lower temperature. 

Understanding the magnetism and the exchange mechanism in Fe~16~N~2~ is extremely important.  In this paper, we have studied the magnetic exchange interactions and the Curie temperatures in Fe~16~N~2~ and V-doped Fe~16~N~2~. The magnetic exchange constants were obtained via a first-principles-based method and the Curie temperatures were calculated via both MFA and MC simulation.  In this paper, we also address the long-standing question about the appearance of anomalously large magnetic moments.  We organize the remainder of our paper as follows: In Section II, we discuss the magnetic exchange interactions and the Curie temperatures.  In the following Section III, we discuss the possible origin of the giant magnetic moments in the system.  We also introduce another potential application of this material in Section IV: Fe~16~N~2~ could be used as a spin injector, in addition to its use as a permanent magnet.  Finally, in Section V and VI respectively, our conclusions and the computational methods are presented. 

Results and Discussion {#Sec2}
======================

Fe~16~N~2~ has a body-centered tetragonal structure with space group I~4~/*mmm* (number 139).  The unit cell contains two N atoms and sixteen Fe atoms.  The three Wyckoff positions that are occupied by the Fe atoms are 8h (eight atoms), 4e (four atoms), and 4d (four atoms).  Our PBE- optimized lattice constants, namely, a = b = 5.68 Å, c = 6.22 Å, and position parameters, namely, x~8*h*~ = 0.243 and z~4*e*~ = 0.294, well agree with the values that were measured by Jack *et al*.^[@CR5]^ (a = b = 5.72 Å and c = 6.29 Å, x~8*h*~ = 0.242, and z~4*e*~ = 0.293).  For comparison, we have considered two additional variations of PBE, namely, PBE-sol^[@CR12]^ and revised PBE^[@CR13]^, which yielded similar results: The average magnetic moments per Fe are 2.43 *μ*~*B*~ (PBE), 2.45 *μ*~*B*~ (PBE-sol) and 2.48 *μ*~*B*~ (rev PBE).  As the magnetic moments only differ in the second decimal place, we performed other calculations with the PBE functional only. 

Exchange constants and curie temperature {#Sec3}
----------------------------------------

In Fig.  [1](#Fig1){ref-type="fig"}, we show the calculated exchange constants that were obtained via the SPR-KKR method using the PBE functional.  Although the trend of the results is similar to that of the results that were obtained by Ke *et al*.^[@CR11]^ via the TB-LMTO method under the local density approximation (LDA) and the GW method, there are substantial differences.  The largest among all interactions is between the 8h-4d sites (30.89 meV) at a distance of 2.54 Å, as shown in Fig.  [1(a)](#Fig1){ref-type="fig"}), while for Ke *et al*.  it was between the 8h-4e sites.  The next-strongest interaction from our calculation is between the 8h-4e sites (20.15 meV) (Fig.  [1(a)](#Fig1){ref-type="fig"}) at a distance of 2.43 Å between two octahedral Fe atoms.  The third-strongest interaction (15.34 meV) is between the octahedral 8h and 4e sites at a distance of 2.67 Å. Due to structural considerations, this interaction should be mediated through a 90° Fe8h-N-Fe4e superexchange interaction.  Among the inter-sublattice interactions, the 4d-4e interaction is the weakest (0.74 meV) (Fig.  [1(b)](#Fig1){ref-type="fig"}).  According to Fig.  [1(a)](#Fig1){ref-type="fig"}, the 8h-8h interactions (between Fe atoms that are residing at different Wyckoff positions) are antiferromagnetic (−2.97 meV) at the smallest distance (2.76 Å).  These interactions may correspond to a 90° Fe8h-N-Fe8h superexchange interaction.  The strongest Fe8h-Fe8h interaction is ferromagnetic (3.58 meV) at 3.9 Å and it is due to a 180° Fe8h-N-Fe8h superexchange interaction.  All the strong magnetic interactions that involve either 8h-4d or 8h-4e sites are nearly along the \[111\]-direction, similar to bcc Fe; however, due to the structural distortion (in the case of Fe~16~N~2~), these lattice vectors are slightly deviated.  For example, the vector that connects the 8h and 4d sites is \[0.25, 0.24, 0.27\].  The interactions between the 4d and 4e sites are nearly along the \[011\]-direction.  For bcc Fe, the most prominent interaction along the \[011\]-direction is anti-ferromagnetic, as calculated by Pajda *et al*.^[@CR14]^; however, in the present case, it is still ferromagnetic, although very weak.  This suggests that Fe~16~N~2~ is a much stronger ferromagnet than bcc Fe. Figure 1Exchange parameters for (**a**) Fe8h, (**b**) Fe4e, and (**c**) Fe4d as centers. (** d**) The calculated Curie temperatures for various values of the radius of the sphere, within which all interactions were considered. 

Among the intra-sublattice interactions (between the Fe atoms that are at the same Wyckoff positions), the 4e-4e interaction is the strongest (12.39 meV), which corresponds to the interaction between two Fe4e atoms in two neighboring octahedra.  The next-strongest 4e-4e interaction is between two epical Fe atoms in the Fe-N octahedra and, therefore, is mediated by a 180° superexchange interaction.  The strongest 4d-4d interaction is 5.3 meV at a pair distance of 3.1 Å and, therefore, is almost negligible. 

In Fig.  [1(d)](#Fig1){ref-type="fig"}, we show the variation of the mean-field Curie temperature as a function of the cluster radius, which is denoted as \|**r**\|.  Even for a small radius of approximately \|**r**\| = 0.5 Å, T~*C*~ is much higher than the experimental value (813 K); hence, the mean-field picture does not explain the T~*C*~ value of this material.  MFA overestimates T~*C*~ due to neglect of the thermal fluctuations of magnetization.  However, in the present case, MFA overestimates T~*C*~ by almost 50%, which is large.  Such overestimation of T~*C*~ of Fe~16~N~2~ by MFA was also observed by Ke *et al*.;  however, they attributed this to the lack of a proper experimental measurement of T~*C*~. Such consistent huge overestimation of T~*C*~ by MFA reflects a low-dimensional behavior of the spin system compared to a pure 3D one.  The origin of such quasi-2D-like behavior could be the poor connectivity between spins that differ in terms of direction.  There is almost no exchange between spins at the 4e and 4d sites and a very weak 4d-4d interaction is observed in Fig.  [1](#Fig1){ref-type="fig"}. 

To obtain more accurate results, we performed MC simulation using the metropolis method^[@CR15]^. In Fig.  [2](#Fig2){ref-type="fig"}, we present our results for the MC simulations.  To obtain the value of T~*C*~, we fit the temperature-dependent magnetization data to the function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$M(T)={(1-\frac{T}{{T}_{C}})}^{\beta }$$\end{document}$, which yields a T~*C*~ value of approximately 765 K according to Fig.  [2(a)](#Fig2){ref-type="fig"}, which is much closer to the experimental value compared to the MFA results that were obtained by us and those that were obtained by Ke *et al*.  For comparison, we also performed MC simulation for the V-doped Fe~16~N~2~. We demonstrated in our previous study^[@CR10]^ that Fe~16~N~2~ can be stabilized by doping V at site 8h; hence, we performed the MC calculations with exchange constants that were obtained via the SPR-KKR method by replacing one Fe-8h with V. Figure [2(b)](#Fig2){ref-type="fig"} shows the obtained results.  The T~*C*~ value that was obtained after fitting is approximately 640 K, which is less than that for bulk Fe~16~N~2~, but still above room temperature. Figure 2Monte Carlo simulation of temperature-dependent magnetization.  The top panel (a) corresponds to pure Fe~16~N~2~ and the bottom panel (b) to V-doped Fe~16~N~2~.

Understanding the *giant magnetic* Moment {#Sec4}
=========================================

We will discuss the origin of the giant magnetic moment in this section.  As already discussed in the introduction, the experimental findings about such moments are not consistent.  Therefore, the reliability of such reports is uncertain.  In the following, we analyze the possibility of a *giant* moment being present in this system.  The following two factors that might influence the saturation magnetization are discussed here: (1) extrinsic effects such as strain, which we analyze from the fixed magnetic moment calculations, and (2) intrinsic effects such as the existence of localized or semi-localized d-states within the octahedral region and exchange interactions that favor the high-spin configuration of the 3d-states, which we investigate based on the magnetic susceptibility and the nature of the exchange interactions. 

Structural sensitivity and magnetism: constrained magnetic calculations {#Sec5}
-----------------------------------------------------------------------

To elucidate the relation between the structure and magnetism, we performed constrained moment calculation with the PBE exchange correlation function.  The magnetic moment of the cell was fixed such that each Fe site had a magnetic moment of approximately 3 *μ*~*B*~. Under this constraint, we relaxed the cell geometry completely, which yielded lattice parameters a = b = 5.84 Å and c = 6.34 Å. The c/a ratio remains close 1.10 (as in the unconstrained case); hence, the large magnetic moment does not result in any additional structural anisotropy, such as additional tetragonality.  The space-group (No.  139) remains almost the same after the structural optimization with giant moments.  The only difference appears in the position parameters with z~4*e*~ = 0.3 (0.293 for unconstrained case) and x~8*h*~ = 0.240 (0.242 in the unconstrained case).  This results in a slight increase in the distance between N and Fe at the 8h and 4e sites.  Next, we performed a calculation that involves the optimization of positions while the volume of the cell is fixed to the volume that was obtained via the constrained magnetic calculation that is described above.  We allowed the atomic spin moments to change with the electronic steps in this case.  This resulted a magnetic moment of 2.52 *μ*~*B*~ per Fe atom, which is slightly larger than the values that are reported in Table [1](#Tab1){ref-type="table"}.  This calculation suggests a weak correlation between the structure and the magnetic moment in Fe~16~N~2~. Therefore, it is highly unlikely for such large moments to appear due to strain effects that arise from substrates or via other means.  Therefore, the origin of giant magnetic moment should be related to intrinsic physics such as partial localization of d-electrons in the Fe-N octahedra or the nature of the exchange interaction promoting the maximization of the local spins through Hund's coupling.  The partial localization of d-electrons in Fe~16~N~2~ was demonstrated by Ji *et al*.^[@CR9]^; however, they demonstrated it in terms of the difference in charge density between inside and outside the Fe~6~N octahedra. Table 1Lattice constants and magnetic moments in three types of GGA calculations. GGAa (Å)c (Å)m~8*h*~ (*μ*~*B*~)m~4*e*~ (*μ*~*B*~)m~4*d*~ (*μ*~*B*~)x~8*h*~z~4*e*~PBE5.686.222.362.182.820.2420.293PBE-sol5.756.282.402.172.840.2420.294Revised PBE5.756.282.422.222.860.2430.294

Intrinsic effects: partial localization of d-electrons in Fe~6~N octahedra and high spin states being favored by exchange interactions {#Sec6}
--------------------------------------------------------------------------------------------------------------------------------------

In Fig.  [3](#Fig3){ref-type="fig"}, we show the calculated magnetic susceptibility as a function of the temperature.  The susceptibility shows a sharp peak at T = T~*C*~, which indicates the presence of localized states in this material.  For itinerant systems, the magnetic susceptibility usually shows a broad peak below T = T~*C*~. Therefore, it is important to investigate the possibility of localized electronic states that contribute to the magnetic moments. Figure 3Magnetic susceptibility of bulk Fe~16~N~2~ that was obtained from MC simulations. 

In Fig.  [4](#Fig4){ref-type="fig"}, we analyze the magnetic exchange interactions (obtained via SPR-KKR) that involve 8h sites within the framework of the tight-binding model.  Within the tight-binding model, the intersite hopping integral, which is denoted as *t*~*dd*~, between the Fe sites that involve only d-orbitals can be expressed as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${t}_{dd} \sim {(r/a\mathrm{0)}}^{-5}$$\end{document}$^[@CR16]^. First, we discuss the strongest among all the interactions: the interaction between the Fe-3d states at the 8h and 4e sites.  In Fig.  [4(a)](#Fig4){ref-type="fig"}, for the 8h-4d interactions, we fit the exchange constants with two itinerant-type interactions: the double-exchange and RKKY interactions.  Typically, the double-exchange interaction is directly proportional to *t*~*dd*~, namely, *r*/*a*0^−5^, while the RKKY interaction is proportional to (*r*/*a*0)^−3^. According to the figure, *J*~8*h*−4*d*~ accords extremely well with the double-exchange type.  Due to its double-exchange nature, the *J*~8*h*−4*d*~ interaction is strongly ferromagnetic.  However, for the 8h-4e interactions (Fig.  [4(b)](#Fig4){ref-type="fig"}), a linear combination of double-exchange and superexchange interactions of the form *J* = *α*(*r*/*a*0)^−5^ + *β*(*r*/*a*0)^−10^ accords well (as the superexchange interaction has $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${t}_{dd}^{2}$$\end{document}$-dependence), where *α* and *β* are the fitting parameters.  Finally, for 8h-8h interactions, we fit the exchange constants with four types of interactions: (1) linear combinations of double-exchange and superexchange interactions (red curve), (2) double-exchange interactions (blue curve), (3) RKKY interactions (green curve), and (4) combinations of RKKY and superexchange interactions (yellow curve).  According to the figure, the behavior of the data most strongly resembles an interaction that is represented by case (4), namely, a combination of double-exchange and superexchange interactions.  The substantial tendency of the exchange interactions within the octahedra toward superexchange interactions suggests partial weak localization of electron states at 8h-sites. Figure 4Numerically obtained magnetic exchange interactions that involve the octahedral 8h site: (**a**) 8h-4d, (**b**) 8h-4e, and (**c**) 8h-8h interactions.  These values are fitted using various models, as discussed in Section IV. 

The strongest magnetic interaction in this material is between 8h and 4d sites (30.89 eV per pair) and it is of the double-exchange type.  This was first predicted by Sakuma *et al*.^[@CR17]^. Our conclusion regarding the high moment is as follows: Due to the presence of N, there is a moderate correlation effect in the octahedral region that localizes partially the d-states within this region (especially the 8h states), along with a strong tendency of double-exchange interactions to prefer high-spin states^[@CR18]^.

To determine the extent to which the partial localization of the octahedral d-states may affect the magnetic behavior or the saturation moment of Fe~16~N~2~, we computed the magnetic moments per unit cell via the DFT + U method^[@CR19]^ with various values of U, as shown in Fig.  [5](#Fig5){ref-type="fig"}.  Switching on the localization effect at a moderate level (3--5 eV) gives rise to a magnetic moment that is as high as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$ \sim 2.8{\mu }_{B}$$\end{document}$. For typical Fe-based insulating systems, the typical choice for U is approximately 6 eV. As Fe~16~N~2~ is a metal, the U value should be smaller than 6 eV.Figure 5Magetic moments per Fe for various U values at the 8h site. 

Fe~16~N~2~ as a Spin Injector {#Sec7}
=============================

The bulk spin polarization that we obtained for the PBE-optimized Fe~16~N~2~ is approximately 60%, which is large compared to the elemental metallic ferromagnets, such as Fe (45%), Co (42%) and Ni (27%)^[@CR20]^. Such large spin polarization falls in between the half-metallic Heusler alloys (\~100%), which were proposed to be the optimal spin-injectors for the semiconductors^[@CR21]^, and the elemental ferromagnets.  However, until now, the Heusler alloys were not highly successful spin injectors due to their large lattice mismatch with the semiconductors^[@CR22],[@CR23]^. Fe~16~N~2~ can be easily grown on III--V semiconductors such as GaAs or InGaAs (historically, the giant moment was observed in epitaxially grown single crystals of Fe~16~N~2~ on either InGaAs (001) or GaAs(001) film).  As Fe~16~N~2~ has both large spin polarization and satisfactory lattice matching with semiconductors, we examined its application as a spin injector as follows.  We performed calculations with 4 monolayers (ML) of Fe~16~N~2~ (001) on 4 ML of InGaAs (001) and 8 ML of InGaAs (001).  The lattice parameters of bulk InGaAs were obtained.  The optimized structure was used to construct the metal/semiconductor (001) interfaces, as shown in Fig.  [6](#Fig6){ref-type="fig"}.  The magnetic moments per Fe for configurations (a), (b) and (c) are 2.48 *μ*~*B*~, 2.51 *μ*~*B*~ and 2.52 *μ*~*B*~, respectively.  The magnitude of the moment depends on the nature of the interface.  The magnetic moment is higher for interfaces that contain either Ga or In.  These moments are not *giant*; however, they are comparable to the maximum value of the saturation moment that is predicted by the Slater-Pauling curve^[@CR24]^. Then, we investigate the possibility of spin injection from the ferromagnet (FM) to the semiconductor (SC) in such systems.  As an example, we have computed the layer-resolved spin polarization, namely,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${P}_{L}=\frac{{D}_{L}^{\uparrow }({E}_{F})-{D}_{L}^{\downarrow }({E}_{F})}{{D}_{L}^{\uparrow }({E}_{F})+{D}_{L}^{\downarrow }({E}_{F})},$$\end{document}$$where L is the layer index and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${D}_{L}^{\uparrow (\downarrow )}({E}_{F})$$\end{document}$ is the density of states at the Fermi energy for the majority (minority) of spin carriers for interface (b) in Fig.  [6](#Fig6){ref-type="fig"}.  The results are shown inFig.  [7](#Fig7){ref-type="fig"}, which presents the layer-dependent magnetic moments in the top panels.  On the Fe~16~N~2~ side, all four layers have moments that are close to 2.5 *μ*~*B*~ (the average is approximately 2.51 *μ*~*B*~), while the spin polarization (bottom panel) decreases as we approach the interface.  However, on the semiconductor side (the 5^*th*^ to 8^*th*^ layers), the spin polarization is high and is non-negligible, even up to the last layer.  The average spin injection efficiency over all four layers is approximately 30%, which is higher than (almost double) what is observed for Fe/GaAs heterostructures^[@CR25]^. This implies a long spin diffusion length in Fe~16~N~2~/InGaAs systems.  The spin polarization changes sign across the interface; hence, the spin imbalance in the semiconductor side is due to the minority spin. Figure 6Fe~16~N~2~(001)/InGaAs(001) thin films with 4 monolayers of metal and semiconductor for three configurations. Figure 7(Top panel) Layer-resolved magnetic moments per atom for the Fe~16~N~2~(001)/InGaAs(001) thin-film configuration (**b**) of Fig.  [6](#Fig6){ref-type="fig"}. ( Bottom panel) The corresponding layer-resolved spin polarization. 

Fe~16~N~2~ seems to be a highly effective spin injector for III--V semiconducting systems such as InGaAs.  As shown in our previous study^[@CR10]^, small V doping stabilizes Fe~16~N~2~; therefore, its interface with III--V semiconductors might be important for spin-optoelectronic or spintronic^[@CR26],[@CR27]^ applications. 

Conclusions {#Sec8}
===========

We have studied the exchange interactions and magnetic properties of Fe~16~N~2~ and V-doped Fe~16~N~2~ via first-principles-based methods.  The magnetism and T~*C*~ are dominated by the exchange interactions between 8h-4d and 8h-4e sites.  By fitting the exchange parameters that were obtained via the first-principles-based method to various theoretical models of exchange interactions, we demonstrate that in addition to the double-exchange interaction, a substantial role is played by the superexchange interaction.  Although the double-exchange interaction being the strongest favors a high-spin state, the presence of the superexchange interaction among the Fe sites within the Fe~6~N octahedral region indicates the existence of localized states, which may also contribute to the occurrence of a giant moment in this material.  In addition, due to its high spin polarization and fantastic lattice matching with the III--V semiconductors, Fe~16~N~2~ can also find important spintronic applications in addition to the permanent magnet application. 

Computational Method {#Sec9}
====================

The electronic structure calculations were performed via first-principles methods within the framework of density functional theory (DFT) with the Perdew-Burke-Ernzerhof exchange correlation energy functional^[@CR28]^ based on a generalized gradient approximation.  We used a projector-augmented wave method as implemented in the Vienna *ab-initio* simulation package (*VASP*)^[@CR29]^. Kohn-Sham wave functions of the valence electrons were expanded in a plane-wave basis with an energy cut-off of 450 eV. The augmentation charge cut-off was set to 627.1 eV. The Brillouin zone sampling was conducted using a Monkhorst Pack grid of 7 × 7 × 7 k-points.  Ionic relaxation was performed using the conjugate-gradient method until forces were reduced to within 0.01 eV/Angstrom.  The calculated structural parameters, such as the lattice constant, positions, and magnetic moments, are employed to obtain the exchange parameters using Korringa-Kohn-Rostoker (KKR) multiple-scattering theory.  We have used the Munich *SPR-KKR* package^[@CR30]^ to obtain the exchange parameters via a real-space approach that was formulated by Liechtenstein *et al*.^[@CR31]^, which maps the total energy of the system to an effective Hamiltonian that is expressed as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$ {\mathcal H} =-\sum _{i\ne j}\,{J}_{ij}{{\bf{e}}}_{i}{{\bf{e}}}_{j}$$\end{document}$$where **e**~*i*~ represents the vector of the local magnetic moment of the i^*th*^ site and *J*~*ij*~ denotes the exchange constant between sites *i* and *j*.  The angular momentum expansion of the basis functions was calculated up to l = 3.  The k-space integration was performed using a grid of 280 k-points in the irreducible part of the Brillouin zone.  We have used 30 complex energy points to perform the integration over the Green's function.  The Curie temperature within MFA was calculated by solving the following coupled equations:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\langle {{\bf{e}}}^{\alpha }\rangle =\frac{2}{3{K}_{B}T}\sum _{\beta }\,{J}_{0}^{\alpha \beta }\langle {{\bf{e}}}^{\beta }\rangle $$\end{document}$$where 〈**e**^*α*^〉 is the average z-component of the magnetic moment in sublattice *α* and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${J}_{0}^{\alpha \beta }={\sum }_{{\rm{r}}}\,{J}_{0{\rm{r}}}^{\alpha \beta }$$\end{document}$. The summation in **r** is performed up to \|**r**\|/*a* = 3.0.  To obtain an accurate value of T~*C*~, we conducted Monte Carlo simulations with the *VAMPIRE* atomistic spin dynamic program^[@CR32]^. We used 67200 atoms with periodic boundary conditions.  The DFT-calculated magnetic moments (from *VASP*) and exchange constants (from *SPR-KKR*) were used as inputs.  The Hamiltonian of the system contained the following exchange and anisotropy terms:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$ {\mathcal H} =-\sum _{ij}\,{J}_{ij}{{\bf{e}}}_{i}{{\bf{e}}}_{j}-{K}_{u}\sum _{i}{({{\bf{e}}}_{i}\cdot {\bf{e}})}^{2}$$\end{document}$$where **e** is the unit vector along the direction of the easy axis and K~*u*~ is the uniaxial anisotropy constant.  The values of K~*u*~ for pure and doped Fe~16~N~2~ were obtained from our previous paper^[@CR10]^.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 

This work was supported by NRF grants that were funded by MSIP, Korea (No.  2009-0082471 and No.  2014R1A2A2A04003865), the Convergence Agenda Program (CAP) of the Korea Research Council of Fundamental Science and Technology (KRCF) and the GKP (Global Knowledge Platform) project of the Ministry of Science, ICT and Future Planning. 

S.B. and S.C.L. conceived the idea.  S.B. performed all the calculations and wrote the paper.  S.C.L. made comments. 

The authors declare no competing interests.