3D Heisenberg universality in the van der Waals antiferromagnet NiPS3

Van der Waals (vdW) magnetic materials have emerged as an excellent platform for investigating low-dimensional magnetism, owing to their unique structural characteristics and tunable magnetic properties^1,2,3. Their remarkable sensitivity to external stimuli, such as electric or magnetic fields, strain, and temperature, makes them highly desirable for applications in next-generation spintronics and quantum computing devices⁴. Furthermore, they offer an ideal platform for investigating fundamental magnetic Hamiltonians in real physical systems⁵. Some of these model Hamiltonians can be solved analytically and have led to many discoveries that have guided the field of low-dimensional magnetism, such as the Berezinskii–Kosterlitz–Thouless (BKT) transition^6,7 in the Heisenberg XY model. As a result, vdW magnetic materials have garnered tremendous attention from both theorists and experimentalists, who seek to explore their fundamental physics and exploit their technological potential^3,8.

The Transition Metal Chalchogenophosphates (TMCPs), include a class of isostructural vdW antiferromagnets (AFMs) which are easy to synthesize and chemically modify, making them exceptional for exploring low dimensional magnetic properties⁹. In addition, the various compounds within this family all exhibit honeycomb magnetic structures, which can often be modeled by principal spin-Hamiltonians^3,5. For instance, the honeycomb arrangement of Ni ions in NiPS₃ bears semblance to graphene, and its layers are similarly easy to exfoliate. Similarly to graphene, NiPS₃ has been established as a rich system for exploring strongly correlated dynamics^{3,10,11,12,13,14,15,16}. However, there are still many open questions related to whether or not the long-range magnetism in NiPS₃ is intrinsically different in its bulk or monolayer forms. It was recently shown that the “zig-zag” AFM order persists when the thickness is reduced down to two layers but disappears altogether in the monolayer limit¹⁶, indicating that pure 2D magnetism may not be stable in this system. Furthermore, the question of the spin-dimensionality in NiPS₃ is still debated at this time, with some studies reporting it as best being described by a two-component spin (XY-like)^12,14,17,18 and others reporting it as a three-component spin (Heisenberg-like)^5,13,19,20. Without a well-established form and set of parameters for the spin-Hamiltonian, some aspects of the magnetic dynamics in NiPS₃ remain to be fully described.

A key current hurdle in studying the detailed magnetic properties of vdW crystals is the scarcity of techniques that can directly characterize the long-range magnetic order within the 2D layer with atomic spatial resolution¹. Historically, neutron scattering, which can be used to analyze the magnetic structure and excitations in crystals, has been instrumental in pioneering our understanding of NiPS₃ and related compounds^{17,19,20,21,22,23}, but these are typically conducted on powder samples or large batches of crystal samples fused together. Single vdW crystals are usually small and fragile, which pose inherent limitations for measurements that require large crystals. In this context, synchrotron and X-ray free-electron laser facilities present an opportunity to provide another perspective for analyzing vdW magnets. While typically it requires sufficient surface polishing because of the smaller penetration depth compared to neutrons, X-ray beam spot sizes however can be focused to a few microns such that tiny crystals can be studied. This combined with a strong X-ray cross-section allows thin vdW crystals to be analyzed without averaging over large sample volumes. Additionally, by utilizing their extreme brightness and tunable energy, the ability to probe electronic and magnetic states with elemental specificity is also advantageous. Thus resonant X-ray diffraction is particularly promising for the study of vdW systems, yet remains largely unexplored.

In this article, we address these challenges by directly studying the properties of the long-range magnetic order in NiPS₃ using resonant X-ray diffraction from a single crystal. We explicitly probe the d-states of Ni with an incident X-ray beam energy tuned to the L-edge (2p → 3d) resonance on a system that can not be polished easily due to the delicate nature of these magnetic vdW layers. The in-plane antiferromagnetic order was revealed by orienting the scattering plane of the X-ray beam parallel to the vdW layers, with the beam incident on the natural crystal surface along the edge of the layers. Using theoretical calculations and measuring the critical scaling of the AFM so-called “zig-zag” sublattice order parameter, we find that the magnetic structure surprisingly scales as the 3D Heisenberg universality class. This is in contradiction to the many 2D universality classes conventionally assigned to the van der Waals magnet family^4,24. This points to the importance of exploring the nature of these magnetic interactions, such as through magnetic fluctuation studies, where the dimensionality largely determines the role of thermal fluctuations.

