# Dimensionality dependent electronic structure of the exfoliated van der
Waals antiferromagnet NiPS3

M. F. DiScala  Department of Physics, Brown University, Providence, RI 02912
D. Staros  Department of Chemistry, Brown University, Providence, RI 02912  A.
de la Torre  Department of Physics, Brown University, Providence, RI 02912  A.
Lopez  Department of Physics, Brown University, Providence, RI 02912  D. Wong
Department of Dynamics and Transport in Quantum Materials, Helmholtz-Zentrum
Berlin für Materialen und Energie, Albert-Einstein-Strasse 15, 12489 Berlin,
Germany  C. Schulz  Department of Dynamics and Transport in Quantum Materials,
Helmholtz-Zentrum Berlin für Materialen und Energie, Albert-Einstein-Strasse
15, 12489 Berlin, Germany  M. Bartkowiak  Department of Dynamics and Transport
in Quantum Materials, Helmholtz-Zentrum Berlin für Materialen und Energie,
Albert-Einstein-Strasse 15, 12489 Berlin, Germany  B. Rubenstein  Department
of Chemistry, Brown University, Providence, RI 02912  K. W. Plumb  Department
of Physics, Brown University, Providence, RI 02912  Authors to whom
correspondence should be addressed: Brenda Rubenstein,
<EMAIL_ADDRESS>and Kemp Plumb<EMAIL_ADDRESS>
###### Abstract

Resonant Inelastic X-ray Scattering (RIXS) was used to measure the local
electronic structure in few-layer exfoliated flakes of the van der Waals
antiferromagnet NiPS3. The resulting spectra show a systematic softening and
broadening of NiS6 multiplet excitations with decreasing layer count from the
bulk to three atomic layers (3L). These trends are driven by a decrease in the
transition metal-ligand and ligand-ligand hopping integrals, and in the
charge-transfer energy: $\Delta$ = 0.60 eV in the bulk and 0.22 eV in 3L
NiPS3. Relevant intralayer magnetic exchange integrals computed from the
electronic parameters exhibit a systematic decrease in the average interaction
strength with thickness and place 2D NiPS3 close to the phase boundary between
stripy and spiral antiferromagnetic order, which may explain the apparent
vanishing of long-range order in the 2D limit. This study explicitly
demonstrates the influence of _inter_ layer electronic interactions on _intra_
layer ones in insulating magnets. As a consequence, the magnetic Hamiltonian
in few-layer insulating magnets can be significantly different from that in
the bulk.

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, 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].

Figure 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.

NiPS3 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 NiPS3 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 NiPS3.

In this letter, we show that Ni-S electronic energy scales are strongly
altered by dimensionality in NiPS3 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
NiS6 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
NiPS3 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 NiPS3, it’s effect will be important to the properties of a broad
class of few-layer insulating vdW magnets.

Figure 2: Calculated energy levels of the NiS6 multiplet ligand-field model as
a function of $\Delta$ (a) compared to experimental bulk spectrum (b). Fixed
model parameters are listed in Table 1 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_{\mathrm{loss}}$ value of fitted peaks in
rightmost panel. Vertical dashed line indicates best fit at $\Delta$ = 0.60
eV.

Single crystals of NiPS3 were grown via standard vapor transport methods [1,
2]. Bulk NiPS3 was exfoliated in air using conventional scotch-tape methods
[3] and deposited either onto a blank SiO2 substrate or onto a SiO2 substrate
pre-treated with a patterned Cu grid. Samples deposited onto blank SiO2
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\theta\\!=\\!90^{\circ}$ was used with an
$\approx$ 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
NiPS3 spectra were collected with an identical scattering geometry, but with
an $\approx$ 177 meV overall energy resolution FWHM.

Figs. 1(a) & (b) show representative RIXS spectra for bulk and three-layer
(3L) NiPS3, respectively. Spectra were collected at the peak of the Ni L3-edge
XAS $E_{i}\\!=\\!853$ eV, corresponding to 2$p_{3/2}$ to 3$d$ electronic
transitions. We concentrate on the low energy region $E_{\mathrm{loss}}$ = 0.2
$\rightarrow$ 2.15 eV that contains excitations within the NiS6 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 NiPS3 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 NiS6 octahedra. The center of mass positions of the two bulk peaks
are assigned to the $t_{2g}\rightarrow 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}T_{2g}$ $\rightarrow$ ${}^{3}A_{1g}$ $+$
${}^{3}E_{g}$ (peaks B & C) and ${}^{3}T_{1g}$ $\rightarrow$ ${}^{3}A_{2g}$
$+$ ${}^{3}E_{g}$ (peaks E & F), in agreement with Raman and optical
measurements [21, 22]. Access to spin-flip ($\Delta S\neq 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}\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 NiPS3.
We verify these peak assignments through the application of the NiS6 multiplet
ligand-field model described below. We note that since our incident energy was
tuned to the peak of the Ni L3-edge XAS, our measurements were not sensitive
to the sharp 1.47 eV peak reported in Ref. [23].

Figure 3: (a) & (b) Normalized Ni L3 RIXS spectra for bulk (a) and 3L (c)
NiPS3 with AIM model from parameters in Table 1 and Ref.[19]. (c) & (d)
Simulated RIXS map as a function of incident energy from resonance ($\Delta
E_{i}^{Res}$) and $E_{\mathrm{loss}}$. $\Delta/T_{pd}$ and $\Delta/U$
parameterize the hybridization and charge transfer characters respectively.

