Robust Ferromagnetism in Hexagonal Honeycomb Transition Metal Nitride Monolayer

Two-dimensional intrinsic magnetic materials with high Curie temperature are promising candidates for next-generation spintronic devices. In this work, we design two kinds of two-dimensional transition metal nitrides, VN2 and FeN2, both with a hexagonal honeycomb lattice. Based on the formation energy, and phonon spectra calculations as well as the molecular dynamics simulations, their structural stability is demonstrated. Then, we determine the ferromagnetic ground states of VN2 and FeN2 monolayers through the energy calculations, and the Curie temperatures of 222 K and 238 K are estimated by solving the Heisenberg model using the Monte Carlo simulation method. Hence, the VN2 and FeN2 monolayers are demonstrated to be new two-dimensional ferromagnetic materials with high temperature ferromagnetism or large-gap half-metallicity.


Introduction
Because of the novel mechanical, thermal, electronic, and magnetic properties associated with the reduced dimensionality, two-dimensional (2D) materials have attracted a lot of research interest in physics and materials.Since graphene was discovered in 2004 [1], various 2D materials have been fabricated in experiments, such as boron nitride [2], silicene [3], borophene [4], and transition metal dichalcogenides [5].Recently, more and more attention has been paid to 2D transition metal nitrides, because they have good mechanical and electronic properties, include flexibility and suitable energy band gap [6][7][8][9][10][11].
Some MXenes [12] are typical 2D transition metal nitrides, which commonly appear as M 2 N, M 3 N 2 , and M 4 N 2 (M is a transition metal atom) [13].In the M 2 N (M = Ti, Cr, and Zr) monolayers, one layer of N atoms is sandwiched between two layers of M atoms.If the two elements of the M 2 N layer are exchanged, is the new structure of the MN 2 layer stable?After the Surface modification with the OH moiety or O atom, the nonmagnetic Cr 2 N monolayer can become the ferromagnetic Cr 2 N(OH) 2 or Cr 2 NO 2 monolayer [14,15].In them, the Cr layer is between the N and O layers.Inspired by the emergence of ferromagnetism in Cr 2 N(OH) 2 and Cr 2 NO 2 , a question is raised.If one layer of metal atoms is intercalated between two layers of N atoms, namely the MN 2 structure mentioned above, would the ferromagnetism emerge?Although the initial idea to build the 2D ferromagnetic monolayer is simple, there are two facts to motivate us to further explore the 2D ferromagnetism in the MN 2 systems.In the previous research on the 2D transition metal nitrides, MoN 2 , ReN 2 , and TaN 2 have been predicted and synthesized [10,16,17].From this, we infer that the MN 2 compounds are probably also stable when M is 3d transition metal atom.On the other hand, the M atoms in the MN 2 structure form a triangle lattice, which is more preferred for ferromagnetic order than antiferromagnetic order because the antiferromagnetic interactions are usually frustrated in a triangle lattice.
The Mermin-Wagner theorem states there is no long-range 2D magnetic order at any finite temperature in an isotropic spin orientation case [18].This theorem can lead to the misleading impression that it is difficult to prepare magnetic two-dimensional materials experimentally and then hinder the experimental development of the magnetic 2D materials.In fact, the long-range 2D magnetic order can exist in a 2D materials when the isotropic spin orientation is broken.Until 2017, the 2D ferromagnetism in CrI 3 monolayer with the Curie temperature of 45 K was discovered [19], which is the first experimental realization of the magnetic long-range order in a 2D material.The milestone discovery has stimulated great research interest in the ferromagnetism in 2D materials.Subsequently, 2D Cr 2 Ge 2 Te 6 [20] was synthesized, and the ferromagnetism with the Curie temperature of about 30 K was measured by scanning magneto-optic Kerr microscopy.In the 2D VSe 2 layer synthesized by Bonilla et al. [21], the room-temperature ferromagnetism was observed.Recently, in atomically thin Cr 2 Te 3 , a conspicuous ferromagnetic transition from Stoner to Heisenberg-type is directly observed, indicating that the Heisenberg model is suitable to describe the magnetic interactions in 2D magnets with the dual nature of localized and itinerant ferromagnetism [22].Anomalous Hall effect was confirmed by magnetotransport measurements in single-sheet 2D Fe 3 GaTe 2 crystals with room-temperature ferromagnetism [23].In the CoPS 3 , FePS 3 , and NiPS 3 2D layers with an antiferromagnetic zigzag order, the magneto-elastic response was investigated and spontaneous anisotropic magnetostriction was found [24].Magnetic 2D materials offer an emerging platform for fundamental studies of magnetism in the 2D limit.In addition, the Curie temperature is a vital parameter of ferromagnetic 2D materials, and the low Curie temperature is the main obstacle limiting the application of ferromagnetic 2D materials in electronic devices.Up to now, ferromagnetic 2D materials with high Curie temperatures are scarce.Therefore, it is an urgent and necessary task to explore new two-dimensional ferromagnetic materials for both fundamental physics studies and application research on functional materials.
In this work, we propose two new ferromagnetic 2D transition metal nitrides, VN 2 and FeN 2 monolayers, and systematically investigate their dynamical, electronic, and magnetic properties.Based on the Heisenberg model and the Monte Carlo simulations, the Curie temperatures are estimated, and the high-temperature ferromagnetism in the 2D FeN 2 and VN 2 layers is demonstrated.

