Field-free coexistence of skyrmions and anti-skyrmions induced by higher-order interactions and biaxial strain in the NiI₂ monolayer

Abstract

The spin spiral (SS) state in monolayer (ML) NiI₂ presents a promising avenue for the exploration of two-dimensional multiferroics, but the mechanisms underlying this SS state and its effective modulation remain insufficiently understood. In this study, we employ first-principles calculations to reveal the major features and physical mechanisms governing the magnetism of ML-NiI₂ under in-plane strain, explicitly considering the effects of higher-order interactions (HOIs), which have been overlooked in prior research. Our findings identify the magnetic ground state of ML-NiI₂ as an SS state characterized by a propagation vector q = 0.23 and demonstrate that in-plane biaxial strain serves as an effective method for tuning the period and stability of the SS state, the frustrated ratio, magnetic anisotropy energy and the strength of HOIs. Notably, the strength of HOIs significantly increases under higher strain, particularly for the four-site-four-spin exchange interaction, which is substantial and crucial to the properties of ML-NiI₂. Additionally, through atomistic spin dynamics simulations, the coexistence of skyrmions and anti-skyrmions with diameters of at least 2 nm is achieved in ML-NiI₂ via in-plane biaxial strain, without the necessity for an external magnetic field, which is more conducive to spintronic applications. The presented results establish ML-NiI₂ as a prospective candidate for next-generation spintronic devices and contribute to a deeper understanding of HOIs in two-dimensional magnets.

---

New concepts

Our work demonstrates the field-free stabilization of 2-nm skyrmions and anti-skyrmions in a centrosymmetric monolayer magnet (NiI₂) through synergistic strain engineering and higher-order interactions (HOIs). In stark contrast to conventional skyrmions stabilized by Dzyaloshinskii–Moriya interaction (DMI) in non-centrosymmetric systems, the stabilization here is achieved in the absence of DMI and is dominated by frustrated exchange and four-site-four-spin interactions (F). Critically, biaxial strain (8%) amplifies HOIs by 200% and tunes magnetic frustration, breaking prior limitations where such textures required large DMI or micron-scale sizes. This contrasts with existing studies on centrosymmetric magnets, which overlooked HOIs and achieved only larger (>5 nm) or field-dependent skyrmions. Our findings reveal that strain-activated HOIs mimic DMI-like canting, a mechanism generalizable to other 2D frustrated magnets. This introduces a paradigm for designing ultra-dense spintronic materials by leveraging intrinsic spin interactions and strain, rather than relying on structural asymmetry or external fields. This work advances materials science by (1) establishing HOIs as critical drivers of nanoscale topology in centrosymmetric systems, (2) providing a strain-mediated pathway to stabilize and miniaturize spin textures, and (3) bridging the gap between frustrated magnetism and practical device engineering.

1. Introduction

Since the first experimental observation of magnetic skyrmion lattices in 2009,¹ noncollinear magnetic spin textures have attracted considerable interest due to their unique physical properties and potential applications in spintronics. Among them, magnetic skyrmions, as noncollinear magnetic vortex structures, exhibit remarkable stability against perturbations and can be driven by ultralow currents, making them promising candidates for high-density information storage and processing devices.^2–9

While most early studies focused on skyrmions stabilized by Dzyaloshinskii–Moriya interaction (DMI) in non-centrosymmetric systems,^10–12 recent attention has shifted toward centrosymmetric magnets, where alternative mechanisms such as frustrated exchange and higher-order interactions can promote topological spin textures.^13,14 In particular, skyrmions originating from symmetric exchange interactions can reach sizes as small as 1 nm, far smaller than those stabilized by DMI or dipolar interactions, offering superior integration density and enhanced topological transport properties.^15–19 Noncollinear magnetism, encompassing spin canting,²⁰ vortices,²¹ cycloidal textures,²² and spin spirals,^23,24 etc., has been an eye-catching field and garnered significant attention for its application potential in innovative nonvolatile spintronic devices with high storage density and low energy dissipation.

Among various centrosymmetric platforms, bulk van der Waals magnet nickel diiodide (NiI₂) has emerged as a key material due to its rich magnetic phase diagram and strong spin–lattice coupling, especially the spin spiral ordering has been observed as the magnetic ground state experimentally.^25–27 Its spiral magnetic ground state, characterized by a proper-screw helimagnetic (HM) state with a propagation vector q ≈ (0.138, 0, 1.457),²⁵ disrupts inversion symmetry and induces spontaneous electric polarization in the plane,²⁸ classifying it as a type-II multiferroic material.

