Elucidating the Role of Dimensionality on the Electronic Structure of the Van der Waals Antiferromagnet NiPS3

Abstract

The sustained interest in investigating magnetism in the 2D limit of insulating antiferromagnets is driven by the possibilities of discovering, or engineering, novel magnetic phases through layer stacking. However, due to the difficulty of directly measuring magnetic interactions in 2D antiferromagnets, it is not yet understood how intralayer magnetic interactions in insulating, strongly correlated, materials can be modified through layer proximity. Herein, the impact of reduced dimensionality in the model van der Waals antiferromagnet NiPS₃ is explored by measuring electronic excitations in exfoliated samples using Resonant Inelastic X‐ray Scattering (RIXS). The resulting spectra shows systematic broadening of NiS₆ multiplet excitations with decreasing layer count from bulk down to three atomic layers (3L). It is shown that these trends originate from a decrease in transition metal‐ligand and ligand–ligand hopping integrals, and by charge‐transfer energy evolving from Δ = 0.83 eV in the bulk to 0.37 eV in 3L NiPS₃. Relevant intralayer magnetic exchange integrals computed from the electronic parameters exhibit a decrease in the average interaction strength with thickness. This study underscores the influence of interlayer electronic interactions on intralayer ones in insulating magnets, indicating that magnetic Hamiltonians in few‐layer insulating magnets can greatly deviate from their bulk counterparts.

Introduction

In 2D magnets, enhanced fluctuations and lattice connectivity strike a balance from which collective states unobtainable in three dimensions may emerge. The ability to prepare isolated monolayers of van der Waals (vdW) magnets has enabled access to new magnetic phases and tests of fundamental theorems of magnetism,^[1–6] and opens up possibilities for controlling or engineering unconventional states through stacking.^[7] Much work has concentrated on vdW materials with a net ferromagnetism in the 2D limit.^[5,8–10] However, antiferromagnets may offer more possibilities to explore complex magnetic order, topological spin textures, or quantum spin-liquids that arise from frustrated interactions and which are stabilized in 2D.^[11–15]

NiPS₃ stands as one of the few exfoliatable materials that exhibits both antiferromagnetic order and strong correlations.^[16,17] Recent Raman scattering measurements suggest that the magnetic order in NiPS₃ is highly sensitive to dimensionality and find that long-range order vanishes in the monolayer limit in favor of a fluctuating magnetic phase.^[3] Based on the magnetic Hamiltonian that was determined by inelastic neutron scattering on bulk samples,^[18] the thickness-dependent Raman data were associated with the proliferation of vorticies through a Berezinskii–Kosterlitz–Thouless phase transition in the 2D material. However, this explanation assumes that the few-layer magnetic Hamiltonian is identical to that in the bulk. More direct experimental access to the electronic energy scales and magnetic interactions is thus necessary to resolve the nature of the magnetic state in exfoliated NiPS₃.

In this work, we demonstrate that Ni-S electronic energy scales are strongly altered by dimensionality in NiPS₃ leading to a few-layer magnetic Hamiltonian that differs from that of the bulk. We use resonant Inelastic X-ray Scattering (RIXS) to directly interrogate correlated electronic states in exfoliated flakes of NiPS₃ and reveal a systematic softening and broadening of NiS₆ multiplet excitations with decreasing thickness that is reproduced by a multiplet ligand-field model. Overall, we find that decreased hopping integrals and charge transfer energy in 2D result in a more covalent character for Ni d-orbitals (Figure 1a). We additionally compute the relevant magnetic exchange integrals and find a decrease in the second- and third-nearest neighbor magnetic interaction strengths, and an increase in the first-nearest neighbor interaction strength. This change moves NiPS₃ closer to the boundary between the stripy antiferromagnetic and spiral ordered phases of the honeycomb antiferromagnet. Lastly, we show how the change of electronic energy scales in thinner samples occurs due to decreased electronic vdW delocalization across layers in the 2D limit. Since this mechanism is not necessarily specific to NiPS₃, its effect will be important to the properties of a broad class of strongly correlated few-layer vdW magnets.

[IMAGE OMITTED. SEE PDF]

Experimental Results