Atomic Structure
We design three structures of MN 2 (M = V and Fe) with the M layer intercalated between two N layers, shown in Figure 1a-c.The first structure is a hexagonal (H) phase with the symmetry of the P 6m2 (No. 187) space group, in which the M atom is coordinated by six N atoms and located at the center of a regular triangular prism of six nitrogen atoms.This is the hexagonal honeycomb geometry, similar to the well-known graphene, boron nitride, and MoS 2 monolayers.The second structure is a trigonal (T) phase with the symmetry of the P 3m1 (No. 164) space group, in which the M atom is in an octahedral configuration bonded to six N atoms.The third structure is a tetragonal phase (Tetra) with the symmetry of the P 4m2 (No. 115) space group, in which the M atom is at the center of a tetrahedron formed by four N atoms.Table 1 presents the energies of the MN 2 (M = V and Fe) monolayers in the hexagonal, trigonal, and tetragonal phases.As shown, the hexagonal phase has the lowest energy and is determined to be the most likely structural phase of the 2D MN 2 (M = V and Fe) compounds.In the hexagonal MN 2 monolayer, the N-N distances are 1.568 Å and 1.336 Å for VN 2 and FeN 2 , respectively, which are close to the length of the N-N single bond, 1.46 Å .The formation of the N-N bond is the main reason why the hexagonal structure has a lower energy compared to the trigonal and tetragonal phases.So, the hexagonal VN 2 and FeN 2 monolayers are selected as our research objects, and the mechanical, electronic, and magnetic properties are investigated in the following sections.

Evidence for Structural Stability
As analyzed in the above paragraph, the hexagonal phase of FeN 2 and VN 2 layers is the low-energy structural phase, in which we note that one unit cell only contains three atoms, one Fe (or V atom) and two N atoms.To further inspect the structural stability of hexagonal FeN 2 and VN 2 monolayers, we compute the formation energies with the standard generalized gradient approximation in the Perdew-Burke-Ernzerhof form.The formation energy per atom is expressed as in which E MN 2 is the energy of one unit cell of the MN 2 monolayer, E M is the energy of one metal atom in the bulk metal, E N 2 is the energy of a nitrogen molecule, and the number 3 is the total number of atoms in one unit cell. in which E MN 2 , E M , and E N 2 are the energies of the MN 2 monolayer, bulk metal per atom, and nitrogen molecule, respectively.The formation energy averaged over the number of atoms is −0.113 eV/atom and 0.665 eV/atom for FeN 2 and VN 2 , respectively.Because N 2 molecules have low energy due to their N≡N triple bond, metal nitrides usually have positive formation energies.For example, the formation energies of the nitrides CuN 3 [25], PtN 2 [26] and g-C 3 N 4 [27] that have already been synthesized experimentally are 0.42 eV/atom, 0.55 eV/atom, and 0.35 eV/atom, respectively.The formation energies of the 2D FeN 2 and VN 2 monolayers are comparable to the ones of these nitrides, which indicates that similar to the CuN 3 , PtN 2 , and g-C 3 N 4 , the FeN 2 and VN 2 monolayers can exist under ambient conditions.The phonon spectra is an important means to judge the structural stability of solid materials.Based on the DFTP method, we perform the calculation of the phonon spectra of FeN 2 and VN 2 to inspect its dynamic stability.The highly symmetric points in reciprocal space are Γ (0 0 0), M (0.5 0 0), and K (1/3 1/3 0).If there exists an unstable phonon mode, it is presented as an imaginary frequency.As shown in Figure 2a,c, the absence of imaginary modes in the entire Brillouin zone indicates that the monolayer VN 2 and FeN 2 are dynamically stable.
Next, we perform the first-principles molecular dynamic simulations to examine their thermal stability.The temperature is set to 300 K, and the time length is 5 ps.In Figure 2b,d, the total potential energy of the MN 2 (M = Fe and V) monolayer fluctuates around a certain value, and no distinct drop of energy occurs, which indicates that the monolayer framework is retained at the temperature of 300 K and there is no bond broken.The final configurations of FeN 2 and VN 2 monolayer after 5 ps molecular dynamics simulations are shown in the inserted diagrams in Figure 2b,d, reflecting the robustness of the FeN 2 and VN 2 structures.Therefore, the 2D MN 2 (M = Fe and V) monolayers have good dynamic and thermal stability.

