Magnetoelastic and Magnetoelectric Coupling in Two-Dimensional Nitride MXenes: A Density Functional Theory Study

Two-dimensional multiferroic (2D) materials have garnered significant attention due to their potential in high-density, low-power multistate storage and spintronics applications. MXenes, a class of 2D transition metal carbides and nitrides, were first discovered in 2011, and have become the focus of research in various disciplines. Our study, utilizing first-principles calculations, examines the lattice structures, and electronic and magnetic properties of nitride MXenes with intrinsic band gaps, including V2NF2, V2NO2, Cr2NF2, Mo2NO2, Mo2NF2, and Mn2NO2. These nitride MXenes exhibit orbital ordering, and in some cases the orbital ordering induces magnetoelastic coupling or magnetoelectric coupling. Most notably, Cr2NF2 is a ferroelastic material with a spiral magnetic ordered phase, and the spiral magnetization propagation vector is coupled with the direction of ferroelastic strain. The ferroelectric phase can exist as an excited state in V2NO2, Cr2NF2, and Mo2NF2, with their magnetic order being coupled with polar displacements through orbital ordering. Our results also suggest that similar magnetoelectric coupling effects persist in the Janus MXenes V8N4O7F, Cr8N4F7O, and Mo8N4F7O. Remarkably, different phases of Mo8N4F7O, characterized by orbital ordering rearrangements, can be switched by applying external strain or an external electric field. Overall, our theoretical findings suggest that nitride MXenes hold promise as 2D multiferroic materials.

Considering their rich diversity in structure and composition, MXenes are promising 2D materials to realize multiferroicity. Tahir et al. [19] reported the first observation of ferroelectric and multiferroic properties in a delaminated Ti 3 C 2 T x MXene film. Moreover, several theoretical studies have predicted intrinsic ferroelectricity or magnetism in some MXenes, such as i-MXene (Ta 2/3 Fe 1/3 ) 2 CO 2 as a type-I multiferroic material [20], Hf 2 VC 2 F 2 as a type-II multiferroic material [21], and magnetic Janus MXenes, that can control their magnetic order via an out-of-plane electric field [22]. The recent advances in the synthesis of monoor few-layer MXenes such as molten salt etching methods [23], halogen-based etching methods [24], the intercalation-alloying-expansion-microexplosion mechanism [25], and delamination based on the water freezing and thawing (FAT) method [26] have promised the production of large-lateral-size and high-quality MXenes that go beyond the wellstudied Ti 3 C 2 T x . These new techniques give us the flexibility to modulate the composition of MXenes, including the control of functional group termination, and thus, examine the theoretical prediction of multiferroic properties in some MXenes. Nevertheless, there are only a limited number of reports concerning the fabrication of nitride MXenes, with the examples of V 2 NT x [27], Mo 2 NT x [27], Ti 2 NT x [28], and Ti 4 N 3 T x [29]. These reports primarily emphasize the fabrication processes of nitride MXenes and highlight their enhanced electrical conductivity compared to their carbide MXene counterparts. However, there is a notable absence of experimental studies regarding their magnetic properties and the possibility of multiferroicity within them.
Ab initio density functional theory (DFT) simulations can facilitate further exploration of 2D multiferroics for studying the electronic, magnetic, and mechanical properties of materials at a quantum mechanical level. Kumar et al. [30] carried out a systematic study on the magnetic properties of the nitride MXene M 2 NT 2 and predicted that high-Curie-temperature ferromagnetism exists in multiple nitride MXenes, where M is an early transition metal (Ti, V, Cr, or Mn), N is Nitrogen, and T is a surface terminal group (O, OH, or F). Their results also hinted at the presence of orbital ordering in magnetic nitride MXenes, but did not fully elucidate its nature. Orbital ordering [31] is a common occurrence in oxides and plays important roles in many fascinating physical phenomena such as superconductivity [32][33][34][35], colossal magnetoresistance [36][37][38], and ferroelectric polarization [39]. Furthermore, it has been shown that electric fields can be used to control orbital ordering, Jahn-Tellar distortions [40], and the magnetic anisotropy of multiferroic PTO/LTO superlattices [41].
In this work, DFT simulations were used to investigate the structural, electronic, and magnetic properties of nitride MXenes M 2 NT 2 , with transition metals (M = V, Cr, Mn, or Mo), selected for their magnetic properties, and different surface terminations (T = O or F) chosen as functional groups to passivate the reactive transition metal surface. Our results show that the nitride MXenes V 2 NF 2 , V 2 NO 2 , Cr 2 NF 2 , Mo 2 NO 2 , Mo 2 NF 2 , and Mn 2 NO 2 are semiconductors, with band gaps ranging from 0.20 eV to 2.1 eV, and exhibit intrinsic magnetism and orbital ordering. We found out that orbital ordering plays an important role in mediating magnetoelastic coupling in Cr 2 NF 2 and magnetoelectric coupling in V 2 NO 2 and the Janus MXene Mo 8 N 4 F 7 O.