As for the monolayer (ML) NiI₂, its magnetic properties remain under active debate. Early theoretical studies proposed a ferromagnetic (FM) ground state,^29–32 while later works, accounting more fully for spin–orbit coupling and noncollinear configurations, identified anti-ferromagnetic (AFM)³³ and spin spiral phases^34–36 as energetically favorable. Such a spiral phase arises from the magnetic frustration induced by the competition between FM first-neighbor and AFM third-neighbor magnetic exchange interactions.^36,37 Even under a large external magnetic field, Amoroso et al. predicted a spontaneous chiral magnetic anti-biskyrmion lattice in ML-NiI₂.³⁸ However, significant discrepancies persist in the reported propagation vector q for the spin spiral state, with values ranging from ∼0.137 to 0.250 across different theoretical and experimental studies. For instance, Li et al. determined the spiral period to be λ = 7.3 a₀ (|q| = ∼0.137),³⁴ while Fumega et al. analyzed the HM state with a 7 × 1 × 1 supercell, yielding |q| ∼ 0.143,³⁵ both values near the bulk phase of 0.138. In contrast, Sødequist et al. reported |q| to be 0.14 for ML-NiI₂ and noted a shallow saddle point at |q| = 0.21 with 2 meV energy difference.³⁶ Subsequently, Miao et al. experimentally measured the period of ML-NiI₂ as 1.76 nm, corresponding to λ = 4.54 a₀ (|q| = ∼0.220), and confirmed this result through DFT calculations.³⁹ Other reports include propagation vectors of |q| = −0.250 using a 4 × 8 × 1 supercell⁴⁰ and |q| = 0.220.³³ These inconsistencies highlight the need for further investigation into the fundamental properties of monolayer NiI₂, particularly regarding its potential to host skyrmion states in a centrosymmetric setting.

Critically, prior research has already demonstrated that in-plane strain can substantially tune magnetic couplings and phase stability in monolayer NiI₂.⁴⁰ Building on these insights, we conduct a comprehensive first-principles study of strain-modulated magnetism in NiI₂, with particular emphasis on the role of higher-order interactions (HOIs), which have thus far been overlooked. We not only identify a spin spiral ground state with |q| = 0.23 but also show how biaxial strain can selectively enhance HOIs, leading to the field-free stabilization of 2-nm skyrmions and antiskyrmions. The presented results position ML-NiI₂ as a promising candidate for next-generation spintronic devices and provide theoretical guidance for controlling topological magnetism in two-dimensional van der Waals magnets.

2. Calculation methods

Our first-principles calculations are carried out within the framework of density-functional theory (DFT) using the Vienna ab initio simulation package (VASP)⁴¹ with the projector augmented wave (PAW) method.⁴² The exchange–correlation functional is described by the generalized gradient approximation (GGA) within the Perdew–Burke–Ernzerhof (PBE) functional.⁴³ To improve the description of on-site Coulomb interaction between the localized 3d states of Ni, the GGA + U method⁴⁴ is implemented within the Dudarev scheme, with an effective Hubbard-like term U set to be 4.0 eV, in accordance with prior studies.^30,45 The cut-off energy for plane-wave expansions is set to be 550 eV. The atomic position is fully relaxed with the Monkhorst–Pack k-point mesh of 11 × 11 × 1, until the total energy and Feynman–Hellmann force on each atom are less than 1 × 10^-6 eV per atom and 0.01 eV Å⁻¹, respectively. To ensure high accuracy in calculating the single ion anisotropy (SIA), which involves small energy differences between spin orientations, stricter convergence criteria of 1 × 10^-7 eV for total energy and 0.001 eV Å⁻¹ for Hellmann–Feynman forces are applied. A vacuum space of at least 15 Å is applied between the layer and its adjacent periodic images along the out-of-plane axis to avoid the spurious interactions. The phonon dispersion was obtained by using the Phonopy package⁴⁶ based on density functional perturbation theory (DFPT).

3. Results and discussion

The geometric structure of ML-NiI₂, shown in Fig. 1(a), features a Ni atom triangular lattice sandwiched between two layers of iodine atoms, with each Ni cation octahedrally coordinated by six I anions, forming a 1T phase. The optimized lattice constant for ML-NiI₂ is a₀ = 3.986 Å, consistent with previous studies.^30,45 We start by discussing the energy dispersion E(q) of spin spirals in ML-NiI₂. To identify the stable existence of a noncollinear ground state, we calculate the energy difference between spin spiral and FM states as a function of the magnetic propagation vector q using scalar-relativistic approximation (excluding SOC due to its incompatibility with generalized Bloch theorem (gBT)^36,47). As shown in Fig. 1(b), energy minima appear at |q| = 0.23 along the Γ-K direction, while the AFM and FM states have higher energies of 20.3 and 4.1 meV u.c.⁻¹, respectively. A shallow saddle point at |q| = 0.20 along the Γ-M direction has an energy just 0.35 meV above the minima. These results confirm the spin spiral state's stability over FM and AFM states, supporting previous theoretical work by Ni et al.³³

Fig. 1 (a) Top and side views of the atomic structure of ML-NiI₂. The gray and purple spheres represent Ni and I atoms, respectively. The dotted line areas in the top and side views show the primitive cell. (b) Energy dispersion E(q) of flat spin spirals along the high-symmetry directions of the 2DBZ for the ML-NiI₂. (c) The schematic view of the spin-spiral state in ML-NiI₂ (taking q = 0.2 for example, its period is 5 lattices). Only Ni atoms (brown circles) are shown for visual clarity. Red arrows indicate the directions of the noncollinear magnetic moments. (d) Effects of biaxial strain on the spin-spiral state. (e) The relationship between strain and the ground state magnetic propagation vector q, extracted from the black dash curve in (d).