Two resonant X-ray diffraction experiments were carried out at the Ni L-edge at 853 eV. These two experiments were performed at the Stanford Synchrotron Radiation Light Source (SSRL) and Advanced Light Source (ALS) facilities. Details about these experimental setups, as well as an additional survey of the magnetic and structural reflections at the Ni K-edge (1s → 3d) at 8.332 keV from the Advanced Photon Source (APS), are discussed in Supplementary materials.

Obtaining the (0k0)-reflection from NiPS₃ poses a practical challenge because it requires scattering along the magnetic layer, i.e. making the scattering plane co-planar with the vdW layers (ab-plane). The scattering geometry we employed for all these experiments is shown in Fig. 1a together with a scanning electron microscopy (SEM) image of the crystal Fig. 1b used in the X-ray measurements taken with a small beam size at the ALS facility (see experimental section). This was used to explore the preferred scattering region on the crystal edge, where optimal surface regions were identified with a length-scale several times larger than the incident beam spot size. The SEM image shows the morphology of the crystal as seen from its side-profile, where some layers are delaminated from the rest of the bulk crystal in certain regions along the edge. These regions were avoided by marking the position of region with the highest quality and using a small beam size.

a Magnetic unit cell of NiPS₃ viewed along the c-axis and X-ray scattering setup schematic showing the orientation of the crystal and copper mount with respect to the beam (shown in red). The scattering plane is parallel to the ab-plane. b Scanning electron micrograph of the edge of the NiPS₃ crystal used in the X-ray experiments, viewed from its side profile. The sample is oriented to be viewed along the b-axis in order to locate the optimal location to place the beam to avoid regions where the layers have deformed or delaminated from the bulk crystal. The beam spot size was approximately 8 μm. Inset: the preferred scattering region where the vdW layers exhibit a smooth surface.

Full size image

The atomic structure is shown in Fig. 1a, with the Ni atoms shown in solid red. The (0k0) peaks are forbidden by this crystal structure for odd k, but since the “zig-zag” AFM order (one spin chain with spin pointing along a, and the opposite spin chain pointing in the (bar{a})-direction) possesses a different symmetry than the parent lattice, the magnetic structure can be directly studied at these forbidden reflections without background from non-resonant charge scattering. No scattering was observed at the magnetic reflections off-resonance at base temperature. Tuning the beam energy to the Ni L-edge, an incredibly strong resonance was observed with an energy width of ~2 eV at the same momentum transfer. Figure 2 shows the temperature-dependent resonant magnetic scattering intensity from integrated rocking curves of the (010) reflection. The measurements show a decrease in intensity as the sample is warmed until vanishing at the T_N, corresponding to the antiferromagnetic/paramagnetic phase transition, which was also confirmed with magnetic susceptibility measurements of the same sample (see Fig. 2c).

a, b θ-scan intensities of the magnetic (010) Bragg peak collected at 853 eV (Ni L-edge) for different temperatures from SSRL and ALS lightsources, respectively. Intensity traces are scaled and offset for visualization. Insets show the total intensity at each temperature, normalized to their highest temperature value where the intensity is lowest. c Magnetic susceptibility data using vibrating sample magnetometry with an applied B-field of 1 T ∥ with the a-axis for both field-cooled (FC) and zero-field-cooled (ZFC) measurements. d Integrated intensity of SSRL (shown in blue) and ALS (shown in orange) normalized to their value at T = 0.97 × T_N shows their functional form as they approach the magnetic phase transition. The Néel temperature T_N = 155.5 K is marked with a vertical line. The curves corresponding to the fitted scaling value β = 0.367 ± 0.013 and the theoretical values for the other models for members of this family of β = 0.231 (2D XY) and β = 0.125 (2D Ising) are shown together with the data.

Full size image