In order to facilitate an efficient exploration of the parameter 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\sigma$ and $pd\pi$, S
3p\- orbital level splitting, $T_{pp}$ = $pp\sigma$ $-$ $pp\pi$, cubic crystal
field (10Dq) and trigonal distortion $\delta$, Ni $3d$-$3d$ and $2p$-$3d$ on-
site Coulomb interactions $U_{dd}$ and $U_{pd}\\!=\\!1.2U_{dd}$, and charge
transfer energy $\Delta$. For initial comparisons of this model to our data,
we used physically meaningful parameters for octahedrally coordinated NiS6
[25, 26, 27]. We carried out a search of the parameter space for $\Delta$,
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 [19, 6]. Fig. 2 shows the calculated energy levels for a
NiS6 cluster as a function of the charge transfer energy $\Delta$ for fixed
parameters that give the best agreement between measured and calculated peak
energies [19]. Previous optical and X-ray absorption (XAS) studies classified
NiPS3 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 $\Delta\\!=\\!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 NiPS3 at $\approx$ 38 nm. From bulk to 3L, the
${}^{3}T_{2g}$ modes (${}^{3}A_{1g}$ \+ ${}^{3}E_{g}$) shift slightly upward
in $E_{\mathrm{loss}}$, while ${}^{1}E_{g}$ has the largest upward shift of
$\Delta E_{loss}$ = +149(24) meV; the ${}^{3}T_{1g}$ modes (${}^{3}A_{2g}$ \+
${}^{3}E_{g}$) also shift in $E_{\mathrm{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 NiPS3.

Table 1: Fixed values in eV of hopping integrals extracted from ab initio calculations of the nonmagnetic configuration. Charge transfer energy $\Delta$, and intra-orbital Coulomb repulsion $U=F^{0}_{dd}+\frac{4}{49}(F^{2}_{dd}+F^{4}_{dd})$ extracted from RIXS modeling [7, 30]. |  $pd\sigma$ |  $pd\pi$ |  $pp\sigma$ |  $pp\pi$ |  $T_{pp}$ |  $\Delta$ |  U
---|---|---|---|---|---|---|---
Bulk |  -1.07 |  0.67 |  0.89 |  -0.09 |  0.98 |  0.60 |  8.3
3L |  -0.93 |  0.46 |  0.62 |  -0.01 |  0.63 |  0.22 |  8.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 NiPS3 in both the bulk and monolayer
geometries, generated maximally-localized 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\sigma$ and $pd\pi$, as well as $T_{pp}$ were fixed to those obtained from
ab initio calculations of the nonmagnetic configuration [Table 1] while
$\Delta$ 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_{\mathrm{loss}}$. We can reproduce
the observed broadening and softening of NiS6 multiplet excitations with
thickness by a decrease in charge transfer energy and transition metal-ligand
hopping integrals, as parameterized by $\Delta/T_{pd}$, and $\Delta/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})$ [7, 30]: $\Delta/T_{pd}$ =
0.32 and $\Delta/U$ = 0.07 in bulk, and $\Delta/T_{pd}$ = 0.14 and $\Delta/U$
= 0.03 in 3L.

Figure 4: (a) Calculated ground state ratio of
$\ket{3d^{8}}$/$\ket{3d^{9}\underline{L}^{1}}$ as a function of $T_{pp}$ and
$\Delta/T_{pd}$. Star and hexagon points indicate bulk and 3L values of
$T_{pp}$ and $\Delta/T_{pd}$ respectively. (b) Calculated change in
$J_{3}/J_{1}$ as function of $T_{pp}$ (top-axis) and $\Delta/T_{pd}$ (bottom-
axis) for fixed bulk values of $\Delta/T_{pd}$ and $T_{pp}$, respectively.

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 NiPS3 gets
thinner, metal-ligand $\pi$-hopping is reduced ($pd\pi$ decreases) due to the
removal of stabilizing, $\pi$-like, interlayer vdW interactions. The same
effect also causes $pd\sigma$ and $T_{pp}$ to change significantly because of
the mixed $\sigma$\- and $\pi$-bonding character present in the
$sp_{3}$-hybridized phosphorus atoms that bridge the NiS6 clusters. In the
context of our MLWFs, this is reflected in a change in the largest tight-
binding energies used to solve for $pd\sigma$ and $pd\pi$ [19].

The combination of RIXS measurements and ab initio calculations constrain the
electronic ground state that underlies the magnetic properties of NiPS3. In
Fig. 4(a), 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
$\Delta/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,
$\absolutevalue{\alpha}^{2}/\,\absolutevalue{\beta}^{2}=1.08$, and in 3L,
$\absolutevalue{\alpha}^{2}/\,\absolutevalue{\beta}^{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 $\Delta/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}$.

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
$\ket{3d^{9}\underline{L}^{1}}$ states, and fourth and sixth order terms for
the $\ket{3d^{8}}$ states [27]; a detailed description of these expressions is
given in Ref. [19]. In bulk NiPS3, we find $J_{1}^{B}\\!\sim\\!-4.0$ meV,
$J_{2}^{B}\\!\sim\\!0.25$ meV, and $J_{3}^{B}\\!\sim\\!17$ meV, in excellent
agreement with recently reported values from inelastic neutron scattering [34,
2]. The decrease in $T_{pp}$ and $\Delta/T_{pd}$ leads to an overall
enhancement of $J_{1}^{3L}\\!\sim\\!-4.5$ meV, a vanishing $J^{3L}_{2}$, and
decrease in $J_{3}^{3L}\\!\sim\\!10$ meV. In Fig. 4(b), we summarize the
dependence of $J_{3}/J_{1}$ on $T_{pp}$ and $\Delta/T_{pd}$. We find that
$J_{1}$ is dominated by the $\ket{3d^{9}\underline{L}^{1}}$ state, while
$J_{2}$ and $J_{3}$ 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_{1}$. 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_{3}/J_{1}$ from -4.2
in bulk to -2.2 in 3L, positions 3L NiPS3 closer to a phase boundary between
the stripy AFM phase and a spiral ordered phase [35]. It is likely that, in
the 2D limit, NiPS3 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 NiPS3 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.

##  I Acknowledgments

###### Acknowledgements.

We thank Mark Dean for helpful discussions and comments on this manuscript. We
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. is 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 computational
resources and services at the Center for Computation and Visualization, Brown
University.