We also examine the influence of compressive and tensile biaxial strains, defined as , where positive values indicate tensile strain. Fig. 1(d) displays the energy differences (E_q–E_FM) under varying strains, revealing that the spin spiral state is highly sensitive to strain. Fig. 1(e) shows how the ground state magnetic propagation vector q changes with strain, indicating that compressive strain increases |q| from 0.23 to 0.29 as ε goes from 0 to −8%, shortening the spin spiral period. Conversely, tensile strain decreases |q| to 0.14 at ε = 6%, lengthening the spin spiral period. This behavior shows an approximately linear relationship between strain (within −8% to 6%) and the ground state |q|. Intriguingly, at ε = 8%, we also find a significant decrease in ground state |q| accompanied by a transition from the spin spiral state to the FM state. We compute the phonon dispersions under ±8% biaxial strain to confirm the mechanical stability of ML-NiI₂, which show no imaginary frequencies (see Fig. S2 in the SI). This confirms the structural integrity of the system at these strain levels. Note that such strain levels as high as 12–13% have been successfully demonstrated in 2D van der Waals materials through advanced techniques.^48–50 Therefore, biaxial strain effectively modulates the spin spiral ground state of ML-NiI₂ potentially expanding its application in spintronic devices.

The easy magnetic axis orientation is a crucial factor in understanding the properties of the magnetic states. We estimate that the single ion anisotropy (SIA) of pristine ML-NiI₂ is 0.23 meV as shown in Fig. 2(a), suggestive of a preference for an out-of-plane spin orientation, consistent with other theoretical studies.^33,34 As tensile strain increases, SIA increases steadily, maintaining an out-of-plane direction. In contrast, under compressive strain, the SIA initially decreases to a minimum of −0.30 meV at ε = −4%, causing a shift from out-of-plane to in-plane orientation. With further compressive strain, SIA sharply increases, resulting in a second switch back to out-of-plane orientation at ε below −6%. Notably, SIA peaks at 0.59 meV under −8% strain. These findings show that the orientation and magnitude of magnetization in ML-NiI₂ can be effectively tuned using biaxial strain.

Fig. 2 (a) The relationship between single-ion anisotropy (SIA) energy and biaxial strain ε in ML-NiI₂. (b) Selective atomic-layer-resolved SIA of ML-NiI₂ as a function of biaxial strain. (c) and (d) Orbital-resolved SIA for I-p orbitals in ML-NiI₂ under ε = −8% and −4%. (e) The relationship between strain and the SIA contributed by each I-p orbital in ML-NiI₂.

To further investigate the sharp increase in SIA under compressive strain, we calculate the element-resolved and orbital-resolved SIA by examining the energy differences of SOC matrix elements.⁵¹ The element-resolved SIA in Fig. 2(b) shows that the Ni atom contributes a stable negative value of about 0.2 meV across the strain range, while the contribution from I atoms substantially enhances under both tensile and large compressive strains, demonstrating a pronounced V-shaped variation that aligns with the overall SIA trend. Consequently, the variation of I atoms exerts a predominant influence on the overall SIA of ML-NiI₂, highlighting that heavier non-magnetic atoms contribute more significantly than magnetic atoms, a trend also seen in MnBi₂Te₄ and Mn₂I₃Br₃.^52,53 We also explore the orbital-resolved SIA of I-p orbitals, where the hybridization of occupied and unoccupied states, as described by second-order perturbation theory, underlies the SIA. As depicted in Fig. 2(c) and (d), the hybridization between I-p_y and I-p_z orbitals contributes positively (∼0.67 meV at ε = 0%) to SIA, while the I-p_xp_y hybridized orbital provides a negative contribution of −0.59 meV at ε = 0%, consistent with previous calculated results.^45,54 The relationship between strain and SIA from each I-p orbital shows that I-p_yp_z orbits maintain stable positive contributions. As the strain varies from 0% to −8%, however, the negative SIA stemming from I-p_xp_y hybridized orbitals initially increases, reaching a peak with a value of about −0.70 meV at ε = −4%. Subsequently, the SIA contributed by I-p_xp_y diminishes to −0.54 meV under a compressive strain of ε = −8%, serving as the primary driver for the great improvement of SIA in ML-NiI₂ under the compressive strain.

Next, we discuss how biaxial strain affects magnetic exchange interactions. Fig. 3(a) illustrates the first, second and third nearest neighbor Heisenberg exchange interactions J₁, J₂ and J₃ as functions of strain, based on the spin Hamiltonian equation as shown in the SI. For strain-free ML-NiI₂, J₁, J₂ and J₃ are −3.94 meV, −0.36 meV and 3.28 meV, respectively, aligning well with other theoretical findings.^33,40 It is observed that both J₁ and J₂ exhibit negative values, indicative of FM coupling. While J₃ is positive, signifying AFM coupling. J₁ remains FM throughout the entire strain range, but its strength diminishes. Specifically, the reduction in J₁ under compressive strain (∼2 meV at ε = −8%) significantly exceeds that under tensile strain (∼0.74 meV at ε = 8%). Meanwhile, the change in J₂ is negligible, remaining close to zero across all strain levels. In contrast, J₃ demonstrates more pronounced tunability with a sharp rise of approximately 160.4% with compressive strain increasing to ε = −8%, reaching 8.54 meV. This leads to magnetic frustration from the strong competition of FM (negative J₁′) and AFM (positive J₃′) interactions, as reflected in frustration ratio |J₃|/J₁, which increases with increasing compressive strain, as shown in Fig. 3(b), indicating a substantial enhancement in the antiferromagnetic coupling of the system.