NiPS₃ crystallizes in the monoclinic C2/m space group with Ni (S = 1) atoms on honeycomb lattices in the ab-planes. Each Ni atom is octahedrally coordinated by six S atoms, with P atoms situated between the 2D sheets of NiS₆ that form the vdW gap along the c-axis with an interlayer spacing of ≈ 0.636 nm.^[2,16] NiPS₃ is known to magnetically order at a transition temperature of T_N = 155 K; ferromagnetic zig-zag chains form with moments parallel to the a-axis and antiferromagnetically coupled along the b-axis, with a small out of plane component.^[16,21]

Figure 1c,e show representative RIXS spectra for bulk and three-layer (3L) NiPS₃, respectively. Spectra were collected at the peak of the Ni L₃-edge XAS E_i = 853 eV, corresponding to 2p_3/2 to 3d electronic transitions. We concentrate on the low energy region E_loss = 0.2 → 2.15 eV that contains excitations within the NiS₆ multiplet. The bulk and three-layer (3L) spectra are qualitatively similar except for an overall energy broadening and softening that is readily visible in the 3L data (Figure 1b). While the qualitative similarity between bulk and 3L spectra is consistent with the fact that there are no drastic structural reconstructions upon exfoliation, the apparent broadening and softening indicates a change in the electronic structure of NiPS₃ with thickness.

In order to elucidate the origin of this change, we first concentrate our analysis on the bulk spectra and identify all relevant features. We found that a minimum of six Gaussian modes were required to fit the bulk data, labeled A - F in Figure 1c, with minimum energy widths fixed by the experimental energy resolution (≈ 55meV FWHM). Each of these features can be identified as an excitation within the electronic multiplet of trigonally distorted NiS₆ octahedra (Figure 2a; Figure SV, Supporting Information). The center of mass positions of the two bulk peaks are assigned to the t_2g → e_g (d-d) excitations of ³T_2g and ³T_1g symmetry, respectively, in good agreement with optical measurements.^{[17, 22]} The trigonal distortion introduces a D_3d symmetry which splits ³T_2g → ³A_1g + ³E_g (peaks B & C) and ³T_1g → ³A_2g + ³E_g (peaks E & F), in agreement with Raman and optical measurements.^[22,23] Access to spin-flip (ΔS ≠ 0) excitations in the RIXS cross-section leads us to assign peak D ¹E_g symmetry as the next highest excited state above ³T_2g in a 3d⁸ system. Lastly, peak A is assigned to a charge transfer excitation with $$3d^{9}\underline{L}^{1}$$ character, where $$\underline{L}^n$$ denotes n ligand holes. The 800 meV energy scale of this peak indicates a small charge transfer energy in NiPS₃. We verify these peak assignments through the application of the NiS₆ multiplet ligand-field model described below. We note that since our incident energy was tuned to the peak of the Ni L₃-edge XAS, our measurements were not sensitive to the sharp 1.47 eV peak reported in Ref. [24].

[IMAGE OMITTED. SEE PDF]

We now bring our attention to the spectra collected on the exfoliated samples, and in particular, the 3L sample. Following the fitting procedure from the bulk spectra, we again fit six Gaussian modes; however, Empirically fitting the 3L spectra to a minimum of six Gaussian peaks resulted in two scenarios of equally good fit quality. In order to address the width of peak A, we explore the two fitting scenarios for the 3L spectra. In scenario one, the widths of all peaks were held fixed at the experimental resolution; this fit converged with a systematic softening of all peaks between the bulk and 3L data sets. In scenario two, the FWHM of peak A was allowed to relax, while peaks B - D were fixed to the experimental resolution, and peaks E - F were fixed to the bulk fitted FWHM to experimental resolution ratio; this fit converged with minimal softening of all peaks, but systematic broadening and increased intensity attributed to peak A. We found that peak energies extracted from scenario one could only be reproduced within physically meaningful parameters using a negative charge transfer energy, while scenario two is reproduced with a small positive charge transfer energy.^[25] A negative charge transfer energy for the 3L sample implies a zero-crossing of the charge transfer energy as a function of thickness between bulk and 2D exfoliated samples. We rule out such a transition based on the smooth evolution of thickness dependent RIXS data and Raman spectra.^[23] The full details of each fitting scenario to the 3L spectrum, including details of the bulk spectrum fitting, can be found in the Supporting Information.^[25]