## References

  * Park [2016] J.-G. Park, Opportunities and challenges of 2D magnetic van der Waals materials: Magnetic graphene?, J. Phys.: Condens. Matter 28, 301001 (2016).
  * Kuo  _et al._ [2016] 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, and J.-G. Park, Exfoliation and Raman Spectroscopic Fingerprint of Few-Layer NiPS3 Van der Waals Crystals, Sci Rep 6, 20904 (2016).
  * Kim  _et al._ [2019] K. Kim, S. Y. Lim, J.-U. Lee, S. Lee, T. Y. Kim, K. Park, G. S. Jeon, C.-H. Park, J.-G. Park, and H. Cheong, Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS3, Nat Commun 10, 345 (2019).
  * Sivadas  _et al._ [2018] N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao, Stacking-Dependent Magnetism in Bilayer CrI 3, Nano Lett. 18, 7658 (2018).
  * Gong  _et al._ [2017] 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, and X. Zhang, Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals, Nature 546, 265 (2017).
  * Tang  _et al._ [2023] Y. Tang, K. Su, L. Li, Y. Xu, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, C.-M. Jian, C. Xu, K. F. Mak, and J. Shan, Evidence of frustrated magnetic interactions in a Wigner–Mott insulator, Nat. Nanotechnol. , 1 (2023).
  * Hellman  _et al._ [2017] 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, and B. L. Zink, Interface-induced phenomena in magnetism, Rev. Mod. Phys. 89, 025006 (2017).
  * Huang  _et al._ [2017] 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, and X. Xu, Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit, Nature 546, 270 (2017).
  * Burch  _et al._ [2018] K. S. Burch, D. Mandrus, and J.-G. Park, Magnetism in two-dimensional van der Waals materials, Nature 563, 47 (2018).
  * Gibertini  _et al._ [2019] M. Gibertini, M. Koperski, A. F. Morpurgo, and K. S. Novoselov, Magnetic 2D materials and heterostructures, Nat. Nanotechnol. 14, 408 (2019).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics January Special Issue, 321, 2 (2006).
  * Plumb  _et al._ [2014] K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, $\alpha$ \- RuCl 3 : A spin-orbit assisted Mott insulator on a honeycomb lattice, Phys. Rev. B 90, 041112(R) (2014).
  * Nasu  _et al._ [2016] J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, and R. Moessner, Fermionic response from fractionalization in an insulating two-dimensional magnet, Nature Phys 12, 912 (2016).
  * Lee  _et al._ [2018] K. H. Lee, S. B. Chung, K. Park, and J.-G. Park, Magnonic quantum spin Hall state in the zigzag and stripe phases of the antiferromagnetic honeycomb lattice, Phys. Rev. B 97, 180401(R) (2018).
  * Zhou  _et al._ [2019] B. Zhou, Y. Wang, G. B. Osterhoudt, P. Lampen-Kelley, D. Mandrus, R. He, K. S. Burch, and E. A. Henriksen, Possible structural transformation and enhanced magnetic fluctuations in exfoliated $\alpha$-RuCl3, Journal of Physics and Chemistry of Solids Spin-Orbit Coupled Materials, 128, 291 (2019).
  * Wildes  _et al._ [2015] 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, and T. J. Hicks, Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS 3, Phys. Rev. B 92, 224408 (2015).
  * Kim  _et al._ [2018] 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, and T. W. Noh, Charge-Spin Correlation in van der Waals Antiferromagnet NiPS 3, Phys. Rev. Lett. 120, 136402 (2018).
  * Lançon  _et al._ [2018] D. Lançon, R. A. Ewings, T. Guidi, F. Formisano, and A. R. Wildes, Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS 3, Phys. Rev. B 98, 134414 (2018).
  * [19] See supplementary material 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.
  * Schulz  _et al._ [2020] C. Schulz, K. Lieutenant, J. Xiao, T. Hofmann, D. Wong, and K. Habicht, Characterization of the soft X-ray spectrometer PEAXIS at BESSY II, J Synchrotron Rad 27, 238 (2020).
  * Afanasiev  _et al._ [2021] 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, and A. D. Caviglia, Controlling the anisotropy of a van der Waals antiferromagnet with light, Science Advances 7, eabf3096 (2021).
  * Wang  _et al._ [2022] 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, and X. Ling, Electronic Raman scattering in the 2D antiferromagnet NiPS3, Science Advances 8, eabl7707 (2022).
  * Kang  _et al._ [2020] 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, and J.-G. Park, Coherent many-body exciton in van der Waals antiferromagnet NiPS3, Nature 583, 785 (2020).
  * Haverkort  _et al._ [2012] M. W. Haverkort, M. Zwierzycki, and O. K. Andersen, Multiplet ligand-field theory using Wannier orbitals, Phys. Rev. B 85, 165113 (2012).
  * Bocquet  _et al._ [1992] A. E. Bocquet, T. Mizokawa, T. Saitoh, H. Namatame, and A. Fujimori, Electronic structure of 3 _d_ -transition-metal compounds by analysis of the 2 _p_ core-level photoemission spectra, Phys. Rev. B 46, 3771 (1992).
  * Krishnakumar and Sarma [2003] S. R. Krishnakumar and D. D. Sarma, X-ray photoemission study of NiS 2 - x Se x ( x = 0.0 – 1 . 2 ), Phys. Rev. B 68, 155110 (2003).
  * Takubo  _et al._ [2007] K. Takubo, T. Mizokawa, J.-Y. Son, T. Nambu, S. Nakatsuji, and Y. Maeno, Unusual Superexchange Pathways in an NiS 2 Triangular Lattice with Negative Charge-Transfer Energy, Phys. Rev. Lett. 99, 037203 (2007).
  * Ghiringhelli  _et al._ [2005] G. Ghiringhelli, M. Matsubara, C. Dallera, F. Fracassi, R. Gusmeroli, A. Piazzalunga, A. Tagliaferri, N. B. Brookes, A. Kotani, and L. Braicovich, NiO as a test case for high resolution resonant inelastic soft x-ray scattering, J. Phys.: Condens. Matter 17, 5397 (2005).
  * Yan  _et al._ [2021] M. Yan, Y. Jin, Z. Wu, A. Tsaturyan, A. Makarova, D. Smirnov, E. Voloshina, and Y. Dedkov, Correlations in the Electronic Structure of van der Waals NiPS 3 Crystals: An X-ray Absorption and Resonant Photoelectron Spectroscopy Study, J. Phys. Chem. Lett. 12, 2400 (2021).
  * Shen  _et al._ [2022] 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, and M. P. M. Dean, Role of Oxygen States in the Low Valence Nickelate La$_4$Ni$_3$O$_8$, Phys. Rev. X 12, 011055 (2022), arXiv:2110.08937 [cond-mat] .
  * Wang  _et al._ [2019] Y. Wang, G. Fabbris, M. P. M. Dean, and G. Kotliar, EDRIXS: An open source toolkit for simulating spectra of resonant inelastic x-ray scattering, Computer Physics Communications 243, 151 (2019).
  * Jefferson  _et al._ [1992] J. H. Jefferson, H. Eskes, and L. F. Feiner, Derivation of a single-band model for CuO 2 planes by a cell-perturbation method, Phys. Rev. B 45, 7959 (1992).
  * Eskes and Jefferson [1993] H. Eskes and J. H. Jefferson, Superexchange in the cuprates, Phys. Rev. B 48, 9788 (1993).
  * Wildes  _et al._ [2022] A. R. Wildes, J. R. Stewart, M. D. Le, R. A. Ewings, K. C. Rule, G. Deng, and K. Anand, Magnetic dynamics of $nips_{3}$, Phys. Rev. B 106, 174422 (2022).
  * Fouet  _et al._ [2001] J. Fouet, P. Sindzingre, and C. Lhuillier, An investigation of the quantum J 1 - J 2 - J 3 model on the honeycomb lattice, Eur. Phys. J. B 20, 241 (2001).

Supplementary Information for ”Dimensionality dependent electronic structure
of the exfoliated van der Waals antiferromagnet NiPS3”

##  II Sample preparation

###  II.1 Growth and Fabrication

