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

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 NiPS3 is explored by measuring electronic excitations in exfoliated samples using Resonant Inelastic X‐ray Scattering (RIXS). The resulting spectra shows systematic broadening of NiS6 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 NiPS3. 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.

In two-dimensional 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,[3][4][5][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 two-dimensional (2D) limit [5,[8][9][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 are stabilized in 2D [11][12][13][14][15].
NiPS 3 stands as one of the few exfoliatable materials that exhibits both antiferromagnetic order and strong correlations [1,17].Recent Raman scattering measurements suggest the magnetic order in NiPS 3 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 [2], the thickness-dependent Raman data were associated with the proliferation of vorticies through a Berezinskii-Kosterlitz-Thoules 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 necessary to resolve the nature of the magnetic state in exfoliated NiPS 3 .6 S 3d 9  Ni 3d 8  Ni 3d 9  Ni FIG. 1.(a) Schematic electronic density of states from bulk to 3L NiPS3.Ni L3 RIXS spectra for bulk (b) and 3L (d) NiPS3.Black points are experimental data with errorbars smaller than the symbol size, red lines show Gaussian fit to the data.(c) Gaussian peak position versus layer count.Dashed linear fits highlight the electronic structure change with thickness; gray regions follow the numerical center of mass and full width half max of the overall experimental peaks.
In this letter, we show that Ni-S electronic energy scales are strongly altered by dimensionality in NiPS 3 arXiv:2302.07910v1[cond-mat.str-el]15 Feb 2023 FIG. 2. Calculated energy levels of the NiS6 multiplet ligandfield model as a function of ∆ (a) compared to experimental bulk spectrum (b).Fixed model parameters are listed in Table I and Ref. [19].In (a), black and gray lines show energy levels calculated with and without SOC respectively.Symmetry labels adapted from a calculation without SOC.Horizontal dashed lines and shaded regions show E loss value of fitted peaks in rightmost panel.Vertical dashed line indicates best fit at ∆ = 0.60 eV.
and thus, the few-layer magnetic Hamiltonian differs from that of the bulk.Resonant Inelastic X-ray Scattering (RIXS) on exfoliated flakes reveals a systematic softening and broadening of NiS 6 multiplet excitations with decreasing thickness that is reproduced by a multiplet ligand-field model.Decreased hopping integrals and charge transfer energy in 2D result in a more covalent character [Fig.1(a)].We compute the relevant magnetic exchange integrals and find a systematic 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 3 closer to the boundary between the stripy antiferromagnetic and spiral ordered phases of the honeycomb antiferromagnet.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 specific to NiPS 3 , it's effect will be important to the properties of a broad class of few-layer insulating vdW magnets.
Single crystals of NiPS 3 were grown via standard vapor transport methods [1,2].Bulk NiPS 3 was exfoliated in air using conventional scotch-tape methods [3] and deposited either onto a blank SiO 2 substrate or onto a SiO 2 substrate pre-treated with a patterned Cu grid.Samples deposited onto blank SiO 2 substrates were later patterned with Cu fiducial markers using electron-beam lithography [19].Exfoliated samples were spin coated with a PMMA protective layer and stored in an Ar atmosphere to prevent degradation.The PMMA coatings were removed immediately prior to loading the samples into the RIXS vacuum chamber via washing with acetone and isopropyl alcohol.Room temperature RIXS measurements were carried out on the PEAXIS beamline at BESSY II [20].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 minute segments to minimize sample exposure to the X-ray beam.Bulk NiPS 3 spectra were collected with an identical scattering geometry, but with an ≈ 177 meV overall energy resolution FWHM.
Figs. 1(a) & (b) show representative RIXS spectra for bulk and three-layer (3L) NiPS 3 , respectively.Spectra were collected at the peak of the Ni L 3 -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 6 multiplet.The bulk and 3-layer (3L) spectra are qualitatively similar except for a systematic overall energy broadening and softening that is readily visible in the 3L data [Fig. 1 (b)].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 3 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 six resolution-limited Gaussian modes were required to fit the bulk data, labeled A -F [Fig.1(b)].Each of these features can be identified as an excitation within the electronic mulitplet on the slightly trigonally distorted NiS 6 octahedra.The center of mass positions of the two bulk peaks are assigned to the t 2g → e g (d-d ) excitations of 3 T 2g and 3 T 1g symmetry, respectively, in good agreement with optical measurements [17,21].The trigonal distortion introduces a D 3d symmetry which splits 3 in agreement with Raman and optical measurements [21,22].Access to spin-flip (∆S = 0) excitations in the RIXS cross-section leads us to assign peak D 1 E g symmetry as the next highest excited state above 3 T 2g in a 3d 8 system.Peak A is assigned to a charge transfer excitation with 3d 9 L 1 character, where L n denotes n ligand holes.The 800 meV energy scale of this peak indicates a small charge transfer energy in NiPS 3 .We verify these peak assignments through the application of the NiS 6 multiplet ligand-field model described below.We note that since our incident energy was tuned to the peak of the Ni L 3 -edge XAS, our measurements were not sensitive to the sharp 1.47 eV peak reported in Ref. [23].
In order to facilitate an efficient exploration of the pa-  ) and E loss .∆/T pd and ∆/U parameterize the hybridization and charge transfer characters respectively.
rameter space, we utilize a basis of symmetry-adapted linear combinations of ligand orbitals [7] within a multiplet ligand-field model.Our 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 ∆.For initial comparisons of this model to our data, we used physically meaningful parameters for octahedrally coordinated NiS 6 [25][26][27].We carried out a search of the 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 [6,19].Fig. 2 shows the calculated energy levels for a NiS 6 cluster as a function of the charge transfer energy ∆ for fixed parameters that give the best agreement between measured and calculated peak energies [19].Previous optical and X-ray absorption (XAS) studies classified NiPS 3 as a negative charge transfer insulator [17], while more recent RIXS and XAS measurements indicate a positive charge transfer gap [23,29].We find that a small positive ∆ = 0.60 eV was necessary to give an accurate match to the data.
We now bring our attention the 3L sample.Empirically fitting the 3L spectra to a minimum of six Gaussian peaks resulted in two scenarios of equally good fit quality.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, all peak widths were allowed to relax; this fit converged with minimal softening of all peaks, but systematic broadening and increased spectral weight 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 [19].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 [22].
Fits for scenario two are shown in Fig. 1(d), while Fig. 1(c) 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 3 at ≈ 38 nm.From bulk to 3L, the 3 T 2g modes ( 3 A 1g + 3 E g ) shift slightly upward in E loss , while 1 E g has the largest upward shift of ∆E loss = +149 (24) meV; the 3 T 1g modes ( 3 A 2g + 3 E g ) also shift in E loss .However, both the 3 A 2g and 3 E g modes become mixed with higher energy 1 T 2g ( 1 D) modes, split by D 3d symmetry, in few-layer samples [19].We find a 161 (10)% increase in the FWHM of peak A over the bulk data, suggesting a change in the charge transfer energy in few-layer samples.This observed systematic softening and broadening of excitations signifies an electronic structure intricately connected to sample thickness in NiPS 3 .
Having identified a clear empirical trend, we then turned to ab initio calculations for further insight.We used Density Functional Theory (DFT) to converge the electronic ground state of NiPS 3 in both the bulk and monolayer geometries, generated maximallylocalized Wannier functions (MLWF) which spanned the ground state DFT subspace, and used the corresponding tight-binding energy cross-terms to solve for the hopping integrals within the two-center approximation [19].These calculations reduced ambiguities in our parameter assignments by directly providing physically-grounded constraints on our fits.In Figs. 3 (a) & (b), we show RIXS spectra calculated using the open-source toolkit EDRIXS [31] 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 I] 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 6 multiplet excitations with thickness by a decrease in charge transfer energy and transition metal-ligand hopping integrals, as parameterized by ∆/T pd , and ∆/U , where T pd = − √ 3pdσ and on-site 3d Coulomb repulsion U = F 0 dd + 4 49 (F 2 dd + F 4 dd ) [7,30]: ∆/T pd = 0.32 and ∆/U = 0.07 in bulk, and ∆/T pd = 0.14 and ∆/U = 0.03 in 3L.
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 3 gets thinner, metal-ligand π-hopping is reduced (pdπ decreases) due to the removal of stabilizing, π-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 3 -hybridized phosphorus atoms that bridge the NiS 6 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π [19].
The combination of RIXS measurements and ab initio calculations constrain the electronic ground state that underlies the magnetic properties of NiPS 3 .In Fig. 4(a), we investigate the change in the ground state character of |Ψ g = α 3d 8 + β 3d 9 L 1 + γ 3d 10 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 3d 10 L 2 state and nearly equal populations of the 3d 8 and 3d 9 L 1 states.In bulk, |α| 2 / |β| 2 = 1.08, and in 3L, |α| 2 / |β| 2 = 1.2, implying a small increase in the magnitude of the paramagnetic Ni moment.As shown in Fig. 4(a), 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 3d 9 L 1 character, this is offset by the reduction in the ligand-ligand hybridization T pp .
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 [27,32,33].The first J 1 , second J 2 , and third J 3 nearest neighbor expressions are given by the second order perturbation terms for the 3d 9 L 1 states, and fourth and sixth order terms for the 3d 8 states [27]; a detailed description of these expressions is given in Ref. [19].In bulk NiPS 3 , we find J B 1 ∼ −4.0 meV, J B 2 ∼ 0.25 meV, and J B 3 ∼ 17 meV, in excellent agreement with recently reported values from inelastic neutron scattering [2,34].The decrease in T pp and ∆/T pd leads to an overall enhancement of J 3L 1 ∼ −4.5 meV, a vanishing J 3L 2 , and decrease in J 3L 3 ∼ 10 meV.In Fig. 4(b), we summarize the dependence of J 3 /J 1 on T pp and ∆/T pd .We find that J 1 is dominated by the 3d 9 L 1 state, while J 2 and J 3 are dominated by the 3d 8 state.Thus, the decrease in T pp is directly responsible for an increased 3d 9 L 1 contribution to J 1 .As a consequence of the overall reduction in the average exchange interaction strength, the magnetic transition temperature is expected to be reduced in fewlayer samples compared to bulk samples.Furthermore, the decrease in J 3 /J 1 from -4.2 in bulk to -2.2 in 3L, positions 3L NiPS 3 closer to a phase boundary between the stripy AFM phase and a spiral ordered phase [35].It is likely that, in the 2D limit, NiPS 3 is driven into a highly frustrated regime on this phase boundary.
In summary, we used RIXS to access the electronic ground state properties of an exfoliated, correlated antiferromagnet in the 2D limit.We found that 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 3 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 heterostuctures by tuning interfacial energy scales towards the design of the next generation of 2D strongly correlated magnetic materials.
Single crystal samples of NiPS 3 were grown by vapor transport, following previously published methods [1,2].Furnace temperature settings and duration are shown in Table SI.Bulk NiPS 3 was exfoliated using conventional scotchtape methods [3] and deposited either onto a blank SiO In both cases, Cu was chosen as a material that could provide a fluorescence contrast to SiO 2 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.

Sample Preservation
Sample degradation in air through oxidation is a known issue for exfoliated flakes [3,4].To ensure our samples remained intact, SiO 2 chips were spin-coated with a layer of PMMA immediately after exfoliation and then stored in an Ar glovebox.We tested the viability of this method by preparing two 5L flakes of NiPS 3 , treating one with PMMA while the other was left out in air for one week.In coating showed significant changes and a nearly vanishing XAS intensity.
In addition to a change in XAS spectra, we have also observed a degredation of RIXS spectra as a function of beam exposure time.After 11.5 hours of measurements, we began to observe a monotonic decrease in relative intensity between elastic line and NiS 6 multiplet excitations [Fig.SIII].From T = 0 m to T = 210 m, an overall broadening of the spectra can be seen coupled with a near vanishing of higher energy loss (E loss 2.5 eV) charge transfer excitations.

FITTING PROCEDURE Bulk Fitting and Peak Assignment
We began our analysis by exploring the parameter space of the d 8 Tanabe-Sugano (TS) diagram within a single-ion model with octahedral (O h ) symmetry, and subsequently, trigonal (D 3d ) symmetry [5].We thus minimized the difference between the calculated energies and those of peaks A -F for the bulk spectra.dd .10Dq was initially set to the E loss value of the nominal 10Dq line from the RIXS spectrum [6], with a best fit value of 1.05 eV.Additionally, we found the best agreement with the fitted peaks with a 49% reduction of F 0,2,4 dd from its atomic Hartree-Fock value, and a range of D 3d splitting between δ ∼ -60 and -100 meV.This analysis showed that the peaks between E loss ∼ 1 -1.75 eV can be explained as excitations within the valence band of Ni 2+ in the presence of a trigonally-distorted octrahedral complex.
We now turn our attention to the determination of peak A, which the single-ion model failed to capture.We repeated the above analysis as a function of D 3d CF splitting within the multiplet ligand field (MLFM) model in order to understand the change from Ni 2+ valence band excitations to full NiS 6 multiplet excitations, and thus, elucidate the origin of peak A. plore the full range of free MFLM parameters, including ∆, 10Dq, and the reduction factor to the Slater-Condon parameters F 0,2,4 dd , denoted F dd R. We find best agreement to our data with a small charge transfer energy, ∆ = 0.60 eV, which gives rise to the low energy excitation present in Fig. SV(b), and absent in Fig. SV(a).We thus conclude that the origin of peak A is a low energy charge transfer excitation and label it as 3d 9 L 1 .

3L Fitting
As mentioned in the main text, empirically fitting the 3L RIXS spectra resulted in two scenarios of equally good fit quality with a minimum of six Gaussian peaks.Using physically meaningful parameters, these two scenarios can be summarized by the sign of the charge transfer energy (∆) extracted from modeling.Scenario one, where peak widths were fixed to the experimental resolution, resulted in a negative charge transfer energy, while scenario two, where peak widths were allowed to relax, resulted in a positive charge transfer energy.

Positive Charge Transfer Scenario
Using physically meaningful parameters for octahedrally coordinated NiS 6 , we carried out a search of the parameter space for charge transfer energy ∆, ligandligand and metal-ligand hoppings T pp and pdσ, and trigonal distortion δ, keeping F (G) pd fixed to 80% of their Hartree-Fock values.The reduction factor for F 0,2,4 dd was additionally kept fixed at the fitted bulk value.Fig. SVII shows the calculated energy of all excited states within a NiS 6 cluster as a function of these MLFM parameters shown with best fit values for ∆ and trigonal distortion δ, along with fixed hopping parameters for T pp and pdσ obtained from DFT calculations.

Negative Charge Transfer Scenario
Following the same minimization procedure as the positive charge transfer scenario, Fig. SVIII shows the calculated energy levels and best fit values for ∆ and trigonal distortion δ, along with fixed hopping parameters T pp , pdσ.In addition to a negative charge transfer energy, ∆ = −0.55eV, a reduced cubic crystal field splitting relative to bulk, 10Dq = 0.31 eV, was also necessary to describe this fitting scenario.
In order to facilitate a comprehensive search of the MFLM parameter space within fitting scenario one, we tracked the change in peak A as a function of ∆ vs. (pdσ, T pp , δ), shown as colormaps in Fig. SIX.Within the range of fixed hopping parameters for bulk and 1L NiPS 3 , we find a clear delineation between the two fitting scenarios, where a peak A position of ∼0.6 eV is only achievable with negative charge transfer values.We start with an outline of our approach for obtaining key impurity model parameters from first principles calculations since, to our knowledge, no one has explicitly documented a procedure for obtaining metal-ligand and ligand-ligand LCAO hopping parameters (e.g., pdσ, pdπ, T pp ) from DFT in slightly distorted materials like NiPS 3 .Originally, we aimed to employ a localization technique similar to Ref. 7 wherein we obtain atomic orbital-like Wannier functions which could be used in an LCAO/TB 0.8 0.0 0.8 (eV)  framework to calculate NiPS 3 's multiplet ligand field theory parameters.However, the slight trigonal distortion of NiPS 3 and its more complex coordination environment compared to the ideal NiO octahedron of Ref. 7 motivated our choice to instead use symmetry-adapted, maximally localized Wannier functions (MLWFs) [8,9] since MLWFs are constructed directly from the converged DFT potential and do not require reconstruction of the potential between spherically symmetric local ba-sis functions and interstitial plane wave basis functions as in the N MTO approach [10].Chronologically, our DFT/Wannier90 workflow consisted of: 1.A self-consistent field calculation (described in further detail below) to solve for the set of singleparticle orbitals (SPO's) and energies which converged the ground state charge density of nonmagnetic NiPS 3 at the PBE+U (=4 eV) level of theory, followed by;  3. Calculation of the band structure in the first Brillouin zone of the P 1 bulk and monolayer structures to find the optimal band subspace {ψ nk } from which to obtain the MLWF's w n (r − R), and corresponding unitary transformations U k mn .We chose the 27 isolated bands highlighted in of the S and P p-shaped MLWF's ranged from 2-4 Å, indicating physically reasonable localization.Further, the ratios of imaginary to real components of the MLWF's were 1 × 10 −6 or less, which indicated that the converged MLWF's were good-quality.
Next, we needed to select energy expressions from the overconstrained systems of LCAO equations presented in Ref. 12 to solve for the metal-ligand hopping parameters.The presence of inequivalent sulfur sites and trigonal distortion in NiPS 3 yields many possible pdσ and pdπ values upon solving for hopping parameters based upon different Ni-S pairs.To remove this ambiguity, we computed our pdσ and pdπ values based upon the two largest energies w i (r)| Ĥ|w j (r) between Wannier function pairs based on the intuition that those with the largest energies contribute most significantly to the electronic structure of the material.Table SIII provides the pairs of energies used to solve for the hopping parameters for each Ni-S pair in ML and bulk NiPS 3 where d b is the corresponding bond length ( Å).As an example, solving for pdσ and pdπ for the bulk Ni1-S3 pair entailed solving the system of equations where l 1 , m 1 , and n 1 are the direction cosines between the atomic Ni and S sites.We explicitly obtained symbolic expressions for the orbital combinations not considered in Ref. [12] using the compact closed-form treatment provided in Ref. [13] and thereby verified their cyclic permutational symmetries [12].The ligand-ligand hopping parameters ppσ and ppπ (and thus T pp = ppσ − ppπ) were similarly determined via with the direction cosines l 2 , m 2 , and n 2 between the two S sites.Since only S3 and S4 were nearest neighbors, we report only one set of ppσ, ppπ, and T pp parameters for the bulk and monolayer systems in the main text.Additionally, it is noted that the same procedure using the bulk geometry (Table SV) with an added vertical vacuum of > 25 Å to prevent interaction of images in the out-of-plane direction, 'bulk+vac', yielded monolayer hopping parameters as reported in Table SIV.The resultant average of pdσ and pdπ values are respectively 0.29 eV and 0.06 eV larger in magnitude than those of the PBE-relaxed monolayer geometry 'ML' (Table SV), which reflects the sensitivity of pdσ and corresponding lack of sensitivity of pdπ to a Ni-S bond distance which is 0.06 Å shorter.Fits to the experimental data constrained by the unrelaxed monolayer 'bulk+vac' hopping parameters were not robust whereas those contrained by the PBE-relaxed monolayer 'ML' were, leading us to conclude that the change in RIXS response of NiPS 3 in the 2D limit is accompanied by only a slight increase in the 2D lattice constant relative to the bulk (0.5%), and changes in pdσ by 0.14 eV and pdπ by 0.21 eV.All simulations of structural and electronic properties were performed using DFT as implemented within the Quantum ESPRESSO package.[14] The PBE+U functional with U = 4 eV was selected to model this material since previous reports have shown this functional to yield structural parameters and band gaps that are consistent with experimental results.[15,16] The key features of the real part of the experimental optical conductivity, σ(ω), of bulk (zigzag AFM) NiPS 3 are also reproduced when σ(ω) is calculated from DFT+U using U = 4 eV.[17] Further, a noncollinear Neél antiferromagnetic (AFM) spin configuration was used for our relaxations, DFT parameter convergence, and preliminary band structure convergence calculations despite the experimentally observed zig-zag configuration of bulk NiPS 3 , since experimentally comparable lattice parameters were recovered with Neél ordering.[1] Also, previous work, verified by our own calculations, demonstrates that the noncollinear Neél antiferromagnetic and zigzag magnetic states are virtually degenerate for two-dimensional NiPS 3 [18] and thus yield magnetic moments and band gaps with negligible differences.The final analysis for obtaining MLWF's and hopping parameters was however done for nonmagnetic NiPS 3 using the Neél-obtained plane-wave cutoff and kgrid parameters.