Electronic Structure
GGA + U is a simple method to deal with the highly localized and strongly correlated electron systems in which the Hubbard U parameter has a clear physical significance.To determine the value of Hubbard U for a certain material, M. Cococcioni and S. de Gironcoli [28] developed a linear response approach inclosed in the VASP code.By this method, we obtain the Hubbard U values of 2.2 eV and 3.42 eV for the FeN 2 and VN 2 monolayer.
In the hexagonal MN 2 (M = Fe and V) monolayer, the M atom is in a regular triangular prism and the crystal field of the M atom has the C 3h symmetry.The five 3d suborbitals are split into three groups, namely d x 2 −y 2 /d xy , d xz /d yz , and d z 2 , which are labeled as E ′ , E ′′ , and A ′ in terms of symmetry signals in the character table of C 3h .The magnetic moment of Fe or V atom is from the spin moment of unpaired electrons, which is determined by the difference between the number of electrons occupying the spin-up and spin-down states.In Figure 3, the states below Fermi energy are fully occupied and those electronic states above the Fermi energy are empty.For the density of states, it means the number of states per one energy unit (1 eV).By integrating the density of states from the lowest energy to the Fermi energy, we can get the number of occupied electrons.Consequently, for each atomic orbital, the net spin-polarized charge in this orbital is derived and its contribution to the total moment of Fe or V atom is figured out.For the FeN 2 monolayer, the partial density of states for five 3d suborbitals of the Fe atom is displayed in Figure 3a-c.As shown, these orbitals are strongly spin-polarized.Accordingly, the large moment of 3.2 µ B is formed around Fe atoms.In the two degenerate E ′ orbitals, d x 2 −y 2 and d xy , all the spin-up states and a small number of spin-down states are occupied, and the electrons in the two orbitals give rise to the magnetic moment of about 1.5 µ B .In the two degenerate E ′′ suborbitals, d xz and d yz , the electrons occupied in the spin-up states are slightly more than the electrons in the spin-down states, and they can give rise to the moment of about 0.5 µ B .As for the d z 2 suborbital, the spin-up states are fully occupied and the spin-down states are empty, which lead to 1.0 µ B moment.In addition, there is some spin-polarization in the s and p states of the Fe atom, which has a small contribution to the total moment.On the other hand, the V atom in the VN 2 monolayer, the moment is 1.0 µ B .As seen from Figure 3f-h, most of the moment comes from the dz 2 suborbital and a small part of it is from other d suborbitals, namely 0.67 µ B and 0.33 µ B , respectively.
The energy bands of MN 2 (M = Fe and V) are shown in Figure 4a,b.Their spin-up and spin-down bands have an obvious spin-splitting, and two bands cross the Fermi energy, indicating that the FeN 2 and VN 2 are magnetic metals.For the VN 2 bands, there are only spin-up bands near the Fermi energy, and no spin-down bands appears in the energy range from −0.82 eV to 0.31 eV.So, the 2D VN 2 acts as a conductor to electrons in the spin-up states, but as a semiconductor to the spin-down electrons.The VN 2 monolayer exhibits good half-metallicity with a large spin gap of 1.1 eV, which can provide fully spin-polarized currents and large magnetoresistance in some spintronic devices, such as spin-based transistors, diodes, and spin Seebeck devices.For the partial density of states of FeN 2 monolayer in Figure 3a-e, the spin-down states of Fe d suborbitals emerge near the Fermi energy, and the spin-up states of N 2p orbital also exist in the energy range.For the partial density of states of the VN 2 monolayer in Figure 3f-j