Single crystal samples of NiPS3 were grown by vapor transport, following
previously published methods [1, 2]. Furnace temperature settings and duration
are shown in Table SI.

Table SI: Furnace temperature settings and duration Zone 1 (∘C) |  Zone 2 (∘C) |  Duration (Days)
---|---|---
700 |  750 |  2
670 |  620 |  16

Bulk NiPS3 was exfoliated using conventional scotch-tape methods [3] and
deposited either onto a blank SiO2 substrate or onto a SiO2 substrate pre-
treated with a patterned Copper (Cu) grid. The patterned Cu grid consisted of
100 $\mu m$ x $\mu m$ SiO2 cells separated by 200 $\mu m$ of 50 nm thickness
Cu [Fig. SI(b) & (c)]. The 3L sample was deposited onto a blank SiO2 substrate
and was later patterned with a Cu fiducial marker using electron-beam
lithography [Fig. SI(a)], again with a Cu thickness of 50 nm. In both cases,
Cu was chosen as a material that could provide a fluorescence contrast to SiO2
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.

Figure SI: (a) Cu fiducial marker patterned around pre-exfoliated 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.

###  II.2 Sample Preservation

Figure SII: XAS spectra taken before and after 4 hours of beam exposure for
two $\lesssim$ 5L samples of NiPS3. (a) Sample preserved with PMMA immediately
after exfoliation. (b) Sample left in air for one week prior to measurements.

Sample degradation in air through oxidation is a known issue for exfoliated
flakes [3, 4]. To ensure our samples remained intact, SiO2 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
$\lesssim$5L flakes of NiPS3, treating one with PMMA while the other was left
out in air for one week. In Fig. SII we show XAS measurements taken on these
two samples at the PEAXIS beamline pre and post a four hour beam exposure at
the Ni L3 edge in a specular scattering configuration. The sample coated in
PMMA showed minimal changes to XAS spectra, while the sample without the PMMA
protective 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 NiS6 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_{\mathrm{loss}}$ $\gtrsim$ 2.5 eV)
charge transfer excitations.

##  III Fitting Procedure

###  III.1 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. Fig. SIV shows how the excited state energy levels
for divalent nickel change from their free-ion values as a function of $O_{h}$
and $D_{3d}$ crystal field (CF) splittings compared to $E_{\mathrm{loss}}$
values of fitted peaks. For this calculation, our model included the $O_{h}$
and $D_{3d}$ CF splitting terms, 10Dq and $\delta$, and Slater-Condon
parameters, $F_{dd}^{0,2,4}$. 10Dq was initially set to the
$E_{\mathrm{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_{dd}^{0,2,4}$ from its atomic
Hartree-Fock value, and a range of $D_{3d}$ splitting between $\delta\sim$ -60
and -100 meV. This analysis showed that the peaks between
$E_{\mathrm{loss}}$$~{}\sim$ 1 - 1.75 eV can be explained as excitations
within the valence band of Ni2+ in the presence of a trigonally-distorted
octrahedral complex.

Figure SIII: RIXS spectra at 160K at Ni L3 resonance in 30 minute segments for
an $\sim$ 5L samples of NiPS3. RIXS spectra segments were collected after a
previous 11.5 hours of beam exposure. Figure SIV: Calculated excited state
energy levels for Ni2+ as a function of $O_{h}$ CF splitting 10Dq (left) and
$D_{3d}$ CF splitting $\delta$ (right) in the absence of spin-orbit coupling.
$F_{dd}^{0,2,4}$ was fit to 49% of its atomic Hartree-Fock value. Horizontal
dashed lines and shaded regions show $E_{\mathrm{loss}}$ value of fitted
peaks. Only symmetry labels relevant to the discussion in the main text are
included.

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 Ni2+ valence band excitations to full NiS6
multiplet excitations, and thus, elucidate the origin of peak A. Fig. SV shows
the TS diagrams for these two models. Fig. SV(a) is equivalent to the
rightmost panel in Fig. SIV, presented here including spin-orbit coupling.
Fig. SV(b) shows the MLFM TS diagram using parameters discussed in the main
text. Comparing peaks between $E_{\mathrm{loss}}$$~{}\sim$ 1 - 1.75 eV in
these two TS diagrams shows how the introduction of additional interaction
energies in the MFLM alter the Ni2+ valence band excitations. Including the
analysis in Fig. SV(b), the complete minimization between calculated energies
and peaks A-F is continued in Fig. SVI, where we explore the full range of
_free_ MFLM parameters, including $\Delta$, 10Dq, and the reduction factor to
the Slater-Condon parameters $F_{dd}^{0,2,4}$, denoted $F_{dd}$R. We find best
agreement to our data with a small charge transfer energy, $\Delta$ = 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}\underline{L}^{1}$.

###  III.2 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 ($\Delta$)
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_.

###  III.3 Positive Charge Transfer Scenario

Using physically meaningful parameters for octahedrally coordinated NiS6, we
carried out a search of the parameter space for charge transfer energy
$\Delta$, ligand-ligand and metal-ligand hoppings $T_{pp}$ and $pd\sigma$, and
trigonal distortion $\delta$, keeping $F(G)_{pd}$ fixed to 80% of their
Hartree-Fock values. The reduction factor for $F_{dd}^{0,2,4}$ was
additionally kept fixed at the fitted bulk value. Fig. SVII shows the
calculated energy of all excited states within a NiS6 cluster as a function of
these MLFM parameters shown with best fit values for $\Delta$ and trigonal
distortion $\delta$, along with fixed hopping parameters for $T_{pp}$ and
$pd\sigma$ obtained from DFT calculations.

###  III.4 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 $\Delta$ and trigonal distortion $\delta$, along with fixed hopping
parameters $T_{pp}$, $pd\sigma$. In addition to a negative charge transfer
energy, $\Delta=-0.55$ eV, a reduced cubic crystal field splitting relative to
bulk, 10Dq = 0.31 eV, was also necessary to describe this fitting scenario.

Figure SV: Calculated 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_{\mathrm{loss}}$ value of fitted peaks

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
$\Delta$ vs. ($pd\sigma$, $T_{pp}$, $\delta$), shown as colormaps in Fig. SIX.
Within the range of fixed hopping parameters for bulk and 1L NiPS3, we find a
clear delineation between the two fitting scenarios, where a peak A position
of $\sim$0.6 eV is only achievable with negative charge transfer values.

Figure SVI: Calculated energy of all excited states as a function of 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. Figure 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. Figure 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. Figure SIX: Calculated energy of
$3d^{9}\underline{L}^{1}$ excited state, peak A, as a function of MLFM
parameters. Contour lines indicate energies extracted from the two fitting
routines for the 3L data, vertical dashed lines are parameter values from DFT
for Bulk (violet) and 3L (green). These values are $pd\sigma$ = -0.93, -1.07,
and $T_{pp}$ = 0.98, 0.63 for Bulk and 3L respectively.

##  IV Electronic Structure Details

In this work, we used first principles simulations to accurately determine and
constrain the ligand field theory parameters employed in the Anderson impurity
model (AIM) used to model the experimental RIXS spectra. Our model which has a
total of 14 adjustable parameters: $pd\sigma$, $pd\pi$, $pp\sigma$, $pp\pi$,
charge transfer energy $\Delta$, spin-orbit coupling coefficients
$\zeta_{d,i}$ of the initial state and $\zeta_{d,n}$ of the intermediate
state, Coulomb energies $U_{dd}$ between $d$-orbitals and between $p/d$
orbitals $U_{dp}$, the reduction prefactors to the atomic Slater integrals
$F_{dd}$, $F_{dp}$, $G_{dp}$ obtained from Hartree-Fock calculations of a free
ion, crystal field splitting energy $10Dq$, and trigonal distortion $\delta$.

Table SII: Parameters and values used in MLFM calculations. Parameters not explicitly listed are the Ni 2p-3d on-site Coulomb interaction $U_{dp}$, fixed to 1.2$U_{dd}$, and the 3d SOC for the initial and intermediate RIXS states, fixed to $\zeta_{3d,i}=83$ and $\zeta_{3d,n}=102$ meV respectively. $F_{dd}R$ and $F(G)_{pd}R$ are the Slater integral scalings applied to their atomic Hartree-Fock values. $F_{pd}$ and $G_{pd}$ are fixeed to 80% of their atomic Hartree-Fock values.  Parameter |  Bulk Model |  3L Model
---|---|---
10Dq | 0.43 eV | 0.40 eV
$\delta$ | -0.08 eV | -0.08 eV
Udd | 5.5 eV | 5.5 eV
$\Delta$ | 0.60 eV | 0.22 eV
$pd\sigma$ | -1.07 eV | -0.93 eV
$pd\pi$ | 0.67 eV | 0.46 eV
Tpp | 0.98 eV | 0.63 eV
FddR | 0.85 | 0.85
|  FdpR =
---
GdpR
0.8 | 0.8

###  IV.1 Workflow for Determining Hopping Parameters from First Principles
DFT Calculations

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\sigma$, $pd\pi$, $T_{pp}$) from DFT
in slightly distorted materials like NiPS3. 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 framework to calculate
NiPS3’s multiplet ligand field theory parameters. However, the slight trigonal
distortion of NiPS3 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
basis functions and interstitial plane wave basis functions as in the $N$MTO
approach [10]. Chronologically, our DFT/Wannier90 workflow consisted of:

  1. 1.