Fits for scenario two are shown in Figure 1e, while Figure 1d summarizes the centroids of the fitted peaks for the various sample thicknesses measured. Minimal differences were found between the bulk and 60L, placing a lower limit on bulk behavior for exfoliated NiPS₃ at ≈ 38 nm. From bulk to 3L, peaks B - F vary slightly in E_loss. We find a broadening in the FWHM of peak A by 360(50) meV over the bulk data that suggests a change in the charge transfer energy in few-layer samples. Furthermore, the observed systematic broadening of excitations signifies an electronic structure intricately connected to sample thickness in NiPS₃.

Theoretical Analysis

Having identified a clear empirical trend, we perform a thorough fit of the experimental data for paramagnetic NiPS₃ using a high-fidelity multiplet ligand-field model (MLFM) in the basis of symmetry-adapted linear combinations of ligand orbitals.^[26] Crucially, we use rigorously converged parameter-free solutions of the Schrödinger equation to physically guide our search of the MLFM parameter space. Our MLFM model includes Slater-Condon parameters($$F^{0,2,4}_{dd}$$, $$F^{0,2}_{pd}$$, $$G^{1,3}_{pd}$$), covalent hopping integrals between S 3p- and Ni 3d- orbitals, pdσ and pdπ, S 3p- orbital level splitting, T_pp = ppσ − ppπ, cubic crystal field (10Dq) and trigonal distortion δ, Ni 3d-3d and 2p-3d on-site Coulomb interactions U_dd and U_pd = 1.2U_dd, and charge transfer energy Δ. The values of the transition metal-ligand (pdσ, pdπ) and ligand–ligand (ppσ, ppπ) hopping integrals for bulk and monolayer geometries are calculated independently by self-consistently converging to the electronic ground state of the nonmagnetic configuration of NiPS₃ using Density Functional Theory (DFT). We employ a substantial plane wave energy cutoff of more than 3,400 eV, sample the Brillouin zone such that all energies are converged to within less than 0.001 eV, and use correlation-consistent pseudopotentials obtained from fully many-body ab initio calculations.^[38] We then generate maximally-localized Wannier functions (MLWF) which span a substantial subspace of the rigorously converged DFT ground states, and use the resulting tight-binding inter-orbital expectation values to solve for the hopping integrals within the two-center linear combination of atomic orbitals approximation.^[25]

We carried out a search of the remaining parameter space for Δ, 10Dq, and $$F^{0,2,4}_{dd}$$ by minimizing the difference between calculated energies peaks A - F while keeping F(G)_pd fixed to 80% of their atomic Hartree-Fock values.^[25,27]. For initial comparisons of this model to our data, we used physically meaningful parameters for octahedrally coordinated NiS₆.^[28–30] Figure 2a,c shows the calculated energy levels for a NiS₆ cluster as a function of the charge transfer energy Δ using fixed parameters that give the best agreement between measured and calculated peak energies for bulk and 3L NiPS₃, respectively.^[25] Previous optical and X-ray absorption (XAS) studies classified NiPS₃ as a negative charge transfer insulator,^[17] while more recent RIXS and XAS measurements indicate a positive charge transfer gap.^[24,31] We find that, for bulk NiPS₃, a small positive Δ = 0.83 eV was necessary to give an accurate match to the data; for 3L NiPS₃, the charge transfer gap decreases, and accordingly, the best match to the data is obtained with Δ = 0.37 eV.

In Figure 3a,b, we show RIXS spectra calculated using the open-source toolkit EDRIXS^[32] compared to the experimental data. Intensities were normalized to the nominal 10Dq line. Covalent hopping integrals, pdσ and pdπ, as well as T_pp were fixed to those obtained from ab initio calculations of the nonmagnetic configuration (Table 1) while Δ was allowed to vary. To account for broadening of excitations not captured by our multiplet ligand-field model and facilitate better comparison with experimental data, the calculated spectra were broadened by increasing the final-state lifetime above 2 eV in E_loss. We can reproduce the observed broadening and softening of NiS₆ multiplet excitations in the 3L spectrum by a decrease in charge transfer energy and transition metal-ligand hopping integrals, as parameterized by Δ/T_pd, and Δ/U, where T_pd = $$-\sqrt {3}pd\sigma$$ and on-site 3d Coulomb repulsion $$U = F^0_{dd} + \frac{4}{49}(F^2_{dd} + F^4_{dd})$$^[26,33]: Δ/T_pd = 0.45 and Δ/U = 0.10 in bulk, and Δ/T_pd = 0.23 and Δ/U = 0.04 in 3L.