Ferromagnetism and Curie Temperature
Finally, we demonstrate that the FeN 2 and VN 2 monolayers are high-temperature FM materials.In Figure 5, the magnetic couplings between two neighbor moments and three magnetic orders are sketched.For clarity, the atomic structure is displayed with the wire frame in Figure 5b-e.To determine the magnetic ground state of the FeN 2 and VN 2 monolayers, the spin-polarized calculations are carried out for three magnetic orders, including the FM order, antiferromagnetic I order (AFM-I), antiferromagnetic II order (AFM-II), and antiferromagnetic III order (AFM-III).The energy per formula cell and the magnetic moment around Fe or V atom are listed in Table 2.In order to present more information on magnetic properties, we illustrate the spin-charge density isosurfaces of FeN 2 and VN 2 monolayers in FM and three AFM orders in Figure 6.The isovalues are set to 0.057 and 0.018 e/Å3 for FeN2 and VN2 layers, and these isosurfaces are roughly the same size.The spin density difference, reflected by the two isovalues, is consistent with the difference of magnetic moment, a large moment of 3.3 µ B around Fe atom and a small moment of 1.0 µ B around V atom.The yellow and cyan isosurfaces indicate up-spin and down-spin density surfaces.In the AFM-I and AFM-II phases in Figure 6f,g, the distribution of spin-charge density is slightly different from the ones in FM and AFM-III phases in Figure 6e,h, which results in a small change around 1.0 µ B in the different magnetic orders.The energy differences among FM, AFM-I, AFM-II orders, and AFM-III orders are assumed to only arise from the magnetic interactions.Based on the energy differences, the nearestneighbor coupling J 1 , next-nearest-neighbor coupling J 2 , and next-next-coupling J 3 can be derived according to the following expressions, The method for calculating exchange couplings comes from the Appendix in the Ref. [29], which is widely used in the study of ground state magnetism in the magnetic materials [30,31].The values of J 1 , J 2 , and J 3 are presented in Table 2.In addition, we compute the single-site magnetic anisotropy energy (MAE) in terms of the definition, MAE = E 100 − E 001 , and list the values in Table 2.The negative MAE value means the easy magnetic axis is in the ab plane, namely an in-plane magnetic anisotropy.Table 2.The energies in FM, AFM-I, AFM-II, and AFM-III orders, local magnetic moment, exchange coupling J 1 (meV/S 2 ), J 2 (meV/S 2 ), and J 3 (meV/S 2 ), MAE(µeV/S 2 ) and T c (K) for the monolayer FeN 2 and VN 2 , respectively.In previous studies, the 2D Ising model and the Heisenberg model were usually employed to estimate the critical temperature for magnetic phase transition [9,32,33].For the Ising model in a triangular spin-lattice, the critical temperature is T c = J 1 • 4/ln3 [34].In terms of the formula, the Curie temperatures of FeN 2 and VN 2 monolayers are 830 K and 450 K, respectively, which are overestimated.Ising model is the limit of the Heisenberg model.That is to say, with magnetic anisotropy going to infinity, the Heisenberg model would be reduced to the Ising model.Because of the small magnetic anisotropy in real 2D FM materials, the T c is usually overestimated by the Ising model.Therefore, Heisenberg model is a more precise model and has been successfully used to estimate the T c of CrI 3 and other synthesized 2D FM materials [35,36].In this paper, we solve the Heisenberg model with Monte Carlo method to evaluate the critical temperature.The Hamiltonian of the Heisenberg model is defined as (3) in which j, j ′ , and j ′′ represent the nearest, next-nearest, and third-nearest neighboring sites of the i site in the triangle lattice and α is a coordinate component x, y, or z.A is the single-site magnetic anisotropic energy, which is the energy difference when the spin is along the (1 0 0) and (0 0 1) directions.The variation of average magnetic moment and the susceptibility with temperature are presented in Figure 7a,b, where the susceptibility . For the FeN 2 monolayer and VN 2 monolayers, the Curie temperatures are 222 K and 238 K, respectively.Among 44 kinds of ferromagnetic 2D materials reported in Ref. [35], the highest Curie temperature predicted by the Heisenberg model is 261 K, which indicates that the FeN 2 monolayer and the VN 2 monolayer are high-temperature ferromagnetic 2D materials.

