Crucial role of magnetic dipole and anomalous field dependence in $\mathrm{NiPS_3}$

Magnetic van der Waals (vdW) materials have emerged as fertile platforms for both fundamental physics and device applications [2, 3, 4, 5] [1]. NiPS$_3$ is a particular representative case under active recent investigation [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] as a prototypical antiferromagnet with XY-like (or XXZ) spins, supported by susceptibility and thickness-dependent magnetization measurements [6, 20, 25, 26, 27]. Recently, a series of experiments report the intriguing many-body physics in NiPS$_3$ [6, 7, 8, 9, 10, 11, 12, 19, 23]. In order to scrutinize the characteristic behaviors of this material, new theoretical models or pictures have been suggested such as Hund’s exciton [12] and coherent many-body exciton [6, 7, 8] associated with Zhang-Rice triplet-singlet transition [6, 13, 21, 28].

It is quite certain that these observations are intrinsically attributed to magnetism, but remain not captured by the simple XY model. For example, the sudden change in photoluminescence (PL) peak polarization likely indicates higher-order anisotropy [7, 15]. The critical exponent was found to be larger than that of XY model [29]. Neutron measurements also suggest symmetry-breaking terms within the easy-plane [14, 18, 29] [30]. Optical pump polarization is thought to be strongly associated with in-plane anisotropy [9, 10], and resonant inelastic X-ray scattering (RIXS) data are interpreted as having a correlated excitonic origin [6, 12]. While these experiments demonstrate how the great advancements of modern measurement techniques unveil otherwise hidden critical details of the low-energy excitations and responses, their comprehensive understanding remains far from clear and is still under active debate [6, 7, 12, 13].

A full account of anisotropic interactions beyond the XY model is required.

In principle, such an effective magnetic Hamiltonian beyond the XY model can be constructed by incorporating relativistic corrections of order $1/(mc)^2$ as derived from quantum electrodynamics [31]. For real materials, however, constructing a Hamiltonian that includes terms such as Kitaev-$\Gamma$, Dzyaloshinskii–Moriya (DM), and single-ion anisotropy, and clarifying their roles, remains highly challenging [32, 33, 34, 35]. This requires precise determination of energy scales typically in the $\mu$eV to sub-meV range, which remains demanding even for state-of-the-art theoretical and experimental methods.

To unravel the anisotropic magnetic profile in NiPS$_3$, we conducted torque magnetometry together with relativistic density functional theory (DFT) calculations. The high-field torque response ($H \geq9$ T) reveals an intriguing anomaly. Furthermore, complete three-dimensional (3D) understanding of torque data, combined with the results of relativistic DFT and magnetic force calculation, provides the detailed anisotropic interaction profile. From this we established a generalized spin Hamiltonian incorporating terms down to the $\mu$eV scale, which consistently explains both the previously reported many-body magneto-optical features [7, 15, 21, 22] and torque anomaly associated with an out-of-plane spin-flop transition. Remarkably, the dipolar interaction plays an essential role, being directly responsible for the characteristic torque and PL behaviors by favoring a specific energetic order. Our results establish a unified framework to interpret recent experiments and guide future applications.

NiPS$_3$ consists of edge-shared NiS$_6$ octahedra forming a honeycomb lattice and exhibits a zigzag antiferromagnetic (AFM) order below $T_N = 155$ K with $C_{2/m}$ symmetry due to the stacking as shown in Figure 1(a,Figure 1b) [6, 12, 25, 29, 36, 37] [30]. Previous magnetization studies were limited to $H \leq9$ T [29, 38, 39] [30]. Here, we extended the measurements up to $H = 16$ T, as summarized in Figure 1(c) (see Figure 1(e) for the angle convention). At low fields $H \leq9$ T, the magnetization increases linearly and shows larger values along the $c^*$-axis than along the $b$- and $a$-axes ($M_{\parallel a} < M_{\parallel b} \leq M_{\parallel c^*}$), consistent with previous reports [29, 38, 39] [30]. The in-plane magnetization exhibits easy-axis anisotropy, already deviating from the ideal XY behavior [26, 27], in agreement with the susceptibility data (Fig. S4). At high fields around $H \approx10$ T, a characteristic up-turn feature appears, un-

reported in earlier studies, and is observed only when the field is applied nearly along the $a$-axis. Above $H \geq12$ T, the magnetization curves for the $a$- and $b$-axes eventually merge. Notably, this critical field coincides with the reported $H \approx10.5$ T at which a sudden change in PL polarization occurs [7]. Below we will show that these two experiments share the same origin.