The soft X-ray data confirmed a temperature dependence in agreement with the “zig-zag” magnetic order parameter and transition from an AFM to a paramagnetic phase at T_N= 155.5 K. The X-ray scattering intensity I—here the integration of the θ-rocking curve of the (010) reflection—is proportional to the sublattice magnetization squared, and can be fit to the form, (I propto alpha (-tau)^{2beta}) for the ordered phase (T < T_N), where τ is the reduced temperature (T − T_N)/T_N and α is a constant of proportionality. In order to compare the intensities measured at both instruments, both curves are normalized to their value closest to T = 0.97 × T_N, as shown in Fig. 2d. Note that critical scaling analysis is only meaningful very near to the transition temperature²⁵. The fitted value β = 0.367 ± 0.013 was obtained by fitting both sets of data from the SSRL and ALS measurements together in this region. These dual soft X-ray studies confirmed the magnetic ground state, Néel temperature, and yielded the critical index, which we can use to compare our observations with predictions made from theory. Universality of the critical temperature scaling points to the 3D Heisenberg universality class, which we will explore further with ab initio calculations.

We performed independent first-principles Density Functional Theory (DFT) calculations, and use these results to initiate additional theoretical calculations using Density-Matrix Renormalization Group (DMRG) and Monte Carlo simulations. We outline the results here while the theoretical details can be found in the Methods section. Because there is a large discrepancy in the literature for NiPS₃, we used first-principles DFT²⁶ to extract the magnetic exchange parameters using two different methods, designated “Method A” and “Method B” (as shown in Table 1). DFT was performed here using the Vienna ab-initio simulation package (VASP)^27,28. Second, to calculate the magnetic exchange parameters, we used the single-particle Green’s function constructed from a low-energy tight-binding model by the atom-centered Wannier function using Ni-d and s–p orbitals in the Wannier90 code suite²⁹ (See Methods section for details). The exchange parameters are then used to build an atomistic spin model based on the Heisenberg Hamiltonian in Eq. (1) and simulate the temperature effect on the magnetization following the Monte Carlo Metropolis algorithm³⁰ as implemented in the Vampire software package³¹. This predicted the Néel temperature to be T_N = 153 K and was confirmed to have good agreement with the magnetic susceptibility and X-ray data (Fig. 2). The Monte Carlo temperature dependence of the squared magnetization is shown in Fig. 2d. Finally, to confirm the zig-zag magnetic ordering in the ground state, the exchange parameters were also used to calculate the ground state using the DMRG method. Due to the 3d⁸ Ni atom having two unpaired electrons in the e_g orbitals, we use the spin-1 Heisenberg model with exchanges J₁, J₂, J₃, and J₄ on the honeycomb lattice with single-ion anisotropy A_ii. The definition of this spin-Hamiltonian is,

where the S_i are the quantum spin operators acting on the ith site. J_ij equals J₁ when i and j are the nearest neighbors and equals J₂ and J₃ for second and third nearest neighbors, respectively. For the Monte Carlo simulation of the bulk spin dynamics, the small interlayer exchange J₄ is included. The DMRG result confirms the stability of the magnetic structure in the single layer and that the magnetic exchange parameters yield a zig-zag order for the ground state measured here, with ferromagnetic chains along the a-axis and antiferromagnetic modulation along the b-axis with modulation vector Q_mag = [010] as shown in Fig. 1. The spins are canted slightly above and below the ab-plane.

Although these types of TMCP materials are primarily thought to exhibit 2D magnetism due to their vdW structure, the critical exponent measured here, and confirmed with theory, is surprisingly consistent with the universality class of the 3D Heisenberg model. We have demonstrated that our experimental findings suggest a 3D Heisenberg-like spin Hamiltonian, and that this model predicts a static ground state that is consistent with our data. To further explore this, the atomistic spin model which was used to calculate the temperature effect on the magnetization via Monte Carlo simulation was used for annealing the magnetic state to simulate the temperature dependence of the order parameter. Fitting these calculated values as a function of temperature also shows the order parameter exponent to be β = 0.36, within our error bars as measured through the soft X-ray measurements. This further confirms that our experimental observations and first-principles calculations using the Heisenberg 3D model yield a consistent picture in NiPS₃.

The critical exponent closely matches the three-dimensional Heisenberg universality class (β = 0.366), which is characterized by a three-component order parameter and an O(3) symmetry, usually used to describe isotropic magnets^4,32. This result is intriguing as it is well above the value for the critical exponent assigned to spin-Hamiltonians from 2D universality classes²⁴. This was determined in the same sample for low temperatures¹⁷, and in related systems, i.e. the 2D XY model (β = 0.231) for MnPS₃³³, or the 2D Ising model (β = 0.125) for FePS₃³⁴.