Table 1 Fixed values in eV of hopping integrals extracted from ab initio calculations of the nonmagnetic configuration. Charge transfer energy Δ, and intra-orbital Coulomb repulsion $$U = F^0_{dd} + \frac{4}{49}(F^2_{dd} + F^4_{dd})$$ extracted from RIXS modeling.^[26,33]

	pdσ	pdπ	ppσ	ppπ	Tpp	Δ	U
Bulk	−1.07	0.67	0.89	−0.09	0.98	0.83	8.3
3L	−0.93	0.46	0.62	−0.01	0.63	0.37	8.3

[IMAGE OMITTED. SEE PDF]

We determine that the underlying mechanism responsible for the significant change in the RIXS signal with thickness is predominantly electronic rather than structural in origin, though the lattice constant is slightly overestimated in the PBE-optimized monolayer. We find that as NiPS₃ gets thinner, metal-ligand π-hopping is reduced (pdπ decreases) due to the removal of π-like interlayer vdW interactions. The same effect also causes pdσ and T_pp to change significantly because of the mixed σ- and π-bonding character present in the sp₃-hybridized phosphorus atoms that bridge the NiS₆ clusters. In the context of our MLWFs, this is reflected in a change in the largest tight-binding energies used to solve for pdσ and pdπ (Figure 1b).^[25]

The combination of RIXS measurements and ab initio calculations constrain the electronic ground state that underlies the magnetic properties of NiPS₃. In Figure 4c, we investigate the change in the ground state character of $$\ket {\Psi _g} = \alpha \ket {3d^8} + \beta \ket {3d^9\underline{L}^1} + \gamma \ket {3d^{10}\underline{L}^2}$$ extracted from our multiplet ligand-field model as a function of T_pp and Δ/T_pd. We find a negligible contribution from the $$\ket {3d^{10}\underline{L}^2}$$ state and nearly equal populations of the $$\ket {3d^8}$$ and $$\ket {3d^9\underline{L}^1}$$ states. In bulk, $$\lvert {\alpha }\rvert^2/\,\lvert {\beta }\rvert^2 = 1.18$$, and in 3L, $$\lvert {\alpha }\rvert^2/\,\lvert {\beta }\rvert^2 = 1.3$$, implying a small increase in the magnitude of the paramagnetic Ni moment. As shown in Figure 4c, this minimal change is accounted for by the dependence of the ground state character on both the hybridization Δ/T_pd and T_pp. While an increased transition metal-ligand hybridization tends to enhance the $$\ket {3d^9\underline{L}^1}$$ character, this is offset by the reduction in the ligand–ligand hybridization T_pp.

[IMAGE OMITTED. SEE PDF]

Despite the small change in ground state character, the large change in hybridization and hopping parameters influences the magnetic exchange interactions. Using the parameters obtained from our RIXS measurements and ab initio modeling, we compute the superexchange interactions up to the third-nearest neighbor within a sixth order cell-perturbation.^[30,34,35] The first J₁, second J₂, and third J₃ nearest neighbor expressions are given by the second order perturbation terms for the $$\ket {3d^9\underline{L}^1}$$ states, and fourth and sixth order terms for the $$\ket {3d^8}$$ states;^[30] Figure 4a shows the first-nearest neighbor superexchange pathway for the $$\ket {3d^8}$$ configuration, while a detailed description of all these expressions is given in the Supporting Information.^[25] In bulk NiPS₃, we find $$J_1^{B}\!\approx \!-3.6$$ meV, $$J_2^{B}\!\approx \!0.17$$ meV, and $$J_3^{B}\!\approx \!12$$ meV, in excellent agreement with recently reported values from inelastic neutron scattering.^[18,21] The decrease in T_pp and Δ/T_pd leads to an overall enhancement of $$J_1^{3L}\!\approx \!-4.2$$ meV, a vanishing $$J^{3L}_{2}$$, and decrease in $$J_3^{3L}\!\approx \!7.3$$ meV. In Figure 4b, we summarize the dependence of J₃/J₁ on T_pp and Δ/T_pd. We find that J₁ is dominated by the $$\ket {3d^9\underline{L}^1}$$ state, while J₂ and J₃ are dominated by the $$\ket {3d^8}$$ state. Thus, the decrease in T_pp is directly responsible for an increased $$\ket {3d^9\underline{L}^1}$$ contribution to J₁. As a consequence of the overall reduction in the average exchange interaction strength, the magnetic transition temperature is expected to be reduced in few-layer samples compared to bulk samples. Furthermore, the decrease in J₃/J₁ from -3.3 in bulk to -1.7 in 3L, positions 3L NiPS₃ closer to a phase boundary between the stripy AFM phase and a spiral ordered phase.^[36] It is likely that, in the 2D limit, NiPS₃ is driven into a highly frustrated regime on this phase boundary.