Electronic Structure Analysis
To ensure that our DFT calculations were fully converged, first, we converged the Monkhorst-Pack k-point mesh and plane wave energy cutoff parameters for the PBE-relaxed bulk and monolayer structures (details may be found in SI ) with the Neél-antiferromagnetic configuration, using the PBE+U family of functionals.[19,20] A k-point mesh of 6 × 6 × 6 for the bulk and 4 × 4 × 2 for the monolayer was sufficient to converge their respective total energies to within 5×10 −5 Ry.The plane-wave cutoff energy was converged to within chemical accuracy on the converged k-grid at a high cutoff of 600 Ry (bulk shown in Due to the steep computational expense expected to accompany non-self consistent field and band structure calculations of NiPS 3 with plane wave cutoffs of 600 Ry, we also performed tests to explicitly see how the band structure of bulk Neél AFM NiPS 3 with the Materials Project (C2/m) structure is affected when it is calculated using partially converged plane wave cutoff energies.As illustrated in Figure SXIII, a band structure obtained with a plane wave cutoff of 250 Ry is essentially identical to the E cut = 600 Ry-obtained band structure.We thus deemed the subspace spanned by the wavefunctions obtained at E cut = 250 Ry sufficient for obtaining accurate Wannier90 tight-binding parameters, and all of our Wannier90 analysis is done on wavefunctions obtained at this level of theory.

Obtaining Bulk and Monolayer Crystal Structures
To obtain relaxed bulk and monolayer NiPS 3 crystal structures (see Figure SXII and Table SV), the bulk structure with experimentally-determined lattice param- eters having a monoclinic (C2/m) symmetry and a 10atom unit cell consisting of two NiPS 3 formula units was first obtained from the Materials Project (MP) website [22].The lattice vectors of the bulk structure were constrained to the experimental neutron diffraction values of crystalline NiPS 3 reported by Ref. 1 and atomic coordinates were relaxed using the Broyden-Fletcher-Goldfarb-Shanno quasi-Newton algorithm until ionic forces were smaller than 10 −4 a.u.This yielded the bulk structure having P 1 symmetry on which our DFT and Wannierization was performed (Figure SXIIa).One monolayer of NiPS 3 was cleaved from the starting MP structure and a vacuum of > 30 Å was added to prevent interlayer interactions between images for the monolayer simulations.Then, to inform our choice of exchange-correlation functional for monolayer variablecell relaxation since no experimental monolayer lattice parameters exist, we performed variable-cell relaxations of bulk NiPS 3 in the collinear Neél AFM magnetic configuration starting from the MP structure using the LDA and PBE functionals with and without Hubbard U = 2.0 eV, which yielded the lattice parameters in Table SVI and which are displayed alongside the experimental lattice parameters of paramagnetic bulk NiPS 3 obtained by Ref. 1.All bulk structural relaxations were performed such that all lattice vectors and angles could vary under only the constraint of C2/m cell symmetry; the PBE+Uconverged k-point meshes were employed for all cell relaxations.Evidently, the cell relaxation using the PBE functional yielded the closest lattice parameters to experiment; this thus provided our rationale for using the PBE-optimized monolayer NiPS 3 crystal structure (having P 1 symmetry) for our monolayer electronic structure and MLWF analysis (see Figure SXIIb).The exact ground state wavefunction within an NiS 6 cluster is given by |Ψ g = α 3d 8 + β 3d 9 L 1 + γ 3d 10 L 2 .Due to the magnitude of the trigonal field obtained from modeling (δ = −80 meV), we simplify the calculation by approximating Ni-S bonds within an NiS 6 cluster as orthogonal.As mentioned in the main text, the ground state character for bulk and 3L NiPS 3 finds a nearly equal mixture of 3d 8 and 3d 9 L 1 (|α| 2 ≈ |β| 2 ), with a negligibly small 3d 10 L 2 character.Thus, we independently apply the cell-perturbation to these two states.For fully occupied (inactive) t 2g orbitals, 3d 8 and 3d 9 L 1 are given by 3d 8 = d x 2 −y 2 d 3z 2 −r 2 x ), and p i x,y,z denote holes in the S 3p orbitals.Utilizing the hopping parameters obtained from ab initio calculations, we evaluate the twocenter atomic overlap integrals within the Slater-Koster scheme of linear combinations of atomic orbitals.The superexchange for 3d 8 and 3d 9 L 1 states are defined as J α i and J β i , where i ∈ {1, 2, 3}.The 1 NN superexchange is given by the fourth order perturbation term in 3d 8 , and second order term in