This result shows NiPS₃ can be described by a three-dimensional and three-component spin-1 Heisenberg model and a single-ion anisotropy with a weak interlayer exchange J₄. This was hinted at in the neutron diffraction work mentioned above which also analyzed the temperature dependence of the magnetic order parameter near the transition and found a value of β = 0.3, implying a crossover to 3D-like behavior¹⁷. While the value of β = 0.3 for high-T puts the spin-Hamiltonian in the 3-dimensional regime, the question of the precise form taken among 3D Ising (β = 0.326), 3D XY (β = 0.348) and 3D Heisenberg (β = 0.366), which are distinguished by their number of relevant spin components (n = 1, n = 2, and n = 3 respectively), was left unanswered. The data fitted from our experiments directly confirm the 3D Heisenberg behavior from the critical exponent. This has been suggested using other methods, such as a recent study by Scheie et al, which mapped the magnetic Brillouin zone using Inelastic Neutron Scattering (INS) which observed a weak inter-plane magnon dispersion corresponding to a ferromagnetic exchange of −0.38 meV²⁰. Taken together, these studies provide complimentary insights into the magnetic properties of NiPS₃, and share similar qualitative conclusions about the dimensionality of the spin dynamics.

We further find that the AFM Heisenberg exchange, in particular the third-nearest neighbor exchange J₃ dominates the Hamiltonian. Though the interlayer ferromagnetic exchange J₄ is crucial for introducing the 3D character, it is much smaller in comparison to the in-plane interactions and insufficient for creating long-range correlations along the c-axis. Hence, NiPS₃ can be thought of as an easy-plane antiferromagnet, with magnetic vdW layers which are weakly coupled ferromagnetically. Further exploration will be needed to determine an exact description of the Hamiltonian, where details are emerging in the lower energy magnon excitations, such as possible bi-axial characteristics and other effects that cannot be explained by linear spin wave theory^20,35. New tools are being developed in both the X-ray and neutron scattering frontiers to elucidate these small energetic differences at the fringe of current experimental capabilities^36,37.

The main motivation here however was to use resonant X-rays to probe the magnetic structure as well as its length scale to provide direct evidence of the AFM state. Measuring the in-plane AFM order using X-rays poses significant challenges due to the inherent limitations of accessing in-plane features of thin vdW crystals, which is typically attempted in a grazing incidence from the plane orthogonal to the c-axis direction. However, our study demonstrates that it is indeed feasible to probe the AFM order in NiPS₃ using both soft (and hard) X-rays by scattering resonant X-rays directly from the edge of the layer. We show clean growth is enough to provide a smooth region of the surface to be accessed, without resorting to polishing or cleaving the sample, which would result in either powderization or destroying the material. This work opens up new possibilities for investigating AFM materials using more advanced X-ray scattering techniques, such as pump-probe studies at X-ray free-electron lasers³⁶, where ultra-fast magnetism can be explored, or fast fluctuations using advanced pulse structure methods³⁸ combined with coherent X-rays, such as in X-ray Photon Fluctuation Spectroscopy³⁹. By successfully measuring the AFM order directly with X-rays, we provide experimental evidence that expands the applicability of X-ray methods to a wide range of low-dimensional magnetic systems.

In conclusion, our study demonstrates the use of resonant X-ray diffraction to measure the AFM order directly along the vdW layer in NiPS₃, using single crystals, opening new avenues for investigation in these types of quasi-2D magnetic systems. The magnetic reflections were strong enough to be observed for two different types of experiments, at soft X-ray energies using the L-edge in Ni, where magnetic signal is usually large for transition metals, as well as for higher energies using resonance at the K-edge (see Methods Section). These results were verified at three different X-ray synchrotron light sources. Our observations confirm the “zig-zag” AFM structure previously measured by neutron diffraction, and is further supported by our DMRG ground state calculations for a spin-1 Heisenberg Hamiltonian on a honeycomb lattice. This theoretical calculation was made possible by first-principles DFT, used to extract the parameters needed for the ground state. Additional theoretical evidence was reported in solving the Landau–Lifshitz–Gilbert spin dynamics equation through Monte Carlo simulations, which could reproduce the Néel temperature and the critical exponent within the errors of our measurement. Through this experimental and theoretical evidence, we find the vdW antiferromagnet NiPS₃—conventionally considered a 2D XY magnet—is unexpectedly more closely described by the 3D Heisenberg universality class. In other words, the magnetism can not be thought of as a two-component order parameter, nor can the interlayer interaction be neglected. It will be interesting to explore this finding further by studying other members of the TMCP family with X-rays, described by different Hamiltonians, to see if the 3D nature is more general or a special case for the largely in-plane spin-1 system. Moreover, because fluctuations dominate the behavior in 2D, it would be enlightening to explore the fluctuations directly of these magnetic vdW systems as has been done in other systems⁴⁰, to see if they map on to the requisite universality classes.