Computational Methods
The calculations are implemented in the Vienna Ab-initio Simulation Package (VASP) [37,38], and the generalized gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof (PBE) and the projected augmented wave (PAW) pseudopotential are used [39,40].The GGA+U method [41] is also used to improve the description of the electronic correlation of the Fe and V 3d electrons.The criterions of 10 −6 eV and 0.01 eV/Å are set for the total energy and atomic force convergence.The vacuum layer, plane-wave cutoff energy, and k-point mesh are 20 Å, 600 eV, and 39 × 39 × 1, respectively.When the interval between two k-points is set to 0.01 × 2π•Å −1 in the Brillouin Zone, the number of k-points is 39 along a or b axis.The phonon dispersion is calculated by the density functional perturbation theory (DFPT) method as implemented in the Phonopy code [42].The first-principles molecular dynamics simulation in the mole volume-temperature (NVT) ensemble of a 5 × 5 × 1 supercell lasts for 5 ps with a time step of 1 fs at 300 K [43].
The temperature of the phase transition is evaluated by solving the Heisenberg model with the Monte Carlo method [44].Monte Carlo simulations are performed with a 50 × 50 × 1 spin-lattice.For a given temperature, the spin system evolves into a thermal equilibrium state, and the statistical observable such as energy and magnetic moment are directly obtained and other physical quantities can be derived.

Discussion and Conclusions
It is noted that the space group of VN 2 and FeN 2 monolayers is P-6m2, which is a non-centrosymmetric group.In the last decade, topological spin textures such as skyrmions have attracted extensive attention because they have promising applications in novel spintronic devices.Some non-centrosymmetric layered materials have been demonstrated to possess Dzyaloshinskii-Moriya (DM) interaction, which can stabilize the topological spin textures [45,46].For example, Cr-intercalated hexagonal TaS 2 with the non-centrosymmetry of P6 3 22 space group and Fe 2.8 GeTe 2 with the non-centrosymmetry of P6 3 mc group are recently reported.The VN 2 and FeN 2 monolayers predicted in this work have similar hexagonal lattices and non-centrosymmetry.Hence, the topological spin textures can be expected to emerge in the MN 2 (M = V, Fe, and other 3d metals) monolayers.
In summary, by the first-principles calculations, we propose two new 2D transition metal nitrides MN 2 (M = Fe and V).The stability has been demonstrated by the calculations of phonon spectra, molecular dynamics simulation, and formation energy.The electronic structures showed that VN 2 monolayers is FM half-metals with a large half-metallic gap.More importantly, the high-temperature ferromagnetism in the FeN 2 and VN 2 monolayers are determined through solving the Heisenberg model using the Monte Carlo simulation method.Our work provides a new route to designing two-dimensional materials with intrinsic magnetism.

Figure 1 .
Figure 1.Top view, oblique view, and side view of the MN 2 (M = V and Fe) structures in the hexagonal (a), trigonal (b), and tetragonal (c) phases.

Figure 2 .
Figure 2. (a,c), Phonon spectra of FeN 2 and VN 2 monolayers in the ferromagnetic ordering phase.(b,d), Total energy evolution with respect to time at 300 K in the molecular dynamics simulations of FeN 2 and VN 2 monolayers.The insets are the top and side views of their final configurations after 5 ps molecular dynamics simulations at the temperatures of 300 K.
, there are only the spin-up states of V d suborbitals and N 2p states near the Fermi energy.This results in the half-metallicity appearing in the VN 2 monolayer instead of the FeN 2 monolayer.

Figure 3 .
Figure 3. Partial density of states of Fe 3d, V 3d, and N 2p suborbitals in the FeN 2 and VN 2 monolayers.(a-c) d yz /d xz , d xy /d x 2 −y 2 , and d z 2 suborbitals of Fe atom in the FeN 2 monolayer.(d,e) p x /p y and p z suborbitals of N atom in the FeN 2 monolayer.(f-h) d yz /d xz , d xy /d x 2 −y 2 , and d z 2 of V atom in VN 2 monolayer.(i,j) p x /p y and p z suborbitals of N atom in the VN 2 monolayer.The positive and negative values correpond to the spin-up and spin-down channels.

Figure 4 .
Figure 4. Band structure of FeN 2 (a) and VN 2 (b) monolayers.Spin-up and spin-down bands are displayed as blue and red lines, respectively.

Figure 5 .
Figure 5. (a) The nearest, next-nearest, and next-next-nearest neighbored exchange interaction J 1 , J 2 , and J 3 , marked in the MN 2 (M = Fe or V) layer with the double headed arrows.(b-e) FM, AFM-I, AFM-II, and AFM-III orders.The atomic structures are displayed with the wire frame in (b-e).

Figure 7 .
Figure 7. Average magnetic moment M and susceptibility χ as functions of temperature.The vertical dashed lines indicate the phase transition temperature.