Density Functional Theory
Density functional theory (DFT) calculations were performed using the Vienna ab initio simulation package (VASP) [42] with projector-augmented wave (PAW) pseudopotentials [43] and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [44]. Electron correlations in 3d and 4d electrons were treated using the rotational invariant on-site Coulomb parameters in Dudarev's formulation [45]. Herein, the effective U values for V [46,47], Cr [46,48], Mn [46], and Mo [49][50][51] atoms are set to be 3 eV. A vacuum space of 20Å was used to prevent the interactions between the adjacent images along the vertical direction. Γ-centered 12 × 12 × 1 Monkhorst-Park k-point meshes were employed for structural optimizations and calculations of electronic properties, with the cutoff energy of the plane-wave basis set to 600 eV. To determine the ground state configuration, we performed structural optimization on the in-plane lattice constants and internal atomic coordinates of a 2 × 2 supercell with 5 different initial magnetic configurations, including four antiferromagnetic phases (AFM1, AFM2, AFM3, and AFM4 in Figure 1) and a ferromagnetic phase (FM), until the Hellman-Feynman force on each atom was less than 5 × 10 −3 eV/Å. The nudged elastic band method [52] was used to determine the paths and energy barriers of ferroelectric switching, and ferroelectric polarization was calculated using the Berry phase method [53]. The ISOTROPY tool [54] was used to aid with the symmetry and space group analysis.

LKAG Formalism and Monte Carlo Simulations
To calculate the isotropic exchange parameters, J ij , of two-dimensional nitride MXenes, we first derive the tight binding Hamiltonian based on the maximally localized Wannier function (MLWF) with the Wannier90 package [55,56], and then calculate the J ij by treating local spin rotation as a perturbation, using LKAG formalism [57] as implemented in the TB2J package [58]. These values of J ij were later used as the input parameters for the Monte Carlo (MC) simulations to verify the magnetic ground states and estimate the transition temperatures, using the Metropolis algorithm as implemented in the code SPIRIT [59]. The simulations adopted the Heisenberg spin model with a magnetic anisotropy term: with classical spin S i,j of length one at atomic sites i, j; single-ion magnetic anisotropy constant, A; component of the spin along the easy-magnetization axis, S z i . A 2D supercell of 50 × 100 (20,000 magnetic atoms) with a periodic boundary condition was adopted for the simulations, and for each temperature 2 × 10 5 steps were taken for thermalization and another 2 × 10 5 steps were performed for sampling.