Fig. 3 The evolutions of (a) Heisenberg exchange interactions J and (b) frustration ratio |J₃|/J₁ in ML-NiI₂ under biaxial strain. J₁′, J₂′ and J₃′ are the modified Heisenberg exchange parameters. (c) The value of bond angles and bond lengths corresponding to the exchange paths between the first three nearest neighbors. The insets of (c) show the possible paths for J₁, J₂ and J₃. (d) The variation of calculated higher-order interactions (HOIs) with strain. The inset of (d) represents the schematic diagram of different nearest neighboring J₁, J₂, J₃ and HOIs. B, Y and F represent biquadratic interaction, three-site-four-spin interaction and four-site-four-spin interaction, respectively.

To understand the microphysical mechanisms behind Heisenberg exchange interactions under strain, we further explore the variations of super exchange and direct exchange interactions. Direct overlap of 3d orbitals between neighbouring Ni atoms can result in AFM-preferred coupling, with its strength highly dependent on the distances between magnetic atoms. In contrast, super-exchange involves indirect interaction via intermediary I-p orbitals, for which the strength can be roughly estimated by the bond angle ∠Ni–I–Ni following the semiempirical Goodenough–Kanamori–Anderson (GKA) rule.^55–57 According to GKA rules, the 90° bond angle tends to exhibit FM coupling through super-exchange interaction. Typically, the competition between AFM-preferred direct exchange and FM-preferred super-exchange determines the final magnetic state of the system. Taking J₁ as an example, in the absence of strain, the approximately 90° of ∠Ni–I–Ni (∼91.97°) ensures an FM state through super exchange. However, the application of strain, whether tensile or compressive, causes the bond angle θ₁ to deviate from 90°, weakening FM strength. Additionally, compressive strain reduces the 〈Ni–Ni〉 distance d₁, substantially enhancing AFM-preferred direct exchange coupling. While in tensile strain, the increase of d₁ decreases the AFM strength of direct exchange, thereby enabling the system to retain FM characteristics to the greatest extent possible. As a result, the FM strength of J₁ under compressive strain decreases more rapidly compared to its tensile counterpart.

The exchange mechanisms of J₂ and J₃ are more complicated than that of J₁, involving multiple exchange pathways described as super-super-exchange (SSE) coupling via the p orbitals of two I anions, forming an extended Ni–I⋯I–Ni SSE interaction. As mentioned by Goodenough,⁵⁸ these extended super exchange interaction follows the same principles as the conventional super exchange interaction. Fig. 3(c) illustrates that each SSE path of J₃ features two Ni–I⋯I bond angles of approximately 135° (ranging from 133.1° to 138.4°), resulting in J₃ being expected to exhibit robust AFM coupling. Conversely, all possible SSE paths of J₂ involve Ni–I⋯I bond angles approaching 90° (θ₂ and θ₄), making J₂ significantly weaker than J₃. Moreover, both I anions in all SSE paths of J₃ lie in the same plane, enhancing I–I hybridization compared to J₂, thus making J₃ considerably larger in magnitude than J₁ and J₂. The negative integrated crystal Hamilton population (–ICOHP), calculated using the LOBSTER code,⁵⁹ further analyses the super exchange interaction intensity, where higher –ICOHP values signify a stronger interaction. For J₃ without strain, the average −ICOHP is 0.379, increasing to 0.650 under ε = −8% strain, and dropping to −0.032 under ε = 8% strain, showing a trend consistent with J₃ (see the SI for more details). Moreover, to ensure that the applied biaxial strain does not accidentally induce a breaking of inversion symmetry, we calculate the DMI for the strained systems as well. Our results confirm that the DMI remains zero under all strain conditions considered in this study, preserving the centrosymmetric character of the monolayer (see the SI for more details).

Beyond classic Heisenberg interaction, higher-order interactions (HOIs) from multiple-electron hopping across 2–4 atomic sites can also influence the magnetic properties of system. Previous studies have often overlooked HOIs, but it has been shown that the higher-order biquadratic exchange interaction B has a strong impact on magnetic order of 2D magnets,^33,60,61 and it is suggested to be non-negligible in low temperature magnets. For ML-NiI₂, the calculated B value is 0.52 meV, close to the value of 0.59 meV reported in ref. 33 and comparable to the Heisenberg exchange J₂. Except for B, more complicated HOIs like three-site-four-spin exchange Y and four-site-four-spin exchange F were also considered through multi-Q states calculations (details are given in the SI). The calculated parameters for HOIs are illustrated in Fig. 3(d). In the strain-free scenario, values of Y and F are −0.01 and −0.56 meV, respectively. Under strain, B and Y show slight variations, while F changes significantly. Incorporating HOIs slightly modifies the Heisenberg exchange interaction, denoted as J₁′, J₂′ and J₃′ (see the SI for more details). As shown in Fig. 3(a), the modified parameters display minimal deviation from the originals, and their trend under strain is basically the same.