Sample preparation

Single crystals of NiPS₃ were used for the series of experiments in this study. They were grown with optimized growth parameters by HQ Graphene using chemical vapor transport. The crystals naturally grow as hexagonal platelets. From a typical batch, the thickest crystals were selected with the most well-defined edges for the resonant X-ray experiments. Energy-dispersive X-ray spectroscopy was then used to confirm the composition and the a and b edges were identified by orienting the crystals in a Bruker D8 Venture single-crystal X-ray diffractometer.

X-ray beamline experiments

All experiments were conducted with a single NiPS₃ crystal mounted on a copper sample holder using silver paint. Different single-crystal specimens were used in each experiment, however their shape and dimensions are approximately the same. The thickness of the edge in the c-direction is roughly 100 μm and the mass of the crystals is ~2.5 mg. Mechanically, the crystals are soft and easily deformed. The quality of the edge surface is especially critical for soft X-ray experiments due to the increased surface sensitivity at lower X-ray photon energies. Additionally, the magnetic characteristics of vdW crystals are known to vary dramatically with lattice strain and stacking arrangement^1,4. With this in mind, the crystals were surveyed with SEM at the Stanford Nano Shared Facilities laboratory in order to locate the regions along the edge with minimal surface disorder.

Two resonant X-ray diffraction experiments were carried out using soft X-rays at the Ni L₃-edge (853 eV). The first was at Beamline 13-3 of the Stanford Synchrotron Radiation Light Source (SSRL), and the second was at Beamline 7.0.1.1 of the Advanced Light Source (ALS). Resonant X-ray diffraction was also measured at the Ni K-edge (8.332 keV) at Beamline 4-ID-D of the APS using a four-circle diffractometer. Vibrating Sample Magnetometry (VSM) was used to characterize the temperature dependence of the magnetic susceptibility and confirm AFM magnetic ordering in the sample. No differences between the field-cooled and zero-field-cooled were observed, and the magnetic susceptibility was in good agreement with previously reported measurements¹⁷.

The soft X-ray data from beamline 13-3 of the SSRL was collected at the Ni L-edge (853 eV) in a vertical scattering geometry as illustrated in Fig. 1a using a 4-circle diffractometer. The sample was mounted to a copper puck using colloidal silver paint. The energy was selected using a spherical grating monochromater. The sample was cooled using a liquid helium cryostat and sustained under ultra-high vacuum. The beam spot size at the sample location was ~250 microns. X-ray scattering intensity was measured using a Hamamatsu GaAsP photodiode.

The soft X-ray data from beamline 7.0.1.1. of the ALS was collected at the Ni L-edge (853 eV) in a vertical scattering geometry as illustrated in Fig. 1a. The energy was selected using a spherical grating monochromater. The sample was cooled using a liquid nitrogen cryostat and sustained under high vacuum. The sample was mounted to a copper puck using colloidal silver paint. X-ray scattering intensity was measured using a Si photodiode. A pinhole was used to achieve a small beam spot with size ~8 microns.

The hard X-ray data (shown in Supplementary materials) at the Ni K-edge (8.332 keV) was collected in a vertical scattering geometry with horizontal polarization as illustrated in Fig. 1a using a 4-circle diffractometer. The energy was selected using a double crystal Si (111) monochromator, which was detuned by 2 arcsec in order to reduce the intensity to 87% of the maximum. The beam shape was measured to ~0.18 mm in the horizontal direction and 0.11 mm in the vertical. The sample was mounted to a copper puck using colloidal silver paint and was cooled with a liquid helium cryostat. An energy-dispersive Silicon drift detector was used to measure the elastic X-ray scattering intensity. This removes the Kα fluorescence but not the Kβ, which can be seen in Supplementary Fig. 3 as background in the fixed-Q energy scan measurements.

Density functional theory