A self-consistent field calculation (described in further detail below) to
solve for the set of single-particle orbitals (SPO’s) and energies which
converged the ground state charge density of nonmagnetic NiPS3 at the
PBE+$U$(=4 eV) level of theory, followed by;

  2. 2.

A non-self consistent mapping of the SPO’s onto a dense $6\times 6\times 6$
Monkhorst pack k-grid;

  3. 3.

Calculation of the band structure in the first Brillouin zone of the $P1$ bulk
and monolayer structures to find the optimal band subspace
$\left\\{\psi_{n\textbf{k}}\right\\}$ from which to obtain the MLWF’s
$w_{n}(\textbf{r}-\textbf{R})$,

$w_{n}(\textbf{r}-\textbf{R})=\int_{BZ}\sum_{m}U^{\textbf{k}}_{mn}\psi_{m\textbf{k}}(\textbf{r})e^{-i\textbf{k}\cdot\textbf{r}}d\textbf{k}$
(1)

and corresponding unitary transformations $U_{mn}^{\textbf{k}}$. We chose the
27 isolated bands highlighted in Figure SX as the subspace from which to
obtain our 25 MLWF’s since they contained almost all of the Ni $d$ and S $p$
partial density of states. Further details regarding the choice of NiPS3
crystal structures and the DFT methodology are respectively provided in SI
IV.2.2 and SI IV.2.1;

Figure 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.

  4. 4.

Calculation of the initial overlaps of the Bloch-periodic parts of the SPO’s;

  5. 5.

Conversion of the SPO’s to Wannier90-readable formats using pw2wannier90 [11];

  6. 6.

Determination of the MLWF’s according to Equation 1 [9]; and lastly

  7. 7.

Calculation of metal-ligand and ligand-ligand hopping parameters as described
hereafter.

The Wannier90-obtained tight-binding Hamiltonian was a 25 $\times$ 25 matrix
consisting of energies $E_{m,n}=\left\langle
w_{m}(\textbf{r})|\hat{H}_{1}|w_{n}(\textbf{r})\right\rangle$, where $m$ and
$n$ respectively denote row and column indices and $w_{m}(\textbf{r})$ is the
$m^{th}$ MLWF in the primitive cell. Our MLWF basis included the $d$-orbitals
of both Ni atoms, the $p$-orbitals of each of the three distinct sulfurs
coordinated to Ni1 (see Figure SXII), and the $p$-orbitals of the two P atoms
in our simulation cell, for completeness. The radial spreads of the Ni
$d$-shaped MLWF’s ranged from 0.6-1.1 Å and those 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
$\times$ 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 NiPS3 yields many possible $pd\sigma$ and $pd\pi$ values upon solving for
hopping parameters based upon different Ni-S pairs. To remove this ambiguity,
we computed our $pd\sigma$ and $pd\pi$ values based upon the two largest
energies $\left\langle
w_{i}(\textbf{r})|\hat{H}|w_{j}(\textbf{r})\right\rangle$ 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 NiPS3 where $d_{b}$ is the
corresponding bond length (Å).

Table SIII: Predicted NiPS3 metal-ligand hopping parameters in eV and their corresponding largest energies. Sys | pair |  $E_{i}$, $E_{j}$ |  $pd\sigma$ |  $pd\pi$ |  $d_{b}$ (Å)
---|---|---|---|---|---
bulk | Ni1-S3 |  $E_{x,xz}$, $E_{y,xy}$ |  -1.07 |  0.81 | 2.384
bulk | Ni1-S4 |  $E_{z,z2}$, $E_{y,z2}$ |  -0.92 |  0.77 | 2.379
bulk | Ni1-S6 |  $E_{x,x2y2}$, $E_{z,z2}$ |  -1.23 |  0.43 | 2.400
ML | Ni1-S3 |  $E_{y,xy}$, $E_{x,z2}$ |  -1.01 |  0.57 | 2.453
ML | Ni1-S4 |  $E_{x,xz}$, $E_{z,x2y2}$ |  -1.08 |  0.32 | 2.452
ML | Ni1-S6 |  $E_{x,x2y2}$, $E_{z,z2}$ |  -0.71 |  0.50 | 2.453