Orbital Ordering, Electronic, and Magnetic Properties
The crystal symmetry analysis using ISOTROPY revealed that the space groups of the nitride MXenes V 2 NO 2 , CrNF 2 , MoNF 2 , and MnNO 2 are monoclinic P2 1 /m (Figure 2b), while the space groups of V 2 NF 2 and MoNO 2 are triclinic P1 (Figure 2c). Group theory analysis [60] using ISODISTORT shows that the distortion of the high-symmetry parent structure of space group P3m1 (Figure 2a) into the P2 1 /m phase can be decomposed into three modes: the interplane antipolar displacements of T along the c-axis, with irreducible representation Γ + 1 (a, 0, 0) ( Figure 2d); the strain and interplane antipolar displacements of M and T along the a-axis, with irreducible representation Γ + 3 (a, (a) . The algorithm used to obtain these irreducible representations is summarized in Dorian and Harolds' work [60].
The distortions of the P1 structure can be decomposed into four modes: (Figure 2g), and M − 1 (0, a, 0) mode ( Figure 2h). The Γ + 1 (a, 0, 0) mode and M − 2 (0, a, 0) mode are identical to their counterparts in the distortion of the P1 structure, while the Γ + 3 (0, a, 0) mode's antipolar displacements and strain are along the a -axis and this strain along the a -axis deforms the structure to a non-orthogonal cell. Other than that, there are additional intraplane antipolar displacements of M, N, and T along the b-axis with irreducible representation M − 1 (0, a, 0). The M − 1 (0, a, 0) intraplane displacements cause every metal-ligand bond in the octahedral complex to have different lengths.
In order to better understand the electronic properties of the nitride MXenes, the ligand's p-orbitals and the transition metal's d-orbitals' resolved band structures were plotted along high-symmetry paths in the Brillouin zone (Figure 4a-f). The results revealed that V 2 NF 2 , V 2 NO 2 , Cr 2 NF 2 , Mo 2 NO 2 , and Mo 2 NF 2 exhibited the characteristics of a Mott transition, i.e., the valence band maximum (VBM) and conduction band minimum (CBM) were found to be primarily comprised of d electron bands. The band gap of these Mott insulating nitride MXenes is approximately 1 eV if the t 2g orbitals are partially occupied, and approximately 2 eV for those with fully occupied t 2g orbitals. In contrast, the band structure of Mn 2 NO 2 shows a pd charge-transfer type of band gap, with the VBM being a mixture of e g orbitals from Mn and p orbitals from O, and the CBM dominated by e g and t 2g orbitals. Next, we analyze the nature of the orbital ordering by computing the magnetic properties. Table 1 summarizes the magnetic moments of metal ions (in nitride MXenes) at sublattice 1 and sublattice 2 sites, along with their ground state magnetic orders based on DFT calculations. The magnetic moments of the metal ions are calculated using Wannier functions based on the Wannier90 results of their respective magnetic ground states; as expected, the magnitudes of the magnetic moments of the metal ions are larger (smaller) at sublattice 1 (sublattice 2) as the number of occupied d orbitals is higher (lower). Since these transition metal complexes are in distorted octahedral or distorted tetrahedral symmetry, the d orbitals cannot be classified as pure e g or t 2g orbitals, but rather a mixture of the two. For instance, the d-orbital occupancy of Cr 2+ at the sublattice 1 of Cr 2 NF 2 (Figure 3a) was found to be about three electrons per site, predominantly comprising a mixture of d xy , d yz , and d zx orbitals (Figure 3c). On the other hand, for Cr 3+ at sublattice 2 (Figure 3b), the degeneracy of e g orbitals was lifted, with the d z 2 orbital largely occupying the CBM; while the electronic states in the energy range of approximately 0.5 to 1 eV below the Fermi level were mainly occupied by d xy , d zx , and d yz orbitals (Figure 3d).
majority-spin minority-spin t 2g orbital of metals e g orbital of metals s, p orbitals of ligands    (Figure 5a), with the magnetization at the atomic site k of sublattice i in layer j, with R k(i,j) approximated by (20) • and q ≈ ( π a , π b , 0), with a and b representing the lengths of the lattice vector along the − → a -and − → b -axes, respectively. The initial phase of magnetization, φ i,j is summarized in Figure 5c. It is worth noting that the long-range magnetic ordering in V 2 NF 2 , Mo 2 NO 2 , and Mo 2 NF 2 is predicted to prevail above room temperature. Such room temperature magnetic behavior in nitride MXenes was also predicted by Kumar et al. [30], where Cr 2 NF 2 and Mn 2 NT 2 (T = F, O, and OH) are reported to show ferromagnetism, with Curie temperatures well above room temperature; pristine and La-doped Ti 3 C 2 MXenes are reported to have ferromagnetism and antiferromagnetism coexisting at temperatures as high as 300 K [61]. Currently, most of the experimentally discovered 2D ferro-/antiferromagnets have their Curie/Néel temperatures limited below room temperature [2,3,[5][6][7][8][9][10][11][12][13], with only VSe 2 [4] and MnSe x [7] showing strong ferromagnetism above 300 K. Our results and previous studies suggest that the MXenes are an important family of 2D materials to search for a new generation of room temperature ferro-/antiferromagnets.  Figure 5. Spin textures of (a) Mn 2 NO 2 at 0 K and (b) Cr 2 NF 2 at 0 K generated using a Monte Carlo simulation in a 50 × 100 lattice with periodic boundaries. The color map is the magnitude of spin S i polarized along the out-of-plane axis; the green (black) arrows in sublattice 1 (2) represent the in-plane components of spin S i ; the cyan and blue spheres correspond to metal ions in sublattices 1 and 2, respectively. Spin components of (c) Mn 2 NO 2 projected to a-b plane and (d) Cr 2 NF 2 projected to a-c plane; φ i,j represents the initial phase of the spin components in sublattice i in layer j.
For Cr 2 NF 2 , the Monte Carlo simulation at 0 K predicted a magnetic configuration (Figure 5b) that is similar to an out-of-plane Néel-type spin spiral, with the magnetization M k(i,j) = M i (cos(q · R k(i,j) + φ i,j ), 0, sin(q · R k(i,j) + φ i,j )), where q = (0.2 × 2π a , 0, 0), where a represents the length of the lattice vector along the − → a -axis. The formation of an intralayer spin spiral magnetization in Cr 2 NF 2 is induced by frustrated isotropic couplings between Cr atoms, which can be described by a simplified model involving four values of the interatomic exchange coupling parameter J (Figure 6a-d). The strong antiferromagnetic couplings (J 1 = −55 meV and J 2 = −52 meV) between interlayer Cr atoms explain three key features of the Monte Carlo simulation results: (i) nearly antiparallel magnetization between Cr atoms from successive layers; (ii) nearly parallel magnetization between intralayer Cr atoms along the b-axis; and (iii) a small initial phase difference (∆φ ≈ 30 • ) (Figure 5d) between intralayer Cr atoms within the same unit cell. The frustrated isotropic couplings are attributed to the competing interactions between intralayer ferromagnetic coupling J 3 (18 meV) and interlayer ferromagnetic coupling J 4 (28 meV), which favor ferromagnetic and antiferromagnetic couplings of intralayer Cr atoms along the a-axis and b-axis, respectively. Since magnetization of intralayer Cr atoms within the same unit cell are highly in phase, and the wavevector q is parallel to the a-axis, the value of q can be approximated using the classical J N -J NN model for one-dimensional frustrated ferromagnets, i.e., q a = |J N |/(4J NN ) [62]. Using J NN = 2 × 18 meV and J NN = 28 meV, we obtain q a = 0.19, which is similar to the prediction of the Monte Carlo simulation (q a = 0.20). Interestingly, although the interacting ion pairs of J 5 are separated by comparable distances with similar atomic environments when compared with those of J 1 and J 2 , J 5 has a much weaker antiferromagnetic coupling (−4.4 meV) relative to J 1 and J 2 . This unusual weak coupling is caused by orbital ordering in Cr 2 NF 2 , and is explained in detail in the next paragraph.
To understand the origin of antiferromagnetic and ferromagnetic interactions among Cr ions, we explain the magnetic interactions using two processes: the antiferromagnetic direct exchange process and the ferromagnetic double exchange process (Figure 6a-d). For J 1 and J 2 , the interatomic distances between interacting ions are 2.772 Å and 3.149 Å, respectively, and the major contribution to the antiferromagnetic coupling originates from the direct exchange process between the Cr ions; the contribution from double exchange to J 1 is negligible since the overlap between p-ligands and d-metals is lacking due to the non-occupancy of d x 2 −y 2 orbitals in metal ions; on the other hand, the distance between the connecting ligand and Cr ions in J 2 coupling is too large, making double exchange processes unfavorable.
For the ferromagnetic J 3 interaction, a double exchange process arises as the electron hops between the occupied d yz orbital in Cr 3+ and unoccupied d z 2 orbital in Cr 2+ , mediated by a non-magnetic N ion. While for the double exchange process in J 4 's interaction, an electron from the occupied d z 2 in Cr 2+ couples with the unoccupied d z 2 in Cr 3+ , through a non-magnetic N ion. Finally, a strong direct exchange process is expected in J 5 , as the separation between the interacting ions is only 2.946 Å. Nevertheless, a ferromagnetic double exchange process also arises from the electron coupling between the occupied d yz orbital in Cr 3+ and unoccupied d z 2 orbital in Cr 3+ at the opposite site. The competing interaction between the double exchange process and direct exchange process results in the weak coupling strength of J 5 (−4.4 meV).

Magnetoelastic Coupling in Cr 2 NF 2
The material Cr 2 NF 2 possesses a significant ferroelastic strain (10.9%), i.e., it has three equally stable orientation variants in the lattice structure that can be switched by applying an external stress. The ferroelastic strain is defined as (( b a − 1) × 100%), with lattice parameters a = 6.03 Å and b = 6.689 Å, calculated along the three diagonal lines ( a 1 , a 2 , and a 3 ) of a deformed hexagon (Figure 7a). The shorter lattice parameter, a, is laying along a 1 , and the propagation vector of the spin spiral magnetization is perpendicular with respect to a 1 , showing coupling between the ferroelastic strain and magnetic order. In addition to the previously mentioned orientation state (Figure 7a), there are two equivalent orientation states that can be achieved in Cr 2 NF 2 , by switching the values of a 2 (Figure 7b) or a 3 (Figure 7c) to 6.03 Å. The transitions between the two orientation states are analyzed using the nudged elastic band (NEB) method and result in an activation energy barrier of 56 meV per atom. In comparison to other reported 2D ferroelastic materials, Cr 2 NF 2 has a low ferroelastic strain relative to BP 5 (41.4%) [63], borophane (42%) [64], and phosphorene (37.9%) [65], but is similar to AgF 2 (13.5%) [66] and t-YN (14.4%) [67]. The activation energy barrier of Cr 2 NF 2 is also moderate, comparable to AgF 2 (51 meV) [66] and t-YN (33 meV) [67], but lower than that of BP 5 (320 meV) [63], borophane (100 meV) [64], and phosphorene (200 meV) [65]. These results suggest that the moderate activation energy barrier of Cr 2 NF 2 allows for experimental manipulation of its ferroelasticity while maintaining stability, which could facilitate the control of its spin spiral magnetization.

Orbital-Ordering-Induced Magnetoelectric Coupling
It is noteworthy that structural optimization has revealed the existence of stabilized ferroelectric (FE) phases in the materials V 2 NO 2 , Cr 2 NF 2 , and Mo 2 NF 2 . These FE phases are characterized by energies that are 31 meV, 70 meV, and 47 meV higher than the energies of their corresponding ground states, respectively.
As an illustration, we use V 2 NO 2 to highlight the role of orbital ordering in magnetoelectric coupling among the aforementioned MXenes. Figure 8a,b present contour plots of the spin density of V 2 NO 2 in the antiferroelectric (AFE) phase (ground state) and FE phase (excited state), respectively. Both phases exhibit AFM3 orders. The transition from the AFE phase to the FE phase results in a significant increase in the average metal-ligand bond length of the V L ion (circled in black in Figure 8a,b) in the lower layer and a decrease in the average metal-ligand bond length of the V U ion (circled in red in Figure 8a,b) in the upper layer. This in turn leads to the transformation of V L and V U ions from V 4+ and V 3+ to V 3+ and V 4+ , respectively. As a consequence, the distribution of V 4+ and V 3+ ions becomes asymmetrical, with the ratio of V 4+ to V 3+ being 1:3 and 3:1 in the upper and lower layers of V 2 NO 2 , respectively.
To further verify the stability of the AFM3 order of V 2 NO 2 , we compared the energies of the FE phase at different magnetic orders and performed Monte Carlo simulations. The results confirmed that the AFM3 order remains the most stable magnetic configuration. The energy required to switch from the AFE to the FE phase is estimated to be 51 meV per V atom using the nudged elastic band (NEB) method. The Berry phase method calculates the polarization along the a-and c-axes to be 6.9 µC/cm 2 and 2.3 µC/cm 2 , respectively, expressed in bulk form by considering its thickness. Since the magnetic moment of V 3+ (≈2 µ B ) is larger than that of V 4+ (≈1 µ B ), the asymmetrical distribution of V 4+ and V 3+ ions in the FE phase results in a net average magnetic moment of 0.25 µ B per V atom, demonstrating the coupling between the electrical and magnetic degrees of freedom.
Nevertheless, the excitation energy (32 meV) of the ferroelectric phase of V 2 NO 2 is likely too high to be stabilized experimentally. One potential method to induce spontaneous ferroelectric polarization is to fabricate Janus MXenes, which are MXene structures with different compositions of the upper and lower functional groups. The absence of inversion symmetry in Janus MXenes leads to a large dipole moment. However, this type of ferroelectric polarization usually cannot be switched by an external electric field, which limits its application in nanoelectronic devices.
In this study, we demonstrate that ferroelectric polarization and magnetization in the Janus MXene Mo 8 N 4 F 7 O can be altered by an external electric field or external strain. We found out that the Janus MXenes V 8 N 4 O 7 F, Cr 8 N 4 F 7 O, and Mo 8 N 4 F 7 O have excited states characterized by orbital ordering rearrangements, and the excitation energies of these Janus MXenes (18.5 meV, 18.6 meV, and 6.6 meV, respectively) are much lower than those of their non-Janus counterparts, i.e., V 2 NO 2 , Cr 2 NF 2 , and Mo 2 NF 2 . Interesting results are found in the transition of Mo 8 N 4 F 7 O from the ground state phase (P1 phase) to the excited phase (N1 phase). Our results reveal that Mo 8 N 4 F 7 O is a bipolar magnetic semiconductor (BMS) in both the P1 and N1 phases. BMSs are characterized by a VBM and CBM that are fully spin polarized in opposite spin channels (Figure 9a,b), offering an ideal platform to manipulate the polarization of spin current using gate voltage [52,68,69]. In the P1 (N1) phase, the VBM is dominated by electronic contributions from lower (upper)-layer Mo ions, while the CBM is dominated by electronic contributions from upper (lower)-layer Mo ions, with a band gap of 1.4 (0.4) eV. The transition from the P1 to N1 phase also results in an increase in the ferroelastic strain, ( b a − 1) × 100% (Figure 9c), from 6.0% to 13.1%, with a P1 = 6.57 Å and b P1 = 6.96 Å, while a N1 = 6.38 Å and b N1 = 7.22 Å. By linearly interpolating a and b between (a P1 ,b P1 ) and (a N1 ,b N1 ), we are able to further reduce the excitation energy until the N1 and P1 phases become nearly energetically degenerate (around a ferroelastic strain of 10.2%). To illustrate the orbital ordering rearrangements upon the transition, we show the (contour plot of) spin density of the P1 phase ( Figure 9c) and highlight the changes upon transition (Figure 9d,e). In the P1 phase, the Mo atoms in the lower layer (75% F and 25% O as surface functional groups) have a higher average magnetic moment than the Mo atoms in the upper layer (with only F as a surface functional group), thus give a net magnetic moment of +0.25 µ B per Mo. Upon transition to the N1 phase, the bond length between Mo U (circled in red in Figure 9c) and the nearest N ion along the negative z-direction (of local coordinate system for d-orbitals) increases significantly, thereby increasing the occupancy of the d z 2 orbital and magnetic moment at the Mo U site; since the N ion now receives less electrons from the Mo U site, Mo L ions (circled in black in Figure 9c) donate additional electrons to the N ion, resulting in a slight decrease in occupancy of d xy orbitals and the magnetic moments at the Mo L site. As a result, the transition from the P1 phase to N1 phase changes the average magnetic moment from +0.25 µ B per Mo to −0.25 µ B per Mo.
Because of the broken inversion symmetry in Mo 8 N 4 F 7 O and the difference in electronegativity of O and F, we expect an intrinsic out-of-plane electric field directed from the lower layer to the upper layer. This prediction is confirmed by our calculations of electronic energies under an applied electric field (Figure 9f): both P1 and N1 phases have their electronic energies decrease (increase), with the applied electric field directing upward (downward). However, the electronic energy of the N1 phase decreases more significantly with an upward applied electric field, and increases less significantly with a downward electric field, when compared with the energy of the P1 phase under the same conditions, implying that the overall strength of the intrinsic electric field is higher in the N1 phase relative to the P1 phase. From Figure 9c, we can see that in the P1 phase, the total ionic charge of Mo ions from the upper layer is higher than from those in the lower layer; such an asymmetrical distribution of low and high ionic charge on the Mo ions results in an intrinsic out-of-plane electric field that is opposite to the intrinsic electric field due to the asymmetrical surface functional groups, thus reducing the overall strength of the intrinsic electric field. On the other hand, the intrinsic electric field, due to both the asymmetrical distribution of low/high-ionic-charge Mo ions and of surface functional groups, is pointing upwards, thus enhancing the overall strength of the intrinsic electric field. The results show that at a ferroelastic strain of 10.2%, the energies of Mo 8 N 4 F 7 O in the P1 and N1 phases become degenerate, where switching of the phase can be achieved by application of an external electric field, making non-volatile electric control of spin polarization possible. It should be noted that the origin of magnetoelectric coupling in our study is due to orbital ordering, which is similar to the previously reported multiferroicity in LuFe 2 O 4 [70]. In LuFe 2 O 4 , the non-centrosymmetric distribution of Fe 2+ and Fe 3+ would give rise to spontaneous ferroelectricity and ferrimagnetism. To quantify the coupling between the magnetic order, electronic polarization, ionic displacement, and lattice distortion, a detailed study using a Landau-like model [71,72] is required in the future.

Conclusions
In conclusion, using first-principles calculations, we have shown that multiple nitride MXenes (V 2 NO 2 , V 2 NF 2 , Cr 2 NF 2 , Mo 2 NO 2 , Mo 2 NF 2 , Mn 2 NO 2 ) are magnetic semiconductors, with band gaps ranging from 0.20 eV to 2.1 eV, and exhibit orbital ordering. Among them, Cr 2 NF 2 exhibits spin spiral magnetization induced by orbital ordering, where the propagation vector is coupled with the ferroelastic strain of the structure, suggesting the magnetic order can be switched by applying external stress. We also find that the Janus MXene Mo 8 N 4 F 7 O exhibits two distinct ferroelectric phases marked by different orbital orderings, where the energy difference between the two phases can be reduced by applying an external strain or external electric field. Remarkably, the rearrangement of the orbital ordering during transition would also reverse the direction of net magnetism, illustrating the coupling between the strain tensor, electrical polarization, and magnetization. Our results show that orbital ordering can play an important role in coupling the electronic, magnetic, and mechanical properties of 2D nitride MXenes, making nitride MXenes important candidates for potential application in spintronic devices.