Conclusion

In summary, we used RIXS to access the electronic ground state properties of an exfoliated, correlated antiferromagnet in the 2D limit. We found that the electronic energy scales associated with Ni-S hybridization, and consequently the magnetic exchange interactions, are altered in a non-trivial way though the modification of interlayer energy scales upon exfoliation of NiPS₃ despite minimal structural changes. Our findings demonstrate that magnetic exchange parameters determined from measurements on bulk materials are not applicable in the 2D limit, as interlayer interactions, absent in 2D, affect intralayer ones. The underlying electronic mechanism we have identified points to the possibility of controlling magnetic interactions in strongly correlated van der Waals heterostructures by tuning interfacial energy scales toward the design of the next generation of 2D strongly correlated magnetic materials.

Experimental Section

Crystal Growth

Single crystal samples of NiPS₃ were grown by vapor transport, following previously published methods.^[16,18] Stoichiometric quantities of metallic Nickel (99.994% purity), crystalline Phorphorus (99.999% purity) and Sulfur (99.999% purity), were added to a quartz tube in an argon-filled glovebox. The quartz tube was then evacuated, sealed, and placed in a two-zone furnace and heated to 700 °C/750 °C over 6 h; the furnace was held at this temperature gradient for 2 days, and cooled to 670 °C/620 °C over 8 h and held for an additional 16 days. After a total of 18 days in the furnace, samples were cooled to room temperature over 8 h. The crystals formed were shiny gray metallic and had hexagonal motif characteristic of NiPS₃.

Exfoliated Sample Preparation

Bulk NiPS₃ was exfoliated using conventional scotch-tape methods^[2] and deposited either onto a blank SiO₂ substrate or onto a SiO₂ substrate pre-treated with a patterned Copper (Cu) grid. Sample layer count was determined using a combination of atomic-force microscope thickness measurements and optical contrast with SiO₂ substrate. See Supporting Information and Figure SII (Supporting Information) for more details.^[25] The patterned Cu grid consisted of 100 µ m × µ m SiO₂ cells separated by 200 µ m of 50 nm thickness Cu (Figure SIb,c, Supporting Information).^[25] The 3L sample was deposited onto a blank SiO₂ substrate and was later patterned with a Cu fiducial marker using electron-beam lithography (Figure SIa, Supporting Information), again with a Cu thickness of 50 nm.^[25] In both cases, Cu was chosen as a material that could provide a fluorescence contrast to SiO₂ in the soft X-ray regime. This fluorescence contrast proved invaluable in locating small samples whose signals were weak under an X-ray beam. The Cu grid proved a useful method for locating sample(s) as a unique grid scheme could be defined for each chip if the orientation of each chip remained consistent; however, the fiducial marker had the advantage of of being visible by eye, resulting in unequivocal sample location and removing the requirement of a grid scheme. Exfoliated samples were spin coated with a PMMA protective layer and stored in an Ar atmosphere to prevent degradation.

RIXS Measurements

Room temperature RIXS measurements on exfoliated flakes were carried out on the PEAXIS beamline at BESSY II.^[37] The PMMA coatings were removed immediately prior to loading the samples into the RIXS vacuum chamber via washing with acetone and isopropyl alcohol. All samples chosen for measurement had lateral dimensions larger than, or comparable to, the 15 (H) µ m x 4 (V) µ m x-ray beam spot size at the PEAXIS beamline. A horizontal scattering geometry of 2θ = 90° was used with an ≈ 235 meV energy resolution (full width at half max, FWHM) using linear horizontal polarization and specular geometry. Spectra were collected in 30 min segments to minimize sample exposure to the X-ray beam. It was noted here that the 7 L spectra only was collected with marginally better instrumental resolution of ≈ 190 meV. Bulk NiPS₃ RIXS measurements were carried out on the SIX beamline at NSLS II with identical scattering geometry, but with an ≈ 55 meV energy resolution FWHM. The main resonance peak at the Ni L₃-edge was chosen via X-ray absorption spectroscopy (XAS) preformed on-site, prior to each RIXS measurements.