FIG. 3 .
FIG. 3. (a) & (b) Normalized Ni L3 RIXS spectra for bulk (a) and 3L (c) NiPS3 with AIM model from parameters in Table I and Ref.[19].(c) & (d) Simulated RIXS map as a function of incident energy from resonance (∆E Res 2 substrate or onto a SiO 2 substrate pre-treated with a patterned Copper (Cu) grid.The patterned Cu grid consisted of 100 µm x µm SiO 2 cells separated by 200 µm of 50 nm thickness Cu [Fig.SI(b) & (c)].The 3L sample was deposited onto a blank SiO 2 substrate and was later patterned with a Cu fiducial marker using electronbeam lithography [Fig.SI(a)], again with a Cu thickness of 50 nm.
FIG. SI.(a) Cu fiducial marker patterned around preexfoliated 3L flake of NiPS3.Optical contrast of x100 magnification optical image was crucial in determining sample thickness.(b) Sample image of exfoliated flakes on SiO2 substrate pre-treated with patterned Cu grid.(c) x100 magnification optical image of Cu grid cell depicting measured 7L sample and dimensions.
FIG. SII.XAS spectra taken before and after 4 hours of beam exposure for two 5L samples of NiPS3.(a) Sample preserved with PMMA immediately after exfoliation.(b) Sample left in air for one week prior to measurements.change from their free-ion values as a function of O h and D 3d crystal field (CF) splittings compared to E loss values of fitted peaks.For this calculation, our model included the O h and D 3d CF splitting terms, 10Dq and δ, and Slater-Condon parameters, F 0,2,4dd .10Dq was initially set to the E loss value of the nominal 10Dq line from the RIXS spectrum[6], with a best fit value of 1.05 eV.Additionally, we found the best agreement with the fitted peaks with a 49% reduction of F 0,2,4 FIG. SIII.RIXS spectra at 160K at Ni L3 resonance in 30 minute segments for an ∼ 5L samples of NiPS3.RIXS spectra segments were collected after a previous 11.5 hours of beam exposure.
FIG. SIV.Calculated excited state energy levels for Ni 2+ as a function of O h CF splitting 10Dq (left) and D 3d CF splitting δ (right) in the absence of spin-orbit coupling.F 0,2,4 dd was fit to 49% of its atomic Hartree-Fock value.Horizontal dashed lines and shaded regions show E loss value of fitted peaks.Only symmetry labels relevant to the discussion in the main text are included.
FIG. SV. energy levels of excited states as a function of D 3d CF splitting within a (a) single-ion model and (b) multiplet ligand-field model.Black and gray lines show energy levels calculated with and without SOC respectively.Horizontal dashed lines and shaded regions show E loss value of fitted peaks

pd
FIG. SVII.Calculated energy of all excited states as a function of MLFM parameters within fitting scenario two in comparison two 3L RIXS spectra.Horizontal dashed lines with respective shaded regions show the centroid position of fitted peaks.Vertical dashed lines represent the best fit values for each parameter.
FIG. SVIII.Calculated energy of all excited states as a function of MLFM parameters within fitting scenario one in comparison two 3L RIXS spectra.Horizontal dashed lines with respective shaded regions show the centroid position of fitted peaks.Vertical dashed lines represent the best fit values for each parameter.

2 .
A non-self consistent mapping of the SPO's onto a dense 6 × 6 × 6 Monkhorst pack k-grid; FIG. SX.Bulk NiPS3 band structure with the band subspace selected in this work outlined in red, and the corresponding partial density of states indicating the Ni d-orbital (purple), P p-orbital (green), and S p-orbital characters (blue) depicted on the right.
FIG. SXI.PBE+U (= 2 eV)-calculated total energies (8 × 8 × 8) of bulk coordinate-relaxed NiPS3 constrained to the experimental cell parameters of Ref. 1 as a function of the plane wave energy cutoff E P W cut .
FIG. SXII.Top views of the crystal structures of (a) monolayer NiPS3 optimized with PBE and (b) bulk NiPS3 constrained to experimental lattice vectors[1] with atomic coordinates using PBE+(U =4 eV).

TABLE I .
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 30] 49 (F 2 dd + F 4 dd ) extracted from RIXS modeling[7,30].FIG. 4. (a) Calculated ground state ratio of 3d 8 / 3d 9 L 1 as a function of Tpp and ∆/T pd .Star and hexagon points indicate bulk and 3L values of Tpp and ∆/T pd respectively. (b)Calculated change in J3/J1 as function of Tpp (top-axis) and ∆/T pd (bottom-axis) for fixed bulk values of ∆/T pd and Tpp, respectively.

TABLE SI .
Furnace temperature settings and duration FIG. SVI.Calculated energy of all excited states as a function MLFM parameters in comparison to bulk RIXS spectra.Horizontal dashed lines with respective shaded regions show the centroid position of fitted peaks.Vertical dashed lines represent the best fit values for each parameter.

TABLE SIV .
Predicted "bulk" with vacuum NiPS3 metalligand hopping parameters in eV and their corresponding largest energies.

TABLE SV .
Relaxed NiPS3 crystal structure parameters.

TABLE SVI .
Bulk NiPS3 crystal structure parameters after full-cell relaxation with different DFT functionals, Vxc, compared with those from experiment.Where a U is specified in the table, we employed U = 4 eV.