It should be noted that additional factors, like the anisotropic Kitaev interaction, have been proposed to influence the ground state in few-layer NiI₂ films;^34,62 however, recent theoretical studies have shown that its strength is significantly reduced in the monolayer, being approximately three times weaker than in the bilayer and bulk systems.⁶³ Furthermore, atomic-scale experiments⁶⁴ and theoretical analyses⁶⁵ concur that a minimal classical Heisenberg J₁–J₃ model, rather than the Kitaev interaction, dominates and adequately describes the magnetic ground state in the ultrathin limit.

Frustrated magnets are commonly considered as a fertile ground for novel magnetic phenomena, with magnetic frustration playing a key role in spin textures by adding internal degrees of freedom and enhancing morphologies diversity.^66,67 In the strongly frustrated ML-NiI₂ system, HOIs can compete with or offset the Heisenberg exchange interaction and SIA, creating complex magnetic spin textures. To visualize these spin textures and examine the influence of HOIs under various biaxial strains, we performed atomistic spin dynamics simulations with parameters derived from the first-principles calculation. As shown in Fig. 4(c), a spin spiral state, also known as distorted spin spiral⁶⁸ or a multi-band spin helix,⁶⁹ emerges even without strain or a magnetic field. The simulated spin textures are sensitive to biaxial strain. Upon the application of compressive strain, the magnetic domain wall becomes thinner (see Fig. 4(a) and (b)), resulting from the significant increase of frustration ratio |J₃|/J₁ as shown in Fig. 3(b). At ε = −4%, an in-plane spin spiral state is observed due to the in-plane magnetic anisotropy, while greater tensile strain gradually thickens the domain wall and leads to a longer-period spin spiral as shown in Fig. 4(d) and (e), consistent with the propagation vector q in Fig. 1(e). Notably, when the strain reaches ε = 8%, coexistent magnetic field-free skyrmions and anti-skyrmions amidst wormlike domains can be found (Fig. 4(e)), which allows for more advanced logical operations in spintronic devices. Next, we include the HOI terms in the Hamiltonian and neglect the smaller three-site-four-spin interaction Y in the subsequent simulations given its insignificance compared to the biquadratic interaction B and four-site-four-spin exchange F. Simulation results displayed in Fig. 4(h) and (m) indicate that HOIs have minimal impact under strain-free conditions. However, at ε = 8% tensile strain, biquadratic interaction B removes worm-like domain walls, leading to larger magnetic domain clumps with the localized FM state (Fig. 4(j)), aligning with the findings shown in Fig. 1(d). Further introduction of F interaction results in small-size isolated skyrmions combined with clumps of large magnetic domain (Fig. 4(o)), where the diameter of skyrmions can be as small as only ∼2 nm, potentially enhancing the integration rate in spintronics. This phenomenon can be attributed to the fact that F can produce Dzyaloshinskii–Moriya interaction-like effects between neighboring spins, canting the magnetic moment of Ni atoms and moderately increasing topological charge density.^17,61,70 Additional finite-temperature spin dynamics simulations (see the SI for more details) confirm that skyrmions and antiskyrmions become metastable under thermal fluctuations and disappear near the transition temperature, consistent with the low Néel temperature. While in the case of ε = −8% with the equivalent large value of HOIs (see Fig. 4(k)), the spin texture splits into multiple domains in various directions, significantly enhancing antiferromagnetism.

Fig. 4 Magnetic field-free spin textures of ML-NiI₂ under different biaxial strains. The upper panels (a)–(e) only contain Heisenberg exchange J, the middle panels (f)–(j) show the interplay between J and the biquadratic interaction B, and the lower panels (k)–(o) show the interplay between J, B and four-site-four-spin exchange F. The insets are enlarged images of the yellow frame with dimensions of 10 nm × 10 nm. The arrows indicate the orientation of in-plane spin component, and color bar indicates the out-of-plane spin component.

A vertical magnetic field H with a gradient of 0.1T is applied to study the spin texture under varying strains and magnetic interactions, as shown in Fig. 5. For moderate strains of −4% to 4%, the spin spiral states stably persist under magnetic field, as illustrated in Fig. 5(a) and (b). The spin textures exhibit minimal alterations compared to those at zero magnetic field displayed in Fig. 4(l) and (m), corroborating with previous results.^68,69 At a larger compressive strain of ε = −8%, the inclusion of HOIs and magnetic field accentuates anti-ferromagnetism nature of the system (see Fig. 5(a)). Notably, at a tensile strain of ε = 8%, the spin texture is significantly influenced by the external magnetic field and HOIs. Focusing solely on the Heisenberg interaction J, as the magnetic field H increases from 0 to 1.5T, worm-like domain walls progressively collapse, creating a mixed state of skyrmions and anti-skyrmions with shorter domains, as shown in Fig. 5(d). Further escalating the magnetic field to H = 2.2T causes total disappearance of the domain walls, leading to a phase with the coexistence of skyrmions and anti-skyrmions (Fig. 5(g)). Obviously, the enhancement of magnetic field is in favor of the generation of small-sized spin texture. When B interactions are taken into account, isolated anti-skyrmions appear even at a low magnetic field of only 0.1T (see Fig. 5(e)). More appealingly, the mixed state of small-sized skyrmions and anti-skyrmions can be identified upon increasing the magnetic field to H = 0.4T (Fig. 5(h)), while in the presence of all HOI interactions, a smaller Bloch-type skyrmion with a diameter of only ∼1.5 nm forms at H = 0.1T (Fig. 5(f)), and increasing the magnetic field enhances the density of skyrmions and anti-skyrmions (Fig. 5(i)). We also test the Kitaev interaction and found that while it slightly modulates skyrmion size and density, the essential stabilization of ultra-small skyrmions and antiskyrmions remains dominated by strain-enhanced higher-order interactions (see the SI for more details).