As an example, solving for $pd\sigma$ and $pd\pi$ for the bulk Ni1-S3 pair
entailed solving the system of equations

$\displaystyle
E^{Ni-S}_{x,xz}=\sqrt{3}l_{1}^{2}n_{1}(pd\sigma)+n_{1}(1-2l_{1}^{2})(pd\pi)$
(2) $\displaystyle
E^{Ni-S}_{y,xy}=\sqrt{3}m_{1}^{2}l_{1}(pd\sigma)+l_{1}(1-2m_{1}^{2})(pd\pi)$
(3)

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\sigma$
and $pp\pi$ (and thus $T_{pp}=pp\sigma-pp\pi$) were similarly determined via

$\displaystyle E^{S-S}_{x,x}=l_{2}^{2}(pp\sigma)+(1-l_{2}^{2})(pp\pi)$ (4)
$\displaystyle E^{S-S}_{x,y}=l_{2}m_{2}(pp\sigma)-l_{2}m_{2}(pp\pi)$ (5)

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\sigma$, $pp\pi$, 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\sigma$ and
$pd\pi$ 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\sigma$ and corresponding lack of sensitivity of $pd\pi$
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 NiPS3 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\sigma$ by 0.14 eV and $pd\pi$
by 0.21 eV.

Table SIV: Predicted “bulk” with vacuum NiPS3 metal-ligand hopping parameters in eV and their corresponding largest energies. Sys | pair |  $E_{i}$, $E_{j}$ |  $pd\sigma$ |  $pd\pi$ |  $d_{b}$ (Å)
---|---|---|---|---|---
bulk+vac | Ni1-S3 |  $E_{x,xz}$, $E_{y,xy}$ |  -1.16 |  0.64 | 2.384
bulk+vac | Ni1-S4 |  $E_{z,z2}$, $E_{z,yz}$ |  -0.96 |  0.62 | 2.379
bulk+vac | Ni1-S6 |  $E_{x,x2y2}$, $E_{z,z2}$ |  -1.97 |  0.30 | 2.400

###  IV.2 Density Functional Theory Methodology

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, $\sigma(\omega)$,
of bulk (zigzag AFM) NiPS3 are also reproduced when $\sigma(\omega)$ 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 NiPS3, 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 NiPS3[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
NiPS3 using the Neél-obtained plane-wave cutoff and k-grid parameters.

####  IV.2.1 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
IV.2.2) with the Neél-antiferromagnetic configuration, using the PBE+$U$
family of functionals.[19, 20] A k-point mesh of $6\times 6\times 6$ for the
bulk and $4\times 4\times 2$ for the monolayer was sufficient to converge
their respective total energies to within 5$\times 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 Fig. SXI) despite the use of
ultra-soft ccECP’s designed for plane-wave calculations.[21]

Figure SXI: PBE+$U(=2$ eV$)$-calculated total energies ($8\times 8\times 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_{cut}^{PW}$.

Due to the steep computational expense expected to accompany non-self
consistent field and band structure calculations of NiPS3 with plane wave
cutoffs of 600 Ry, we also performed tests to explicitly see how the band
structure of bulk Neél AFM NiPS3 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.

####  IV.2.2 Obtaining Bulk and Monolayer Crystal Structures

Figure 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).

To obtain relaxed bulk and monolayer NiPS3 crystal structures (see Figure SXII
and Table SV), the bulk structure with experimentally-determined lattice
parameters having a monoclinic ($C2/m$) symmetry and a 10-atom unit cell
consisting of two NiPS3 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
NiPS3 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 $P1$
symmetry on which our DFT and Wannierization was performed (Figure SXIIa).

Table SV: Relaxed NiPS3 crystal structure parameters. |  $a$ (Å) |  $b$ (Å) |  $d_{Ni-Ni}$ (Å) | $\overline{d}_{Ni-S}$ |  $d_{P-P}$ (Å)
---|---|---|---|---|---
bulk | 5.8114 | 10.0640 | 3.3532 | 2.39(1) | 2.1848
ML | 5.8507 | 10.1347 | 3.3781 | 2.452(1) | 2.1927

One monolayer of NiPS3 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 variable-cell relaxation since no experimental
monolayer lattice parameters exist, we performed variable-cell relaxations of
bulk NiPS3 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
NiPS3 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+$U$-converged 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 NiPS3 crystal structure
(having $P1$ symmetry) for our monolayer electronic structure and MLWF
analysis (see Figure SXIIb).

Table SVI: Bulk NiPS3 crystal structure parameters after full-cell relaxation with different DFT functionals, $V_{xc}$, compared with those from experiment. Where a $U$ is specified in the table, we employed $U=4$ eV. | $V_{xc}$ |  $a$ (Å) |  $b$ (Å) |  $c$ (Å) |  $\beta$ (°)
---|---|---|---|---|---
Exp.[1] | – | 5.811 | 10.064 | 6.896 | 106.22
DFT | PBE | 5.855 | 10.140 | 7.149 | 105.81
DFT | PBE+$U$ | 5.906 | 10.231 | 7.223 | 105.79
DFT | LDA | 5.727 | 9.929 | 6.262 | 107.64
DFT | LDA+$U$ | 5.734 | 9.934 | 6.408 | 107.25
Figure SXIII: Band structures of Neél-AFM bulk NiPS3 for PW cutoffs of (a) 250
Ry, and (b) 600 Ry (fully converged). Crystal structure was vc-relaxed with
PBE+$U(=4$ eV)

##  V Superexchange Expressions

The exact ground state wavefunction within an NiS6 cluster is given by
$\ket{\Psi_{g}}=\alpha\ket{3d^{8}}+\beta\ket{3d^{9}\underline{L}^{1}}+\gamma\ket{3d^{10}\underline{L}^{2}}$.
Due to the magnitude of the trigonal field obtained from modeling
($\delta=-80$ meV), we simplify the calculation by approximating Ni-S bonds
within an NiS6 cluster as orthogonal. As mentioned in the main text, the
ground state character for bulk and 3L NiPS3 finds a nearly equal mixture of
$\ket{3d^{8}}$ and $\ket{3d^{9}\underline{L}^{1}}$
($|\alpha|^{2}\approx|\beta|^{2}$), with a negligibly small
$\ket{3d^{10}\underline{L}^{2}}$ character. Thus, we independently apply the
cell-perturbation to these two states. For fully occupied (inactive) $t_{2g}$
orbitals, $\ket{3d^{8}}$ and $\ket{3d^{9}\underline{L}^{1}}$ are given by