To map out the detailed magnetic Hamiltonian profile, we performed the torque magnetometry measurements as depicted in Figure 1(f). Below $T_N$, the torque is expressed as $\vec{\tau} = (\vec{\chi}\cdot\vec{H})\times\vec{H}$. Owing to the mirror symmetry, the off-diagonal components of susceptibility ($\chi_{ab}$, $\chi_{bc}$) vanish [40]. Under an in-plane rotating field $\vec{H}=H(\cos\varphi,\sin\varphi,0)$, the torque thus simplifies to describe the angular dependence in the linear-response regime:

Figure 2(a) shows that the measured torque $\tau_{ab}$ indeed follows this $\sin2\varphi$ dependence well in the low field $H \leq5$ T. This simple trigonometric form of data already shows the need of in-plane magnetic anisotropy terms beyond the ideal XY spin model, as discussed above with magnetization data.

An intriguing new feature starts to emerge as the field strength increases. As shown in Figure 2(a), a characteristic “shoulder” develops at $\varphi\approx180^\circ$ when the field strength is greater than $H \approx7$–$9$ T. It eventually evolves into a well-distinguished second peak at $H \geq12$ T. Obviously, this new shape cannot be captured by $\sin2\varphi$. Hypothesizing that it originates from an intrinsic material characteristic, we consider a general form of spin Hamiltonian as follows:

Here, $\hat{S}_i$ and $\hat{S}_j$ are the unit spin vectors at the local sites $i$ and $j$, respectively. $\mathbf{J}_{ij}$ denotes the exchange coupling tensor that includes both isotropic and anisotropic exchange interactions such as Heisenberg, DM, and Kitaev-$\Gamma$ interactions [32, 33, 34, 35, 41, 42]. The second term represents the single-ion anisotropy, expressed as $\sum_i K_i(\hat{e}_i\cdot\hat{S}_i)^2$ for a uniaxial case with strength $K_i$ along the local axis $\hat{e}_i$. Both exchange and single-ion anisotropies originate from spin–orbit coupling, and therefore one needs to resort to relativistic DFT. In the current study, they were obtained respectively from the magnetic force linear response theory (MFT) and the DFT total energy calculation. The third term $D(\hat{S}_i,\hat{S}_j)$ corresponds to the magnetic dipole–dipole interaction and is given by

where $\mu_0$ is the vacuum permeability and $m_i,m_j$ denote the magnetic moments at each site. The torque data can be compared by taking the external field into account; $\mathcal{H}_{\mathrm{ext}}=-\sum_i\vec{m}_i\cdot\vec{H}$. Once the parameters of each term are given, the corresponding total energy $E(\{\hat{S}_i\};\vec{H})$ depends on the moment directions, and one can obtain the torque by minimizing the energy:

where $\vec{M}_{\mathrm{total}}$ is the induced magnetization. Below, we try to construct the Hamiltonian and systematically investigate each term. Finally, Equation (4) gives rise to the result shown in Figure 2(b) in excellent agreement with the experiment.

We now examine the contribution of each term. The first term of Eq. (2) is calculated using MFT combined with relativistic non-collinear DFT [41, 43, 44]. Previous studies treated the exchange couplings as isotropic Heisenberg ($J_{xx}=J_{yy}=J_{zz}$) as a common assumption [14, 18, 45] [30]. In contrast, our non-collinear MFT makes no a priori assumptions and directly evaluates the full exchange tensor $\mathbf{J}_{ij}$ for each Ni–Ni and Ni–S interaction. For comparison with earlier works, we present the isotropic average $J=(J_{xx}+J_{yy}+J_{zz})/3$, summarized in Figure 3(e) (see Table S1). The third-nearest intra-layer Ni–Ni interaction at 6.7 Å (see also in Fig. S6) is found to be the dominant coupling, stabilizing the known zigzag ground state [29, 45, 46] [30]. In addition, our calculations reveal previously unexplored Ni–S interactions and long-range profiles. Off-diagonal components of $\mathbf{J}_{ij}$, including antisymmetric DM and symmetric Kitaev-$\Gamma$ terms, are also finite with a magnitude of approximately 0.01 meV (Fig. S5).