Fig. 5 Modulation of spin textures by external magnetic field and HOIs. Spin textures under (a) ε = 0% and 0.5T magnetic field, (b) ε = −4% and H = 0.3T magnetic field, and (c) ε = −8% and H = 0.5T magnetic field, in which all the exchange interactions are considered. (d) and (g) The spin textures under ε = 8% tensile strain and magnetic field H of 1.5T and 2.2T, in which only Heisenberg exchange J interaction is considered. (e) and (h) Spin textures under ε = 8% and magnetic field of H = 0.1T and 0.4T, in which the interplay between J and B is considered. (f) and (i) Spin textures under ε = 8% and magnetic field of 0.1T and 0.3T, in which all the exchange interactions are involved. The insets are enlarged images of the yellow frame with dimensions of 15 nm × 15 nm and 7 nm × 7 nm.

To sum up, the introduction of HOIs enables the system to generate smaller and more complex spin textures at lower magnetic fields, providing valuable insights for the future application of spintronic devices based on skyrmions. As compared with the previously proposed anisotropic exchange pathway,³⁸ our work does not contradict these previous findings but rather reveals an alternative and complementary mechanism for stabilizing topological magnetism that significantly expands the mechanistic understanding in monolayer NiI₂. This HOIs-based mechanism differs from the anisotropic exchange fundamentally in three key aspects. First, the origin of non-collinearity stems from symmetric HOIs rather than two-site anisotropic exchange, specifically the strain-enhanced four-site-four-spin (F). Second, the resulting spin textures differ from each other. Our approach leads to the coexistence of fundamental skyrmions and antiskyrmions instead of the anti-biskyrmion lattice. These elementary textures offer a simpler topological profile and are potentially more suitable for spintronic encoding. Most importantly, the operational requirements are readily achievable. Unlike previous approaches that rely on enormous magnetic fields (∼100 T), our strain-mediated mechanism stabilizes textures at zero field and enables nucleation and control with minimal magnetic fields as low as 0.1–0.4 T. This orders-of-magnitude reduction in field strength significantly enhances practical feasibility for device integration.

4. Conclusions

In summary, we have proposed a novel theoretical mechanism for stabilizing topological spin textures in centrosymmetric monolayer NiI₂ through a combination of first-principles calculations and atomistic spin dynamics simulations. Our study theoretically identifies a spin spiral ground state with propagation vector q = 0.23, which exhibits robust stability under moderate biaxial strain and small external magnetic fields. The spin spiral period and phase stability can be effectively tuned via strain engineering, with a notable transition to a ferromagnetic state under 8% tensile strain. Furthermore, the single-ion anisotropy is significantly enhanced under both compressive and tensile strains, primarily governed by modifications in the p_xp_y hybridized orbitals of iodine atoms. A key emphasis of our work lies in highlighting the crucial role of higher-order interactions (HOIs), particularly the four-site-four-spin exchange F, which has been overlooked previously. Our simulations propose that the synergetic effect between HOIs and biaxial strain leads to the emergence of diverse spin textures, including the coexistence of skyrmions and antiskyrmions with diameters around 2 nm, without requiring an external magnetic field. The strain-tunable frustration ratio and external magnetic fields play critical roles in modulating the density and stability of these topological structures. Our findings provide a deeper understanding of spin textures mediated by magnetic frustration and HOIs, offering promising theoretical insights for the design of ultra-dense and low-power spintronic devices.

Author contributions

H. Y. X. and L. Q. conceived and supervised the project. Z. B. W. and W. G. H. performed the theoretical calculations and corresponding analysis, and wrote the manuscript. All authors contributed to the interpretation of the results and revised the manuscript.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: the effect of strain on the atomic structure of monolayer NiI₂; phonon dispersions for the monolayer NiI₂ under 0%, +8%, and −8% biaxial strain; calculation methods; calculation of the spin Hamiltonian; single-ion anisotropy (SIA) energy and orbital-resolved SIA; choice of Hubbard-U value; higher-order exchange interactions from DFT; specific values of magnetic exchange parameters; estimation of Dzyaloshinskii–Moriya interaction (DMI); negative integrated crystal Hamilton population (–ICOHP) for J₃; magnetic transition temperature and thermal stability; finite-size convergence tests; and influence of the Kitaev interaction on spin textures. See DOI: https://doi.org/10.1039/d5nh00706b.