Acknowledgements

The authors thank Mark Dean for helpful discussions and comments on this manuscript. The authors also thank Naiyuan J. Zhang and Erin Morissette for their guidance and consultation on pattern fabrication. M.F.D. and K.W.P. were supported by the National Science Foundation under grant no. OMA-1936221. A.D.L.T. was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award Number DE-SC0021. D.S. and B.R. were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and the Center for Predictive Simulation of Functional Materials, while A.L. was supported by the Brown University Diversity Fellowship. RIXS measurements were carried out at the U41-PEAXIS beamline at the BESSY II electron storage ring operated by the Helmholtz-Zentrum Berlin für Materialien und Energie. This research was conducted using the computational resources of the National Energy Research Scientific Computing Center under Contract No. DE- AC0205CH11231, which is a U.S. Department of Energy Office of Science User Facility, as well as resources and services at the Center for Computation and Visualization, Brown University.

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data that support the findings of this study are openly available in Zenodo at , reference number 10125572.

References

J.‐G. Park, J. Phys.: Condens. Matter 2016, 28, [eLocator: 301001].

C.‐T. Kuo, M. Neumann, K. Balamurugan, H. J. Park, S. Kang, H. W. Shiu, J. H. Kang, B. H. Hong, M. Han, T. W. Noh, J.‐G. Park, Sci. Rep. 2016, 6, [eLocator: 20904].

K. Kim, S. Y. Lim, J.‐U. Lee, S. Lee, T. Y. Kim, K. Park, G. S. Jeon, C.‐H. Park, J.‐G. Park, H. Cheong, Nat. Commun. 2019, 10, 345.

N. Sivadas, S. Okamoto, X. Xu, Craig. J. Fennie, D. Xiao, Nano Lett. 2018, 18, 7658.

C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, X. Zhang, Nature 2017, 546, 265.

Y. Tang, K. Su, L. Li, Y. Xu, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, C.‐M. Jian, C. Xu, K. F. Mak, J. Shan, Nat. Nanotechnol. 2023, 18, 233.

F. Hellman, A. Hoffmann, Y. Tserkovnyak, G. S. D. Beach, E. E. Fullerton, C. Leighton, A. H. MacDonald, D. C. Ralph, D. A. Arena, H. A. Dürr, P. Fischer, J. Grollier, J. P. Heremans, T. Jungwirth, A. V. Kimel, B. Koopmans, I. N. Krivorotov, S. J. May, A. K. Petford‐Long, J. M. Rondinelli, N. Samarth, I. K. Schuller, A. N. Slavin, M. D. Stiles, O. Tchernyshyov, A. Thiaville, B. L. Zink, Rev. Mod. Phys. 2017, 89, [eLocator: 025006].

B. Huang, G. Clark, E. Navarro‐Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo‐Herrero, X. Xu, Nature 2017, 546, 270.

K. S. Burch, D. Mandrus, J.‐G. Park, Nature 2018, 563, 47.

M. Gibertini, M. Koperski, A. F. Morpurgo, K. S. Novoselov, Nat. Nanotechnol. 2019, 14, 408.

A. Kitaev, Annals of Physics 2006, 321, 2.

K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.‐Y. Kee, Y.‐J. Kim, Phys. Rev. B 2014, 90, [eLocator: 041112(R)].

J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, R. Moessner, Nature Phys 2016, 12, 912.

K. H. Lee, S. B. Chung, K. Park, J.‐G. Park, Phys. Rev. B 2018, 97, [eLocator: 180401].

B. Zhou, Y. Wang, G. B. Osterhoudt, P. Lampen‐Kelley, D. Mandrus, R. He, K. S. Burch, E. A. Henriksen, J. Phys. Chem. Solids 2019, 128, 291.

A. R. Wildes, V. Simonet, E. Ressouche, G. J. McIntyre, M. Avdeev, E. Suard, S. A. J. Kimber, D. Lançon, G. Pepe, B. Moubaraki, T. J. Hicks, Phys. Rev. B 2015, 92, [eLocator: 224408].

S. Y. Kim, T. Y. Kim, L. J. Sandilands, S. Sinn, M.‐C. Lee, J. Son, S. Lee, K.‐Y. Choi, W. Kim, B.‐G. Park, C. Jeon, H.‐D. Kim, C.‐H. Park, J.‐G. Park, S. J. Moon, T. W. Noh, Phys. Rev. Lett. 2018, 120, [eLocator: 136402].