The magneto-crystalline anisotropy (MCA), corresponding to the first two terms of Eq. (2), is directly evaluated from relativistic DFT total energies. The results for rotating spin directions are shown in Figure 3(a,Figure 3b). As expected, the MCA is small, with a maximum of 8.5 $\mu$eV/Ni for the $c^*a$ rotation. The total energy profile identifies the $ab$-plane as the easy plane and the $c$-axis as the hard axis. A subtle feature appears for the $c^*a$ rotation, where the energy minimum ($\theta\approx82^\circ$) slightly deviates from $\theta=90^\circ$, reflecting the broken mirror and residual $C_2$ symmetry [40]. The $ab$- and $bc^*$-plane retain mirror symmetry with the sinusoidal dependence. Despite the small energy scale, the DFT results capture the essential symmetry features of the anisotropy, supporting their reliability.

Here, we pay special attention to the energetics revealed by MCA. The calculated energy profile distinguishes the $a$-axis from the $b$-axis as expected from the low-field torque data (Figure 2(a)). However, the high-field features exhibiting ‘shoulders’ cannot be reconciled with the energetic order $E_{\parallel b}<E_{\parallel a}<E_{\parallel c^*}$ (Figure 3(a,Figure 3b)). The torque computed from the corresponding Hamiltonian (first two terms in Eq. (2)) also fails to reproduce the experimental profile (Fig. S6(a)), highlighting the missing contribution of the dipolar interaction.

Magnetic dipole–dipole interaction is typically very weak in real materials and often neglected [47, 48, 49], though they can contribute to anisotropy in large moment systems [50, 51]. Incorporating this effect quantum-mechanically into the DFT scheme is known to be challenging despite recent progress [52, 53]. In NiPS$_3$, with $L=0$, we include this term classically using the Ni ($\sim1.5\ \mu_B$) and S ($\sim0.1\ \mu_B$) spin moments obtained from relativistic DFT. Figures 3(c,d) show the calculated dipolar energy with respect to the in-plane and out-of-plane spin rotation, respectively. Note that the energy minimum is found at $\theta\approx90^\circ$ (in $ca$-plane), consistent with the experimental easy-axis [54]. This dipolar interaction yields $E_{\parallel a}<E_{\parallel c^*}<E_{\parallel b}$.

In Figure 3(f), we summarize the magnetic total energy profile as a 3D contour map encompassing all anisotropy terms in Eq. (2). It is noted that the overall energetic order is now $E_{\parallel a}<E_{\parallel c^*}<E_{\parallel b}$ as indicated by the three dashed lines on the surface and also in Fig. S8. Importantly, this order matches the dipolar energy trend shown in Figure 3(c,Figure 3d), demonstrating that the dipolar interaction plays a pivotal role in determining the magnetic anisotropy of this material.

Using this final form of Hamiltonian, we recalculate torque as a function of field direction and strength. As mentioned above, the result shown in Figure 2(b) is in excellent agreement with the experiment in Figure 3(a). Hereby, we conclude that the dipole interaction is the key intrinsic source of magnetic anisotropy responsible for the observed high-field torque data. It is remarkable that the dipolar contribution plays such an important role that a qualitatively different signature is detected in the experiment.

Interestingly, the observed energy hierarchy, $E_{\parallel a}<E_{\parallel c^{*}}<E_{\parallel b}$, deviates from the easy-plane expectations of the XXZ model ($E_{\parallel a}\leq E_{\parallel b}<E_{\parallel c^{*}}$). XY or XXZ models assume spins confined within the $ab$-plane, allowing only weak in-plane anisotropy. In contrast, our results indicate that the spin degrees of freedom can no longer be considered strictly in-plane, as will become more evident in subsequent discussion of the spin-flop transition. Model simulations, performed with various sets of arbitrarily chosen parameters, consistently reproduce the qualitative shoulder feature as long as the energy hierarchy $E_{\parallel a}<E_{\parallel c^{*}}<E_{\parallel b}$ is retained (Fig. S9). In realistic parameter regimes, the dipolar interaction stabilizes this feature.

Let us now further scrutinize experimental data based on the Hamiltonian constructed above. Figure 4(a) shows the calculated torque as a function of field applied nearly along the $a$-axis ($\varphi=12^\circ$). While $\tau_{ab}$ gradually increases along with the field strength, it becomes notably steeper at around $H\approx8$ T up to the critical field $H_{\mathrm{SF}}=10.2$ T, and suddenly drops at $H\geq H_{\mathrm{SF}}$, indicating a spin-flop transition. It is clearly demonstrated with snapshot pictures in Figure 4(b) with moments (red and blue arrows) and torque (violet arrow). Below $H_{\mathrm{SF}}$, spins lie within the $ab$-plane (i.e., $\theta\approx90^\circ$; see also Figure 1(c)). As the field strength increases and reaches $H\approx H_{\mathrm{SF}}$, they stand up, while maintaining the same AFM spin order, and