$\begin{split}&\ket{3d^{8}}=\ket{d_{x^{2}-y^{2}}d_{3z^{2}-r^{2}}}\\\
&\ket{3d^{9}\underline{L}^{1}}=\frac{1}{\sqrt{2}}(\ket{d_{x^{2}-y^{2}}L_{3z^{2}-r^{2}}}+\ket{L_{x^{2}-y^{2}}d_{x^{2}-y^{2}}}\end{split}$
(6)

where $d_{x^{2}-y^{2}}$ and $d_{3z^{2}-r^{2}}$ denote holes in the two Ni 3d
$e_{g}$ orbitals, and $L_{x^{2}-y^{2}}$ and $L_{3z^{2}-r^{2}}$ denote holes in
the S $3p$ molecular orbitals. These molecular orbitals are given by
$L_{x^{2}-y^{2}}=$ $\frac{1}{2}$($p^{1}_{x}$ $-$ $p^{2}_{y}$ $+$ $p^{3}_{y}$
$-$ $p^{4}_{x}$) and $L_{3z^{2}-r^{2}}=$ $\frac{1}{\sqrt{3}}$($p^{5}_{z}$ $-$
$p^{6}_{z}$) $+$ $\frac{1}{2\sqrt{3}}($$p^{1}_{x}$ $+$ $p^{2}_{y}$ $-$
$p^{3}_{y}$ $-$ $p^{4}_{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 two-center atomic overlap integrals within the
Slater-Koster scheme of linear combinations of atomic orbitals. The
superexchange for $\ket{3d^{8}}$ and $\ket{3d^{9}\underline{L}^{1}}$ states
are defined as $J_{i}^{\alpha}$ and $J_{i}^{\beta}$, where i $\in\\{1,2,3\\}$.

The 1NN superexchange is given by the fourth order perturbation term in
$\ket{3d^{8}}$, and second order term in $\ket{3d^{9}\underline{L}^{1}}$.
These superexchange expressions are

$\begin{split}&J_{1}^{\alpha}=\frac{-\alpha^{2}J_{H}^{S}(pd\sigma)^{4}}{\Delta^{\prime\prime}(\Delta^{\prime\prime}-J_{H}^{S}/2)(\Delta^{\prime}+U_{pd})^{2}}\\\
&J_{1}^{\beta}=\frac{\beta^{2}}{2U_{p}}\biggl{[}\frac{1}{16}[(pp\sigma)+(pp\pi)]^{2}+\frac{1}{144}[(pp\sigma)\\\
&\qquad\qquad\qquad\qquad\qquad\quad+9(pp\pi)]^{2}-\frac{J_{H}^{S}}{9}\biggr{]}\end{split}$
(7)

where $\Delta^{\prime}=\Delta+2U_{pd}$,
$\Delta^{\prime\prime}=\Delta^{\prime}+U_{p}/2$, and $J_{H}^{S}$ is the Hund’s
coupling between S 3p holes on a shared S site, assumed to be $\sim 0.7$ eV
[23]. $U_{p}$ is the on-site Coulomb interaction for S 3p orbitals and is set
to 4 eV [24]. Note, here $U_{pd}=1$ eV and is the Coulomb interaction term
between S 3p and Ni 3d and is not the same as the $U_{pd}$ shown in the main
text. For bulk NiPS3, we find $J_{1}^{\alpha}$ = -1.8 meV and $J_{1}^{\beta}$
= -2.2 meV, then $J_{1}^{\alpha}$ = -1.6 meV and $J_{1}^{\beta}$ = -2.9 meV in
3L.

The 2NN superexchange is now given by the sixth order perturbation term in
$\ket{3d^{8}}$, and remains the second order perturbation term in
$\ket{3d^{9}\underline{L}^{1}}$. These expressions are

$\begin{split}&J_{2}^{\alpha}=\alpha^{2}\Gamma\;(pp\pi)^{2}(pd\sigma)^{4}\\\
&J_{2}^{\beta}=\frac{\beta^{2}}{18U_{p}}(pp\pi)^{2}\\\
&\Gamma=\frac{2(\Delta^{\prime}-U_{pd})^{2}+\Delta^{\prime
2}}{2(\Delta^{\prime}-U_{pd})^{4}(\Delta^{\prime\prime}-U_{pd})\Delta^{\prime
2}}\end{split}$ (8)

These give $J_{2}^{\alpha}$ = 0.2 meV and $J_{2}^{\beta}$ = 0.052 meV for
bulk, while $J_{2}^{\alpha}$ = 0.0021 meV and $J_{2}^{\beta}$ = 0.0003 meV for
3L NiPS3. Lastly, the 3NN superexchange is again given by the sixth order
perturbation term in $\ket{3d^{8}}$ and second order in
$\ket{3d^{9}\underline{L}^{1}}$. We thus have

$\begin{split}&J_{3}^{\alpha}=2\alpha^{2}\Gamma\;(pd\sigma)^{4}\biggl{[}\frac{1}{4}[(pp\sigma)-(pp\pi)]^{2}+(pp\pi)^{2}\biggr{]}\\\
&J_{3}^{\beta}=\frac{\beta^{2}}{24U_{p}}[(pp\sigma)-(pp\pi)]^{2}\end{split}$
(9)

For bulk NiPS3, these expressions give $J_{3}^{\alpha}$ = 12 meV and
$J_{3}^{\beta}$ = 4.6 meV, while for 3L we find $J_{3}^{\alpha}$ = 8.6 meV and
$J_{3}^{\beta}$ = 1.8 meV.

## References

  * Wildes  _et al._ [2015] 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, and T. J. Hicks, “Magnetic structure of the quasi-two-dimensional antiferromagnet NiPS 3,” Phys. Rev. B 92, 180401(R) (2015).
  * Lançon  _et al._ [2018] D. Lançon, R. A. Ewings, T. Guidi, F. Formisano, and A. R. Wildes, “Magnetic exchange parameters and anisotropy of the quasi-two-dimensional antiferromagnet NiPS 3,” Phys. Rev. B 98, 134414 (2018).
  * Kuo  _et al._ [2016] Cheng-Tai Kuo, Michael Neumann, Karuppannan Balamurugan, Hyun Ju Park, Soonmin Kang, Hung Wei Shiu, Jin Hyoun Kang, Byung Hee Hong, Moonsup Han, Tae Won Noh, and Je-Geun Park, “Exfoliation and Raman Spectroscopic Fingerprint of Few-Layer NiPS3 Van der Waals Crystals,” Sci Rep 6, 20904 (2016).
  * Lu  _et al._ [2020] Haoyong Lu, Wei Wang, Yun Liu, Limin Chen, Qiyun Xie, Handi Yin, Guofeng Cheng, and Liang He, “Exfoliation, lattice vibration and air stability characterization of antiferromagnetic van der Waals NiPS3 nanosheets,” Applied Surface Science 504, 144405 (2020).
  * Tanabe and Sugano [1954] Yukito Tanabe and Satoru Sugano, “On the Absorption Spectra of Complex Ions. I,” J. Phys. Soc. Jpn. 9, 753–766 (1954).
  * Ghiringhelli  _et al._ [2005] G Ghiringhelli, M Matsubara, C Dallera, F Fracassi, R Gusmeroli, A Piazzalunga, A Tagliaferri, N B Brookes, A Kotani, and L Braicovich, “NiO as a test case for high resolution resonant inelastic soft x-ray scattering,” J. Phys.: Condens. Matter 17, 5397–5412 (2005).
  * Haverkort  _et al._ [2012] M. W. Haverkort, M. Zwierzycki, and O. K. Andersen, “Multiplet ligand-field theory using Wannier orbitals,” Phys. Rev. B 85, 165113 (2012).
  * Sakuma [2013] R. Sakuma, “Symmetry-adapted Wannier functions in the maximal localization procedure,” Phys. Rev. B 87, 235109 (2013).
  * Marzari and Vanderbilt [1997] Nicola Marzari and David Vanderbilt, “Maximally localized generalized Wannier functions for composite energy bands,” Phys. Rev. B 56, 12847–12865 (1997).
  * Andersen and Saha-Dasgupta [2000] O. K. Andersen and T. Saha-Dasgupta, “Muffin-tin orbitals of arbitrary order,” Phys. Rev. B 62, R16219–R16222 (2000).
  * Pizzi  _et al._ [2020] Giovanni Pizzi, Valerio Vitale, Ryotaro Arita, Stefan Blügel, Frank Freimuth, Guillaume Géranton, Marco Gibertini, Dominik Gresch, Charles Johnson, Takashi Koretsune, Julen Ibañez-Azpiroz, Hyungjun Lee, Jae-Mo Lihm, Daniel Marchand, Antimo Marrazzo, Yuriy Mokrousov, Jamal I Mustafa, Yoshiro Nohara, Yusuke Nomura, Lorenzo Paulatto, Samuel Poncé, Thomas Ponweiser, Junfeng Qiao, Florian Thöle, Stepan S Tsirkin, Małgorzata Wierzbowska, Nicola Marzari, David Vanderbilt, Ivo Souza, Arash A Mostofi, and Jonathan R Yates, “Wannier90 as a community code: New features and applications,” J. Phys.: Condens. Matter 32, 165902 (2020).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, “Simplified LCAO Method for the Periodic Potential Problem,” Phys. Rev. 94, 1498–1524 (1954).
  * Podolskiy and Vogl [2004] A. V. Podolskiy and P. Vogl, “Compact expression for the angular dependence of tight-binding Hamiltonian matrix elements,” Phys. Rev. B 69, 233101 (2004).
  * Giannozzi  _et al._ [2009] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, Andrea Dal Corso, Stefano de Gironcoli, Stefano Fabris, Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, Anton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitsonen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch, “QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials,” J. Phys.: Condens. Matter 21, 395502 (2009).
  * Gu  _et al._ [2019] Yuhao Gu, Qiang Zhang, Congcong Le, Yinxiang Li, Tao Xiang, and Jiangping Hu, “Ni-based transition metal trichalcogenide monolayer: A strongly correlated quadruple-layer graphene,” Phys. Rev. B 100, 165405 (2019).
  * Kim  _et al._ [2019] Heung-Sik Kim, Kristjan Haule, and David Vanderbilt, “Mott Metal-Insulator Transitions in Pressurized Layered Trichalcogenides,” Phys. Rev. Lett. 123, 236401 (2019).
  * Kim  _et al._ [2018] So Yeun Kim, Tae Yun Kim, Luke J. Sandilands, Soobin Sinn, Min-Cheol Lee, Jaeseok Son, Sungmin Lee, Ki-Young Choi, Wondong Kim, Byeong-Gyu Park, C. Jeon, Hyeong-Do Kim, Cheol-Hwan Park, Je-Geun Park, S. J. Moon, and T. W. Noh, “Charge-Spin Correlation in van der Waals Antiferromagnet NiPS 3,” Phys. Rev. Lett. 120, 136402 (2018).
  * Lane and Zhu [2020] Christopher Lane and Jian-Xin Zhu, “Thickness dependence of electronic structure and optical properties of a correlated van der Waals antiferromagnetic NiPS 3 thin film,” Phys. Rev. B 102, 075124 (2020).
  * Anisimov  _et al._ [1991] Vladimir I. Anisimov, Jan Zaanen, and Ole K. Andersen, “Band theory and Mott insulators: Hubbard _U_ instead of Stoner _I_ ,” Phys. Rev. B 44, 943–954 (1991).
  * Ernzerhof and Scuseria [1999] Matthias Ernzerhof and Gustavo E. Scuseria, “Assessment of the Perdew–Burke–Ernzerhof exchange-correlation functional,” The Journal of Chemical Physics 110, 5029–5036 (1999).
  * Annaberdiyev  _et al._ [2018] Abdulgani Annaberdiyev, Guangming Wang, Cody A. Melton, M. Chandler Bennett, Luke Shulenburger, and Lubos Mitas, “A new generation of effective core potentials from correlated calculations: 3d transition metal series,” The Journal of Chemical Physics 149, 134108 (2018).
  * Jain  _et al._ [2013] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, and Kristin A. Persson, “Commentary: The Materials Project: A materials genome approach to accelerating materials innovation,” APL Materials 1, 011002 (2013).
  * Autieri  _et al._ [2022] Carmine Autieri, Giuseppe Cuono, Canio Noce, Milosz Rybak, Kamila M. Kotur, Cliò Efthimia Agrapidis, Krzysztof Wohlfeld, and Magdalena Birowska, “Limited Ferromagnetic Interactions in Monolayers of MPS 3 (M = Mn and Ni),” J. Phys. Chem. C 126, 6791–6802 (2022).
  * Babuka  _et al._ [2017] T. Babuka, K. Glukhov, Y. Vysochanskii, and M. Makowska-Janusik, “New insight into strong correlated states realised in a ferroelectric and paraelectric chalcogenide Sn 2 P 2 S 6 crystal,” RSC Adv. 7, 27770–27779 (2017).