To obtain insight into the magnetic structure of NiPS₃, we performed first-principles calculations within the framework of the density functional theory (DFT)²⁶ using VASP^27,28. Ground-state electronic structure was obtained using the projector-augmented-wave pseudo-potential. Electron exchange-correlation effects were treated by using the generalized gradient approximation⁴¹ with Perdew–Burke–Ernzerhof (PBE) parametrization. Strong correlation effects on Ni-d orbitals were accounted for by adding an effective onsite Hubbard potential (U_eff)^42,43. Throughout our calculations, we used U_eff = 4.0 eV. The weak interlayer vdW interactions were treated using the DFT-D3 correction of Grimme with zero-damping function⁴⁴. An energy cut-off of 350 eV was used for the plane-wave basis set, and Brillouin zone integrations were performed with an 11 × 9 × 11 Γ-centered k − mesh⁴⁵. Total energies were converged to 10⁻⁵ eV. Experimental lattice parameters were used, but ionic positions were optimized until the residual force on each ion was less than 10⁻² eV/Å and the stress tensors became negligible.

Two methods were used to calculate magnetic exchange parameters. (1) Single-particle Green’s function with rigid spin-rotation as a perturbation after Liechtenstein, Katsnelson, Antropov, and Gubanov (LKAG)⁴⁶ as implemented in TB2J code⁴⁷. To obtain the single-particle Green’s function, we constructed a low-energy tight-binding model using Ni-d and s–p orbitals in the Wannier90 code suite²⁹. (2) The four-state method⁴⁸ in which the atom pair involving the magnetic interactions is isolated in a (2sqrt{2}times sqrt{2}times 2) super-cell.

Ground state calculations

We used the complex number DMRG method to confirm the magnetic order of the ground state observed by magnetic X-ray diffraction⁴⁹ using the spin-1 Heisenberg model with J₁, J₂, and J₃ on the single-layered honeycomb lattice with single-ion anisotropy A_ii as defined in Eq. (1).

The coupling constants obtained from the four-state energy mapping method⁴⁸ are: J₁ = 3.81 meV, J₂ = 0.15 meV, and J₃ = −14.6 meV. The minus sign of J₃ implies that the superexchange between two spins is antiferromagnetic; otherwise, it is ferromagnetic.

We also used the four-state energy mapping method to compute the Single-Ion Anisotropy (SIA) and show that the diagonal part of the SIA Hamiltonian has a small non-zero value, i.e. ({A}_{ii}^{zz}{({S}_{i}^{z})}^{2}) with ({A}_{ii}^{zz}=-0.1) meV. In DMRG calculations, we used J₁ as the energy unit and redefined J₂₍₃₎/J₁ and ({A}_{ii}^{zz}/{J}_{1}) for simplicity. The lattice geometry used was cylindrical with two basis vectors e₁ = (a, 0) and ({{bf{e}}}_{2}=({bf{a}}/2,sqrt{3}{bf{b}}/2)) for each unit cell with periodic and open boundary conditions along the e₂ and e₁ directions, respectively, with corresponding width N_y and length N_x, where N_y and N_x are the numbers of spins along the e₂ and e₁ directions. Note that each unit cell contains two spins. This paper focuses primarily on the four-leg cylinder with N_y = 4 and N_x = 12 because the DMRG method is more precise for computing quasi-one dimensional systems with finite cylinder widths but a relatively longer cylinder length. By keeping up to 2187 states in each DMRG sweep with a typical truncation error ϵ ~ 10⁻⁵, we obtained a variational ground state with magnetic ordering. The spin symmetry-breaking directions of each spin have three components, (langle {S}_{i}^{x}rangle), (langle {S}_{i}^{y}rangle) and (langle {S}_{i}^{z}rangle). Without SIA, all spins form the zig-zag order, but each spin is canted with more pronounced (langle {S}_{i}^{x}rangle) and (langle {S}_{i}^{y}rangle) components and a small out-of-plane (langle {S}_{i}^{z}rangle) component. A finite SIA further suppresses (langle {S}_{i}^{z}rangle) and drives the magnetic ordering into a perfect zig-zag order within the ab-plane. Conclusively, from our simulations, even without the interlayer coupling, the zig-zag ordered state can be stabilized by the intralayer Heisenberg spin exchange J₁, J₂, and J₃ with ferromagnetic J₁ and J₂ plus more significant antiferromagnetic J₃. Introduction of SIA further suppresses the out of ab-plane spin component, allowing the 2D zig-zag order to form.