D. Lançon, R. A. Ewings, T. Guidi, F. Formisano, A. R. Wildes, Phys. Rev. B 2018, 98, [eLocator: 134414].

K. Momma, F. Izumi, J. Appl. Cryst. 2011, 44, 1272.

B. O. Community, Blender ‐ a 3D modelling and rendering package, Blender Foundation, Stichting Blender Foundation, Amsterdam 2018.

A. R. Wildes, J. R. Stewart, M. D. Le, R. A. Ewings, K. C. Rule, G. Deng, K. Anand, Phys. Rev. B 2022, 106, [eLocator: 174422].

D. Afanasiev, J. R. Hortensius, M. Matthiesen, S. Mañas‐Valero, M. Šiškins, M. Lee, E. Lesne, H. S. J. van der Zant, P. G. Steeneken, B. A. Ivanov, E. Coronado, A. D. Caviglia, Sci. Adv. 2021, 7, [eLocator: eabf3096].

X. Wang, J. Cao, H. Li, Z. Lu, A. Cohen, A. Haldar, H. Kitadai, Q. Tan, K. S. Burch, D. Smirnov, W. Xu, S. Sharifzadeh, L. Liang, X. Ling, Sci. Adv. 2022, 8, [eLocator: eabl7707].

S. Kang, K. Kim, B. H. Kim, J. Kim, K. I. Sim, J.‐U. Lee, S. Lee, K. Park, S. Yun, T. Kim, A. Nag, A. Walters, M. Garcia‐Fernandez, J. Li, L. Chapon, K.‐J. Zhou, Y.‐W. Son, J. H. Kim, H. Cheong, J.‐G. Park, Nature 2020, 583, 785.

See supporting information for details on sample preparation and preservation, additional details about the exact diagonalization calculations as well as further information about the analysis and fitting of the rixs spectra, details of the ab initio calculation methodology in obtaining mlwfs, and tm‐l and l‐l hopping parameters, and detailed description of the superexchange expressions.

M. W. Haverkort, M. Zwierzycki, O. K. Andersen, Phys. Rev. B 2012, 85, [eLocator: 165113].

G. Ghiringhelli, M. Matsubara, C. Dallera, F. Fracassi, R. Gusmeroli, A. Piazzalunga, A. Tagliaferri, N. B. Brookes, A. Kotani, L. Braicovich, J. Phys.: Condens. Matter 2005, 17, 5397.

A. E. Bocquet, T. Mizokawa, T. Saitoh, H. Namatame, A. Fujimori, Phys. Rev. B 1992, 46, 3771.

S. R. Krishnakumar, D. D. Sarma, Phys. Rev. B 2003, 68, [eLocator: 155110].

K. Takubo, T. Mizokawa, J.‐Y. Son, T. Nambu, S. Nakatsuji, Y. Maeno, Phys. Rev. Lett. 2007, 99, [eLocator: 037203].

M. Yan, Y. Jin, Z. Wu, A. Tsaturyan, A. Makarova, D. Smirnov, E. Voloshina, Y. Dedkov, J. Phys. Chem. Lett. 2021, 12, 2400.

Y. Wang, G. Fabbris, M. P. M. Dean, G. Kotliar, Comput. Phys. Commun. 2019, 243, 151.

Y. Shen, J. Sears, G. Fabbris, J. Li, J. Pelliciari, I. Jarrige, X. He, I. Bozovic, M. Mitrano, J. Zhang, J. F. Mitchell, A. S. Botana, V. Bisogni, M. R. Norman, S. Johnston, M. P. M. Dean, Phys. Rev. X 2022, 12, [eLocator: 011055].

J. H. Jefferson, H. Eskes, L. F. Feiner, Phys. Rev. B 1992, 45, 7959.

H. Eskes, J. H. Jefferson, Phys. Rev. B 1993, 48, 9788.

J. Fouet, P. Sindzingre, C. Lhuillier, Eur. Phys. J. B 2001, 20, 241.

C. Schulz, K. Lieutenant, J. Xiao, T. Hofmann, D. Wong, K. Habicht, J. Synchrotron Rad. 2020, 27, 238.

A. Annaberdiyev, G. Wang, C. A. Melton, M. C. Bennett, L. Shulenburger, L. Mitas, The Journal of Chemical Physics, 2018 149. [DOI: https://dx.doi.org/10.1063/1.5040472]