make a sizable change in directions (12 T and 14 T). At the same time, the out-of-plane torque component is significantly reduced, consistent with the measured torque and responsible for the characteristic ‘shoulder’ features in Figure 2. For the more clear presentation, see Movie S1.

Finally, we consider the anisotropic PL response recently reported by Wang et al. [7], which has been proposed for coherent many-body exciton manipulation, opto-spintronic devices, and magnetic quantum information [6, 7, 9, 10, 21, 22, 55]. However, understanding this phenomenon remains unclear due to the lack of detailed theoretical analysis. It is clear that the observation indicates anisotropy beyond the conventional XY-type interactions [7, 9, 10, 14, 15]. Interestingly, the critical field for the anomalous PL response ($H_{\mathrm{PL}}=10.5$ T [15, 21, 22, 56]) is close to those observed in our magnetization ($H_{\mathrm{MM}\sim10}$ T) and torque ($H_{\tau}=7$–9 T) measurements.

Based on the Hamiltonian, we calculated the field-dependent response of the Néel vector, $\vec{L}=\vec{S}_{1}-\vec{S}_{2}$ (Figure 4(c)), a key quantity in the previous reports [7, 9]. The Néel vector carries crucial information for understanding the excitonic behavior. The spin-flop transition involves a sudden reorientation of $\vec{L}$. As shown in Figure 4(d), the out-of-plane angle $\theta$ rapidly changes from $90^\circ$ to $0^\circ$ around $H=10$–11 T, while the in-plane angle $\psi$ evolves continuously up to $45^\circ$. This reorientation of $\vec{L}$ accounts for the sudden change in PL, with the dipolar interaction playing a crucial role in shaping this characteristic response. We note that spin-flop is the key to understanding not only the torque anomaly but also the high-field magnetization behavior presented in Figure 1(d). Figure 1(e) presents the calculated magnetization based on the generalized Hamiltonian. It is reasonably well compared with experimental data by capturing the rapid up-turn along the $a$-axis and the linear behavior (no up-turn) along the $b$- and $c$-axes. The analysis shows that applying the high enough field along the $a$-axis flops the spins to align along the $c$-axis (as in Figure 4(b)), causing the rapid increase in magnetization.

Some discrepancies remain between experiment and calculation. The calculation gives the identical $b$- and $c^{*}$-axis magnetizations, whereas experimentally the $c$-axis magnetization is slightly larger. Similar to its sister compound FePS3 [54, 57], the weak vdW interaction allows the formation of 120°-rotated twinned domains along the $c^{*}$-axis [14]. As shown in the Supplemental Materials, coexistence of such twinned domains reproduces not only the weak linear increase of magnetization along the $a$-axis for $H < H_{\mathrm{SF}}$, resulting in $M_{\parallel b} < M_{\parallel c^{*}}$, but also the finite torque values at $\varphi= 0^\circ$ and 180°, consistent with experiment (see Figs. S10 and S11).

To unravel the unusual magnetic behaviors of the vdW magnet NiPS3, we combined torque magnetometry with first-principles calculations. The measured torque reveals anomalous high-field patterns. We analyzed these results by establishing the low-energy anisotropic spin Hamiltonian, incorporating Heisenberg, Kitaev-$\Gamma$, DM, single-ion anisotropy, and dipole interactions. The specific energetic order $E_{\parallel a} < E_{\parallel c^{*}} < E_{\parallel b}$ is crucial for explaining the characteristic shoulder pattern in the torque data, with the dipolar interaction playing a key role in stabilizing it. Furthermore, we identify the spin-flop transition within this energetic hierarchy as the unifying mechanism underlying multiple anomalous responses, including high-field magnetization as well as torque and the anisotropic PL reported previously. Our findings provide key insights into understanding magnetic anisotropy in vdW magnets NiPS3, highlighting the critical role of dipole interactions in shaping high-field responses. Given the fundamental importance of these interactions, our results offer a robust theoretical framework for understanding magnetic anisotropy, guiding future studies and potential applications to vdW magnets.

Acknowledgments

M. J. H. was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (Grant No. RS-2025-00559042). This work was supported by the Samsung Science & Technology Foundation (Grant No. SSTF-BA2101-05). Y. J. was supported by a National Research Foundation of Korea (NRF) grant (RS-2024-00337640, RS-2023-00301914). A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-2128556 and the State of Florida. This work was also supported by the Leading Researcher Program of the National Research Foundation of Korea (Grant No. RS-2020-NR049405).