Acknowledgements

This work was support by National Natural Science Foundation of China (Grant No. 52525208, U1930120) and Sichuan Science and Technology Program (Grant No. 2026NSFSCZY0033). The theoretical calculations were performed using the supercomputer resources at TianHe-1 located at the National Supercomputer Center in Tianjin.

References

1. S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii and P. Böni, Science, 2009, 323(5916), 915–919.
2. J. Tang, Y. Wu, W. Wang, L. Kong, B. Lv, W. Wei, J. Zang, M. Tian and H. Du, Nat. Nanotechnol., 2021, 16(10), 1086–1091.
3. S. Vélez, S. Ruiz-Gómez, J. Schaab, E. Gradauskaite, M. S. Wörnle, P. Welter, B. J. Jacot, C. L. Degen, M. Trassin, M. Fiebig and P. Gambardella, Nat. Nanotechnol., 2022, 17(8), 834–841.
4. S. Luo and L. You, APL Mater., 2021, 9(5), 050901.
5. Y. Zhou, Natl. Sci. Rev., 2018, 6(2), 210–212.
6. H. Du, D. Liang, C. Jin, L. Kong, M. J. Stolt, W. Ning, J. Yang, Y. Xing, J. Wang, R. Che and J. Zang, Nat. Commun., 2015, 6(6), 7637.
7. S. S. Parkin, M. Hayashi and L. Thomas, Science, 2008, 320(5873), 190.
8. X. Zhang, E. Motohiko and Y. Zhou, Sci. Rep., 2015, 5(1), 9400.
9. X. Zhang, W. Cai, X. Zhang, Z. Wang, Z. Li, Y. Zhang, K. Cao, N. Lei, W. Kang, Y. Zhang and H. Yu, Nanotechnology, 2016, 28(8), 08LT02.
10. X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui and Y. Tokura, Nat. Mater., 2011, 10(2), 106–109.
11. K. Shibata, X. Z. Yu, T. Hara, D. Morikawa, N. Kanazawa, K. Kimoto, S. Ishiwata, Y. Matsui and Y. Tokura, Nat. Nanotechnol., 2013, 8(10), 723–728.
12. Z. Shen, C. Song, Y. Xue, Z. Wu, J. Wang and Z. Zhong, Phys. Rev. B, 2022, 106(9), 094403.
13. H. Yoshimochi, R. Takagi, J. Ju, N. D. Khanh, H. Saito, H. Sagayama, H. Nakao, S. Itoh, Y. Tokura, T. Arima, S. Hayami, T. Nakajima and S. Seki, Nat. Phys., 2024, 20(6), 1001–1008.
14. A. Fert, V. Cros and J. Sampaio, Nat. Nanotechnol., 2013, 8(3), 152–156.
15. J. Brede, N. Atodiresei, V. Caciuc, M. Bazarnik, A. Al-Zubi, S. Blügel and R. Wiesendanger, Nat. Nanotechnol., 2014, 9(12), 1018–1023.
16. T. Okubo, S. Chung and H. Kawamura, Phys. Rev. Lett., 2012, 108(1), 017206.
17. S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer and S. Blügel, Nat. Phys., 2011, 7(9), 713–718.
18. T. A. D. S. Garel, Phys. Rev. B:Condens. Matter Mater. Phys., 1982, 26(1), 325–329.
19. N. Nagaosa and Y. Tokura, Nat. Nanotechnol., 2013, 8(12), 899–911.
20. L. Janssen, E. C. Andrade and M. Vojta, Phys. Rev. Lett., 2016, 117(27), 277202.
21. Y. Kato, S. Hayami and Y. Motome, Phys. Rev. B, 2021, 104(22), 224405.
22. K. Du, F.-T. Huang, K. Gamage, J. Yang, M. Mostovoy and S.-W. Cheong, Adv. Mater., 2023, 35(39), 2303750.
23. F. J. dos Santos, M. dos Santos Dias and S. Lounis, Phys. Rev. B, 2020, 102(10), 104436.
24. Y. Tokura and S. Seki, Adv. Mater., 2010, 22(14), 1554–1565.
25. S. R. Kuindersma, J. P. Sanchez and C. Haas, Physica B + C, 1981, 111(2), 231–248.
26. Q. Song, C. A. Occhialini, E. Ergeçen, B. Ilyas, D. Amoroso, P. Barone, J. Kapeghian, K. Watanabe, T. Taniguchi, A. S. Botana, S. Picozzi, N. Gedik and R. Comin, Nature, 2022, 602(7898), 601–605.
27. H. Ju, Y. Lee, K.-T. Kim, I. H. Choi, C. J. Roh, S. Son, P. Park, J. H. Kim, T. S. Jung, J. H. Kim, K. H. Kim, J.-G. Park and J. S. Lee, Nano Lett., 2021, 21(12), 5126–5132.
28. T. Kurumaji, S. Seki, S. Ishiwata, H. Murakawa, Y. Kaneko and Y. Tokura, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87(1), 014429.
29. A. S. Botana and M. R. Norman, Phys. Rev. Mater., 2019, 3(4), 044001.
30. M. Lu, Q. Yao, C. Xiao, C. Huang and E. Kan, ACS Omega, 2019, 4(3), 5714–5721.
31. H. Han, H. Zheng, Q. Wang and Y. Yan, Phys. Chem. Chem. Phys., 2020, 22(46), 26917–26922.
32. V. V. Kulish and W. Huang, J. Mater. Chem. C, 2017, 5(34), 8734–8741.
33. J. Y. Ni, X. Y. Li, D. Amoroso, X. He, J. S. Feng, E. J. Kan, S. Picozzi and H. J. Xiang, Phys. Rev. Lett., 2021, 127(24), 247204.
34. X. Li, C. Xu, B. Liu, X. Li, L. Bellaiche and H. Xiang, Phys. Rev. Lett., 2023, 131(3), 036701.
35. A. O. Fumega and J. L. Lado, 2D Mater., 2022, 9(2), 025010.
36. J. Sødequist and T. Olsen, 2D Mater., 2023, 10(3), 035016.
37. A. Ghojavand, M. Soenen, N. Rezaei, M. Alaei, C. Sevik and M. V. Milošević, Phys. Rev. Mater., 2024, 8(11), 114401.
38. D. Amoroso, P. Barone and S. Picozzi, Nat. Commun., 2020, 11(1), 5784.
39. M. Miao; N. Liu; W. Zhang; D. Wang; W. Ji and A. Y.-S. Fu, arXiv, preprint, arXiv: 2309.16526, 2023 DOI:10.48550/arXiv.2309.16526.
40. X.-s Ni, D.-X. Yao and K. Cao, J. Magn. Magn. Mater., 2025, 613, 172661.
41. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54(16), 11169–11186.
42. P. E. Blochl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50(24), 17953–17979.
43. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77(18), 3865–3868.
44. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 57(3), 1505–1509.
45. Q. Cui, J. Liang, B. Yang, Z. Wang, P. Li, P. Cui and H. Yang, Phys. Rev. B, 2020, 101(21), 214439.
46. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B:Condens. Matter Mater. Phys., 2008, 78(13), 134106.
47. A. Fert, M. Chshiev, A. Thiaville and H. Yang, J. Phys. Soc. Jpn., 2023, 92(8), 081001.
48. C. Di Giorgio, E. Blundo, G. Pettinari, M. Felici, Y. Lu, A. M. Cucolo, A. Polimeni and F. Bobba, Adv. Mater. Interfaces, 2020, 7(23), 2001024.
49. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321(5887), 385–388.
50. E. Blundo, C. Di Giorgio, G. Pettinari, T. Yildirim, M. Felici, Y. Lu, F. Bobba and A. Polimeni, Adv. Mater. Interfaces, 2020, 7(17), 2000621.
51. D.-s Wang, R. Wu and A. J. Freeman, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47(22), 14932–14947.
52. Z. Wu, Z. Shen, Y. Xue and C. Song, Phys. Rev. Mater., 2022, 6(1), 014011.
53. Z. Shen, Y. Xue, Z. Wu and C. Song, Phys. Chem. Chem. Phys., 2022, 24(45), 27612–27618.
54. X. J. Dong, M. J. Ren and C. W. Zhang, Phys. Chem. Chem. Phys., 2022, 24(36), 21631–21637.
55. J. B. Goodenough, Phys. Rev., 1955, 100(2), 564–573.
56. J. Kanamori, J. Phys. Chem. Solids, 1959, 10(2), 87–98.
57. P. W. Anderson, Phys. Rev., 1959, 115(1), 2–13.
58. J. B. Goodenough and F. A. Cotton Interscience-Wiley, New York, 1963.
59. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37(11), 1030–1035.
60. A. Kartsev, M. Augustin, R. F. L. Evans, K. S. Novoselov and E. J. G. Santos, npj Comput. Mater., 2020, 6(1), 150.
61. P. Li, D. Yu, J. Liang, Y. Ga and H. Yang, Phys. Rev. B, 2023, 107(5), 054408.
62. A. Cong and K. Shen, Phys. Rev. B, 2024, 109(22), 224419.
63. N. Liu, C. Wang, C. Yan, C. Xu, J. Hu, Y. Zhang and W. Ji, Phys. Rev. B, 2024, 109(19), 195422.
64. M. Amini, A. O. Fumega, H. González-Herrero, V. Vaňo, S. Kezilebieke, J. L. Lado and P. Liljeroth, Adv. Mater., 2024, 36(18), 2311342.
65. T. V. C. Antao, J. L. Lado and A. O. Fumega, Nano Lett., 2024, 24(49), 15767–15773.
66. A. O. Leonov and M. Mostovoy, Nat. Commun., 2015, 6(1), 8275.
67. X. Yao and S. Dong, Phys. Rev. B, 2022, 105(1), 014444.
68. T. T. J. Mutter, A. O. Leonov and K. Inoue, Phys. Rev. B, 2019, 100(6), 060407.
69. H. Y. Yuan, O. Gomonay and M. Kläui, Phys. Rev. B, 2017, 96(13), 134415.
70. M. Hoffmann and S. Blügel, Phys. Rev. B, 2020, 101(2), 024418.

---